diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 96a10f616a..f1b81b7c00 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3295,7 +3295,7 @@ Enables the "new" SILHS importance sampling scheme with prescribed probabilities -Determine starting SILHS first sampling level (k_lh_start) based on maximum within-cloud rcm. If false, and if l_random_k_lh_start is also false, then the SILHS first sampling level is the column maximum of liquid cloud water. +Determine starting SILHS first sampling level (k_lh_start) based on maximum within-cloud rcm. If false, and if l_random_k_lh_start is also false, then the SILHS first sampling level is the column maximum of liquid cloud water. Flag to use shear in the calculation of Richardson number. -Default: .false. +Default: .false. -Flag to run CLUBB with E3SM settings. +Flag to run CLUBB with E3SM settings. Default: .true. @@ -4039,34 +4039,34 @@ Default: .false. Flag to use Total Kenetic Energy (TKE) in eddy diffusion for wp2 and wp3. -Default: .false. +Default: .false. -Flag to use Total Kenetic Energy (TKE) formulation for wp3 pr_turb (turbulent +Flag to use Total Kenetic Energy (TKE) formulation for wp3 pr_turb (turbulent production) term. -Default: .false. +Default: .false. Flag to use smooth Heaviside 'Peskin' in computation of invrs_tau. -Default: .false. +Default: .false. -This flag determines whether we want to use an upwind differencing approximation -rather than a centered differencing for turbulent advection terms. It affects +This flag determines whether we want to use an upwind differencing approximation +rather than a centered differencing for turbulent advection terms. It affects wpxp only. Default: .false. -This flag determines whether we want to use an upwind differencing approximation -rather than a centered differencing for turbulent advection terms. It affects +This flag determines whether we want to use an upwind differencing approximation +rather than a centered differencing for turbulent advection terms. It affects xpyp only. Default: .false. @@ -5551,6 +5551,79 @@ Force scam to use the lat lon fields specified in the scam namelist not what is Default: FALSE + +SCAM to calculate or read tendencies from a global analysis/dycore +Default: FALSE + + + +Use 1st order upwind for analysis tendencies (instead of 2nd order space centered) +Default: FALSE + + + +Use scam state as center column in stencil for analysis advective tendencies +Default: FALSE + + + +Use scam state as center column in stencil for analysis advective tendencies +Default: FALSE + + + +Use scam state as center column in stencil for analysis advective tendencies +Default: FALSE + + + +Use scam state as center column in stencil for analysis advective tendencies +Default: FALSE + + + +Force scam to use tendencies directly from dycore or analysis (not recalculated) +Default: FALSE + + + +Force scam to use omega directly from dycore or analysis (not recalculated) +Default: FALSE + + + +Interpolate analysis fields to constant pressure surfaces +Default: FALSE + + + + +template for analysis forcing dataset. +Default: set by build-namelist. + + + +templatefull path for analysis forcing dataset. +Default: set by build-namelist. + + + +Force scam to compute large-scale forcing from renalysis or 3D model output +Default: FALSE + + = ptop_ref .and. ptop_ref > 1.e-5_r8 ) then ! around waccm top, below top of waccmx cnst_fixed_ubc(1) = .true. else if ( 1.e1_r8 > ptop_ref .and. ptop_ref > 1.e-2_r8 ) then ! well above top of cam and below top of waccm - call endrun('chem_init: do not know how to set water vapor upper boundary when model top is near mesopause') + if(.not.(single_column)) then + call endrun('chem_init: do not know how to set water vapor upper boundary when model top is near mesopause') + else + cnst_fixed_ubc(1) = .true. + endif endif if ( masterproc ) write(iulog,*) 'chem_init: addfld done' diff --git a/src/chemistry/mozart/upper_bc.F90 b/src/chemistry/mozart/upper_bc.F90 index 71a4a65b0c..8b076623e7 100644 --- a/src/chemistry/mozart/upper_bc.F90 +++ b/src/chemistry/mozart/upper_bc.F90 @@ -157,11 +157,16 @@ subroutine ubc_init() use mo_snoe, only: snoe_inti use mo_msis_ubc, only: msis_ubc_inti use constituents,only: cnst_get_ind + use scamMod,only: single_column !---------------------------Local workspace----------------------------- logical :: zonal_avg !----------------------------------------------------------------------- - apply_upper_bc = ptop_ref<1._r8 ! Pa + if (single_column) then + apply_upper_bc = .FALSE. + else + apply_upper_bc = ptop_ref<1._r8 ! Pa + endif if (.not.apply_upper_bc) return diff --git a/src/control/history_scam.F90 b/src/control/history_scam.F90 index 2c81ce1a78..5088de5bce 100644 --- a/src/control/history_scam.F90 +++ b/src/control/history_scam.F90 @@ -45,7 +45,10 @@ subroutine scm_intht() call addfld ('UDIFF', (/ 'lev' /), 'A', 'K','difference from observed u wind', gridname='gauss_grid') call addfld ('VDIFF', (/ 'lev' /), 'A', 'K','difference from observed v wind', gridname='gauss_grid') - call addfld ('TOBS', (/ 'lev' /), 'A', 'K','observed temp') + call addfld ('TOBS', (/ 'lev' /), 'A', 'K','observed temp', gridname='gauss_grid') + call addfld ('UOBS', (/ 'lev' /), 'A', 'm/s','observed zonal wind', gridname='gauss_grid') + call addfld ('VOBS', (/ 'lev' /), 'A', 'm/s','observed meridional wind', gridname='gauss_grid') + call addfld ('QDIFF', (/ 'lev' /), 'A', 'kg/kg','difference from observed water', gridname='gauss_grid') call addfld ('QOBS', (/ 'lev' /), 'A', 'kg/kg','observed water', gridname='physgrid') @@ -100,6 +103,58 @@ subroutine scm_intht() call addfld ('NLTEN_PHYS', (/ 'lev' /), 'I','#/kg/s', 'NL vertical advective forcing', gridname='gauss_grid' ) call addfld ('NITEN_PHYS', (/ 'lev' /), 'I','#/kg/s', 'NI vertical advective forcing', gridname='gauss_grid' ) + call addfld ('U_IOP', (/ 'lev' /), 'I', 'm/s', 'Zonal Wind from IOP ', gridname='gauss_grid' ) + call addfld ('V_IOP', (/ 'lev' /), 'I', 'm/s', 'Mer. Wind from IOP ', gridname='gauss_grid' ) + call addfld ('OMEGA_IOP', (/ 'lev' /), 'I', 'Pa/s', 'Vertical velocity (from IOP) ', gridname='gauss_grid' ) + call addfld ('OMEGA_ANA', (/ 'lev' /), 'I', 'Pa/s', 'Vertical velocity (analysis) ', gridname='gauss_grid' ) + call addfld ('ETAD_ANA', (/ 'lev' /), 'I', 'Pa/s', 'Eta_dot (analysis) ', gridname='gauss_grid' ) + call addfld ('ZETA_ANA', (/ 'lev' /), 'I', '1/s', 'Rel. Vorticity (analysis) ', gridname='gauss_grid' ) + call addfld ('T_ANA', (/ 'lev' /), 'I', 'K', 'Temperature (analysis) ', gridname='gauss_grid' ) + call addfld ('Q_ANA', (/ 'lev' /), 'I', 'g/g', 'Spec. humidity (analysis) ', gridname='gauss_grid' ) + call addfld ('U_ANA', (/ 'lev' /), 'I', 'm/s', 'Zonal wind (analysis) ', gridname='gauss_grid' ) + call addfld ('V_ANA', (/ 'lev' /), 'I', 'm/s', 'Mer. Wind (analysis) ', gridname='gauss_grid' ) + call addfld ('TV_ANA', (/ 'lev' /), 'I', 'K', 'Temperature (analysis) ', gridname='gauss_grid' ) + call addfld ('TTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('QTEN_TOTDYN_ANA', (/ 'lev' /), 'I', 'kg/kg/s', 'tracer tendency (analysis)', gridname='gauss_grid' ) + + call addfld ('UTEN_TOTDYN_ANAR', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_TOTDYN_ANAR', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + + call addfld ('TTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('OMEGA_DYCORE_ANA', (/ 'lev' /), 'I', 'Pa/s','Pressure tendency/velocity (analysis)', gridname='gauss_grid' ) + call addfld ('OMEGA_RECALC_ANA', (/ 'lev' /), 'I', 'Pa/s','Pressure tendency/velocity (analysis)', gridname='gauss_grid' ) + + call addfld ('UTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_PRG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_PHIG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_KEG_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_VORT_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_PFRC_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_VADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_HADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('UTEN_CORIOL', (/ 'lev' /), 'I', 'm/s/s', 'Zonal wind tendency (analysis)', gridname='gauss_grid' ) + + + call addfld ('VTEN_DYCORE_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_VORT_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_PFRC_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_VADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_HADV_ANA', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + call addfld ('VTEN_CORIOL', (/ 'lev' /), 'I', 'm/s/s', 'Meridional wind tendency (analysis)', gridname='gauss_grid' ) + + call addfld ('TTEN_VADV_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('TTEN_HADV_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('TTEN_COMP_ANA', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('TTEN_COMP_IOP', (/ 'lev' /), 'I', 'K/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + + call addfld ('QTEN_VADV_ANA', (/ 'lev' /), 'I', '1/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + call addfld ('QTEN_HADV_ANA', (/ 'lev' /), 'I', '1/s', 'Temperature tendency (analysis)', gridname='gauss_grid' ) + + call addfld ('PS_ANA', horiz_only, 'A', 'Pa','Ana. Surface Pressure', gridname='gauss_grid') + call addfld ('PHIS_ANA', horiz_only, 'A', 'm2/s2','Ana. Surface geopotential', gridname='gauss_grid') + end subroutine scm_intht !####################################################################### diff --git a/src/control/scamMod.F90 b/src/control/scamMod.F90 index b18169b340..97b62325c6 100644 --- a/src/control/scamMod.F90 +++ b/src/control/scamMod.F90 @@ -76,6 +76,20 @@ module scamMod character*(max_path_len), public :: lsmsurffile character*(max_path_len), public :: lsminifile +logical, public :: use_scm_ana_frc = .false. +character*(max_path_len), public :: scm_ana_frc_file_template +character*(max_path_len), public :: scm_ana_frc_path + +logical, public :: scm_ana_x_plevels = .false. +logical, public :: scm_ana_direct_omega = .false. +logical, public :: scm_ana_direct_ttend = .false. +logical, public :: scm_ana_t_react = .false. +logical, public :: scm_ana_q_react = .false. +logical, public :: scm_ana_u_react = .false. +logical, public :: scm_ana_v_react = .false. +logical, public :: scm_ana_upwind = .false. +logical, public :: scm_use_ana_iop = .false. + ! note that scm_zadv_q is set to slt to be consistent with CAM BFB testing @@ -250,7 +264,11 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in) scm_cambfb_mode,scm_crm_mode,scm_zadv_uv,scm_zadv_T,scm_zadv_q,& scm_use_obs_T, scm_use_obs_uv, scm_use_obs_qv, & scm_relax_linear, scm_relax_tau_top_sec, & - scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, scm_backfill_iop_w_init + scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, scm_backfill_iop_w_init, & + use_scm_ana_frc, scm_ana_frc_path, scm_ana_frc_file_template, & + scm_ana_x_plevels, scm_ana_direct_omega, & + scm_ana_t_react, scm_ana_q_react, scm_ana_u_react, scm_ana_v_react, & + scm_ana_upwind, scm_ana_direct_ttend, scm_use_ana_iop single_column=single_column_in @@ -306,6 +324,7 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in) use_camiop = .false. endif + ! If we are not forcing the lat and lon from the namelist use the closest lat and lon that is found in the IOP file. if (.not.scm_force_latlon) then call shr_scam_GetCloseLatLon( ncid, scmlat, scmlon, ioplat, ioplon, latidx, lonidx ) @@ -316,7 +335,8 @@ subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in) scmlat = ioplat scmlon = ioplon end if - + + if (masterproc) then write (iulog,*) 'Single Column Model Options: ' write (iulog,*) '=============================' diff --git a/src/dynamics/eul/dyn_comp.F90 b/src/dynamics/eul/dyn_comp.F90 index 442c9f3228..9450a73aa9 100644 --- a/src/dynamics/eul/dyn_comp.F90 +++ b/src/dynamics/eul/dyn_comp.F90 @@ -842,8 +842,8 @@ subroutine process_inidat(fieldname, m_cnst, fh) ret = pio_inq_varid(fh, cnst_name(m_cnst), varid) ret = pio_get_att(fh, varid, 'units', trunits) if (trunits(1:5) .ne. 'KG/KG' .and. trunits(1:5) .ne. 'kg/kg') then - call endrun(sub//': ERROR: Units for tracer ' & - //trim(cnst_name(m_cnst))//' must be in KG/KG') + !call endrun(sub//': ERROR: Units for tracer ' & + ! //trim(cnst_name(m_cnst))//' must be in KG/KG') end if else if (.not. analytic_ic_active()) then diff --git a/src/dynamics/eul/get_ana_dynfrc_4scam.F90 b/src/dynamics/eul/get_ana_dynfrc_4scam.F90 new file mode 100644 index 0000000000..44dc6dee88 --- /dev/null +++ b/src/dynamics/eul/get_ana_dynfrc_4scam.F90 @@ -0,0 +1,1557 @@ +module get_ana_dynfrc_4scam + + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8, & + cs=>SHR_KIND_CS,cl=>SHR_KIND_CL + use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters + pi => shr_const_pi , & + OOmega => shr_const_omega , & + rdair => shr_const_rdair , & + cpair => shr_const_cpdair + + use scamMod, only: use_scm_ana_frc, & + scm_ana_frc_path, & + scm_ana_frc_file_template, & + scm_ana_x_plevels, & + scm_ana_direct_omega, & + scm_ana_t_react, & + scm_ana_q_react, & + scm_ana_u_react, & + scm_ana_v_react, & + scm_ana_upwind, & + scm_ana_direct_ttend + + + + ! shr_const_mod is in ${CESMROOT}/cime/src/share/util/ + + implicit none + private + save + + public get_ana_dynfrc_fv +! +! Private module data +! + + real(r8) , save , allocatable :: T_1(:,:,:) , U_1(:,:,:), V_1(:,:,:), Q_1(:,:,:),PS_1(:,:),PHIS_1(:,:) + real(r8) , save , allocatable :: T_2(:,:,:) , U_2(:,:,:), V_2(:,:,:), Q_2(:,:,:),PS_2(:,:),PHIS_2(:,:) + real(r8) , save , allocatable :: UTCORE_1(:,:,:) , UTCORE_2(:,:,:) + real(r8) , save , allocatable :: VTCORE_1(:,:,:) , VTCORE_2(:,:,:) + real(r8) , save , allocatable :: TTCORE_1(:,:,:) , TTCORE_2(:,:,:) + real(r8) , save , allocatable :: OGCORE_1(:,:,:) , OGCORE_2(:,:,:) + real(r8) , save , allocatable :: lat_ana(:),lon_ana(:),lev_ana(:) + integer , save :: nlev_ana, nlon_ana, nlat_ana + + real(r8) , save , allocatable :: To_1(:,:,:) , Uo_1(:,:,:), Vo_1(:,:,:), Qo_1(:,:,:),PSo_1(:,:),PHISo_1(:,:) + real(r8) , save , allocatable :: To_2(:,:,:) , Uo_2(:,:,:), Vo_2(:,:,:), Qo_2(:,:,:),PSo_2(:,:),PHISo_2(:,:) + real(r8) , save , allocatable :: UTCOREo_1(:,:,:) , UTCOREo_2(:,:,:), UTCOREo_X(:,:,:) + real(r8) , save , allocatable :: VTCOREo_1(:,:,:) , VTCOREo_2(:,:,:), VTCOREo_X(:,:,:) + real(r8) , save , allocatable :: TTCOREo_1(:,:,:) , TTCOREo_2(:,:,:), TTCOREo_X(:,:,:) + real(r8) , save , allocatable :: OGCOREo_1(:,:,:) , OGCOREo_2(:,:,:), OGCOREo_X(:,:,:) + + + + real(r8) , save , allocatable :: ETAD_X(:,:,:) , OMG_X(:,:,:) + real(r8) , save , allocatable :: ZETA_X(:) + real(r8) , save , allocatable :: KEh_X(:,:,:) + real(r8) , save , allocatable :: Tv_X(:,:,:) + + real(r8) , save , allocatable :: pke_X(:,:,:),pko_X(:,:,:),phik_X(:,:,:),Thv_X(:,:,:) + real(r8) , save , allocatable :: ple_X(:,:,:) , plo_X(:,:,:), phi_X(:,:,:) + + real(r8) , save , allocatable :: To_X(:,:,:) , Uo_X(:,:,:), Vo_X(:,:,:), Qo_X(:,:,:),PSo_X(:,:),PHISo_X(:,:) + + +!======================================================================= +contains +!======================================================================= + +subroutine get_ana_dynfrc_fv ( scmlon, scmlat , & + omega_ana, etad_ana, zeta_ana, & + t_ana , tv_ana , & + q_ana , & + u_ana , & + v_ana , & + ps_ana , & + phis_ana , & + uten_hadv_ana , & + vten_hadv_ana , & + uten_pfrc_ana , & + vten_pfrc_ana , & + uten_vort_ana , & + vten_vort_ana , & + qten_hadv_ana , & + tten_hadv_ana , & + uten_vadv_ana , & + vten_vadv_ana , & + tten_vadv_ana , & + qten_vadv_ana , & + tten_comp_ana , & + uten_keg_ana , & + uten_phig_ana , & + uten_prg_ana , & + uten_dycore_ana , & + vten_dycore_ana , & + tten_dycore_ana , & + omega_dycore_ana , & + omega_recalc_ana , & + u_scm, v_scm, t_scm, q_scm, & + u_ana_diag, v_ana_diag, t_ana_diag, q_ana_diag ) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! US and VS are input (D-grid velocities) +!-------------------------------------------- +! ub(i,j,L)= 0.5*(us(i-1,j,L) + us(i,j,L)) +! vb(i,j,L)= 0.5*(vs(i,j,L) + vs(i,j+1,L)) +! +! uc(i,j,L)= 0.5*(ub(i,j,L) + ub(i,j-1,L)) +! vc(i,j,L)= 0.5*(vb(i,j-1,L) + vb(i+1,j-1,L)) +!--------------------------------------------- +! Grid arrangement in FV latlon h,i-files +!--------------------------------------------- +! J=NY +! ... +! +! ub,vb(I,J) us(I,J),vc(I,J+1) +! +! +! vs(I,J),uc(I,J) ua,va,T,p(I,J) vs(I+1,J),uc(I+1,J) +! +! +! vc(I,J) +! +! +! ua,va,T,p(I,J-1) +! +! ... +! J=1 ... +!---------------------------------------------- + + use pmgrid, only : plev, plat, plevp, plon + use hycoef, only: hyai, hybi, ps0, hyam, hybm + use filenames, only: interpret_filename_spec + use time_manager, only: timemgr_time_ge,timemgr_time_inc,get_curr_date,get_step_size,is_first_step + use netcdf + use cam_abortutils, only: endrun + use ref_pres, only: pref_mid ! In Pascal + + real(r8), intent(in) :: scmlon, scmlat + real(r8), intent(out) :: omega_ana( plev ) + real(r8), intent(out) :: etad_ana(plev) + real(r8), intent(out) :: t_ana(plev) , tv_ana(plev) + real(r8), intent(out) :: zeta_ana(plev) + real(r8), intent(out) :: u_ana(plev) + real(r8), intent(out) :: v_ana(plev) + real(r8), intent(out) :: q_ana(plev) + real(r8), intent(out) :: ps_ana, phis_ana + real(r8), intent(out) :: uten_hadv_ana( plev ) + real(r8), intent(out) :: vten_hadv_ana( plev ) + real(r8), intent(out) :: uten_pfrc_ana( plev ) + real(r8), intent(out) :: vten_pfrc_ana( plev ) + real(r8), intent(out) :: qten_hadv_ana( plev ) + real(r8), intent(out) :: tten_hadv_ana( plev ) + real(r8), intent(out) :: qten_vadv_ana( plev ) + real(r8), intent(out) :: tten_vadv_ana( plev ) + real(r8), intent(out) :: uten_vadv_ana( plev ) + real(r8), intent(out) :: vten_vadv_ana( plev ) + + real(r8), intent(out) :: tten_comp_ana( plev ) + + real(r8), intent(out) :: uten_keg_ana( plev ) + real(r8), intent(out) :: uten_prg_ana( plev ) + real(r8), intent(out) :: uten_phig_ana( plev ) + real(r8), intent(out) :: uten_vort_ana( plev ) + real(r8), intent(out) :: vten_vort_ana( plev ) + real(r8), intent(out) :: uten_dycore_ana( plev ) + real(r8), intent(out) :: vten_dycore_ana( plev ) + real(r8), intent(out) :: tten_dycore_ana( plev ) + real(r8), intent(out) :: omega_recalc_ana( plev ) + real(r8), intent(out) :: omega_dycore_ana( plev ) + + real(r8), intent(in) :: u_scm(plev) + real(r8), intent(in) :: v_scm(plev) + real(r8), intent(in) :: t_scm(plev) + real(r8), intent(in) :: q_scm(plev) + + real(r8), intent(out) :: u_ana_diag(plev) + real(r8), intent(out) :: v_ana_diag(plev) + real(r8), intent(out) :: t_ana_diag(plev) + real(r8), intent(out) :: q_ana_diag(plev) + + integer, save :: iax, jax + integer, save :: Read_year2, Read_month2, Read_day2, Read_sec2, Read_YMD2 + integer, save :: nlev_alc, nlon_alc, nlat_alc + + !!logical , parameter :: l_vectinv = .FALSE. + !!real(r8) :: tv_ana(plev) + real(r8) :: rho_ana( plev ), plo_ana(plev) + + + + real(r8) :: scmlonx + + real(r8) :: ana_wgt1 , ana_wgt2 , dx0, dy, darea + real(r8) :: tend_wgt1 , tend_wgt2 + + integer :: nx, ny,i,j,k,L,LM, iav(1),jav(1),iac,jac + + real(r8) , allocatable :: rlats(:),rlons(:) + real(r8) :: zeta(plev),absvo(plev) + ! Horz. gradient profiles (1=X, 2=Y) + real(r8) :: kehg_ana(plev,2),kehg_X(plev,2) + real(r8) :: phig_ana(plev,2),phig_X(plev,2) + real(r8) :: plog_ana(plev,2),plog_X(plev,2) + real(r8) :: teg_ana(plev,2), teg_X(plev,2) + real(r8) :: qg_ana(plev,2), qg_X(plev,2) + real(r8) :: ug_ana(plev,2), ug_X(plev,2) + real(r8) :: vg_ana(plev,2), vg_X(plev,2) + real(r8) :: lin_pfc_ana(plev,2) , lin_pfc_X(plev,2) + + real(r8) :: omega_ana_x(plev) + real(r8) :: alpha_react(plev) + + real(r8) :: lat_alc(3) , lon_alc(3) + real(r8) :: aalc(3,3,plev) + + + character(len=CL):: Ana_File_Template,Ana_file1,Ana_file2,Ana_Path + + + integer :: dyn_year,dyn_month,dyn_day,dyn_sec,year,month,day,sec + integer :: dyn_step,ymd1,ymd2,curr_sec,next_sec,curr_year,curr_month,curr_day,curr_ymd + + integer :: analysis_step + integer :: ana_year1, ana_month1, ana_day1, ana_sec1 + integer :: ana_year2, ana_month2, ana_day2, ana_sec2 + + logical :: l_Read_next_Ana, Alarm_Read_ana, Alarm_Bump_ana, initialize + + character(len=19) :: subname='get_ana_dynfrc_fv: ' + + write(iulog,*) subname//" version 07 of get_ana_dynfrc_4scam ... " + + + Alarm_Read_Ana = .FALSE. + Alarm_Bump_Ana = .FALSE. + + if ( scmlon < 0 ) then + scmlonx = scmlon + 360._r8 + else + scmlonx = scmlon + end if + + ! Default to 6 hour steps between ana + analysis_step = 6 * 3600 + + + Ana_path = trim(scm_ana_frc_path) + Ana_File_Template = trim(Ana_path)//trim(scm_ana_frc_file_template) + + + call get_curr_date(Year,Month,Day,Sec) + + curr_ymd = (Year*10000) + (Month*100) + Day + curr_sec = Sec + + ana_sec1 = ( Sec / analysis_step ) * analysis_step + ana_day1 = Day + ana_month1 = Month + ana_year1 = Year + + YMD1=(Ana_Year1*10000) + (Ana_Month1*100) + Ana_Day1 + + + call timemgr_time_inc(YMD1,Ana_Sec1, & + YMD2,Ana_Sec2,Analysis_Step,0,0) + + Ana_Year2 = YMD2 / 10000 + Ana_Month2 = (YMD2 - Ana_Year2*10000)/100 + Ana_Day2 = YMD2 - Ana_Year2*10000 - Ana_Month2*100 + + Ana_File1 = interpret_filename_spec(Ana_File_Template , & + yr_spec=Ana_Year1 , & + mon_spec=Ana_Month1, & + day_spec=Ana_Day1 , & + sec_spec=Ana_Sec1 ) + + Ana_File2 = interpret_filename_spec(Ana_File_Template , & + yr_spec=Ana_Year2 , & + mon_spec=Ana_Month2, & + day_spec=Ana_Day2 , & + sec_spec=Ana_Sec2 ) + + + l_Read_next_Ana = .FALSE. + ! On first time step, read in 2 analysis files + if (is_first_step().and.masterproc) then + write(iulog,*) subname//" It's now (First time step):" , curr_YMD, curr_sec + write(iulog,*) "Read Initial ana files " + write(iulog,*) Ana_file1 + write(iulog,*) Ana_file2 + Alarm_Read_Ana = .TRUE. + Alarm_Bump_Ana = .FALSE. + else + ! On subsequent steps test to see if "Curr" date is later or same as "Read". + ! If it is, then l_read_next_ana=.TRUE. + call timemgr_time_ge(Read_ymd2, Read_Sec2, curr_YMD, curr_Sec, l_Read_next_ana ) + endif + + if (l_Read_next_Ana) then + Alarm_Read_Ana = .TRUE. + Alarm_Bump_Ana = .TRUE. + endif + + ! Aloocate space for analysis fields. + ! Read in both Initial Analysis files. Nothing to bump yet + if ( (Alarm_Read_Ana ) .AND. .NOT.(Alarm_Bump_Ana) ) then + initialize=.TRUE. + call read_netcdf_ana_fv_ini ( Ana_File1, nlon_ana, nlat_ana, nlev_ana ,iax, jax ) + + if ( plev /= nlev_ana) then + call endrun ("SCAM plev NE nlev_ana") + end if + + ! Full global fields + allocate( lat_ana(nlat_ana) , lon_ana(nlon_ana), lev_ana(nlev_ana) ) + allocate( U_1(nlon_ana, nlat_ana, nlev_ana), V_1(nlon_ana, nlat_ana, nlev_ana), T_1(nlon_ana, nlat_ana, nlev_ana), & + Q_1(nlon_ana, nlat_ana, nlev_ana), PS_1 (nlon_ana, nlat_ana ), PHIS_1 (nlon_ana, nlat_ana ) ) + allocate( U_2(nlon_ana, nlat_ana, nlev_ana), V_2(nlon_ana, nlat_ana, nlev_ana), T_2(nlon_ana, nlat_ana, nlev_ana), & + Q_2(nlon_ana, nlat_ana, nlev_ana), PS_2 (nlon_ana, nlat_ana ), PHIS_2 (nlon_ana, nlat_ana ) ) + + allocate( UTCORE_1(nlon_ana, nlat_ana, nlev_ana), UTCORE_2(nlon_ana, nlat_ana, nlev_ana) ) + allocate( VTCORE_1(nlon_ana, nlat_ana, nlev_ana), VTCORE_2(nlon_ana, nlat_ana, nlev_ana) ) + allocate( TTCORE_1(nlon_ana, nlat_ana, nlev_ana), TTCORE_2(nlon_ana, nlat_ana, nlev_ana) ) + allocate( OGCORE_1(nlon_ana, nlat_ana, nlev_ana), OGCORE_2(nlon_ana, nlat_ana, nlev_ana) ) + + ! SCM "patches" + nlon_alc=3 + nlat_alc=3 + nlev_alc=nlev_ana + + + + ! Patches of full global fields + allocate( Uo_1(nlon_alc, nlat_alc, nlev_alc), Vo_1(nlon_alc, nlat_alc, nlev_alc), To_1(nlon_alc, nlat_alc, nlev_alc), & + Qo_1(nlon_alc, nlat_alc, nlev_alc), PSo_1 (nlon_alc, nlat_alc ), PHISo_1 (nlon_alc, nlat_alc ) ) + allocate( Uo_2(nlon_alc, nlat_alc, nlev_alc), Vo_2(nlon_alc, nlat_alc, nlev_alc), To_2(nlon_alc, nlat_alc, nlev_alc), & + Qo_2(nlon_alc, nlat_alc, nlev_alc), PSo_2 (nlon_alc, nlat_alc ), PHISo_2 (nlon_alc, nlat_alc ) ) + + allocate( UTCOREo_1(nlon_alc, nlat_alc, nlev_alc), UTCOREo_2(nlon_alc, nlat_alc, nlev_alc), UTCOREo_X(nlon_alc, nlat_alc, nlev_alc) ) + allocate( VTCOREo_1(nlon_alc, nlat_alc, nlev_alc), VTCOREo_2(nlon_alc, nlat_alc, nlev_alc), VTCOREo_X(nlon_alc, nlat_alc, nlev_alc) ) + allocate( TTCOREo_1(nlon_alc, nlat_alc, nlev_alc), TTCOREo_2(nlon_alc, nlat_alc, nlev_alc), TTCOREo_X(nlon_alc, nlat_alc, nlev_alc) ) + allocate( OGCOREo_1(nlon_alc, nlat_alc, nlev_alc), OGCOREo_2(nlon_alc, nlat_alc, nlev_alc), OGCOREo_X(nlon_alc, nlat_alc, nlev_alc) ) + + allocate( Uo_X(nlon_alc, nlat_alc, nlev_alc), Vo_X(nlon_alc, nlat_alc, nlev_alc), To_X(nlon_alc, nlat_alc, nlev_alc), & + Qo_X(nlon_alc, nlat_alc, nlev_alc), PSo_X (nlon_alc, nlat_alc ), PHISo_X (nlon_alc, nlat_alc ) ) + allocate( ETAD_X(nlon_alc,nlat_alc,nlev_alc) ) + allocate( OMG_X(nlon_alc,nlat_alc,nlev_alc) ) + allocate( ple_X(nlon_alc, nlat_alc, nlev_alc+1), plo_X(nlon_alc, nlat_alc, nlev_alc), phi_X(nlon_alc, nlat_alc, nlev_alc+1) ) + allocate( pke_X(nlon_alc, nlat_alc, nlev_alc+1), pko_X(nlon_alc, nlat_alc, nlev_alc), phik_X(nlon_alc, nlat_alc, nlev_alc+1) ) + allocate( THv_X(nlon_alc, nlat_alc, nlev_alc ) ) + allocate( zeta_X(nlev_alc) ) + allocate( KEh_X(nlon_alc, nlat_alc, nlev_alc ) ) + allocate( Tv_X(nlon_alc, nlat_alc, nlev_alc ) ) + + call read_netcdf_ana_fv ( Ana_File1, nlon_ana, nlat_ana, nlev_ana, & + U_1, V_1, & + T_1, Q_1, PS_1, PHIS_1, & + lon_ana, lat_ana, lev_ana & + , utcore_1, vtcore_1, ttcore_1, ogcore_1 & + ) + + call read_netcdf_ana_fv ( Ana_File2, nlon_ana, nlat_ana, nlev_ana, & + U_2, V_2, & + T_2, Q_2, PS_2, PHIS_2, & + lon_ana, lat_ana, lev_ana & + , utcore_2, vtcore_2, ttcore_2, ogcore_2 & + ) + + ! Make patches + Uo_1 = U_1(iax-1:iax+1,jax-1:jax+1,:) + Vo_1 = V_1(iax-1:iax+1,jax-1:jax+1,:) + To_1 = T_1(iax-1:iax+1,jax-1:jax+1,:) + Qo_1 = Q_1(iax-1:iax+1,jax-1:jax+1,:) + PSo_1 = PS_1(iax-1:iax+1,jax-1:jax+1 ) + PHISo_1 = PHIS_1(iax-1:iax+1,jax-1:jax+1 ) + UTCOREo_1 = UTCORE_1(iax-1:iax+1,jax-1:jax+1,:) + VTCOREo_1 = VTCORE_1(iax-1:iax+1,jax-1:jax+1,:) + TTCOREo_1 = TTCORE_1(iax-1:iax+1,jax-1:jax+1,:) + OGCOREo_1 = OGCORE_1(iax-1:iax+1,jax-1:jax+1,:) + + Uo_2 = U_2(iax-1:iax+1,jax-1:jax+1,:) + Vo_2 = V_2(iax-1:iax+1,jax-1:jax+1,:) + To_2 = T_2(iax-1:iax+1,jax-1:jax+1,:) + Qo_2 = Q_2(iax-1:iax+1,jax-1:jax+1,:) + PSo_2 = PS_2(iax-1:iax+1,jax-1:jax+1 ) + PHISo_2 = PHIS_2(iax-1:iax+1,jax-1:jax+1 ) + UTCOREo_2 = UTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + VTCOREo_2 = VTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + TTCOREo_2 = TTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + OGCOREo_2 = OGCORE_2(iax-1:iax+1,jax-1:jax+1,:) + + + ! Mark Ana date as read + Read_year2 = Ana_year2 + Read_month2 = Ana_month2 + Read_day2 = Ana_day2 + Read_sec2 = Ana_sec2 + Read_YMD2 =(Ana_Year2*10000) + (Ana_Month2*100) + Ana_Day2 + + end if + + ! Bump second analysis to first postion, and read in next analysis + if ( (Alarm_Read_Ana ) .AND. (Alarm_Bump_Ana) ) then + + Uo_1 = Uo_2 + Vo_1 = Vo_2 + To_1 = To_2 + Qo_1 = Qo_2 + PSo_1 = PSo_2 + PHISo_1 = PHISo_2 + UTCOREo_1 = UTCOREo_2 + VTCOREo_1 = VTCOREo_2 + TTCOREo_1 = TTCOREo_2 + OGCOREo_1 = OGCOREo_2 + + call read_netcdf_ana_fv ( Ana_File2, nlon_ana, nlat_ana, nlev_ana, & + U_2, V_2, & + T_2, Q_2, PS_2, PHIS_2, & + lon_ana, lat_ana, lev_ana & + , utcore_2, vtcore_2, ttcore_2, ogcore_2 & + ) + + + ! Make patches + Uo_2 = U_2(iax-1:iax+1,jax-1:jax+1,:) + Vo_2 = V_2(iax-1:iax+1,jax-1:jax+1,:) + To_2 = T_2(iax-1:iax+1,jax-1:jax+1,:) + Qo_2 = Q_2(iax-1:iax+1,jax-1:jax+1,:) + PSo_2 = PS_2(iax-1:iax+1,jax-1:jax+1 ) + PHISo_2 = PHIS_2(iax-1:iax+1,jax-1:jax+1 ) + UTCOREo_2 = UTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + VTCOREo_2 = VTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + TTCOREo_2 = TTCORE_2(iax-1:iax+1,jax-1:jax+1,:) + OGCOREo_2 = OGCORE_2(iax-1:iax+1,jax-1:jax+1,:) + + + ! Mark Ana date as read + Read_year2 = Ana_year2 + Read_month2 = Ana_month2 + Read_day2 = Ana_day2 + Read_sec2 = Ana_sec2 + Read_YMD2=(Ana_Year2*10000) + (Ana_Month2*100) + Ana_Day2 + end if + + Alarm_Read_Ana = .FALSE. + Alarm_Bump_Ana = .FALSE. + + + + +#if 1 + call dynfrc_timewgts( & + (/ Ana_Year1, Ana_Month1, Ana_day1, Ana_sec1 /) , & + (/ Ana_Year2, Ana_Month2, Ana_day2, Ana_sec2 /) , & + ana_wgt1 , ana_wgt2 ) +#else + ana_wgt1 = 0._r8 ! 0=all weight on t+1 + ana_wgt2 = 1._r8 - ana_wgt1 +#endif + ! This seems correct if time-tendencies are averaged over interval + tend_wgt1 = 0._r8 + tend_wgt2 = 1._r8 - tend_wgt1 + + + if (masterproc) write(iulog,*) subname//" Ana forcing time wgts ",ana_wgt1,ana_wgt2 + + iac=2 + jac=2 + + + + Uo_X = ana_wgt1 * Uo_1 + ana_wgt2 * Uo_2 + Vo_X = ana_wgt1 * Vo_1 + ana_wgt2 * Vo_2 + To_X = ana_wgt1 * To_1 + ana_wgt2 * To_2 + Qo_X = ana_wgt1 * Qo_1 + ana_wgt2 * Qo_2 + PSo_X = ana_wgt1 * PSo_1 + ana_wgt2 * PSo_2 + PHISo_X = ana_wgt1 * PHISo_1 + ana_wgt2 * PHISo_2 + OGCOREo_X = ana_wgt1 * OGCOREo_1 + ana_wgt2 * OGCOREo_2 + + UTCOREo_X = tend_wgt1 * UTCOREo_1 + tend_wgt2 * UTCOREo_2 + VTCOREo_X = tend_wgt1 * VTCOREo_1 + tend_wgt2 * VTCOREo_2 + TTCOREo_X = tend_wgt1 * TTCOREo_1 + tend_wgt2 * TTCOREo_2 + + lon_alc = lon_ana(iax-1:iax+1) + lat_alc = lat_ana(jax-1:jax+1) + + if(masterproc) write(iulog,*) subname//" SCM lon lat: ",scmlonx,scmlat + if(masterproc) write(iulog,*) subname//" Closest Ana lon lat: ",lon_ana( iax ) , lat_ana( jax ) + + + ! Save off analysis fields for diagnostics and + ! other purposes + T_ana_diag(:) = To_X( iac, jac, :) + Q_ana_diag(:) = Qo_X( iac, jac, :) + U_ana_diag(:) = Uo_X( iac, jac, :) + V_ana_diag(:) = Vo_X( iac, jac, :) + + !================================================ + ! Patch in SCM profiles here if wanted. + ! This acts as "dynamical nudging", since + ! horizontal advective tendencies will become + ! stronger if SCM state drifts away from re-ana. + ! Note, this will only be effective w/ upwind + ! scheme, since 2nd order cntrd skips over central + ! point in stencil. + !---- + ! For stability it turns out may be good to scale + ! with pressure so that high-velocity strato winds + ! don't lead to CFL violations. So, as a bad, dirty, + ! dirty short term solution, weight "reaction" by + ! pref_mid. Clearly, better soln would be to + ! sub-step this part of the dynamics as is done + ! for the other "dycores". + !================================================= + ! Calculate "reaction coefficient" + !--------------------------------- + alpha_react(:)=1.0_r8 !1._r8 + + ! Adjust central profiles in stencils + !------------------------------------ + if (scm_ana_t_react) then + To_X( iac, jac, :) = alpha_react(:) * T_scm(:) & + + ( 1._r8-alpha_react(:) ) * To_X( iac, jac, :) + if(masterproc) write(iulog,*) subname//" REACTING to SCM T-state ..... " + else + if(masterproc) write(iulog,*) subname//" No reaction to SCM T-state ..... " + endif + if (scm_ana_q_react) then + Qo_X( iac, jac, :) = alpha_react(:) * Q_scm(:) & + + ( 1._r8-alpha_react(:) ) * Qo_X( iac, jac, :) + if(masterproc) write(iulog,*) subname//" REACTING to SCM Q-state ..... " + else + if(masterproc) write(iulog,*) subname//" No reaction to SCM Q-state ..... " + endif + if (scm_ana_u_react) then + Uo_X( iac, jac, :) = alpha_react(:) * U_scm(:) & + + ( 1._r8-alpha_react(:) ) * Uo_X( iac, jac, :) + if(masterproc) write(iulog,*) subname//" REACTING to SCM U-state ..... " + else + if(masterproc) write(iulog,*) subname//" No reaction to SCM U-state ..... " + endif + if (scm_ana_v_react) then + Vo_X( iac, jac, :) = alpha_react(:) * V_scm(:) & + + ( 1._r8-alpha_react(:) ) * Vo_X( iac, jac, :) + if(masterproc) write(iulog,*) subname//" REACTING to SCM V-state ..... " + else + if(masterproc) write(iulog,*) subname//" No reaction to SCM V-state ..... " + endif + + + + !========================================= + + call virtual_t( nlon_alc,nlat_alc,nlev_alc, & + To_X , Qo_X , Tv_X ) + + call makepr_fv( nlon_alc,nlat_alc,nlev_alc, & + tv_X , pso_X , phiso_X , & + plo_X, ple_X, phi_X ) + call etadot_fv ( nlon_alc , nlat_alc , nlev_alc , lon_alc , lat_alc , & + uo_X , & + vo_X , & + plo_X, ple_X , etad_X , omg_X ) + call zeta_fv( nlon_alc,nlat_alc,nlev_alc, & + lon_alc ,lat_alc , & + uo_X , vo_X , zeta_X ) + + call makepk_fv( nlon_alc,nlat_alc,nlev_alc, & + To_X , Qo_X , & + pso_X , phiso_X , & + pko_X, pke_X, phik_X, thv_X ) + + KEh_X = 0.5 * ( Uo_X**2 + Vo_X**2 ) + + + if (scm_ana_x_plevels) then + call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, & + iac, jac, uo_X , plo_X ) + call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, & + iac, jac, vo_X , plo_X ) + call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, & + iac, jac, to_X , plo_X ) + call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, & + iac, jac, qo_X , plo_X ) + call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc, & + iac, jac, tv_X , plo_X ) + !Retain p-frc calculation on eta??? + !call patch_eta_x_plv ( nlon_alc , nlat_alc , nlev_alc+1, & + ! iac, jac, phi_X , ple_X ) + if(masterproc) write(iulog,*) subname//" calcs on PRESSURE levels " + else + if(masterproc) write(iulog,*) subname//" calcs on ETA levels " + end if + + + zeta_ana = zeta_X + omega_recalc_ana = omg_X( iac,jac,:) + etad_ana = etad_X( iac,jac,:) + plo_ana = plo_X( iac,jac,:) + t_ana = To_X( iac,jac,:) + tv_ana = Tv_X( iac,jac,:) + q_ana = Qo_X( iac,jac,:) + ps_ana = PSo_X( iac,jac ) + phis_ana = PHISo_X( iac,jac ) + + u_ana = Uo_X( iac,jac,:) + v_ana = Vo_X( iac,jac,:) + + rho_ana = plo_ana / ( Rdair * tv_ana ) + + uten_dycore_ana = UTCOREo_X( iac,jac,:) + vten_dycore_ana = VTCOREo_X( iac,jac,:) + tten_dycore_ana = TTCOREo_X( iac,jac,:) + omega_dycore_ana = OGCOREo_X( iac,jac,:) + + + ! Horz. gradient calcs + + kehg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, KEh_X ) + + ! T_x, T_y should be straight T (not virtual) + teg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Tv_X ) !test 05-31-21 + + qg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Qo_X ) + + ug_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Uo_X ) + + vg_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, Vo_X ) + + aalc = 0.5*( PHI_X( :, :, 2:nlev_alc+1) + PHI_X(: , : ,1:nlev_alc) ) + phig_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, aalc ) + + plog_X = grad_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, plo_X(:,:,1:nlev_alc) ) + + + lin_pfc_X = lin_pfc_fv( nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, ple_X, phi_X ) + + kehg_ana = kehg_X + plog_ana = plog_X + phig_ana = phig_X + teg_ana = teg_X + qg_ana = qg_X + ug_ana = ug_X + vg_ana = vg_X + lin_pfc_ana = lin_pfc_X + + !put together pieces for u*grad(u) form of U and V adv tendencies + + if ( scm_ana_upwind .OR. scm_ana_u_react ) then + uten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Uo_X ) + else + uten_hadv_ana = -u_ana * ug_ana(:,1) - v_ana * ug_ana(:,2) + end if + if ( scm_ana_upwind .OR. scm_ana_v_react ) then + vten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Vo_X ) + else + vten_hadv_ana = -u_ana * vg_ana(:,1) - v_ana * vg_ana(:,2) + end if + + ! Coriolis terms + !====================================== + absvo = 2._r8 * OOmega * sin( lat_ana(jax) * PI/180._r8 ) + !Allow Coriolis to react to SCM winds + uten_vort_ana = absvo * v_ana + vten_vort_ana = -absvo * u_ana + ! ----- Diags for VI form (0-out) + uten_keg_ana = 0._r8 ! fill with 0 + + if (.FALSE.) then ! No horz. p-gradient in p-coords + uten_pfrc_ana = - phig_ana(:,1) + vten_pfrc_ana = - phig_ana(:,2) + else + !put together pieces for Pressure and Phi gradient tencency terms + uten_pfrc_ana = -(1._r8/rho_ana) * plog_ana(:,1) - phig_ana(:,1) + vten_pfrc_ana = -(1._r8/rho_ana) * plog_ana(:,2) - phig_ana(:,2) + end if + + + if ( scm_ana_upwind .OR. scm_ana_t_react ) then + tten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Tv_X ) + else + tten_hadv_ana = -u_ana * teg_ana(:,1) - v_ana * teg_ana(:,2) ! should be straight T (not virtual) + end if + if ( scm_ana_upwind .OR. scm_ana_q_react ) then + qten_hadv_ana = upwind_hadv(nlon_alc,nlat_alc,nlev_alc,iac,jac,lon_alc,lat_alc, u_ana, v_ana, Qo_X ) + else + qten_hadv_ana = -u_ana * qg_ana(:,1) - v_ana * qg_ana(:,2) + end if + + if (.not.(scm_ana_direct_omega)) then + omega_ana = omega_recalc_ana ! use reconstructed omega + if(masterproc) write(iulog,*) subname//" Omega recalc from ana U,V etc." + else + omega_ana = omega_dycore_ana ! use direct omega from dycore/ana + if(masterproc) write(iulog,*) subname//" Omega direct from ana" + end if + + + if (.not.(scm_ana_x_plevels)) then + !Tendencies due to vertical advection (etadot * D_eta ... ) + uten_vadv_ana = vadv_fv( nlev_alc, etad_ana, u_ana ) + vten_vadv_ana = vadv_fv( nlev_alc, etad_ana, v_ana ) + tten_vadv_ana = vadv_fv( nlev_alc, etad_ana, tv_ana ) ! should be straight T (not virtual) + qten_vadv_ana = vadv_fv( nlev_alc, etad_ana, q_ana ) + else + !Tendencies due to vertical advection (Omega * D_p ... ) + uten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, u_ana ) + vten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, v_ana ) + tten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, t_ana ) ! should be straight T (not virtual) + qten_vadv_ana = vadv_fv_press( nlev_alc, omega_ana, plo_ana, q_ana ) + end if + + tten_comp_ana = (1./cpair)*( omega_ana / rho_ana ) + + !DIags for pressure/geop grad forces + uten_phig_ana = - phig_ana(:,1) + uten_prg_ana = - (1._r8/rho_ana) * plog_ana(:,1) + + end subroutine get_ana_dynfrc_fv + +!----------------------------------------------------- +! Stuff ... useful ojala +!----------------------------------------------------- + !------------------------- + function vadv_fv( nlev, etad, aa ) result( tend ) + use hycoef, only: hyai, hybi, ps0, hyam, hybm + integer, intent(in) :: nlev + real(r8), intent(in) :: etad(nlev) , aa(nlev) + real(r8) :: tend(nlev) + real(r8) :: eta(nlev) + integer :: L + + eta = hybm+hyam + + do L=2,nlev-1 + tend(L) = etad(L)* ( aa(L+1) - aa(L-1) ) / ( eta(L+1) - eta(L-1) ) + end do + L=1 + tend(L) = etad(L)* ( aa(L+1) - aa(L) ) / ( eta(L+1) - eta(L) ) + L=nlev + tend(L) = etad(L)* ( aa(L) - aa(L-1) ) / ( eta(L) - eta(L-1) ) + + tend = -1._r8*tend ! for RHS consistency + + end function vadv_fv +!--------------------------- + !------------------------- + function vadv_fv_press( nlev, omega, plo, aa ) result( tend ) + integer, intent(in) :: nlev + real(r8), intent(in) :: omega(nlev) , aa(nlev),plo(nlev) + real(r8) :: tend(nlev) + integer :: L + + do L=2,nlev-1 + tend(L) = omega(L)* ( aa(L+1) - aa(L-1) ) / ( plo(L+1) - plo(L-1) ) + end do + L=1 + tend(L) = omega(L)* ( aa(L+1) - aa(L) ) / ( plo(L+1) - plo(L) ) + L=nlev + tend(L) = omega(L)* ( aa(L) - aa(L-1) ) / ( plo(L) - plo(L-1) ) + + tend = -1._r8*tend ! for RHS consistency + + end function vadv_fv_press +!--------------------------- + function lin_pfc_fv( nlon,nlat,nlev,iax,jax,lons,lats, pre, phi ) result( pfc ) + + integer, intent(in) :: nlon,nlat,nlev,iax,jax + real(r8), intent(in) :: pre(nlon,nlat,nlev+1),phi(nlon,nlat,nlev+1) + real(r8), intent(in) :: lats(nlat),lons(nlon) + real(r8) :: pfc(nlev,2) + real(r8) :: pfxW(nlev) , pfxE(nlev) + real(r8) :: pfyS(nlev) , pfyN(nlev) + real(r8) :: rlats(nlat),rlons(nlon),dx,dy,ds + real(r8) :: pr1,pr2,pr3,pr4, ph1,ph2,ph3,ph4 + integer :: L , igg + + ! Begin + rlons(:) = lons(:) * PI/180._r8 + rlats(:) = lats(:) * PI/180._r8 + + dx=( rlons(2)-rlons(1) ) * Rearth + dy=( rlats(2)-rlats(1) ) * Rearth + + ds = MAX( dx*cos(rlats(jax)) , .1_r8 ) + igg = iax + do L=1,nlev + pr1 = pre(igg-1,jax,L+1) + pr2 = pre(igg ,jax,L+1) + pr3 = pre(igg ,jax,L ) + pr4 = pre(igg-1,jax,L ) + ph1 = phi(igg-1,jax,L+1) + ph2 = phi(igg ,jax,L+1) + ph3 = phi(igg ,jax,L ) + ph4 = phi(igg-1,jax,L ) + pfxW(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) ) + end do + igg = iax +1 + do L=1,nlev + pr1 = pre(igg-1,jax,L+1) + pr2 = pre(igg ,jax,L+1) + pr3 = pre(igg ,jax,L ) + pr4 = pre(igg-1,jax,L ) + ph1 = phi(igg-1,jax,L+1) + ph2 = phi(igg ,jax,L+1) + ph3 = phi(igg ,jax,L ) + ph4 = phi(igg-1,jax,L ) + pfxE(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) ) + end do + ds = dy + igg = jax + do L=1,nlev + pr1 = pre(iax,igg-1,L+1) + pr2 = pre(iax,igg ,L+1) + pr3 = pre(iax,igg ,L ) + pr4 = pre(iax,igg-1,L ) + ph1 = phi(iax,igg-1,L+1) + ph2 = phi(iax,igg ,L+1) + ph3 = phi(iax,igg ,L ) + ph4 = phi(iax,igg-1,L ) + pfyS(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) ) + end do + igg = jax +1 + do L=1,nlev + pr1 = pre(iax,igg-1,L+1) + pr2 = pre(iax,igg ,L+1) + pr3 = pre(iax,igg ,L ) + pr4 = pre(iax,igg-1,L ) + ph1 = phi(iax,igg-1,L+1) + ph2 = phi(iax,igg ,L+1) + ph3 = phi(iax,igg ,L ) + ph4 = phi(iax,igg-1,L ) + pfyN(L) = ( (pr2-pr4)*(ph1-ph3) + (pr1-pr3)*(ph4-ph2) ) /( ds * ( (pr2-pr4) + (pr1-pr3) ) ) + end do + + + do L=1,nlev + pfc(L,1) = 0.5_r8*( pfxW(L) + pfxE(L) ) + pfc(L,2) = 0.5_r8*( pfyS(L) + pfyN(L) ) + end do + + + + end function lin_pfc_fv + !------------------------- + function grad_fv( nlon,nlat,nlev,iax,jax,lons,lats, aa ) result( ga ) + + integer, intent(in) :: nlon,nlat,nlev,iax,jax + real(r8), intent(in) :: aa(nlon,nlat,nlev) + real(r8), intent(in) :: lats(nlat),lons(nlon) + real(r8) :: ga(nlev,2) + real(r8) :: rlats(nlat),rlons(nlon),dx,dy + integer :: L + + ! Begin + rlons(:) = lons(:) * PI/180._r8 + rlats(:) = lats(:) * PI/180._r8 + + dx=( rlons(2)-rlons(1) ) * Rearth + dy=( rlats(2)-rlats(1) ) * Rearth + + do L=1,nlev + ga(L,1) = (aa(iax+1,jax,L) - aa(iax-1,jax,L))/( 2._r8*dx*cos(rlats(jax)) + 0.1_r8 ) + ga(L,2) = (aa(iax,jax+1,L) - aa(iax,jax-1,L))/( 2._r8*dy ) + end do + + + + end function grad_fv + !------------------------- + function upwind_hadv( nlon,nlat,nlev,iax,jax,lons,lats,u,v, aa ) result( hadv_tend ) + !use shr_kind_mod, only: r8 => shr_kind_r8 + !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters + ! pi => shr_const_pi , & + ! omega => shr_const_omega + + integer, intent(in) :: nlon,nlat,nlev,iax,jax + real(r8), intent(in) :: aa(nlon,nlat,nlev) + real(r8), intent(in) :: lats(nlat),lons(nlon),u(nlev),v(nlev) + real(r8) :: hadv_tend(nlev) + real(r8) :: rlats(nlat),rlons(nlon),dx,dy,xten(nlev),yten(nlev) + integer :: L + + ! Begin + rlons(:) = lons(:) * PI/180._r8 + rlats(:) = lats(:) * PI/180._r8 + + dx=( rlons(2)-rlons(1) ) * Rearth + dy=( rlats(2)-rlats(1) ) * Rearth + + do L=1,nlev + if ( u(L) >= 0._r8 ) then + xten(L) = u(L) * ( aa(iax,jax,L) - aa(iax-1,jax,L))/( dx*cos(rlats(jax)) + 0.1_r8 ) + else + xten(L) = u(L) * ( aa(iax+1,jax,L) - aa(iax,jax,L))/( dx*cos(rlats(jax)) + 0.1_r8 ) + end if + end do + do L=1,nlev + if ( v(L) >= 0._r8 ) then + yten(L) = v(L) * ( aa(iax,jax,L) - aa(iax,jax-1,L))/( dy ) + else + yten(L) = v(L) * ( aa(iax,jax+1,L) - aa(iax,jax,L))/( dy ) + end if + end do + + hadv_tend(:) = -1._r8 * ( xten(:) + yten(:) ) + + + end function upwind_hadv +!========================================= + subroutine makepk_fv( nlon,nlat,nlev, t, q, ps, phis, pko, pke, phi, th ) + use hycoef, only: hyai, hybi, ps0, hyam, hybm + integer, intent(in) :: nlon,nlat,nlev + real(r8), intent(in) :: t(nlon,nlat,nlev),q(nlon,nlat,nlev),ps(nlon,nlat),phis(nlon,nlat) + real(r8), intent(out) :: pko(nlon,nlat,nlev),th(nlon,nlat,nlev),pke(nlon,nlat,nlev+1), phi(nlon,nlat,nlev+1) + real(r8) :: ple(nlon,nlat,nlev+1),plo(nlon,nlat,nlev),rv(nlon,nlat,nlev) + real(r8) :: kappa, p00 + integer :: L + + do L=1,nlev+1 + ple(:,:,L) = hyai(L)*ps0 + hybi(L)*ps(:,:) + end do + do L=1,nlev + plo(:,:,L) = hyam(L)*ps0 + hybm(L)*ps(:,:) + end do + + kappa=rdair/cpair + + pko = plo**kappa + pke = ple**kappa + + p00 = 100000._r8 + th = ( ( p00 / plo)**kappa ) * t + + rv = 1._r8/(1._r8 - q) - 1._r8 + th = th*(1._r8 + 0.61_r8 * rv ) + + phi(:,:,nlev+1) = phis(:,:) + do L=nlev,1,-1 + phi(:,:,L) = phi(:,:,L+1) - ( CpAir * Th(:,:,L) ) * ( pke(:,:,L) - pke(:,:,L+1) ) / (p00**kappa ) + end do + + + end subroutine makepk_fv + +!============================================================================= + subroutine makepr_fv( nlon,nlat,nlev, t, ps, phis, plo, ple, phi ) + use hycoef, only: hyai, hybi, ps0, hyam, hybm + use shr_const_mod, only: rdair => shr_const_rdair + integer, intent(in) :: nlon,nlat,nlev + real(r8), intent(in) :: t(nlon,nlat,nlev),ps(nlon,nlat),phis(nlon,nlat) + real(r8), intent(out) :: plo(nlon,nlat,nlev), ple(nlon,nlat,nlev+1), phi(nlon,nlat,nlev+1) + real(r8) :: lnple(nlon,nlat,nlev+1) + integer :: L + + do L=1,nlev+1 + ple(:,:,L) = hyai(L)*ps0 + hybi(L)*ps(:,:) + end do + do L=1,nlev + plo(:,:,L) = hyam(L)*ps0 + hybm(L)*ps(:,:) + end do + + lnple = log( ple ) + phi(:,:,nlev+1) = phis(:,:) + do L=nlev,1,-1 + phi(:,:,L) = phi(:,:,L+1) - (RdAir * T(:,:,L) ) * ( lnple(:,:,L) - lnple(:,:,L+1) ) + !phi(:,:,L) = phi(:,:,L+1) - (RdAir * T(:,:,L) / plo(:,:,L) ) * ( ple(:,:,L) - ple(:,:,L+1) ) + end do + + end subroutine makepr_fv + +!============================================================================= + subroutine virtual_t( nlon,nlat,nlev, t, q, tv ) + use hycoef, only: hyai, hybi, ps0, hyam, hybm + use shr_const_mod, only: rdair => shr_const_rdair + integer, intent(in) :: nlon,nlat,nlev + real(r8), intent(in) :: t(nlon,nlat,nlev),q(nlon,nlat,nlev) + real(r8), intent(out) :: tv(nlon,nlat,nlev) + real(r8) :: rv(nlon,nlat,nlev) + integer :: L + + + rv = 1._r8/(1._r8 - q) - 1._r8 + tv = t*(1._r8 + 0.61_r8 * rv ) + + + end subroutine virtual_t + + !------------------------- + subroutine zeta_fv( nlon,nlat,nlev,lons,lats, u,v, zeta ) + !use shr_kind_mod, only: r8 => shr_kind_r8 + !use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters + ! pi => shr_const_pi , & + ! omega => shr_const_omega + + integer, intent(in) :: nlon,nlat,nlev + real(r8), intent(in) :: u(nlon,nlat,nlev),v(nlon,nlat,nlev) + real(r8), intent(out) :: zeta(nlev) + real(r8), intent(in) :: lats(nlat),lons(nlon) + real(r8) :: rlats(nlat),rlons(nlon) + real(r8) :: dy,dx0,dx,darea,voo,voo2 + + integer :: iap,jap,iam,jam,i,j,L,iax,jax + + iax=2 + jax=2 + + rlons(:) = lons(:) * PI/180._r8 + rlats(:) = lats(:) * PI/180._r8 + + dx0 = rearth* ( rlons(2)-rlons(1) ) + dy = rearth* ( rlats(2)-rlats(1) ) + + darea = dy*dx0*cos( rlats(jax) ) + + + do L =1,nlev + zeta(L) = & + ( V(iax+1,jax, L) - V(iax-1,jax,L) ) / ( 2._r8*dx0*cos( rlats(jax) ) ) & + - ( U(iax,jax+1, L) - U(iax,jax-1,L) ) / ( 2._r8*dy ) + end do + + + end subroutine zeta_fv +!================================================================ + subroutine etadot_fv ( nlon, nlat, nlev, lons, lats, u, v, plo, ple, etadot , omega ) + use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8 + use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters + pi => shr_const_pi + use hycoef, only: hyai, hybi, ps0, hyam, hybm + + integer, intent(in) :: nlon,nlat,nlev + real(r8), intent(in) :: lons(nlon),lats(nlat) + real(r8), intent(in) :: u(nlon,nlat,nlev) , v(nlon,nlat,nlev) , plo( nlon,nlat,nlev) , ple( nlon,nlat,nlev+1) + real(r8), intent(out) :: etadot( nlon,nlat,nlev) ,omega(nlon,nlat,nlev) + !real(r8), intent(in) :: uc(:,:,:) , vc(:,:,:) , ple(:,:,:) + + ! Local variables + real(r8),allocatable :: div(:,:,:) + real(r8),allocatable :: mass(:,:,:), fuc(:,:,:),fvc(:,:,:) + real(r8) :: rlats(nlat), rlons(nlon), rcos1, eta(nlev+1) , dx,dy! radians + real(r8), allocatable :: etadot_t1(:,:), etadot_t2(:,:,:) + integer :: i,j,L,im1,jm1,ip1,jp1 + real :: uc_ijL , vc_ijL + + allocate ( div(nlon,nlat,nlev) ) + allocate ( mass(nlon,nlat,nlev), fuc(nlon,nlat,nlev),fvc(nlon,nlat+1,nlev) ) + allocate ( etadot_t1(nlon,nlat), etadot_t2(nlon,nlat,nlev) ) + + div = 0._r8 + fuc = 0._r8 + fvc = 0._r8 + mass = 0._r8 + etadot = 0._r8 + etadot_t1 = 0._r8 + etadot_t2 = 0._r8 + + rlons(:) = lons(:) * PI/180._r8 + rlats(:) = lats(:) * PI/180._r8 + + do L=1,nlev+1 + eta(L) = hyai(L) + hybi(L) ! 1._r8*L/(nlev+1) + end do + do L=1,nlev + mass(:,:,L) = ( ple(:,:,L+1)-ple(:,:,L) )/( eta(L+1)-eta(L) ) + end do + + ! calculate mass fluxes at gridbox edges, using upwind algorithm + do L=1,nlev + do j=1,nlat + do i=2,nlon + im1=i-1 + !if ( i == 1) im1=nlon + uc_ijL = 0.5_r8*( u(im1,j,L) + u(i,j,L) ) + if ( uc_ijL < 0._r8 ) fuc(i,j,L)= uc_ijL * mass(i,j,L) + if ( uc_ijL >= 0._r8 ) fuc(i,j,L)= uc_ijL * mass(im1,j,L) + end do + end do + end do + ! Note: cos(lat) term incorporated into fluxes + do L=1,nlev + do j=2,nlat + do i=1,nlon + jm1=j-1 + vc_ijL = 0.5_r8 * ( v(i,jm1,L)+v(i,j,L) ) + if ( vc_ijL < 0._r8 ) fvc(i,j,L)= vc_ijL * mass(i,j,L) *cos( rlats(j) ) + if ( vc_ijL >= 0._r8 ) fvc(i,j,L)= vc_ijL * mass(i,jm1,L) *cos( rlats(jm1) ) + end do + end do + end do + + + ! now calculate HORZ divergence of (FUC,FVC). Note coslat term already + ! incorporated in FVC. + do L=1,nlev + do j=1,nlat-1 + do i=1,nlon-1 + ip1=i+1 + jp1=j+1 + rcos1 = 1._r8 /( Rearth*cos( rlats(j) ) ) + div(i,j,L) = rcos1 * ( FUC(ip1,j,L)-FUC(i,j,L) ) / (rlons(ip1)-rlons(i) ) & + + rcos1 * ( FVC(i,jp1,L)-FVC(i,j,L) ) / (rlats(jp1)-rlats(j) ) + end do + end do + end do + + + etadot_t1(:,:)=0._r8 + etadot_t2(:,:,:)=0._r8 + do L=1,nlev + etadot_t1(:,:) = etadot_t1(:,:) + div(:,:,L)*(eta(L+1)-eta(L)) + end do + do L=2,nlev + etadot_t2(:,:,L) = etadot_t2(:,:,L-1) + div(:,:,L)*(eta(L+1)-eta(L)) + end do + do L=1,nlev + etadot(:,:,L) = ( hybm(L)*etadot_t1(:,:) - etadot_t2(:,:,L) ) / mass(:,:,L) + end do + + dx=( rlons(2)-rlons(1) ) * Rearth + dy=( rlats(2)-rlats(1) ) * Rearth + omega = 0._r8 + + do L=1,nlev + do j=2,nlat-1 + do i=2,nlon-1 + omega(i,j,L) = u(i,j,L) * (plo(i+1,j,L)-plo(i-1,j,L))/( 2._r8*dx*cos(rlats(j)) + 0.1_r8 ) & + + v(i,j,L) * (plo(i,j+1,L)-plo(i,j-1,L))/( 2._r8*dy ) & + - etadot_t2(i,j,L) + end do + end do + end do + + +end subroutine etadot_fv + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Reading netcdf files +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!================================================================ + subroutine read_netcdf_ana_fv_ini( anal_file , nlon, nlat, nlev,lonidx,latidx ) + ! + ! READ_NETCDF_ANAL_INI: + ! Open the given analyses data file. Query dimesnisons. + ! Close. + !=============================================================== + use cam_abortutils, only : endrun + use netcdf + use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8 + use hycoef, only: hyai, hybi, ps0, hyam, hybm + use shr_const_mod, only: rdair => shr_const_rdair + use scammod, only: scmlon,scmlat + use shr_scam_mod, only: shr_scam_getCloseLatLon ! Standardized system subroutines + + !------------- + character(len=*),intent(in):: anal_file + + integer, intent(out) :: nlon,nlat,nlev,latidx,lonidx + + ! Local values + !------------- + integer :: ncid,varid,istat + integer :: ilat,ilon,ilev + integer :: i,j,L + + real(r8) :: closelon,closelat + + logical :: l_have_us , l_have_vs + + character(len=24) :: subname='read_netcdf_ana_fv_ini: ' + + l_have_us = .FALSE. + l_have_vs = .FALSE. + + ! masterporc does all of the work here + !----------------------------------------- + if(masterproc) then + + ! Open the given file + !----------------------- + istat=nf90_open(trim(anal_file),NF90_NOWRITE,ncid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//'NF90_OPEN: failed for file ',trim(anal_file) + write(iulog,*) nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + ! Read in Dimensions + !-------------------- + istat=nf90_inq_dimid(ncid,'lon',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_inquire_dimension(ncid,varid,len=nlon) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_dimid(ncid,'lat',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_inquire_dimension(ncid,varid,len=nlat) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_dimid(ncid,'lev',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_inquire_dimension(ncid,varid,len=nlev) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + call shr_scam_getCloseLatLon(ncid ,scmlat,scmlon,closelat,closelon,latidx,lonidx) + + ! Close the analyses file and exit + !----------------------- + istat=nf90_close(ncid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_EUL') + endif + + endif ! (masterproc) then + + + end subroutine read_netcdf_ana_fv_ini + +!================================================================ + subroutine read_netcdf_ana_fv( anal_file , nlon, nlat, nlev, & + u, v, & + t, q, ps, phis, & + lons, lats, levs & + , utcore, vtcore, ttcore, ogcore & + ) + ! + ! READ_NETCDF_ANAL : + ! Open the given analyses data file, read in + ! U,V,T,Q, and PS values as well as Lons, Lats. + !=============================================================== + use cam_abortutils, only : endrun + use netcdf + use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8 + use hycoef, only: hyai, hybi, ps0, hyam, hybm + use shr_const_mod, only: rdair => shr_const_rdair + ! Arguments + !------------- + character(len=*),intent(in):: anal_file + + integer, intent(in ) :: nlon,nlat,nlev + real(r8), intent(out) :: U(nlon,nlat,nlev), V(nlon,nlat,nlev) + real(r8), intent(out) :: T(nlon,nlat,nlev), Q(nlon,nlat,nlev) + real(r8), intent(out) :: PS(nlon,nlat), PHIS(nlon,nlat) + !real(r8), intent(out) :: PHI(nlon,nlat,nlev+1),PLE(nlon,nlat,nlev+1),PLO(nlon,nlat,nlev) + real(r8), intent(out) :: Lats(nlat),Lons(nlon),Levs(nlev) + + real(r8), intent(out) :: UTCORE(nlon,nlat,nlev), VTCORE(nlon,nlat,nlev), TTCORE(nlon,nlat,nlev) + real(r8), intent(out) :: OGCORE(nlon,nlat,nlev) + + ! Local values + !------------- + integer :: ncid,varid,istat + integer :: ilat,ilon,ilev + integer :: i,j,L + + logical :: l_have_us , l_have_vs + + character(len=20) :: subname='read_netcdf_ana_fv: ' + + l_have_us = .FALSE. + l_have_vs = .FALSE. + + ! masterporc does all of the work here + !----------------------------------------- + if(masterproc) then + + ! Open the given file + !----------------------- + istat=nf90_open(trim(anal_file),NF90_NOWRITE,ncid) + if(istat.ne.NF90_NOERR) then + write(iulog,*) subname//'NF90_OPEN: failed for file ',trim(anal_file) + write(iulog,*) nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + end if + end if + + + + if(masterproc) then + + istat=nf90_inq_varid(ncid,'lon',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid,Lons) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'lat',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid,Lats) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'lev',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid,Levs) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + + ! Read in, transpose lat/lev indices, + ! and scatter data arrays + !---------------------------------- + ! First block reads U + !---------------------------------- + istat=nf90_inq_varid(ncid,'U',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid, U ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'V',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid, V ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'T',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid, T ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'Q',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid, Q ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'PS',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + istat=nf90_get_var(ncid,varid,PS ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'PHIS',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_SE') + endif + istat=nf90_get_var(ncid,varid,PHIS ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + endif + + istat=nf90_inq_varid(ncid,'UTEND_CORE',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//"No UTEND_CORE on file: " + write(iulog,*) trim(anal_file) + utcore(:,:,:)=-9999._r8 + else + istat=nf90_get_var(ncid,varid,utcore ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + end if + end if + + istat=nf90_inq_varid(ncid,'VTEND_CORE',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//"No VTEND_CORE on file: " + write(iulog,*) trim(anal_file) + vtcore(:,:,:)=-9999._r8 + else + istat=nf90_get_var(ncid,varid,vtcore ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + end if + end if + + istat=nf90_inq_varid(ncid,'DTCORE',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//"No TTEND_CORE on file: " + write(iulog,*) trim(anal_file) + ttcore(:,:,:)= -9999._r8 + else + istat=nf90_get_var(ncid,varid,ttcore ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + end if + end if + + istat=nf90_inq_varid(ncid,'OMEGA',varid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//"No OMEGA (core) on file: " + write(iulog,*) trim(anal_file) + ogcore(:,:,:)=-9999._r8 + else + istat=nf90_get_var(ncid,varid,ogcore ) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_FV') + end if + end if + + + ! Close the analysis file + !----------------------- + istat=nf90_close(ncid) + if(istat /= NF90_NOERR) then + write(iulog,*) subname//nf90_strerror(istat) + call endrun ('UPDATE_ANALYSES_EUL') + endif + + end if ! (masterproc) then + !------------ + + return + end subroutine read_netcdf_ana_fv +!================================================================ +!================================================================ + subroutine dynfrc_timewgts ( & + ana_prev_date, ana_next_date, & + wgt1 , wgt2 ) + + + use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8 + use ESMF + use time_manager, only:timemgr_time_ge,timemgr_time_inc,get_curr_date,get_step_size + + integer, intent(in) :: ana_prev_date(4), ana_next_date(4) + real(r8) , intent(out) :: wgt1,wgt2 + + type(ESMF_Time) :: Date1,Date2,Date0 + type(ESMF_TimeInterval) :: DateDiff2,DateDiff0,DateDiff, AnaDiff + integer :: DeltaT0, DeltaT2 , YMD, Year,Month,Day,Sec, Ana_interval, rc + + call get_curr_date(Year,Month,Day,Sec) + YMD=(Year*10000) + (Month*100) + Day + + call ESMF_TimeSet(Date0,YY=Ana_prev_date(1), MM=Ana_prev_date(2) , & + DD= Ana_prev_date(3) , S= Ana_prev_date(4) ) + call ESMF_TimeSet(Date1,YY=Year,MM=Month,DD=Day,S=Sec) + + call ESMF_TimeSet(Date2,YY=Ana_next_date(1), MM=Ana_next_date(2) , & + DD= Ana_next_date(3) , S= Ana_next_date(4) ) + AnaDiff =Date2-Date0 + call ESMF_TimeIntervalGet(AnaDiff,S=Ana_interval ,rc=rc) + + DateDiff2 =Date2-Date1 + call ESMF_TimeIntervalGet(DateDiff2,S=DeltaT2,rc=rc) + DateDiff0 =Date1-Date0 + call ESMF_TimeIntervalGet(DateDiff0,S=DeltaT0,rc=rc) + + wgt1 = 1._r8 - ( 1._r8 * DeltaT0 ) / Ana_interval + wgt2 = 1._r8 - ( 1._r8 * DeltaT2 ) / Ana_interval + +end subroutine dynfrc_timewgts + + +!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine patch_eta_x_plv ( nx , ny, nL,ix, jx, aa, plo ) + integer, intent(in) :: nx,ny,nl,ix,jx + real(r8), intent(in) :: plo(nx,ny,nL) + real(r8), intent(inout) :: aa(nx,ny,nL) + + real(r8) :: plx(nL),plq(nL),aax(nL),aaq(nL),aat(nx,ny,nL) + real(r8) :: dp,dpk,dpk1,wtk,wtk1 + integer :: i,j,L,k + + + plx(:) = plo(ix,jx,:) ! target pressures + + do j=1,ny + do i=1,nx + plq(:) = plo(i,j,:) + aaq(:) = aa(i,j,:) + do L=1,nl + do k=2,nl + if ( ( plx(L) <= plq(k) ).AND.(plx(L) > plq(k-1) ) ) then + dp = plq(k)-plq(k-1) + dpk1 = plx(L)-plq(k-1) + dpk = plq(k)-plx(L) + wtk1 = 1._r8 - dpk1 / dp + wtk = 1._r8 - dpk / dp + aax(L) = wtk * aaq(k) + wtk1 * aaq(k-1) + end if + end do + if ( plx(L) <= plq(1) ) aax(L)=aaq(1) + if ( plx(L) > plq(NL) ) aax(L)=aaq(NL) + end do + + aat(i,j,:)=aax(:) + end do + end do + + aa=aat + + end subroutine patch_eta_x_plv + + +end module get_ana_dynfrc_4scam diff --git a/src/dynamics/eul/scmforecast.F90 b/src/dynamics/eul/scmforecast.F90 index f9c0cbc6a8..83fb528a16 100644 --- a/src/dynamics/eul/scmforecast.F90 +++ b/src/dynamics/eul/scmforecast.F90 @@ -9,7 +9,7 @@ module scmforecast use spmd_utils, only: masterproc use cam_logfile, only: iulog use cam_control_mod, only: adiabatic - + use get_ana_dynfrc_4scam, only: get_ana_dynfrc_fv implicit none private save @@ -59,10 +59,16 @@ subroutine forecast( lat , nlon , ztodt , & scm_relax_tau_sec,scm_relax_tau_top_sec,scm_relax_top_p, & scm_relaxation,scm_use_obs_qv,scm_use_obs_t,scm_use_obs_uv,scm_zadv_q,scm_zadv_t, & scm_zadv_uv,tdiff,tobs,uobs,use_3dfrc,use_camiop,vertdivq, & - vertdivt,vertdivu,vertdivv,vobs,wfld,qinitobs,scm_relax_fincl + vertdivt,vertdivu,vertdivv,vobs,wfld,qinitobs,scm_relax_fincl, & + scmlon,scmlat, & + scm_ana_direct_ttend, & + scm_use_ana_iop use time_manager, only : get_curr_calday, get_nstep, get_step_size, is_first_step use cam_abortutils, only : endrun use string_utils, only: to_upper + use shr_const_mod, only: rearth => shr_const_rearth , & ! =6.37122e6_R8 meters + pi => shr_const_pi , & + OOmega => shr_const_omega implicit none @@ -71,6 +77,7 @@ subroutine forecast( lat , nlon , ztodt , & ! ---------------------- ! character(len=*), parameter :: subname = "forecast" + real(r8),parameter :: hugebad=9.99e12_r8 ! --------------------------------------------------- ! ! x = t, u, v, q ! @@ -83,16 +90,16 @@ subroutine forecast( lat , nlon , ztodt , & integer, intent(in) :: nlon real(r8), intent(in) :: ztodt ! Twice time step unless nstep = 0 [ s ] - real(r8), intent(in) :: ps(plon) ! Surface pressure [ Pa ] - real(r8), intent(in) :: psm1(plon) ! Surface pressure [ Pa ] - real(r8), intent(in) :: psm2(plon) ! Surface pressure [ Pa ] + real(r8), intent(inout) :: ps(plon) ! Surface pressure [ Pa ] + real(r8), intent(inout) :: psm1(plon) ! Surface pressure [ Pa ] + real(r8), intent(inout) :: psm2(plon) ! Surface pressure [ Pa ] - real(r8), intent(in) :: t3m1(plev) ! Temperature [ K ] - real(r8), intent(in) :: t3m2(plev) ! Temperature [ K ] - real(r8), intent(in) :: u3m1(plev) ! Zonal wind [ m/s ] - real(r8), intent(in) :: u3m2(plev) ! Zonal wind [ m/s ] - real(r8), intent(in) :: v3m1(plev) ! Meridional wind [ m/s ] - real(r8), intent(in) :: v3m2(plev) ! Meridional wind [ m/s ] + real(r8), intent(inout) :: t3m1(plev) ! Temperature [ K ] + real(r8), intent(inout) :: t3m2(plev) ! Temperature [ K ] + real(r8), intent(inout) :: u3m1(plev) ! Zonal wind [ m/s ] + real(r8), intent(inout) :: u3m2(plev) ! Zonal wind [ m/s ] + real(r8), intent(inout) :: v3m1(plev) ! Meridional wind [ m/s ] + real(r8), intent(inout) :: v3m2(plev) ! Meridional wind [ m/s ] real(r8), intent(inout) :: q3m1(plev,pcnst) ! Tracers [ kg/kg, #/kg ] real(r8), intent(inout) :: q3m2(plev,pcnst) ! Tracers [ kg/kg, #/kg ] @@ -156,6 +163,7 @@ subroutine forecast( lat , nlon , ztodt , & real(r8) vten_zadv(plev) ! Vertical advective forcing of v [ m/s/s ] real(r8) qten_zadv(plev,pcnst) ! Vertical advective forcing of tracers [ #/kg/s, kg/kg/s ] + ! --------------------------- ! ! For 'scm_relaxation' switch ! ! --------------------------- ! @@ -165,11 +173,91 @@ subroutine forecast( lat , nlon , ztodt , & real(r8) relax_u(plev) real(r8) relax_v(plev) real(r8) relax_q(plev,pcnst) - ! +++BPM: allow linear relaxation profile + ! allow linear relaxation profile real(r8) rslope ! [optional] slope for linear relaxation profile real(r8) rycept ! [optional] y-intercept for linear relaxtion profile -!+++ BPM check what we have: + +! ------------------------------------ ! +! Quantities derived from Analyses ! +! ------------------------------------ ! +!======================================! + real(r8) dynfrcp_mom(plev) ! Scaling factor for ana-derived momentum tends + real(r8) dynfrcp_therm(plev) ! Scaling factor for ana-derived thermo tends + logical l_vectinv + real(r8) omega_ana(plev) ! Vertical pressure velocity [ Pa/s ] + real(r8) etad_ana(plev) ! "Eta dot" velocity [ Pa/s ] + real(r8) T_ana(plev), Q_ana(plev) , Tv_ana(plev) ! + real(r8) u_ana(plev), v_ana(plev) ! + real(r8) zeta_ana(plev) ! + real(r8) ps_ana , phis_ana + real(r8) ps_ana_w(plon) , phis_ana_w(plon) + real(r8) T_ana_diag(plev), Q_ana_diag(plev) ! + real(r8) u_ana_diag(plev), v_ana_diag(plev) ! + ! ----------------------------------- ! + ! vertical advective tendencies ! + ! ----------------------------------- ! + real(r8) tten_vadv_ana(plev) ! Vertical advective forcing of t [ K/s ] + real(r8) uten_vadv_ana(plev) ! Vertical advective forcing of u [ m/s/s ] + real(r8) vten_vadv_ana(plev) ! Vertical advective forcing of v [ m/s/s ] + real(r8) qten_vadv_ana(plev) ! Vertical advective forcing of tracers [ #/kg/s, kg/kg/s ] + ! ------------------------------------- ! + ! Horizontal advective/other tendencies ! + ! ------------------------------------- ! + real(r8) uten_hadv_ana(plev) ! of u [ m/s/s ] + real(r8) vten_hadv_ana(plev) ! of v [ m/s/s ] + real(r8) uten_pfrc_ana(plev) ! of u [ m/s/s ] + real(r8) vten_pfrc_ana(plev) ! of v [ m/s/s ] + real(r8) uten_vort_ana(plev) ! of u [ m/s/s ] + real(r8) vten_vort_ana(plev) ! of v [ m/s/s ] + real(r8) tten_hadv_ana(plev) ! of t [ K/s ] + real(r8) qten_hadv_ana(plev) ! of tracers [ #/kg/s, kg/kg/s ] + + !---------------------------------! + ! Adiabatic compression tendency ! + !---------------------------------! + real(r8) tten_comp_ana(plev) ! of t [ K/s ] + + + real(r8) uten_keg_ana(plev) ! of u [ m/s/s ] + real(r8) uten_prg_ana(plev) ! of u [ m/s/s ] + real(r8) uten_phig_ana(plev) ! of u [ m/s/s ] + ! ------------------------------------------ ! + ! Direct dycore or ana tendencies or quants ! + ! Not recalculated. ! + ! (not usually available, ! + ! set=-9999 if missing ) ! + ! ------------------------------------------ ! + real(r8) tten_dycore_ana(plev) ! Total direct Ana forcing of t [ K/s ] + real(r8) vten_dycore_ana(plev) ! Total direct Ana forcing of v [ m/s/s ] + real(r8) uten_dycore_ana(plev) ! Total direct Ana forcing of u [ m/s/s ] + real(r8) omega_dycore_ana(plev) ! Omega direct from Ana/dycore (not recalc) [ Pa/s ] + ! ----------------------------------- ! + ! total recalc. "dycore" tendencies ! + ! ----------------------------------- ! + real(r8) omega_recalc_ana(plev) ! Omega from Ana/dycore (recalculated) [ Pa/s ] + real(r8) tten_totdyn_ana(plev) ! Total Ana forcing of t [ K/s ] + real(r8) uten_totdyn_ana(plev) ! Total Ana forcing of u [ m/s/s ] + real(r8) vten_totdyn_ana(plev) ! Total Ana forcing of v [ m/s/s ] + real(r8) qten_totdyn_ana(plev) ! Total Ana forcing of tracers [ #/kg/s, kg/kg/s ] + real(r8) fcoriol,uten_coriol(plev),vten_coriol(plev) + real(r8) ufcstm2(plev),vfcstm2(plev) + real(r8) ufcor_0(plev),vfcor_0(plev) + real(r8) uten_totdyn_anax(plev) ! Total Ana forcing of u [ m/s/s ] + real(r8) vten_totdyn_anax(plev) ! Total Ana forcing of v [ m/s/s ] + real(r8) tfw0, tfw1, tfw2, tftotw,ztodtn,AA + integer nsubdyn,nt,nstep_curr + + !logical use_ana_iop + logical l_use_reconst_ttend ! use reconstructed T-tendency based on analysis + logical l_use_direct_ttend ! use T-tendency direct from dycore + + + l_use_reconst_ttend = .NOT.( scm_ana_direct_ttend ) + l_use_direct_ttend = .NOT.( l_use_reconst_ttend ) + + + ! check what we have: if (masterproc .and. is_first_step()) write(iulog,*) 'SCAM FORECAST REPORT: ' , & 'have_divq ', have_divq , & 'have_divt ', have_divt , & @@ -211,8 +299,6 @@ subroutine forecast( lat , nlon , ztodt , & 'relaxation ', scm_relaxation , & 'use_3dfrc ', use_3dfrc - !---BPM - ! ---------------------------- ! ! ! @@ -253,8 +339,103 @@ subroutine forecast( lat , nlon , ztodt , & ! = .false. : Use User-generated SCAM IOP file ! ! ------------------------------------------------------- ! - - if( use_camiop ) then + l_vectinv =.FALSE. + + if (scm_use_ana_iop) then + call get_ana_dynfrc_fv ( scmlon, scmlat , & + omega_ana, etad_ana, zeta_ana, & + t_ana , tv_ana , & + q_ana , & + u_ana , & + v_ana , & + ps_ana , & + phis_ana, & + uten_hadv_ana , & + vten_hadv_ana , & + uten_pfrc_ana , & + vten_pfrc_ana , & + uten_vort_ana , & + vten_vort_ana , & + qten_hadv_ana , & + tten_hadv_ana , & + uten_vadv_ana , & + vten_vadv_ana , & + tten_vadv_ana , & + qten_vadv_ana , & + tten_comp_ana , & + uten_keg_ana , & + uten_phig_ana , & + uten_prg_ana , & + uten_dycore_ana , & + vten_dycore_ana , & + tten_dycore_ana , & + omega_dycore_ana , & + omega_recalc_ana , & + u3m2, v3m2, t3m2, q3m2(:,1), & + u_ana_diag, v_ana_diag, t_ana_diag, q_ana_diag ) + else + ! set these to a "bad" value + omega_ana = HugeBad + etad_ana = HugeBad + zeta_ana = HugeBad + t_ana = HugeBad + tv_ana = HugeBad + q_ana = HugeBad + u_ana = HugeBad + v_ana = HugeBad + t_ana_diag = HugeBad + q_ana_diag = HugeBad + u_ana_diag = HugeBad + v_ana_diag = HugeBad + ps_ana = HugeBad + phis_ana = HugeBad + uten_hadv_ana = HugeBad + vten_hadv_ana = HugeBad + uten_pfrc_ana = HugeBad + vten_pfrc_ana = HugeBad + uten_vort_ana = HugeBad + vten_vort_ana = HugeBad + qten_hadv_ana = HugeBad + tten_hadv_ana = HugeBad + uten_vadv_ana = HugeBad + vten_vadv_ana = HugeBad + tten_vadv_ana = HugeBad + qten_vadv_ana = HugeBad + tten_comp_ana = HugeBad + uten_keg_ana = HugeBad + uten_phig_ana = HugeBad + uten_prg_ana = HugeBad + uten_dycore_ana = HugeBad + vten_dycore_ana = HugeBad + tten_dycore_ana = HugeBad + omega_dycore_ana = HugeBad + omega_recalc_ana = HugeBad + endif + + ! -------------------------------------------------------------- ! + ! Re-Calculate midpoint pressure levels if PS_ANA is reasonable ! + ! -------------------------------------------------------------- ! + if (ps_ana < 500000._r8 ) then + psm1=ps_ana + call plevs0( nlon, plon, plev, psm1, pintm1, pmidm1, pdelm1 ) + end if + if(l_vectinv) then + uten_totdyn_ana = uten_hadv_ana + uten_pfrc_ana + uten_vadv_ana + vten_totdyn_ana = vten_hadv_ana + vten_pfrc_ana + vten_vadv_ana + uten_totdyn_anax = uten_hadv_ana + uten_pfrc_ana + uten_vadv_ana + vten_totdyn_anax = vten_hadv_ana + vten_pfrc_ana + vten_vadv_ana + else + uten_totdyn_ana = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana + vten_totdyn_ana = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana + uten_totdyn_anax = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana + vten_totdyn_anax = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana + endif + + tten_totdyn_ana = tten_hadv_ana + tten_vadv_ana + tten_comp_ana + qten_totdyn_ana = qten_hadv_ana + qten_vadv_ana + + +if( use_camiop == .TRUE. ) then do k = 1, plev tfcst(k) = t3m2(k) + ztodt * tten_phys(k) + ztodt * divt3d(k) ufcst(k) = u3m2(k) + ztodt * uten_phys(k) + ztodt * divu3d(k) @@ -269,8 +450,9 @@ subroutine forecast( lat , nlon , ztodt , & enddo enddo - else +else if (use_camiop == .FALSE. ) then + if(scm_use_ana_iop==.FALSE. ) then ! ---------------------------------------------------------------------------- ! ! Compute 'omega'( wfldint ) at the interface from the value at the mid-point. ! ! SCAM-IOP file must provide omega at the mid-point not at the interface. ! @@ -403,19 +585,152 @@ subroutine forecast( lat , nlon , ztodt , & call endrun( subname//':: divq not on the dataset. Unable to forecast Humidity. Stopping') end if - do k = 1, plev - tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) + divt(k) + tten_zadv(k) ) - ufcst(k) = u3m2(k) + ztodt * ( uten_phys(k) + divu(k) + uten_zadv(k) ) - vfcst(k) = v3m2(k) + ztodt * ( vten_phys(k) + divv(k) + vten_zadv(k) ) - do m = 1, pcnst - qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) + divq(k,m) + qten_zadv(k,m) ) + + nstep_curr = get_nstep() + + do k = 1, plev + tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) + divt(k) + tten_zadv(k) ) + ufcst(k) = u3m2(k) + ztodt * ( uten_phys(k) + divu(k) + uten_zadv(k) ) + vfcst(k) = v3m2(k) + ztodt * ( vten_phys(k) + divv(k) + vten_zadv(k) ) + do m = 1, pcnst + qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) + divq(k,m) + qten_zadv(k,m) ) + enddo enddo - enddo + else if ( scm_use_ana_iop ==.TRUE. ) then + !-------------------------------------------------- + ! This block is the ANA/GRIDDED forcing calculation + !-------------------------------------------------- + + nstep_curr = get_nstep() + + ! This initializes SCAM with ANA/GRIDDED + ! fields. + !----------------------------------------- + if (is_first_step()) then + u3m2 = u_ana + v3m2 = v_ana + t3m2 = t_ana + q3m2(:,1) = q_ana + psm2 = ps_ana + u3m1 = u_ana + v3m1 = v_ana + t3m1 = t_ana + q3m1(:,1) = q_ana + psm1 = ps_ana + endif + + + ! ----------------------------------------------------- + ! Applied tendencies are in two + ! categories: 1) physics (includes nudging); + ! and 2) dynamics. Dynamics tendencies are + ! grouped and then scaled by dynfrcp. This is + ! to allow removal of unreliable analysis driven + ! dynamics tendencies above some pressure, + ! typically <~ 10Pa. + !------------------------------------------------------ + dynfrcp_mom(:) = 1._r8 + dynfrcp_therm(:) = 1._r8 + where( pmidm1 < 10._r8) ! changed from 10. Test. + dynfrcp_therm = 0._r8 + dynfrcp_mom = 0._r8 + end where + !------------------------------------------------------ + fcoriol = 2._r8 * OOmega * sin( scmlat * PI/180._r8 ) + uten_coriol = fcoriol * v3m2 + vten_coriol = -fcoriol * u3m2 + nsubdyn = 1 + vfcst = v3m2 + ufcst = u3m2 + ztodtn = ztodt/nsubdyn + + if (l_use_reconst_ttend) then + do nt= 1, nsubdyn + do k = 1, plev + ufcst(k) = ufcst(k) + ztodtn * ( uten_phys(k) & + + dynfrcp_mom(k) * & + ( uten_hadv_ana(k) + uten_vadv_ana(k) & + + uten_vort_ana(k) & + !! + fcoriol * vfcstm2(k) & + + uten_pfrc_ana(k) ) ) + vfcst(k) = vfcst(k) + ztodtn * ( vten_phys(k) & + + dynfrcp_mom(k) * & + ( vten_hadv_ana(k) + vten_vadv_ana(k) & + + vten_vort_ana(k) & + !! - fcoriol * ufcstm2(k) & + + vten_pfrc_ana(k) ) ) + end do + ufcstm2 = ufcst + vfcstm2 = vfcst + end do + end if + + if (l_use_direct_ttend) then + do nt= 1, nsubdyn + do k = 1, plev + ufcst(k) = ufcst(k) + ztodtn * ( uten_phys(k) & + + dynfrcp_mom(k) * uten_dycore_ana(k) ) + vfcst(k) = vfcst(k) + ztodtn * ( vten_phys(k) & + + dynfrcp_mom(k) * vten_dycore_ana(k) ) + end do + ufcstm2 = ufcst + vfcstm2 = vfcst + end do + end if + + ufcor_0 = ufcst + vfcor_0 = vfcst + + uten_totdyn_ana = uten_hadv_ana + uten_vort_ana + uten_pfrc_ana + uten_vadv_ana + vten_totdyn_ana = vten_hadv_ana + vten_vort_ana + vten_pfrc_ana + vten_vadv_ana + + ! This should always be set for diagnostics + ! (08/03/22) + !------------------------------------------- + wfld = omega_ana + + + +!! if (l_use_reconst_ttend) then + do k=1,plev + tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) & + + dynfrcp_therm(k) * & + ( tten_hadv_ana(k) + tten_vadv_ana(k) & + + tten_comp_ana(k) ) ) + end do +!! end if + +#if 0 + if (l_use_direct_ttend) then + do k=1,plev + tfcst(k) = t3m2(k) + ztodt * ( tten_phys(k) & + + dynfrcp_therm(k) * & + ( tten_dycore_ana(k) ) ) + end do + end if +#endif + + do k=1,plev + do m = 1, 1 + qfcst(1,k,m) = q3m2(k,m) + ztodt * ( qten_phys(k,m) & + + dynfrcp_therm(k) * & + ( qten_hadv_ana(k) + qten_vadv_ana(k) ) ) + enddo + enddo + + ps = ps_ana + + + end if ! END use_ana_iop IF block + + ! This code is executed regardless of use_ana_iop value ! ------------------ ! ! Diagnostic Outputs ! ! ------------------ ! - + call outfld( 'TOBS' , tobs, plon, dummy_dyndecomp ) + call outfld( 'UOBS' , uobs, plon, dummy_dyndecomp ) + call outfld( 'VOBS' , vobs, plon, dummy_dyndecomp ) call outfld( 'TTEN_XYADV' , divt, plon, dummy_dyndecomp ) call outfld( 'UTEN_XYADV' , divu, plon, dummy_dyndecomp ) call outfld( 'VTEN_XYADV' , divv, plon, dummy_dyndecomp ) @@ -438,29 +753,31 @@ subroutine forecast( lat , nlon , ztodt , & call outfld( 'UTEN_ZADV' , uten_zadv, plon, dummy_dyndecomp ) call outfld( 'VTEN_ZADV' , vten_zadv, plon, dummy_dyndecomp ) call outfld( 'QVTEN_ZADV' , qten_zadv(:,1), plon, dummy_dyndecomp ) - call outfld( 'TTEN_ZADV' , vertdivt, plon, dummy_dyndecomp ) - call outfld( 'QVTEN_ZADV' , vertdivq(:,1), plon, dummy_dyndecomp ) + !call outfld( 'TTEN_ZADV' , vertdivt, plon, dummy_dyndecomp ) + !call outfld( 'QVTEN_ZADV' , vertdivq(:,1), plon, dummy_dyndecomp ) - call outfld( 'TTEN_PHYS' , tten_phys, plon, dummy ) - call outfld( 'UTEN_PHYS' , uten_phys, plon, dummy ) - call outfld( 'VTEN_PHYS' , vten_phys, plon, dummy ) - call outfld( 'QVTEN_PHYS' , qten_phys(:,1), plon, dummy ) + call outfld( 'TTEN_COMP_IOP', tten_comp_eul, plon, dummy_dyndecomp ) - endif + call outfld( 'TTEN_PHYS' , tten_phys, plon, dummy_dyndecomp ) + call outfld( 'UTEN_PHYS' , uten_phys, plon, dummy_dyndecomp ) + call outfld( 'VTEN_PHYS' , vten_phys, plon, dummy_dyndecomp ) + call outfld( 'QVTEN_PHYS' , qten_phys(:,1), plon, dummy_dyndecomp ) + end if ! END of use_camiop IF BLOCK + + if ( scm_use_ana_iop == .FALSE. ) then ! ---------------------------------------------------------------- ! ! Used the SCAM-IOP-specified state instead of forecasted state ! ! at each time step if specified by the switch. ! ! If SCAM-IOP has 't,u,v,q' profile at a single initial time step. ! - ! ---------------------------------------------------------------- ! - + ! ---------------------------------------------------------------- ! if( scm_use_obs_T .and. have_t ) then do k = 1, plev tfcst(k) = tobs(k) enddo endif - if( scm_use_obs_uv .and. have_u .and. have_v ) then + if( scm_use_obs_uv .and. have_u .and. have_v ) then do k = 1, plev ufcst(k) = uobs(k) vfcst(k) = vobs(k) @@ -483,7 +800,7 @@ subroutine forecast( lat , nlon , ztodt , & relax_u(:) = 0._r8 relax_v(:) = 0._r8 relax_q(:plev,:pcnst) = 0._r8 - ! +++BPM: allow linear relaxation profile + ! allow linear relaxation profile ! scm_relaxation is a logical from scamMod ! scm_relax_tau_top_sec and scm_relax_tau_bot_sec are the relaxation times at top and bottom of layer ! also defined in scamMod @@ -540,7 +857,9 @@ subroutine forecast( lat , nlon , ztodt , & call outfld( 'TRELAX' , relax_T , plon, dummy ) call outfld( 'QRELAX' , relax_q(1:plev,1) , plon, dummy ) call outfld( 'TAURELAX' , rtau , plon, dummy ) - + + end if ! END of 2nd use_ana_iop BLOCK (exec for use_ana_iop=.F.) + ! --------------------------------------------------------- ! ! Assign the final forecasted state to the output variables ! ! --------------------------------------------------------- ! @@ -548,15 +867,82 @@ subroutine forecast( lat , nlon , ztodt , & t3(1:plev) = tfcst(1:plev) u3(1:plev) = ufcst(1:plev) v3(1:plev) = vfcst(1:plev) - q3(1:plev,1:pcnst) = qfcst(1,1:plev,1:pcnst) - + + if ( scm_use_ana_iop == .TRUE. ) then + q3(1:plev,1:1) = qfcst(1,1:plev,1:1) + else + q3(1:plev,1:pcnst) = qfcst(1,1:plev,1:pcnst) + endif + tdiff(1:plev) = t3(1:plev) - tobs(1:plev) qdiff(1:plev) = q3(1:plev,1) - qobs(1:plev) + call outfld( 'QDIFF' , qdiff, plon, dummy_dyndecomp ) call outfld( 'TDIFF' , tdiff, plon, dummy_dyndecomp ) + + call outfld( 'OMEGA_IOP' , wfld, plon, dummy_dyndecomp ) + call outfld( 'OMEGA_ANA' , omega_ana, plon, dummy_dyndecomp ) + call outfld( 'ETAD_ANA' , etad_ana, plon, dummy_dyndecomp ) + call outfld( 'ZETA_ANA' , zeta_ana, plon, dummy_dyndecomp ) + call outfld( 'T_ANA' , T_ana_diag, plon, dummy_dyndecomp ) + call outfld( 'Q_ANA' , Q_ana_diag, plon, dummy_dyndecomp ) + call outfld( 'TV_ANA' , Tv_ana, plon, dummy_dyndecomp ) + call outfld( 'U_ANA' , U_ana_diag, plon, dummy_dyndecomp ) + call outfld( 'V_ANA' , V_ana_diag, plon, dummy_dyndecomp ) + + ps_ana_w(:)=ps_ana + phis_ana_w(:)=phis_ana + + call outfld( 'PS_ANA' , PS_ana_w, plon, dummy_dyndecomp ) + call outfld( 'PHIS_ANA' , PHIS_ana_w, plon, dummy_dyndecomp ) + + call outfld( 'UTEN_CORIOL' , uten_coriol, plon, dummy_dyndecomp ) + call outfld( 'VTEN_CORIOL' , vten_coriol, plon, dummy_dyndecomp ) + + call outfld( 'UTEN_TOTDYN_ANA' , uten_totdyn_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_TOTDYN_ANA' , vten_totdyn_ana, plon, dummy_dyndecomp ) + call outfld( 'TTEN_TOTDYN_ANA' , tten_totdyn_ana, plon, dummy_dyndecomp ) + call outfld( 'QTEN_TOTDYN_ANA' , qten_totdyn_ana, plon, dummy_dyndecomp ) + + call outfld( 'UTEN_TOTDYN_ANAR' , uten_totdyn_anax, plon, dummy_dyndecomp ) + call outfld( 'VTEN_TOTDYN_ANAR' , vten_totdyn_anax, plon, dummy_dyndecomp ) + + call outfld( 'UTEN_DYCORE_ANA' , uten_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_DYCORE_ANA' , vten_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'TTEN_DYCORE_ANA' , tten_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'OMEGA_DYCORE_ANA', omega_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'OMEGA_RECALC_ANA', omega_recalc_ana, plon, dummy_dyndecomp ) + + call outfld( 'UTEN_DYCORE_ANA', uten_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_HADV_ANA' , uten_hadv_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_VADV_ANA' , uten_vadv_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_VORT_ANA' , uten_vort_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_KEG_ANA' , uten_keg_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_PFRC_ANA' , uten_pfrc_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_PRG_ANA' , uten_prg_ana, plon, dummy_dyndecomp ) + call outfld( 'UTEN_PHIG_ANA' , uten_phig_ana, plon, dummy_dyndecomp ) + + call outfld( 'VTEN_DYCORE_ANA', vten_dycore_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_HADV_ANA' , vten_hadv_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_VADV_ANA' , vten_vadv_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_VORT_ANA' , vten_vort_ana, plon, dummy_dyndecomp ) + call outfld( 'VTEN_PFRC_ANA' , vten_pfrc_ana, plon, dummy_dyndecomp ) + + call outfld( 'TTEN_HADV_ANA' , tten_hadv_ana, plon, dummy_dyndecomp ) + call outfld( 'TTEN_VADV_ANA' , tten_vadv_ana, plon, dummy_dyndecomp ) + call outfld( 'TTEN_COMP_ANA' , tten_comp_ana, plon, dummy_dyndecomp ) + + call outfld( 'QTEN_HADV_ANA' , qten_hadv_ana, plon, dummy_dyndecomp ) + call outfld( 'QTEN_VADV_ANA' , qten_vadv_ana, plon, dummy_dyndecomp ) + + if (have_u) call outfld( 'U_IOP' , uobs, plon, dummy_dyndecomp ) + if (have_u) call outfld( 'V_IOP' , vobs, plon, dummy_dyndecomp ) + return end subroutine forecast - end module scmforecast + + +end module scmforecast diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index f83007a62e..69a9419da1 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -4,28 +4,28 @@ module clubb_intr ! Module to interface CAM with Cloud Layers Unified by Bi-normals (CLUBB), developed ! ! by the University of Wisconsin Milwaukee Group (UWM). ! ! ! - ! CLUBB replaces the exisiting turbulence, shallow convection, and macrophysics in CAM5 ! - ! ! + ! CLUBB replaces the exisiting turbulence, shallow convection, and macrophysics in CAM5 ! + ! ! ! Lastly, a implicit diffusion solver is called, and tendencies retrieved by ! ! differencing the diffused and initial states. ! - ! ! + ! ! ! Calling sequence: ! ! ! !---------------------------Code history-------------------------------------------------------------- ! - ! Authors: P. Bogenschutz, C. Craig, A. Gettelman ! - ! Modified by: K Thayer-Calder ! - ! ! + ! Authors: P. Bogenschutz, C. Craig, A. Gettelman ! + ! Modified by: K Thayer-Calder ! + ! ! !----------------------------------------------------------------------------------------------------- ! - use shr_kind_mod, only: r8=>shr_kind_r8 + use shr_kind_mod, only: r8=>shr_kind_r8 use ppgrid, only: pver, pverp, pcols, begchunk, endchunk use phys_control, only: phys_getopts use physconst, only: rairv, cpairv, cpair, gravit, rga, latvap, latice, zvir, rh2o, karman - use spmd_utils, only: masterproc + use spmd_utils, only: masterproc use constituents, only: pcnst, cnst_add use pbl_utils, only: calc_ustar, calc_obklen - use ref_pres, only: top_lev => trop_cloud_top_lev + use ref_pres, only: top_lev => trop_cloud_top_lev use zm_conv_intr, only: zmconv_microp #ifdef CLUBB_SGS use clubb_api_module, only: pdf_parameter, implicit_coefs_terms @@ -44,7 +44,7 @@ module clubb_intr stats_rad_zt(pcols), & ! stats_rad_zt grid stats_rad_zm(pcols), & ! stats_rad_zm grid stats_sfc(pcols) ! stats_sfc - + !$omp threadprivate(stats_zt, stats_zm, stats_rad_zt, stats_rad_zm, stats_sfc) type(grid), target :: dummy_gr @@ -68,7 +68,7 @@ module clubb_intr stats_init_clubb, & stats_zt, stats_zm, stats_sfc, & stats_rad_zt, stats_rad_zm, & - stats_end_timestep_clubb, & + stats_end_timestep_clubb, & #endif clubb_readnl, & clubb_init_cnst, & @@ -94,7 +94,7 @@ module clubb_intr integer, parameter :: & grid_type = 3, & ! The 2 option specifies stretched thermodynamic levels hydromet_dim = 0 ! The hydromet array in SAM-CLUBB is currently 0 elements - + ! Even though sclr_dim is set to 0, the dimension here is set to 1 to prevent compiler errors ! See github ticket larson-group/cam#133 for details real(r8), parameter, dimension(1) :: & @@ -106,28 +106,28 @@ module clubb_intr theta0 = 300._r8, & ! Reference temperature [K] ts_nudge = 86400._r8, & ! Time scale for u/v nudging (not used) [s] p0_clubb = 100000._r8 - - integer, parameter :: & + + integer, parameter :: & sclr_dim = 0 ! Higher-order scalars, set to zero real(r8), parameter :: & wp3_const = 1._r8 ! Constant to add to wp3 when moments are advected - - real(r8), parameter :: & + + real(r8), parameter :: & wpthlp_const = 10.0_r8 ! Constant to add to wpthlp when moments are advected - - real(r8), parameter :: & + + real(r8), parameter :: & wprtp_const = 0.01_r8 ! Constant to add to wprtp when moments are advected - - real(r8), parameter :: & + + real(r8), parameter :: & rtpthlp_const = 0.01_r8 ! Constant to add to rtpthlp when moments are advected - + real(r8), parameter :: unset_r8 = huge(1.0_r8) integer, parameter :: unset_i = huge(1) - ! Commonly used temperature for the melting temp of ice crystals [K] - real(r8), parameter :: meltpt_temp = 268.15_r8 - + ! Commonly used temperature for the melting temp of ice crystals [K] + real(r8), parameter :: meltpt_temp = 268.15_r8 + real(r8) :: clubb_timestep = unset_r8 ! Default CLUBB timestep, unless overwriten by namelist real(r8) :: clubb_rnevap_effic = unset_r8 @@ -183,7 +183,7 @@ module clubb_intr real(r8) :: clubb_detliq_rad = unset_r8 real(r8) :: clubb_detice_rad = unset_r8 real(r8) :: clubb_detphase_lowtemp = unset_r8 - + integer :: & clubb_iiPDF_type, & ! Selected option for the two-component normal ! (double Gaussian) PDF type to use for the w, rt, @@ -192,7 +192,7 @@ module clubb_intr clubb_ipdf_call_placement = unset_i ! Selected option for the placement of the call to ! CLUBB's PDF. - + logical :: & clubb_l_use_precip_frac, & ! Flag to use precipitation fraction in KK microphysics. The ! precipitation fraction is automatically set to 1 when this @@ -255,8 +255,8 @@ module clubb_intr ! that is linearized in terms of wp3. ! (Requires ADG1 PDF and clubb_l_standard_term_ta). clubb_l_godunov_upwind_wpxp_ta, & ! This flag determines whether we want to use an upwind - ! differencing approximation rather than a centered - ! differencing for turbulent advection terms. + ! differencing approximation rather than a centered + ! differencing for turbulent advection terms. ! It affects wpxp only. clubb_l_godunov_upwind_xpyp_ta, & ! This flag determines whether we want to use an upwind ! differencing approximation rather than a centered @@ -299,14 +299,14 @@ module clubb_intr clubb_l_damp_wp3_Skw_squared, & ! Set damping on wp3 to use Skw^2 rather than Skw^4 clubb_l_prescribed_avg_deltaz, & ! used in adj_low_res_nu. If .true., avg_deltaz = deltaz clubb_l_update_pressure ! Flag for having CLUBB update pressure and exner - + ! Constant parameters logical, parameter, private :: & l_implemented = .true., & ! Implemented in a host model (always true) l_host_applies_sfc_fluxes = .false. ! Whether the host model applies the surface fluxes - + logical, parameter, private :: & - apply_to_heat = .false. ! Apply WACCM energy fixer to heat or not (.true. = yes (duh)) + apply_to_heat = .false. ! Apply WACCM energy fixer to heat or not (.true. = yes (duh)) logical :: lq(pcnst) logical :: prog_modal_aero @@ -318,8 +318,8 @@ module clubb_intr integer :: history_budget_histfile_num integer :: edsclr_dim ! Number of scalars to transport in CLUBB integer :: offset - -! define physics buffer indicies here + +! define physics buffer indicies here integer :: & wp2_idx, & ! vertical velocity variances wp3_idx, & ! third moment of vertical velocity @@ -378,8 +378,8 @@ module clubb_intr naai_idx, & ! ice number concentration prer_evap_idx, & ! rain evaporation rate qrl_idx, & ! longwave cooling rate - radf_idx, & - qsatfac_idx, & ! subgrid cloud water saturation scaling factor + radf_idx, & + qsatfac_idx, & ! subgrid cloud water saturation scaling factor ice_supersat_idx, & ! ice cloud fraction for SILHS rcm_idx, & ! Cloud water mixing ratio for SILHS ztodt_idx ! physics timestep for SILHS @@ -399,7 +399,7 @@ module clubb_intr pdf_zm_varnce_w_2_idx, & pdf_zm_mixt_frac_idx - integer, public :: & + integer, public :: & ixthlp2 = 0, & ixwpthlp = 0, & ixwprtp = 0, & @@ -418,7 +418,7 @@ module clubb_intr dnlfzm_idx = -1, & ! ZM detrained convective cloud water num concen. dnifzm_idx = -1 ! ZM detrained convective cloud ice num concen. - ! Output arrays for CLUBB statistics + ! Output arrays for CLUBB statistics real(r8), allocatable, dimension(:,:,:) :: out_zt, out_zm, out_radzt, out_radzm, out_sfc character(len=16) :: eddy_scheme ! Default set in phys_control.F90 @@ -432,14 +432,14 @@ module clubb_intr #ifdef CLUBB_SGS type(pdf_parameter), target, allocatable, public, protected :: & pdf_params_chnk(:) ! PDF parameters (thermo. levs.) [units vary] - + type(pdf_parameter), target, allocatable :: pdf_params_zm_chnk(:) ! PDF parameters on momentum levs. [units vary] - + type(implicit_coefs_terms), target, allocatable :: pdf_implicit_coefs_terms_chnk(:,:) ! PDF impl. coefs. & expl. terms [units vary] #endif contains - + ! =============================================================================== ! ! ! ! =============================================================================== ! @@ -458,12 +458,12 @@ subroutine clubb_register_cam( ) ! Register physics buffer fields and constituents ! !------------------------------------------------ ! - ! Add CLUBB fields to pbuf + ! Add CLUBB fields to pbuf use physics_buffer, only: pbuf_add_field, dtype_r8, dyn_time_lvls use subcol_utils, only: subcol_get_scheme - + call phys_getopts( eddy_scheme_out = eddy_scheme, & - deep_scheme_out = deep_scheme, & + deep_scheme_out = deep_scheme, & history_budget_out = history_budget, & history_budget_histfile_num_out = history_budget_histfile_num ) subcol_scheme = subcol_get_scheme() @@ -478,7 +478,7 @@ subroutine clubb_register_cam( ) cnst_names =(/'THLP2 ','RTP2 ','RTPTHLP','WPTHLP ','WPRTP ','WP2 ','WP3 ','UP2 ','VP2 '/) do_cnst=.true. ! If CLUBB moments are advected, do not output them automatically which is typically done. Some moments - ! need a constant added to them before they are advected, thus this would corrupt the output. + ! need a constant added to them before they are advected, thus this would corrupt the output. ! Users should refer to the "XXXX_CLUBB" (THLP2_CLUBB for instance) output variables for these moments call cnst_add(trim(cnst_names(1)),0._r8,0._r8,0._r8,ixthlp2,longname='second moment vertical velocity',cam_outfld=.false.) call cnst_add(trim(cnst_names(2)),0._r8,0._r8,0._r8,ixrtp2,longname='second moment rtp',cam_outfld=.false.) @@ -508,7 +508,7 @@ subroutine clubb_register_cam( ) call pbuf_add_field('CMELIQ', 'physpkg',dtype_r8, (/pcols,pver/), cmeliq_idx) call pbuf_add_field('QSATFAC', 'physpkg',dtype_r8, (/pcols,pver/), qsatfac_idx) - + call pbuf_add_field('WP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), wp2_idx) call pbuf_add_field('WP3_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), wp3_idx) call pbuf_add_field('WPTHLP_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), wpthlp_idx) @@ -517,7 +517,7 @@ subroutine clubb_register_cam( ) call pbuf_add_field('RTP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), rtp2_idx) call pbuf_add_field('THLP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), thlp2_idx) call pbuf_add_field('UP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), up2_idx) - call pbuf_add_field('VP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), vp2_idx) + call pbuf_add_field('VP2_nadv', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), vp2_idx) call pbuf_add_field('RTP3', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), rtp3_idx) call pbuf_add_field('THLP3', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), thlp3_idx) @@ -563,7 +563,7 @@ subroutine clubb_register_cam( ) call pbuf_add_field('pdf_zm_var_w_2', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), pdf_zm_varnce_w_2_idx) call pbuf_add_field('pdf_zm_mixt_frac', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), pdf_zm_mixt_frac_idx) -#endif +#endif end subroutine clubb_register_cam ! =============================================================================== ! @@ -587,14 +587,14 @@ function clubb_implements_cnst(name) end function clubb_implements_cnst - + ! =============================================================================== ! ! ! ! =============================================================================== ! subroutine clubb_init_cnst(name, latvals, lonvals, mask, q) #ifdef CLUBB_SGS - use clubb_api_module, only: w_tol_sqd, rt_tol, thl_tol + use clubb_api_module, only: w_tol_sqd, rt_tol, thl_tol #endif !----------------------------------------------------------------------- ! @@ -667,7 +667,7 @@ subroutine clubb_init_cnst(name, latvals, lonvals, mask, q) end subroutine clubb_init_cnst - + ! =============================================================================== ! ! ! ! =============================================================================== ! @@ -681,7 +681,7 @@ subroutine clubb_readnl(nlfile) use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_logical, mpi_real8, & mpi_integer use clubb_mf, only: clubb_mf_readnl - + use clubb_api_module, only: & set_default_clubb_config_flags_api, & ! Procedure(s) initialize_clubb_config_flags_type_api, & @@ -694,7 +694,7 @@ subroutine clubb_readnl(nlfile) character(len=*), parameter :: sub = 'clubb_readnl' - logical :: clubb_history = .false., clubb_rad_history = .false. ! Stats enabled (T/F) + logical :: clubb_history = .false., clubb_rad_history = .false. ! Stats enabled (T/F) logical :: clubb_cloudtop_cooling = .false., clubb_rainevap_turb = .false. integer :: iunit, read_status, ierr @@ -728,20 +728,20 @@ subroutine clubb_readnl(nlfile) clubb_l_use_C11_Richardson, clubb_l_use_shear_Richardson, & clubb_l_brunt_vaisala_freq_moist, clubb_l_use_thvm_in_bv_freq, & clubb_l_rcm_supersat_adj, clubb_l_damp_wp3_Skw_squared, & - clubb_l_lmm_stepping, & - clubb_l_e3sm_config, & + clubb_l_lmm_stepping, & + clubb_l_e3sm_config, & clubb_l_use_tke_in_wp3_pr_turb_term, clubb_l_use_tke_in_wp2_wp3_K_dfsn, & clubb_l_smooth_Heaviside_tau_wpxp !----- Begin Code ----- - ! Determine if we want clubb_history to be output + ! Determine if we want clubb_history to be output clubb_history = .false. ! Initialize to false l_stats = .false. ! Initialize to false l_output_rad_files = .false. ! Initialize to false do_cldcool = .false. ! Initialize to false do_rainturb = .false. ! Initialize to false - + ! Initialize namelist variables to clubb defaults call set_default_clubb_config_flags_api( clubb_iiPDF_type, & ! Out clubb_ipdf_call_placement, & ! Out @@ -880,7 +880,7 @@ subroutine clubb_readnl(nlfile) call mpi_bcast(clubb_c_K10, 1, mpi_real8, mstrid, mpicom, ierr) if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_c_K10") call mpi_bcast(clubb_c_K10h, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_c_K10h") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_c_K10h") call mpi_bcast(clubb_beta, 1, mpi_real8, mstrid, mpicom, ierr) if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_beta") call mpi_bcast(clubb_C2rt, 1, mpi_real8, mstrid, mpicom, ierr) @@ -928,21 +928,21 @@ subroutine clubb_readnl(nlfile) call mpi_bcast(clubb_do_energyfix, 1, mpi_logical, mstrid, mpicom, ierr) if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_do_energyfix") call mpi_bcast(clubb_C_invrs_tau_bkgnd, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_bkgnd") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_bkgnd") call mpi_bcast(clubb_C_invrs_tau_sfc, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_sfc") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_sfc") call mpi_bcast(clubb_C_invrs_tau_shear, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_shear") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_shear") call mpi_bcast(clubb_C_invrs_tau_N2, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2") call mpi_bcast(clubb_C_invrs_tau_N2_wp2, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_wp2") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_wp2") call mpi_bcast(clubb_C_invrs_tau_N2_xp2, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_xp2") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_xp2") call mpi_bcast(clubb_C_invrs_tau_N2_wpxp, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_wpxp") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_wpxp") call mpi_bcast(clubb_C_invrs_tau_N2_clear_wp3, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_clear_wp3") + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_C_invrs_tau_N2_clear_wp3") call mpi_bcast(clubb_lmin_coef, 1, mpi_real8, mstrid, mpicom, ierr) if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_lmin_coef") call mpi_bcast(clubb_skw_max_mag, 1, mpi_real8, mstrid, mpicom, ierr) @@ -1017,10 +1017,10 @@ subroutine clubb_readnl(nlfile) ! Overwrite defaults if they are true if (clubb_history) l_stats = .true. - if (clubb_rad_history) l_output_rad_files = .true. + if (clubb_rad_history) l_output_rad_files = .true. if (clubb_cloudtop_cooling) do_cldcool = .true. if (clubb_rainevap_turb) do_rainturb = .true. - + ! Check that all namelists have been set if(clubb_timestep == unset_r8) call endrun(sub//": FATAL: clubb_timestep is not set") if(clubb_rnevap_effic == unset_r8) call endrun(sub//": FATAL:clubb_rnevap_effic is not set") @@ -1078,7 +1078,7 @@ subroutine clubb_readnl(nlfile) if(clubb_detice_rad == unset_r8) call endrun(sub//": FATAL: clubb_detice_rad not set") if(clubb_ipdf_call_placement == unset_i) call endrun(sub//": FATAL: clubb_ipdf_call_placement not set") if(clubb_detphase_lowtemp == unset_r8) call endrun(sub//": FATAL: clubb_detphase_lowtemp not set") - if(clubb_detphase_lowtemp >= meltpt_temp) & + if(clubb_detphase_lowtemp >= meltpt_temp) & call endrun(sub//": ERROR: clubb_detphase_lowtemp must be less than 268.15 K") call initialize_clubb_config_flags_type_api( clubb_iiPDF_type, & ! In @@ -1126,7 +1126,7 @@ subroutine clubb_readnl(nlfile) clubb_l_e3sm_config, & ! In clubb_l_vary_convect_depth, & ! In clubb_l_use_tke_in_wp3_pr_turb_term, & ! In - clubb_l_use_tke_in_wp2_wp3_K_dfsn, & ! In + clubb_l_use_tke_in_wp2_wp3_K_dfsn, & ! In clubb_l_smooth_Heaviside_tau_wpxp, & ! In clubb_config_flags ) ! Out @@ -1215,7 +1215,7 @@ subroutine clubb_ini_cam(pbuf2d) #ifdef CLUBB_SGS real(kind=time_precision) :: dum1, dum2, dum3 - + ! The similar name to clubb_history is unfortunate... logical :: history_amwg, history_clubb @@ -1274,7 +1274,7 @@ subroutine clubb_ini_cam(pbuf2d) pdf_params_chnk(begchunk:endchunk), & pdf_params_zm_chnk(begchunk:endchunk), & pdf_implicit_coefs_terms_chnk(pcols,begchunk:endchunk) ) - + do j = 1, pcols, 1 do l = begchunk, endchunk, 1 call init_pdf_implicit_coefs_terms_api( pverp+1-top_lev, sclr_dim, & @@ -1283,10 +1283,10 @@ subroutine clubb_ini_cam(pbuf2d) enddo ! j = 1, pcols, 1 ! ----------------------------------------------------------------- ! - ! Determine how many constituents CLUBB will transport. Note that - ! CLUBB does not transport aerosol consituents. Therefore, need to + ! Determine how many constituents CLUBB will transport. Note that + ! CLUBB does not transport aerosol consituents. Therefore, need to ! determine how many aerosols constituents there are and subtract that - ! off of pcnst (the total consituents) + ! off of pcnst (the total consituents) ! ----------------------------------------------------------------- ! call phys_getopts(prog_modal_aero_out=prog_modal_aero, & @@ -1294,7 +1294,7 @@ subroutine clubb_ini_cam(pbuf2d) history_clubb_out=history_clubb) ! Select variables to apply tendencies back to CAM - + ! Initialize all consituents to true to start lq(1:pcnst) = .true. edsclr_dim = pcnst @@ -1308,12 +1308,12 @@ subroutine clubb_ini_cam(pbuf2d) if (prog_modal_aero) then ! Turn off modal aerosols and decrement edsclr_dim accordingly call rad_cnst_get_info(0, nmodes=nmodes) - + do m = 1, nmodes call rad_cnst_get_mode_num_idx(m, lptr) lq(lptr)=.false. edsclr_dim = edsclr_dim-1 - + call rad_cnst_get_info(0, m, nspec=nspec) do l = 1, nspec call rad_cnst_get_mam_mmr_idx(m, l, lptr) @@ -1321,7 +1321,7 @@ subroutine clubb_ini_cam(pbuf2d) edsclr_dim = edsclr_dim-1 end do end do - + ! In addition, if running with MAM, droplet number is transported ! in dropmixnuc, therefore we do NOT want CLUBB to apply transport ! tendencies to avoid double counted. Else, we apply tendencies. @@ -1345,7 +1345,7 @@ subroutine clubb_ini_cam(pbuf2d) l_stats_samp = .false. l_grads = .false. - ! Overwrite defaults if needbe + ! Overwrite defaults if needbe if (l_stats) l_stats_samp = .true. ! Define physics buffers indexes @@ -1354,7 +1354,7 @@ subroutine clubb_ini_cam(pbuf2d) ast_idx = pbuf_get_index('AST') ! Stratiform cloud fraction alst_idx = pbuf_get_index('ALST') ! Liquid stratiform cloud fraction aist_idx = pbuf_get_index('AIST') ! Ice stratiform cloud fraction - qlst_idx = pbuf_get_index('QLST') ! Physical in-stratus LWC + qlst_idx = pbuf_get_index('QLST') ! Physical in-stratus LWC qist_idx = pbuf_get_index('QIST') ! Physical in-stratus IWC dp_frac_idx = pbuf_get_index('DP_FRAC') ! Deep convection cloud fraction icwmrdp_idx = pbuf_get_index('ICWMRDP') ! In-cloud deep convective mixing ratio @@ -1385,17 +1385,17 @@ subroutine clubb_ini_cam(pbuf2d) ! ----------------------------------------------------------------- ! ! Define number of tracers for CLUBB to diffuse - ! ----------------------------------------------------------------- ! - + ! ----------------------------------------------------------------- ! + if (clubb_l_do_expldiff_rtm_thlm) then offset = 2 ! diffuse temperature and moisture explicitly - edsclr_dim = edsclr_dim + offset + edsclr_dim = edsclr_dim + offset end if - + ! ----------------------------------------------------------------- ! ! Setup CLUBB core ! ----------------------------------------------------------------- ! - + ! Read in parameters for CLUBB. Just read in default values call set_default_parameters_api( & C1, C1b, C1c, C2rt, C2thl, C2rtthl, & @@ -1449,7 +1449,7 @@ subroutine clubb_ini_cam(pbuf2d) clubb_params ) ! Fill in dummy arrays for height. Note that these are overwrote - ! at every CLUBB step to physical values. + ! at every CLUBB step to physical values. do k=1,nlev+1 zt_g(k) = ((k-1)*1000._r8)-500._r8 ! this is dummy garbage zi_g(k) = (k-1)*1000._r8 ! this is dummy garbage @@ -1504,7 +1504,7 @@ subroutine clubb_ini_cam(pbuf2d) clubb_params(iC_invrs_tau_N2_xp2) = clubb_C_invrs_tau_N2_xp2 clubb_params(iC_invrs_tau_N2_wpxp) = clubb_C_invrs_tau_N2_wpxp clubb_params(iC_invrs_tau_N2_clear_wp3) = clubb_C_invrs_tau_N2_clear_wp3 - + ! Set up CLUBB core. Note that some of these inputs are overwritten ! when clubb_tend_cam is called. The reason is that heights can change ! at each time step, which is why dummy arrays are read in here for heights @@ -1551,12 +1551,12 @@ subroutine clubb_ini_cam(pbuf2d) ! ----------------------------------------------------------------- ! ! Initialize eddy diffusivity module - + ntop_eddy = 1 ! if >1, must be <= nbot_molec nbot_eddy = pver ! currently always pver - + call init_hb_diff( gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, karman, eddy_scheme ) - + ! ----------------------------------------------------------------- ! ! Add output fields for the history files ! ----------------------------------------------------------------- ! @@ -1584,7 +1584,7 @@ subroutine clubb_ini_cam(pbuf2d) call addfld ('WPRCP_CLUBB', (/ 'ilev' /), 'A', 'W/m2', 'Liquid Water Flux') call addfld ('CLOUDFRAC_CLUBB', (/ 'lev' /), 'A', 'fraction', 'Cloud Fraction') call addfld ('RCMINLAYER_CLUBB', (/ 'lev' /), 'A', 'kg/kg', 'Cloud Water in Layer') - call addfld ('CLOUDCOVER_CLUBB', (/ 'lev' /), 'A', 'fraction', 'Cloud Cover') + call addfld ('CLOUDCOVER_CLUBB', (/ 'lev' /), 'A', 'fraction', 'Cloud Cover') call addfld ('WPTHVP_CLUBB', (/ 'ilev' /), 'A', 'W/m2', 'Buoyancy Flux') call addfld ('RVMTEND_CLUBB', (/ 'lev' /), 'A', 'kg/kg /s', 'Water vapor tendency') call addfld ('STEND_CLUBB', (/ 'lev' /), 'A', 'J/(kg s)', 'Static energy tendency') @@ -1593,7 +1593,7 @@ subroutine clubb_ini_cam(pbuf2d) call addfld ('UTEND_CLUBB', (/ 'lev' /), 'A', 'm/s /s', 'U-wind Tendency') call addfld ('VTEND_CLUBB', (/ 'lev' /), 'A', 'm/s /s', 'V-wind Tendency') call addfld ('ZT_CLUBB', (/ 'lev' /), 'A', 'm', 'Thermodynamic Heights') - call addfld ('ZM_CLUBB', (/ 'ilev' /), 'A', 'm', 'Momentum Heights') + call addfld ('ZM_CLUBB', (/ 'ilev' /), 'A', 'm', 'Momentum Heights') call addfld ('UM_CLUBB', (/ 'lev' /), 'A', 'm/s', 'Zonal Wind') call addfld ('VM_CLUBB', (/ 'lev' /), 'A', 'm/s', 'Meridional Wind') call addfld ('WM_ZT_CLUBB', (/ 'lev' /), 'A', 'm/s', 'Vertical Velocity') @@ -1611,8 +1611,8 @@ subroutine clubb_ini_cam(pbuf2d) call addfld ('FQTENDICE', (/ 'lev' /), 'A', 'fraction', 'Frequency of Ice Saturation Adjustment') call addfld ('DPDLFLIQ', (/ 'lev' /), 'A', 'kg/kg/s', 'Detrained liquid water from deep convection') - call addfld ('DPDLFICE', (/ 'lev' /), 'A', 'kg/kg/s', 'Detrained ice from deep convection') - call addfld ('DPDLFT', (/ 'lev' /), 'A', 'K/s', 'T-tendency due to deep convective detrainment') + call addfld ('DPDLFICE', (/ 'lev' /), 'A', 'kg/kg/s', 'Detrained ice from deep convection') + call addfld ('DPDLFT', (/ 'lev' /), 'A', 'K/s', 'T-tendency due to deep convective detrainment') call addfld ('RELVAR', (/ 'lev' /), 'A', '-', 'Relative cloud water variance') call addfld ('CLUBB_GRID_SIZE', horiz_only, 'A', 'm', 'Horizontal grid box size seen by CLUBB') @@ -1652,7 +1652,7 @@ subroutine clubb_ini_cam(pbuf2d) call addfld ( 'edmf_S_AWV' , (/ 'ilev' /), 'A', 'm2/s2' , 'Sum of a_i*w_i*v_i (EDMF)' ) call addfld ( 'edmf_thlflx' , (/ 'ilev' /), 'A', 'W/m2' , 'thl flux (EDMF)' ) call addfld ( 'edmf_qtflx' , (/ 'ilev' /), 'A', 'W/m2' , 'qt flux (EDMF)' ) - end if + end if ! Initialize statistics, below are dummy variables dum1 = 300._r8 @@ -1660,13 +1660,13 @@ subroutine clubb_ini_cam(pbuf2d) dum3 = 300._r8 if (l_stats) then - + do i=1, pcols call stats_init_clubb( .true., dum1, dum2, & nlev+1, nlev+1, nlev+1, dum3, & stats_zt(i), stats_zm(i), stats_sfc(i), & stats_rad_zt(i), stats_rad_zm(i)) - end do + end do allocate(out_zt(pcols,pverp,stats_zt(1)%num_output_fields)) allocate(out_zm(pcols,pverp,stats_zm(1)%num_output_fields)) @@ -1676,12 +1676,12 @@ subroutine clubb_ini_cam(pbuf2d) allocate(out_radzm(pcols,pverp,stats_rad_zm(1)%num_output_fields)) end if - + ! ----------------------------------------------------------------- ! ! Make all of this output default, this is not CLUBB history ! ----------------------------------------------------------------- ! - - if (clubb_do_adv .or. history_clubb) then + + if (clubb_do_adv .or. history_clubb) then call add_default('RELVAR', 1, ' ') call add_default('RHO_CLUBB', 1, ' ') call add_default('UP2_CLUBB', 1, ' ') @@ -1714,14 +1714,14 @@ subroutine clubb_ini_cam(pbuf2d) call add_default('UTEND_CLUBB', 1, ' ') call add_default('VTEND_CLUBB', 1, ' ') call add_default('ZT_CLUBB', 1, ' ') - call add_default('ZM_CLUBB', 1, ' ') + call add_default('ZM_CLUBB', 1, ' ') call add_default('UM_CLUBB', 1, ' ') call add_default('VM_CLUBB', 1, ' ') call add_default('WM_ZT_CLUBB', 1, ' ') call add_default('PBLH', 1, ' ') call add_default('CONCLD', 1, ' ') end if - + if (history_amwg) then call add_default('PBLH', 1, ' ') end if @@ -1750,10 +1750,10 @@ subroutine clubb_ini_cam(pbuf2d) call add_default( 'edmf_qtflx' , 1, ' ') end if - if (history_budget) then + if (history_budget) then call add_default('DPDLFLIQ', history_budget_histfile_num, ' ') call add_default('DPDLFICE', history_budget_histfile_num, ' ') - call add_default('DPDLFT', history_budget_histfile_num, ' ') + call add_default('DPDLFT', history_budget_histfile_num, ' ') call add_default('STEND_CLUBB', history_budget_histfile_num, ' ') call add_default('RCMTEND_CLUBB', history_budget_histfile_num, ' ') call add_default('RIMTEND_CLUBB', history_budget_histfile_num, ' ') @@ -1761,7 +1761,7 @@ subroutine clubb_ini_cam(pbuf2d) call add_default('UTEND_CLUBB', history_budget_histfile_num, ' ') call add_default('VTEND_CLUBB', history_budget_histfile_num, ' ') end if - + ! --------------- ! ! First step? ! @@ -1780,12 +1780,12 @@ subroutine clubb_ini_cam(pbuf2d) call pbuf_set_field(pbuf2d, thlp2_idx, thl_tol**2) call pbuf_set_field(pbuf2d, up2_idx, w_tol_sqd) call pbuf_set_field(pbuf2d, vp2_idx, w_tol_sqd) - + call pbuf_set_field(pbuf2d, rtp3_idx, 0.0_r8) call pbuf_set_field(pbuf2d, thlp3_idx, 0.0_r8) call pbuf_set_field(pbuf2d, up3_idx, 0.0_r8) call pbuf_set_field(pbuf2d, vp3_idx, 0.0_r8) - + call pbuf_set_field(pbuf2d, upwp_idx, 0.0_r8) call pbuf_set_field(pbuf2d, vpwp_idx, 0.0_r8) call pbuf_set_field(pbuf2d, wpthvp_idx, 0.0_r8) @@ -1823,10 +1823,10 @@ subroutine clubb_ini_cam(pbuf2d) call pbuf_set_field(pbuf2d, pdf_zm_mixt_frac_idx, 0.0_r8) end if - + ! The following is physpkg, so it needs to be initialized every time call pbuf_set_field(pbuf2d, fice_idx, 0.0_r8) - + ! --------------- ! ! End ! ! Initialization ! @@ -1834,21 +1834,21 @@ subroutine clubb_ini_cam(pbuf2d) #endif end subroutine clubb_ini_cam - - + + ! =============================================================================== ! ! ! ! =============================================================================== ! subroutine clubb_tend_cam( & state, ptend_all, pbuf, hdtime, & - cmfmc, cam_in, & + cmfmc, cam_in, & macmic_it, cld_macmic_num_steps,dlf, det_s, det_ice) !------------------------------------------------------------------------------- -! Description: Provide tendencies of shallow convection, turbulence, and +! Description: Provide tendencies of shallow convection, turbulence, and ! macrophysics from CLUBB to CAM -! +! ! Author: Cheryl Craig, March 2011 ! Modifications: Pete Bogenschutz, March 2011 and onward ! Origin: Based heavily on UWM clubb_init.F90 @@ -1864,12 +1864,12 @@ subroutine clubb_tend_cam( & use constituents, only: cnst_get_ind, cnst_type use camsrfexch, only: cam_in_t - use time_manager, only: is_first_step + use time_manager, only: is_first_step use cam_abortutils, only: endrun use cam_logfile, only: iulog use tropopause, only: tropopause_findChemTrop use time_manager, only: get_nstep, is_first_restart_step - + #ifdef CLUBB_SGS use hb_diff, only: pblintd use scamMOD, only: single_column,scm_clubb_iop_name @@ -1907,13 +1907,13 @@ subroutine clubb_tend_cam( & use macrop_driver, only: liquid_macro_tend use clubb_mf, only: integrate_mf - + use perf_mod #endif implicit none - + ! --------------- ! ! Input Auguments ! ! --------------- ! @@ -1925,11 +1925,11 @@ subroutine clubb_tend_cam( & real(r8), intent(in) :: cmfmc(pcols,pverp) ! convective mass flux--m sub c [kg/m2/s] integer, intent(in) :: cld_macmic_num_steps ! number of mac-mic iterations integer, intent(in) :: macmic_it ! number of mac-mic iterations - + ! ---------------------- ! ! Input-Output Auguments ! ! ---------------------- ! - + type(physics_buffer_desc), pointer :: pbuf(:) ! ---------------------- ! @@ -1938,11 +1938,11 @@ subroutine clubb_tend_cam( & type(physics_ptend), intent(out) :: ptend_all ! package tendencies - ! These two variables are needed for energy check + ! These two variables are needed for energy check real(r8), intent(out) :: det_s(pcols) ! Integral of detrained static energy from ice real(r8), intent(out) :: det_ice(pcols) ! Integral of detrained ice for energy check - + ! --------------- ! ! Local Variables ! ! --------------- ! @@ -1951,22 +1951,22 @@ subroutine clubb_tend_cam( & type(physics_state) :: state1 ! Local copy of state variable type(physics_ptend) :: ptend_loc ! Local tendency from processes, added up to return as ptend_all - + integer :: i, j, k, t, ixind, nadv integer :: ixcldice, ixcldliq, ixnumliq, ixnumice, ixq integer :: itim_old integer :: ncol, lchnk ! # of columns, and chunk identifier integer :: err_code ! Diagnostic, for if some calculation goes amiss. - integer :: icnt + integer :: icnt logical :: lq2(pcnst) integer :: iter, ierr - + integer :: clubbtop(pcols) - + real(r8) :: frac_limit, ic_limit - real(r8) :: dtime ! CLUBB time step [s] + real(r8) :: dtime ! CLUBB time step [s] real(r8) :: edsclr_in(pcols,pverp+1-top_lev,edsclr_dim) ! Scalars to be diffused through CLUBB [units vary] real(r8) :: wp2_in(pcols,pverp+1-top_lev) ! vertical velocity variance (CLUBB) [m^2/s^2] real(r8) :: wp3_in(pcols,pverp+1-top_lev) ! third moment vertical velocity [m^3/s^3] @@ -1993,7 +1993,7 @@ subroutine clubb_tend_cam( & real(r8) :: um_in(pcols,pverp+1-top_lev) ! meridional wind [m/s] real(r8) :: vm_in(pcols,pverp+1-top_lev) ! zonal wind [m/s] real(r8) :: pre_in(pcols,pverp+1-top_lev) ! input for precip evaporation - real(r8) :: rtp2_mc_out(pcols,pverp+1-top_lev) ! total water tendency from rain evap + real(r8) :: rtp2_mc_out(pcols,pverp+1-top_lev) ! total water tendency from rain evap real(r8) :: thlp2_mc_out(pcols,pverp+1-top_lev) ! thetal tendency from rain evap real(r8) :: wprtp_mc_out(pcols,pverp+1-top_lev) real(r8) :: wpthlp_mc_out(pcols,pverp+1-top_lev) @@ -2025,15 +2025,15 @@ subroutine clubb_tend_cam( & real(r8) :: w_up_in_cloud_out(pcols,pverp+1-top_lev) real(r8) :: zt_g(pcols,pverp+1-top_lev) ! Thermodynamic grid of CLUBB [m] real(r8) :: zi_g(pcols,pverp+1-top_lev) ! Momentum grid of CLUBB [m] - real(r8) :: zt_out(pcols,pverp) ! output for the thermo CLUBB grid [m] + real(r8) :: zt_out(pcols,pverp) ! output for the thermo CLUBB grid [m] real(r8) :: zi_out(pcols,pverp) ! output for momentum CLUBB grid [m] real(r8) :: fcor(pcols) ! Coriolis forcing [s^-1] real(r8) :: sfc_elevation(pcols) ! Elevation of ground [m AMSL] [m] real(r8) :: ubar ! surface wind [m/s] - real(r8) :: ustar ! surface stress [m/s] + real(r8) :: ustar ! surface stress [m/s] real(r8) :: z0 ! roughness height [m] real(r8) :: thlm_forcing(pcols,pverp+1-top_lev) ! theta_l forcing (thermodynamic levels) [K/s] - real(r8) :: rtm_forcing(pcols,pverp+1-top_lev) ! r_t forcing (thermodynamic levels) [(kg/kg)/s] + real(r8) :: rtm_forcing(pcols,pverp+1-top_lev) ! r_t forcing (thermodynamic levels) [(kg/kg)/s] real(r8) :: um_forcing(pcols,pverp+1-top_lev) ! u wind forcing (thermodynamic levels) [m/s/s] real(r8) :: vm_forcing(pcols,pverp+1-top_lev) ! v wind forcing (thermodynamic levels) [m/s/s] real(r8) :: rtm_ref(pcols,pverp+1-top_lev) ! Initial profile of rtm [kg/kg] @@ -2051,7 +2051,7 @@ subroutine clubb_tend_cam( & real(r8) :: wpthlp_sfc(pcols) ! w' theta_l' at surface [(m K)/s] real(r8) :: wprtp_sfc(pcols) ! w' r_t' at surface [(kg m)/( kg s)] real(r8) :: upwp_sfc(pcols) ! u'w' at surface [m^2/s^2] - real(r8) :: vpwp_sfc(pcols) ! v'w' at surface [m^2/s^2] + real(r8) :: vpwp_sfc(pcols) ! v'w' at surface [m^2/s^2] real(r8) :: sclrm_forcing(pcols,pverp+1-top_lev,sclr_dim) ! Passive scalar forcing [{units vary}/s] real(r8) :: wpsclrp_sfc(pcols,sclr_dim) ! Scalar flux at surface [{units vary} m/s] real(r8) :: edsclrm_forcing(pcols,pverp+1-top_lev,edsclr_dim)! Eddy passive scalar forcing [{units vary}/s] @@ -2116,7 +2116,7 @@ subroutine clubb_tend_cam( & real(r8) :: thl2_zt(pcols,pverp+1-top_lev) ! CLUBB Theta-l variance on thermo levs [K^2] real(r8) :: thl2_zt_out(pcols, pverp) ! CLUBB Theta-l variance on thermo levs real(r8) :: wp2_zt(pcols,pverp+1-top_lev) ! CLUBB W variance on theromo levs [m^2/s^2] - real(r8) :: wp2_zt_out(pcols, pverp) + real(r8) :: wp2_zt_out(pcols, pverp) real(r8) :: dlf_liq_out(pcols, pverp) ! Detrained liquid water from ZM [kg/kg/s] real(r8) :: dlf_ice_out(pcols, pverp) ! Detrained ice water from ZM [kg/kg/s] real(r8) :: wm_zt_out(pcols, pverp) ! CLUBB mean W on thermo levs output [m/s] @@ -2151,7 +2151,7 @@ subroutine clubb_tend_cam( & ! --------------- ! ! Pointers ! ! --------------- ! - + real(r8), pointer, dimension(:,:) :: wp2 ! vertical velocity variance [m^2/s^2] real(r8), pointer, dimension(:,:) :: wp3 ! third moment of vertical velocity [m^3/s^3] real(r8), pointer, dimension(:,:) :: wpthlp ! turbulent flux of thetal [m/s K] @@ -2201,16 +2201,16 @@ subroutine clubb_tend_cam( & real(r8), pointer, dimension(:,:) :: qlst ! Physical in-stratus LWC [kg/kg] real(r8), pointer, dimension(:,:) :: qist ! Physical in-stratus IWC [kg/kg] real(r8), pointer, dimension(:,:) :: deepcu ! deep convection cloud fraction [fraction] - real(r8), pointer, dimension(:,:) :: shalcu ! shallow convection cloud fraction [fraction] + real(r8), pointer, dimension(:,:) :: shalcu ! shallow convection cloud fraction [fraction] real(r8), pointer, dimension(:,:) :: khzm ! CLUBB's eddy diffusivity of heat/moisture on momentum (i.e. interface) levels [m^2/s] real(r8), pointer, dimension(:) :: pblh ! planetary boundary layer height [m] real(r8), pointer, dimension(:,:) :: tke ! turbulent kinetic energy [m^2/s^2] real(r8), pointer, dimension(:,:) :: dp_icwmr ! deep convection in cloud mixing ratio [kg/kg] - real(r8), pointer, dimension(:,:) :: ice_supersat_frac ! Cloud fraction of ice clouds (pverp)[fraction] + real(r8), pointer, dimension(:,:) :: ice_supersat_frac ! Cloud fraction of ice clouds (pverp)[fraction] real(r8), pointer, dimension(:,:) :: relvar ! relative cloud water variance [-] real(r8), pointer, dimension(:,:) :: accre_enhan ! accretion enhancement factor [-] real(r8), pointer, dimension(:,:) :: naai - real(r8), pointer, dimension(:,:) :: cmeliq + real(r8), pointer, dimension(:,:) :: cmeliq real(r8), pointer, dimension(:,:) :: cmfmc_sh ! Shallow convective mass flux--m subc (pcols,pverp) [kg/m2/s/] real(r8), pointer, dimension(:,:) :: qsatfac @@ -2287,12 +2287,12 @@ subroutine clubb_tend_cam( & intrinsic :: max character(len=*), parameter :: subr='clubb_tend_cam' - + type(pdf_parameter) :: pdf_params_single_col - + type(grid) :: gr(pcols) integer :: begin_height, end_height - + type(nu_vertical_res_dep) :: nu_vert_res_dep(pcols) ! Vertical resolution dependent nu values real(r8) :: lmin(pcols) @@ -2326,7 +2326,7 @@ subroutine clubb_tend_cam( & frac_limit = 0.01_r8 ic_limit = 1.e-12_r8 - if (clubb_do_adv) then + if (clubb_do_adv) then apply_const = 1._r8 ! Initialize to one, only if CLUBB's moments are advected else apply_const = 0._r8 ! Never want this if CLUBB's moments are not advected @@ -2365,7 +2365,7 @@ subroutine clubb_tend_cam( & ! Determine time step of physics buffer itim_old = pbuf_old_tim_idx() - ! Establish associations between pointers and physics buffer fields + ! Establish associations between pointers and physics buffer fields call pbuf_get_field(pbuf, wp2_idx, wp2, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) call pbuf_get_field(pbuf, wp3_idx, wp3, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) call pbuf_get_field(pbuf, wpthlp_idx, wpthlp, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) @@ -2410,7 +2410,7 @@ subroutine clubb_tend_cam( & call pbuf_get_field(pbuf, rtm_idx, rtm, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) call pbuf_get_field(pbuf, um_idx, um, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) call pbuf_get_field(pbuf, vm_idx, vm, start=(/1,1,itim_old/), kount=(/pcols,pverp,1/)) - + call pbuf_get_field(pbuf, tke_idx, tke) call pbuf_get_field(pbuf, qrl_idx, qrl) call pbuf_get_field(pbuf, radf_idx, radf_clubb) @@ -2444,10 +2444,10 @@ subroutine clubb_tend_cam( & call pbuf_get_field(pbuf, wprtp_mc_zt_idx, wprtp_mc_zt) call pbuf_get_field(pbuf, wpthlp_mc_zt_idx, wpthlp_mc_zt) call pbuf_get_field(pbuf, rtpthlp_mc_zt_idx, rtpthlp_mc_zt) - + ! Allocate arrays in the single column versions of pdf_params call init_pdf_params_api( pverp+1-top_lev, 1, pdf_params_single_col ) - + ! Allocate pdf_params only if they aren't allocated already. if ( .not. allocated(pdf_params_chnk(lchnk)%mixt_frac) ) then call init_pdf_params_api( pverp+1-top_lev, ncol, pdf_params_chnk(lchnk) ) @@ -2456,15 +2456,15 @@ subroutine clubb_tend_cam( & ! Initialize the apply_const variable (note special logic is due to eularian backstepping) if (clubb_do_adv .and. (is_first_step() .or. all(wpthlp(1:ncol,1:pver) == 0._r8))) then - apply_const = 0._r8 ! On first time through do not remove constant - ! from moments since it has not been added yet + apply_const = 0._r8 ! On first time through do not remove constant + ! from moments since it has not been added yet end if ! Set the ztodt timestep in pbuf for SILHS ztodtptr(:) = 1.0_r8*hdtime ! Define the grid box size. CLUBB needs this information to determine what - ! the maximum length scale should be. This depends on the column for + ! the maximum length scale should be. This depends on the column for ! variable mesh grids and lat-lon grids if (single_column) then ! If single column specify grid box size to be something @@ -2472,7 +2472,7 @@ subroutine clubb_tend_cam( & grid_dx(:) = 100000._r8 grid_dy(:) = 100000._r8 else - + call grid_size(state1, grid_dx, grid_dy) end if @@ -2487,11 +2487,11 @@ subroutine clubb_tend_cam( & lq2(1) = .TRUE. lq2(ixcldice) = .TRUE. lq2(ixnumice) = .TRUE. - + latsub = latvap + latice - + call physics_ptend_init(ptend_loc, state%psetcols, 'iceadj', ls=.true., lq=lq2 ) - + stend(:ncol,:)=0._r8 qvtend(:ncol,:)=0._r8 qitend(:ncol,:)=0._r8 @@ -2505,9 +2505,9 @@ subroutine clubb_tend_cam( & ! update local copy of state with the tendencies ptend_loc%q(:ncol,top_lev:pver,1)=qvtend(:ncol,top_lev:pver) - ptend_loc%q(:ncol,top_lev:pver,ixcldice)=qitend(:ncol,top_lev:pver) + ptend_loc%q(:ncol,top_lev:pver,ixcldice)=qitend(:ncol,top_lev:pver) ptend_loc%q(:ncol,top_lev:pver,ixnumice)=initend(:ncol,top_lev:pver) - ptend_loc%s(:ncol,top_lev:pver)=stend(:ncol,top_lev:pver) + ptend_loc%s(:ncol,top_lev:pver)=stend(:ncol,top_lev:pver) ! Add the ice tendency to the output tendency call physics_ptend_sum(ptend_loc, ptend_all, ncol) @@ -2521,61 +2521,61 @@ subroutine clubb_tend_cam( & call outfld( 'QVTENDICE', qvtend, pcols, lchnk ) call outfld( 'QITENDICE', qitend, pcols, lchnk ) call outfld( 'NITENDICE', initend, pcols, lchnk ) - + end if ! Determine CLUBB time step and make it sub-step friendly - ! For now we want CLUBB time step to be 5 min since that is + ! For now we want CLUBB time step to be 5 min since that is ! what has been scientifically validated. However, there are certain - ! instances when a 5 min time step will not be possible (based on + ! instances when a 5 min time step will not be possible (based on ! host model time step or on macro-micro sub-stepping - dtime = clubb_timestep - - ! Now check to see if dtime is greater than the host model - ! (or sub stepped) time step. If it is, then simply - ! set it equal to the host (or sub step) time step. + dtime = clubb_timestep + + ! Now check to see if dtime is greater than the host model + ! (or sub stepped) time step. If it is, then simply + ! set it equal to the host (or sub step) time step. ! This section is mostly to deal with small host model ! time steps (or small sub-steps) - + if (dtime > hdtime) then dtime = hdtime end if - + ! Now check to see if CLUBB time step divides evenly into ! the host model time step. If not, force it to divide evenly. ! We also want it to be 5 minutes or less. This section is ! mainly for host model time steps that are not evenly divisible - ! by 5 minutes - + ! by 5 minutes + if (mod(hdtime,dtime) .ne. 0) then dtime = hdtime/2._r8 - do while (dtime > clubb_timestep) + do while (dtime > clubb_timestep) dtime = dtime/2._r8 end do end if ! If resulting host model time step and CLUBB time step do not divide evenly - ! into each other, have model throw a fit. + ! into each other, have model throw a fit. if (mod(hdtime,dtime) .ne. 0) then call endrun(subr//': CLUBB time step and HOST time step NOT compatible') end if - - ! determine number of timesteps CLUBB core should be advanced, - ! host time step divided by CLUBB time step + + ! determine number of timesteps CLUBB core should be advanced, + ! host time step divided by CLUBB time step nadv = max(hdtime/dtime,1._r8) - + ! Initialize forcings for transported scalars to zero - + sclrm_forcing(:,:,:) = 0._r8 edsclrm_forcing(:,:,:) = 0._r8 sclrm(:,:,:) = 0._r8 - + ! Compute inverse exner function consistent with CLUBB's definition, which uses a constant - ! surface pressure. CAM's exner (in state) does not. Therefore, for consistent - ! treatment with CLUBB code, anytime exner is needed to treat CLUBB variables + ! surface pressure. CAM's exner (in state) does not. Therefore, for consistent + ! treatment with CLUBB code, anytime exner is needed to treat CLUBB variables ! (such as thlm), use "inv_exner_clubb" otherwise use the exner in state do k=1,pver @@ -2589,9 +2589,9 @@ subroutine clubb_tend_cam( & do i=1,ncol inv_exner_clubb_surf(i) = 1._r8/((state1%pmid(i,pver)/p0_clubb)**(rairv(i,pver,lchnk)/cpairv(i,pver,lchnk))) enddo - - ! At each CLUBB call, initialize mean momentum and thermo CLUBB state - ! from the CAM state + + ! At each CLUBB call, initialize mean momentum and thermo CLUBB state + ! from the CAM state do k=1,pver ! loop over levels do i=1,ncol ! loop over columns @@ -2604,11 +2604,11 @@ subroutine clubb_tend_cam( & * inv_exner_clubb(i,k) if (clubb_do_adv) then - if (macmic_it == 1) then + if (macmic_it == 1) then - ! Note that some of the moments below can be positive or negative. - ! Remove a constant that was added to prevent dynamics from clipping - ! them to prevent dynamics from making them positive. + ! Note that some of the moments below can be positive or negative. + ! Remove a constant that was added to prevent dynamics from clipping + ! them to prevent dynamics from making them positive. thlp2(i,k) = state1%q(i,k,ixthlp2) rtp2(i,k) = state1%q(i,k,ixrtp2) rtpthlp(i,k) = state1%q(i,k,ixrtpthlp) - (rtpthlp_const*apply_const) @@ -2623,12 +2623,12 @@ subroutine clubb_tend_cam( & enddo enddo - + if (clubb_do_adv) then - ! If not last step of macmic loop then set apply_const back to - ! zero to prevent output from being corrupted. - if (macmic_it == cld_macmic_num_steps) then - apply_const = 1._r8 + ! If not last step of macmic loop then set apply_const back to + ! zero to prevent output from being corrupted. + if (macmic_it == cld_macmic_num_steps) then + apply_const = 1._r8 else apply_const = 0._r8 end if @@ -2638,8 +2638,8 @@ subroutine clubb_tend_cam( & um(1:ncol,pverp) = state1%u(1:ncol,pver) vm(1:ncol,pverp) = state1%v(1:ncol,pver) thlm(1:ncol,pverp) = thlm(1:ncol,pver) - - if (clubb_do_adv) then + + if (clubb_do_adv) then thlp2(1:ncol,pverp)=thlp2(1:ncol,pver) rtp2(1:ncol,pverp)=rtp2(1:ncol,pver) rtpthlp(1:ncol,pverp)=rtpthlp(1:ncol,pver) @@ -2679,24 +2679,24 @@ subroutine clubb_tend_cam( & s_awv_output(:,:) = 0._r8 mf_thlflx_output(:,:) = 0._r8 mf_qtflx_output(:,:) = 0._r8 - + call t_startf("clubb_tend_cam_i_loop") ! Determine Coriolis force at given latitude. This is never used ! when CLUBB is implemented in a host model, therefore just set ! to zero. - fcor(:) = 0._r8 + fcor(:) = 0._r8 ! Define the CLUBB momentum grid (in height, units of m) do k=1, nlev+1 - do i=1, ncol + do i=1, ncol zi_g(i,k) = state1%zi(i,pverp-k+1)-state1%zi(i,pver+1) - end do + end do end do ! Define the CLUBB thermodynamic grid (in units of m) do k=1, nlev - do i=1, ncol + do i=1, ncol zt_g(i,k+1) = state1%zm(i,pver-k+1)-state1%zi(i,pver+1) end do end do @@ -2706,18 +2706,18 @@ subroutine clubb_tend_cam( & dz_g(i,k) = state1%zi(i,k)-state1%zi(i,k+1) ! compute thickness end do end do - - ! Thermodynamic ghost point is below surface + + ! Thermodynamic ghost point is below surface do i=1, ncol zt_g(i,1) = -1._r8*zt_g(i,2) end do - + do i=1, ncol ! Set the elevation of the surface sfc_elevation(i) = state1%zi(i,pver+1) end do - ! Compute thermodynamic stuff needed for CLUBB on thermo levels. + ! Compute thermodynamic stuff needed for CLUBB on thermo levels. ! Inputs for the momentum levels are set below setup_clubb core do k=1,nlev do i=1, ncol @@ -2731,18 +2731,18 @@ subroutine clubb_tend_cam( & invrs_rho_ds_zt(i,k+1) = 1._r8/(rho_ds_zt(i,k+1)) ! full state (moist) variables - p_in_Pa(i,k+1) = state1%pmid(i,pver-k+1) + p_in_Pa(i,k+1) = state1%pmid(i,pver-k+1) exner(i,k+1) = 1._r8/inv_exner_clubb(i,pver-k+1) thv(i,k+1) = state1%t(i,pver-k+1)*inv_exner_clubb(i,pver-k+1)*(1._r8+zvir*state1%q(i,pver-k+1,ixq) & -state1%q(i,pver-k+1,ixcldliq)) rho_zt(i,k+1) = rga*state1%pdel(i,pver-k+1)/dz_g(i,pver-k+1) - rfrzm(i,k+1) = state1%q(i,pver-k+1,ixcldice) + rfrzm(i,k+1) = state1%q(i,pver-k+1,ixcldice) radf(i,k+1) = radf_clubb(i,pver-k+1) qrl_clubb(i,k+1) = qrl(i,pver-k+1)/(cpairv(i,k,lchnk)*state1%pdel(i,pver-k+1)) end do end do - - ! Compute mean w wind on thermo grid, convert from omega to w + + ! Compute mean w wind on thermo grid, convert from omega to w do k=1,nlev do i=1,ncol wm_zt(i,k+1) = -1._r8*(state1%omega(i,pver-k+1)-state1%omega(i,pver))/(rho_zt(i,k+1)*gravit) @@ -2767,8 +2767,8 @@ subroutine clubb_tend_cam( & qrl_clubb(i,1) = qrl_clubb(i,2) wm_zt(i,1) = wm_zt(i,2) end do - - + + ! ------------------------------------------------- ! ! Begin case specific code for SCAM cases. ! ! This section of code block is NOT called in ! @@ -2786,21 +2786,21 @@ subroutine clubb_tend_cam( & ! Compute surface wind (ubar) ubar = sqrt(um(1,pver)**2+vm(1,pver)**2) if (ubar < 0.25_r8) ubar = 0.25_r8 - + ! Below denotes case specifics for surface momentum ! and thermodynamic fluxes, depending on the case - ! Define ustar (based on case, if not variable) + ! Define ustar (based on case, if not variable) ustar = 0.25_r8 ! Initialize ustar in case no case - + if(trim(scm_clubb_iop_name) == 'BOMEX_5day') then ustar = 0.28_r8 end if - + if(trim(scm_clubb_iop_name) == 'ATEX_48hr') then ustar = 0.30_r8 end if - + if(trim(scm_clubb_iop_name) == 'RICO_3day') then ustar = 0.28_r8 end if @@ -2808,23 +2808,23 @@ subroutine clubb_tend_cam( & if(trim(scm_clubb_iop_name) == 'arm97' .or. trim(scm_clubb_iop_name) == 'gate' .or. & trim(scm_clubb_iop_name) == 'toga' .or. trim(scm_clubb_iop_name) == 'mpace' .or. & trim(scm_clubb_iop_name) == 'ARM_CC') then - + bflx22(1) = (gravit/theta0)*wpthlp_sfc(1) - ustar = diag_ustar(zt_g(1,2),bflx22(1),ubar,zo(1)) + ustar = diag_ustar(zt_g(1,2),bflx22(1),ubar,zo(1)) end if - - ! Compute the surface momentum fluxes, if this is a SCAM simulation + + ! Compute the surface momentum fluxes, if this is a SCAM simulation upwp_sfc(1) = -um(1,pver)*ustar**2/ubar vpwp_sfc(1) = -vm(1,pver)*ustar**2/ubar - + end if - ! Define surface sources for transported variables for diffusion, will + ! Define surface sources for transported variables for diffusion, will ! be zero as these tendencies are done in vertical_diffusion do ixind=1,edsclr_dim do i=1,ncol wpedsclrp_sfc(i,ixind) = 0._r8 - end do + end do end do ! Set stats output and increment equal to CLUBB and host dt @@ -2833,10 +2833,10 @@ subroutine clubb_tend_cam( & stats_nsamp = nint(stats_tsamp/dtime) stats_nout = nint(stats_tout/dtime) - - ! Heights need to be set at each timestep. Therefore, recall - ! setup_grid and setup_parameters for this. - + + ! Heights need to be set at each timestep. Therefore, recall + ! setup_grid and setup_parameters for this. + ! Set-up CLUBB core at each CLUBB call because heights can change ! Important note: do not make any calls that use CLUBB grid-height ! operators (such as zt2zm_api, etc.) until AFTER the @@ -2872,7 +2872,7 @@ subroutine clubb_tend_cam( & vm_ref(:,:) = 0.0_r8 ug(:,:) = 0.0_r8 vg(:,:) = 0.0_r8 - + ! Add forcings for SILHS covariance contributions rtp2_forcing(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, rtp2_mc_zt(1:ncol,:) ) thlp2_forcing(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, thlp2_mc_zt(1:ncol,:) ) @@ -2886,7 +2886,7 @@ subroutine clubb_tend_cam( & wprtp_mc_zt(:,:) = 0.0_r8 wpthlp_mc_zt(:,:) = 0.0_r8 rtpthlp_mc_zt(:,:) = 0.0_r8 - + ! Compute some inputs from the thermodynamic grid ! to the momentum grid @@ -2895,21 +2895,21 @@ subroutine clubb_tend_cam( & invrs_rho_ds_zm(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, invrs_rho_ds_zt(1:ncol,:)) thv_ds_zm(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, thv_ds_zt(1:ncol,:)) wm_zm(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, wm_zt(1:ncol,:)) - + ! Surface fluxes provided by host model - do i=1,ncol + do i=1,ncol wpthlp_sfc(i) = cam_in%shf(i)/(cpairv(i,pver,lchnk)*rho_zt(i,2))! Sensible heat flux wpthlp_sfc(i) = wpthlp_sfc(i)*inv_exner_clubb_surf(i) ! Potential temperature flux wprtp_sfc(i) = cam_in%cflx(i,1)/rho_zt(i,2) ! Moisture flux (check rho) upwp_sfc(i) = cam_in%wsx(i)/rho_zt(i,2) ! Surface meridional momentum flux - vpwp_sfc(i) = cam_in%wsy(i)/rho_zt(i,2) ! Surface zonal momentum flux + vpwp_sfc(i) = cam_in%wsy(i)/rho_zt(i,2) ! Surface zonal momentum flux end do - + ! Need to flip arrays around for CLUBB core do k=1,nlev+1 do i=1,ncol - um_in(i,k) = um(i,pverp-k+1) - vm_in(i,k) = vm(i,pverp-k+1) + um_in(i,k) = um(i,pverp-k+1) + vm_in(i,k) = vm(i,pverp-k+1) upwp_in(i,k) = upwp(i,pverp-k+1) vpwp_in(i,k) = vpwp(i,pverp-k+1) wpthvp_in(i,k) = wpthvp(i,pverp-k+1) @@ -2966,17 +2966,17 @@ subroutine clubb_tend_cam( & do i=1,ncol rcm_inout(i,1) = rcm_inout(i,2) end do - + do k=2,nlev+1 do i=1,ncol pre_in(i,k) = prer_evap(i,pverp-k+1) end do end do - + do i=1,ncol pre_in(i,1) = pre_in(i,2) end do - + ! Initialize these to prevent crashing behavior do k=1,nlev+1 do i=1,ncol @@ -3003,7 +3003,7 @@ subroutine clubb_tend_cam( & end do end do end do - + do ixind=1, hydromet_dim do k=1, nlev+1 do i=1, ncol @@ -3018,7 +3018,7 @@ subroutine clubb_tend_cam( & ! pressure,exner on momentum grid needed for mass flux calc. if (do_clubb_mf) then - + do k=1,pver do i=1,ncol kappa_zt(i,k+1) = (rairv(i,pver-k+1,lchnk)/cpairv(i,pver-k+1,lchnk)) @@ -3026,7 +3026,7 @@ subroutine clubb_tend_cam( & invrs_exner_zt(i,k+1) = inv_exner_clubb(i,pver-k+1) end do end do - + do i=1,ncol kappa_zt(i,1) = kappa_zt(i,2) qc_zt(i,1) = qc_zt(i,2) @@ -3034,21 +3034,21 @@ subroutine clubb_tend_cam( & end do kappa_zm(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, kappa_zt(1:ncol,:)) - + do k=1,pverp do i=1,ncol p_in_Pa_zm(i,k) = state1%pint(i,pverp-k+1) invrs_exner_zm(i,k) = 1._r8/((p_in_Pa_zm(i,k)/p0_clubb)**(kappa_zm(i,k))) end do end do - + end if - - + + if (clubb_do_adv) then if (macmic_it == 1) then - - wp2_in(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, wp2_in(1:ncol,:)) + + wp2_in(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, wp2_in(1:ncol,:)) wpthlp_in(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, wpthlp_in(1:ncol,:)) wprtp_in(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, wprtp_in(1:ncol,:)) up2_in(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, up2_in(1:ncol,:)) @@ -3066,49 +3066,49 @@ subroutine clubb_tend_cam( & vp2_in(i,k) = max(w_tol_sqd,vp2_in(i,k)) end do end do - + end if end if - ! Do the same for tracers + ! Do the same for tracers icnt=0 do ixind=1,pcnst - if (lq(ixind)) then - + if (lq(ixind)) then + icnt = icnt+1 - + do k=1,nlev do i=1,ncol edsclr_in(i,k+1,icnt) = state1%q(i,pver-k+1,ixind) end do end do - + do i=1,ncol edsclr_in(i,1,icnt) = edsclr_in(i,2,icnt) end do - + end if end do - - if (clubb_l_do_expldiff_rtm_thlm) then + + if (clubb_l_do_expldiff_rtm_thlm) then do k=1,nlev do i=1, ncol edsclr_in(i,k+1,icnt+1) = thlm(i,pver-k+1) edsclr_in(i,k+1,icnt+2) = rtm(i,pver-k+1) end do end do - + do i=1, ncol edsclr_in(i,1,icnt+1) = edsclr_in(i,2,icnt+1) - edsclr_in(i,1,icnt+2) = edsclr_in(i,2,icnt+2) + edsclr_in(i,1,icnt+2) = edsclr_in(i,2,icnt+2) end do - + end if do t=1,nadv ! do needed number of "sub" timesteps for each CAM step - + ! Increment the statistics then being stats timestep if (l_stats) then call stats_begin_timestep_api(t, stats_nsamp, stats_nout) @@ -3118,18 +3118,18 @@ subroutine clubb_tend_cam( & !###################### CALL MF DIAGNOSTIC PLUMES ###################### !####################################################################### if (do_clubb_mf) then - + do k=2,pverp do i=1, ncol dzt(i,k) = zi_g(i,k) - zi_g(i,k-1) end do end do - + do i=1, ncol dzt(i,1) = dzt(i,2) invrs_dzt(i,:) = 1._r8/dzt(i,:) end do - + rtm_zm_in(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, rtm_in(1:ncol,:) ) thlm_zm_in(1:ncol,:) = zt2zm_api( pverp+1-top_lev, ncol, gr, thlm_in(1:ncol,:) ) @@ -3158,19 +3158,19 @@ subroutine clubb_tend_cam( & rtm_forcing(i,1) = 0._r8 thlm_forcing(i,1)= 0._r8 end do - + do k=2,pverp do i=1, ncol rtm_forcing(i,k) = rtm_forcing(i,k) - invrs_rho_ds_zt(i,k) * invrs_dzt(i,k) * & ((rho_ds_zm(i,k) * mf_qtflx(i,k)) - (rho_ds_zm(i,k-1) * mf_qtflx(i,k-1))) - + thlm_forcing(i,k) = thlm_forcing(i,k) - invrs_rho_ds_zt(i,k) * invrs_dzt(i,k) * & ((rho_ds_zm(i,k) * mf_thlflx(i,k)) - (rho_ds_zm(i,k-1) * mf_thlflx(i,k-1))) end do end do end if - + ! Arrays are allocated as if they have pcols grid columns, but there can be less. ! Only pass clubb_core the number of columns (ncol) with valid data. ! Advance CLUBB CORE one timestep in the future @@ -3211,8 +3211,8 @@ subroutine clubb_tend_cam( & qclvar_out(:ncol,:), thlprcp_out(:ncol,:), & wprcp_out(:ncol,:), w_up_in_cloud_out(:ncol,:), & rcm_in_layer_out(:ncol,:), cloud_cover_out(:ncol,:), invrs_tau_zm_out(:ncol,:) ) - - + + ! Note that CLUBB does not produce an error code specific to any column, and ! one value only for the entire chunk if ( err_code == clubb_fatal_error ) then @@ -3221,59 +3221,59 @@ subroutine clubb_tend_cam( & write(fstderr,*) "LON: Range:", state1%lon(1), " -- ", state1%lon(ncol) call endrun(subr//': Fatal error in CLUBB library') end if - + if (do_rainturb) then - - do i=1, ncol - rvm_in(i,:) = rtm_in(i,:) - rcm_inout(i,:) + + do i=1, ncol + rvm_in(i,:) = rtm_in(i,:) - rcm_inout(i,:) end do - - do i=1, ncol - + + do i=1, ncol + call copy_multi_pdf_params_to_single( pdf_params_chnk(lchnk), i, & pdf_params_single_col) - - + + call update_xp2_mc_api( gr(i), nlev+1, dtime, cloud_frac_inout(i,:), & rcm_inout(i,:), rvm_in(i,:), thlm_in(i,:), wm_zt(i,:), exner(i,:), pre_in(i,:), pdf_params_single_col, & rtp2_mc_out(i,:), thlp2_mc_out(i,:), & wprtp_mc_out(i,:), wpthlp_mc_out(i,:), & rtpthlp_mc_out(i,:)) end do - - do i=1, ncol + + do i=1, ncol dum1 = (1._r8 - cam_in%landfrac(i)) - ! update turbulent moments based on rain evaporation + ! update turbulent moments based on rain evaporation rtp2_in(i,:) = rtp2_in(i,:) + clubb_rnevap_effic * dum1 * rtp2_mc_out(i,:) * dtime - thlp2_in(i,:) = thlp2_in(i,:) + clubb_rnevap_effic * dum1 * thlp2_mc_out(i,:) * dtime + thlp2_in(i,:) = thlp2_in(i,:) + clubb_rnevap_effic * dum1 * thlp2_mc_out(i,:) * dtime wprtp_in(i,:) = wprtp_in(i,:) + clubb_rnevap_effic * dum1 * wprtp_mc_out(i,:) * dtime wpthlp_in(i,:) = wpthlp_in(i,:) + clubb_rnevap_effic * dum1 * wpthlp_mc_out(i,:) * dtime - + end do - - end if - + + end if + if (do_cldcool) then - + rcm_out_zm(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, rcm_inout(1:ncol,:)) qrl_zm(1:ncol,:) = zt2zm_api(pverp+1-top_lev, ncol, gr, qrl_clubb(1:ncol,:)) thlp2_rad_out(:,:) = 0._r8 - + do i=1, ncol call calculate_thlp2_rad_api(nlev+1, rcm_out_zm(i,:), thlprcp_out(i,:), qrl_zm(i,:), clubb_params, & thlp2_rad_out(i,:)) end do - + do i=1, ncol thlp2_in(i,:) = thlp2_in(i,:) + thlp2_rad_out(i,:) * dtime thlp2_in(i,:) = max(thl_tol**2,thlp2_in(i,:)) end do - + end if - + ! Check to see if stats should be output, here stats are read into ! output arrays to make them conformable to CAM output if (l_stats) then @@ -3286,16 +3286,16 @@ subroutine clubb_tend_cam( & enddo ! end time loop if (clubb_do_adv) then - if (macmic_it == cld_macmic_num_steps) then - - wp2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, wp2_in(1:ncol,:)) + if (macmic_it == cld_macmic_num_steps) then + + wp2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, wp2_in(1:ncol,:)) wpthlp_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, wpthlp_in(1:ncol,:)) wprtp_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, wprtp_in(1:ncol,:)) up2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, up2_in(1:ncol,:)) vp2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, vp2_in(1:ncol,:)) thlp2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, thlp2_in(1:ncol,:)) rtp2_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, rtp2_in(1:ncol,:)) - rtpthlp_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, rtpthlp_in(1:ncol,:)) + rtpthlp_in(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, rtpthlp_in(1:ncol,:)) do k=1,nlev+1 do i=1, ncol @@ -3306,16 +3306,16 @@ subroutine clubb_tend_cam( & vp2_in(i,k) = max(w_tol_sqd, vp2_in(i,k)) end do end do - + end if end if - + ! Convert RTP2 and THLP2 to thermo grid for output rtp2_zt(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, rtp2_in(1:ncol,:)) thl2_zt(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, thlp2_in(1:ncol,:)) wp2_zt(1:ncol,:) = zm2zt_api(pverp+1-top_lev, ncol, gr, wp2_in(1:ncol,:)) - ! Arrays need to be "flipped" to CAM grid + ! Arrays need to be "flipped" to CAM grid do k=1, nlev+1 do i=1, ncol um(i,pverp-k+1) = um_in(i,k) @@ -3371,18 +3371,18 @@ subroutine clubb_tend_cam( & rtp2_zt_out(i,pverp-k+1) = rtp2_zt(i,k) thl2_zt_out(i,pverp-k+1) = thl2_zt(i,k) wp2_zt_out(i,pverp-k+1) = wp2_zt(i,k) - + end do end do do k=1, nlev+1 do i=1, ncol - + mean_rt = pdf_params_chnk(lchnk)%mixt_frac(i,k) & * pdf_params_chnk(lchnk)%rt_1(i,k) & + ( 1.0_r8 - pdf_params_chnk(lchnk)%mixt_frac(i,k) ) & * pdf_params_chnk(lchnk)%rt_2(i,k) - + pdfp_rtp2(i,pverp-k+1) = pdf_params_chnk(lchnk)%mixt_frac(i,k) & * ( ( pdf_params_chnk(lchnk)%rt_1(i,k) - mean_rt )**2 & + pdf_params_chnk(lchnk)%varnce_rt_1(i,k) ) & @@ -3448,18 +3448,18 @@ subroutine clubb_tend_cam( & khzm(i,k) = 0._r8 qclvar(i,k) = 2._r8 end do - end do + end do ! enforce zero tracer tendencies above the top_lev level -- no change icnt=0 do ixind=1,pcnst - if (lq(ixind)) then + if (lq(ixind)) then icnt=icnt+1 - + do i=1, ncol edsclr_out(i,:top_lev-1,icnt) = state1%q(i,:top_lev-1,ixind) end do - + end if end do @@ -3469,61 +3469,60 @@ subroutine clubb_tend_cam( & ! Section below is concentrated on energy fixing for conservation. ! There are two steps to this process. The first is to remove any tendencies - ! CLUBB may have produced above where it is active due to roundoff. + ! CLUBB may have produced above where it is active due to roundoff. ! The second is to provider a fixer because CLUBB and CAM's thermodynamic - ! variables are different. + ! variables are different. ! Initialize clubbtop with the chemistry topopause top, to prevent CLUBB from - ! firing up in the stratosphere + ! firing up in the stratosphere do i=1, ncol clubbtop(i) = troplev(i) do while ((rtp2(i,clubbtop(i)) <= 1.e-15_r8 .and. rcm(i,clubbtop(i)) == 0._r8) .and. clubbtop(i) < pver-1) clubbtop(i) = clubbtop(i) + 1 - end do + end do end do - + ! Compute static energy using CLUBB's variables do k=1,pver do i=1, ncol clubb_s(i,k) = cpairv(i,k,lchnk) * thlm(i,k) / inv_exner_clubb(i,k) & + latvap * rcm(i,k) & + gravit * state1%zm(i,k) + state1%phis(i) - end do - end do - - + end do + end do + ! Compute integrals above layer where CLUBB is active se_upper_a(:) = 0._r8 ! energy in layers above where CLUBB is active AFTER CLUBB is called se_upper_b(:) = 0._r8 ! energy in layers above where CLUBB is active BEFORE CLUBB is called tw_upper_a(:) = 0._r8 ! total water in layers above where CLUBB is active AFTER CLUBB is called tw_upper_b(:) = 0._r8 ! total water in layers above where CLUBB is active BEFORE CLUBB is called - + do i=1, ncol do k=1, clubbtop(i) - + se_upper_a(i) = se_upper_a(i) + (clubb_s(i,k)+0.5_r8*(um(i,k)**2+vm(i,k)**2) & +(latvap+latice)*(rtm(i,k)-rcm(i,k)) & +(latice)*rcm(i,k))*state1%pdel(i,k)/gravit - + se_upper_b(i) = se_upper_b(i) + (state1%s(i,k)+0.5_r8*(state1%u(i,k)**2+state1%v(i,k)**2) & + (latvap+latice)*state1%q(i,k,ixq) & + (latice)*state1%q(i,k,ixcldliq))*state1%pdel(i,k)/gravit - + tw_upper_a(i) = tw_upper_a(i) + rtm(i,k)*state1%pdel(i,k)/gravit - + tw_upper_b(i) = tw_upper_b(i) + (state1%q(i,k,ixq) & +state1%q(i,k,ixcldliq))*state1%pdel(i,k)/gravit end do end do - + ! Compute the disbalance of total energy and water in upper levels, - ! divide by the thickness in the lower atmosphere where we will + ! divide by the thickness in the lower atmosphere where we will ! evenly distribute this disbalance do i=1, ncol se_upper_diss(i) = (se_upper_a(i) - se_upper_b(i))/(state1%pint(i,pverp)-state1%pint(i,clubbtop(i)+1)) tw_upper_diss(i) = (tw_upper_a(i) - tw_upper_b(i))/(state1%pint(i,pverp)-state1%pint(i,clubbtop(i)+1)) end do - + ! Perform a test to see if there will be any negative RTM errors ! in the column. If so, apply the disbalance to the surface do i=1, ncol @@ -3537,50 +3536,50 @@ subroutine clubb_tend_cam( & end do end if end do - + do i=1, ncol - + if (apply_to_surface(i)) then - + tw_upper_diss(i) = (tw_upper_a(i) - tw_upper_b(i))/(state1%pint(i,pverp)-state1%pint(i,pver)) se_upper_diss(i) = (se_upper_a(i) - se_upper_b(i))/(state1%pint(i,pverp)-state1%pint(i,pver)) rtm(i,pver) = rtm(i,pver) + tw_upper_diss(i)*gravit - + if (apply_to_heat) then clubb_s(i,pver) = clubb_s(i,pver) + se_upper_diss(i)*gravit end if - + else - + ! Apply the disbalances above to layers where CLUBB is active do k=clubbtop(i)+1, pver rtm(i,k) = rtm(i,k) + tw_upper_diss(i)*gravit - + if (apply_to_heat) then clubb_s(i,k) = clubb_s(i,k) + se_upper_diss(i)*gravit end if end do - - end if - + + end if + end do - + ! Essentially "zero" out tendencies in the layers above where CLUBB is active do i=1, ncol do k=1, clubbtop(i) if (apply_to_heat) clubb_s(i,k) = state1%s(i,k) rcm(i,k) = state1%q(i,k,ixcldliq) rtm(i,k) = state1%q(i,k,ixq) + rcm(i,k) - end do + end do end do - + ! Compute integrals for static energy, kinetic energy, water vapor, and liquid water ! after CLUBB is called. This is for energy conservation purposes. se_a(:) = 0._r8 ke_a(:) = 0._r8 wv_a(:) = 0._r8 wl_a(:) = 0._r8 - + do k=1,pver do i=1, ncol se_a(i) = se_a(i) + clubb_s(i,k)*state1%pdel(i,k)/gravit @@ -3588,14 +3587,14 @@ subroutine clubb_tend_cam( & wv_a(i) = wv_a(i) + (rtm(i,k)-rcm(i,k))*state1%pdel(i,k)/gravit wl_a(i) = wl_a(i) + (rcm(i,k))*state1%pdel(i,k)/gravit end do - end do - + end do + ! Do the same as above, but for before CLUBB was called. se_b(:) = 0._r8 ke_b(:) = 0._r8 wv_b(:) = 0._r8 - wl_b(:) = 0._r8 - + wl_b(:) = 0._r8 + do k=1, pver do i=1, ncol se_b(i) = se_b(i) + state1%s(i,k)*state1%pdel(i,k)/gravit @@ -3604,23 +3603,23 @@ subroutine clubb_tend_cam( & wl_b(i) = wl_b(i) + state1%q(i,k,ixcldliq)*state1%pdel(i,k)/gravit end do end do - - + + do i=1, ncol ! Based on these integrals, compute the total energy before and after CLUBB call te_a(i) = se_a(i) + ke_a(i) + (latvap+latice) * wv_a(i) + latice * wl_a(i) te_b(i) = se_b(i) + ke_b(i) + (latvap+latice) * wv_b(i) + latice * wl_b(i) - + ! Take into account the surface fluxes of heat and moisture ! Use correct qflux from cam_in, not lhf/latvap as was done previously - te_b(i) = te_b(i) + (cam_in%shf(i)+cam_in%cflx(i,1)*(latvap+latice)) * hdtime + te_b(i) = te_b(i) + (cam_in%shf(i)+cam_in%cflx(i,1)*(latvap+latice)) * hdtime ! Compute the disbalance of total energy, over depth where CLUBB is active se_dis(i) = (te_a(i) - te_b(i))/(state1%pint(i,pverp)-state1%pint(i,clubbtop(i)+1)) end do ! Fix the total energy coming out of CLUBB so it achieves energy conservation. - ! Apply this fixer throughout the column evenly, but only at layers where + ! Apply this fixer throughout the column evenly, but only at layers where ! CLUBB is active. ! ! NOTE: The energy fixer seems to cause the climate to change significantly @@ -3635,12 +3634,12 @@ subroutine clubb_tend_cam( & se_dis(i) = -1._r8*se_dis(i)*gravit/cpairv(i,pver,lchnk) end do end if - + ! Now compute the tendencies of CLUBB to CAM, note that pverp is the ghost point ! for all variables and therefore is never called in this loop rtm_integral_vtend(:) = 0._r8 rtm_integral_ltend(:) = 0._r8 - + do k=1, pver do i=1, ncol @@ -3655,17 +3654,17 @@ subroutine clubb_tend_cam( & end do end do - - + + if (clubb_do_adv) then if (macmic_it == cld_macmic_num_steps) then - + do k=1, pver do i=1, ncol - ! Here add a constant to moments which can be either positive or + ! Here add a constant to moments which can be either positive or ! negative. This is to prevent clipping when dynamics tries to - ! make all constituents positive + ! make all constituents positive wp3(i,k) = wp3(i,k) + wp3_const rtpthlp(i,k) = rtpthlp(i,k) + rtpthlp_const wpthlp(i,k) = wpthlp(i,k) + wpthlp_const @@ -3674,18 +3673,18 @@ subroutine clubb_tend_cam( & ptend_loc%q(i,k,ixthlp2) = (thlp2(i,k) - state1%q(i,k,ixthlp2)) / hdtime ! THLP Variance ptend_loc%q(i,k,ixrtp2) = (rtp2(i,k) - state1%q(i,k,ixrtp2)) / hdtime ! RTP Variance ptend_loc%q(i,k,ixrtpthlp) = (rtpthlp(i,k) - state1%q(i,k,ixrtpthlp)) / hdtime ! RTP THLP covariance - ptend_loc%q(i,k,ixwpthlp) = (wpthlp(i,k) - state1%q(i,k,ixwpthlp)) / hdtime ! WPTHLP + ptend_loc%q(i,k,ixwpthlp) = (wpthlp(i,k) - state1%q(i,k,ixwpthlp)) / hdtime ! WPTHLP ptend_loc%q(i,k,ixwprtp) = (wprtp(i,k) - state1%q(i,k,ixwprtp)) / hdtime ! WPRTP ptend_loc%q(i,k,ixwp2) = (wp2(i,k) - state1%q(i,k,ixwp2)) / hdtime ! WP2 ptend_loc%q(i,k,ixwp3) = (wp3(i,k) - state1%q(i,k,ixwp3)) / hdtime ! WP3 ptend_loc%q(i,k,ixup2) = (up2(i,k) - state1%q(i,k,ixup2)) / hdtime ! UP2 ptend_loc%q(i,k,ixvp2) = (vp2(i,k) - state1%q(i,k,ixvp2)) / hdtime ! VP2 - + end do end do - + else - + do k=1, pver do i=1, ncol ptend_loc%q(i,k,ixthlp2) = 0._r8 @@ -3696,16 +3695,16 @@ subroutine clubb_tend_cam( & ptend_loc%q(i,k,ixwp2) = 0._r8 ptend_loc%q(i,k,ixwp3) = 0._r8 ptend_loc%q(i,k,ixup2) = 0._r8 - ptend_loc%q(i,k,ixvp2) = 0._r8 + ptend_loc%q(i,k,ixvp2) = 0._r8 end do end do - + end if end if - + ! Apply tendencies to ice mixing ratio, liquid and ice number, and aerosol constituents. - ! Loading up this array doesn't mean the tendencies are applied. + ! Loading up this array doesn't mean the tendencies are applied. ! edsclr_out is compressed with just the constituents being used, ptend and state are not compressed icnt=0 do ixind=1,pcnst @@ -3716,17 +3715,17 @@ subroutine clubb_tend_cam( & (ixind /= ixrtpthlp) .and. (ixind /= ixwpthlp) .and.& (ixind /= ixwprtp) .and. (ixind /= ixwp2) .and.& (ixind /= ixwp3) .and. (ixind /= ixup2) .and. (ixind /= ixvp2) ) then - + do k=1, pver do i=1, ncol - ptend_loc%q(i,k,ixind) = (edsclr_out(i,k,icnt)-state1%q(i,k,ixind))/hdtime ! transported constituents + ptend_loc%q(i,k,ixind) = (edsclr_out(i,k,icnt)-state1%q(i,k,ixind))/hdtime ! transported constituents end do end do - + end if end if end do - + call t_stopf("clubb_tend_cam_i_loop") call outfld('KVH_CLUBB', khzm, pcols, lchnk) @@ -3735,7 +3734,7 @@ subroutine clubb_tend_cam( & call outfld('ELEAK_CLUBB', eleak, pcols, lchnk) call outfld('TFIX_CLUBB', se_dis, pcols, lchnk) - ! Add constant to ghost point so that output is not corrupted + ! Add constant to ghost point so that output is not corrupted if (clubb_do_adv) then if (macmic_it == cld_macmic_num_steps) then wp3(:,pverp) = wp3(:,pverp) + wp3_const @@ -3743,7 +3742,7 @@ subroutine clubb_tend_cam( & wpthlp(:,pverp) = wpthlp(:,pverp) + wpthlp_const wprtp(:,pverp) = wprtp(:,pverp) + wprtp_const end if - end if + end if cmeliq(:,:) = ptend_loc%q(:,:,ixcldliq) @@ -3752,44 +3751,44 @@ subroutine clubb_tend_cam( & ! and compute output, etc ! ! ------------------------------------------------- ! - ! Output CLUBB tendencies + ! Output CLUBB tendencies call outfld( 'RVMTEND_CLUBB', ptend_loc%q(:,:,ixq), pcols, lchnk) call outfld( 'RCMTEND_CLUBB', ptend_loc%q(:,:,ixcldliq), pcols, lchnk) call outfld( 'RIMTEND_CLUBB', ptend_loc%q(:,:,ixcldice), pcols, lchnk) call outfld( 'STEND_CLUBB', ptend_loc%s,pcols, lchnk) call outfld( 'UTEND_CLUBB', ptend_loc%u,pcols, lchnk) - call outfld( 'VTEND_CLUBB', ptend_loc%v,pcols, lchnk) + call outfld( 'VTEND_CLUBB', ptend_loc%v,pcols, lchnk) call outfld( 'CMELIQ', cmeliq, pcols, lchnk) call physics_ptend_sum(ptend_loc,ptend_all,ncol) call physics_update(state1,ptend_loc,hdtime) - - ! Due to the order of operation of CLUBB, which closes on liquid first, - ! then advances it's predictive equations second, this can lead to - ! RHliq > 1 directly before microphysics is called. Therefore, we use - ! ice_macro_tend to enforce RHliq <= 1 everywhere before microphysics is called. - + + ! Due to the order of operation of CLUBB, which closes on liquid first, + ! then advances it's predictive equations second, this can lead to + ! RHliq > 1 directly before microphysics is called. Therefore, we use + ! ice_macro_tend to enforce RHliq <= 1 everywhere before microphysics is called. + if (clubb_do_liqsupersat) then - + ! -------------------------------------- ! ! Ice Saturation Adjustment Computation ! ! -------------------------------------- ! - + latsub = latvap + latice lq2(:) = .FALSE. lq2(ixq) = .TRUE. lq2(ixcldliq) = .TRUE. lq2(ixnumliq) = .TRUE. - + call physics_ptend_init(ptend_loc, state%psetcols, 'iceadj', ls=.true., lq=lq2 ) - + stend(:ncol,:)=0._r8 qvtend(:ncol,:)=0._r8 qctend(:ncol,:)=0._r8 inctend(:ncol,:)=0._r8 - + call liquid_macro_tend(npccn(1:ncol,top_lev:pver), state1%t(1:ncol,top_lev:pver), & state1%pmid(1:ncol,top_lev:pver), state1%q(1:ncol,top_lev:pver,ixq), & state1%q(1:ncol,top_lev:pver,ixcldliq), state1%q(1:ncol,top_lev:pver,ixnumliq), & @@ -3801,13 +3800,13 @@ subroutine clubb_tend_cam( & ptend_loc%q(:ncol,top_lev:pver,ixcldliq)=qctend(:ncol,top_lev:pver) ptend_loc%q(:ncol,top_lev:pver,ixnumliq)=inctend(:ncol,top_lev:pver) ptend_loc%s(:ncol,top_lev:pver)=stend(:ncol,top_lev:pver) - + ! Add the ice tendency to the output tendency call physics_ptend_sum(ptend_loc, ptend_all, ncol) - + ! ptend_loc is reset to zero by this call call physics_update(state1, ptend_loc, hdtime) - + ! Write output for tendencies: ! oufld: QVTENDICE,QCTENDICE,NCTENDICE,FQTENDICE temp2d(:ncol,:pver) = stend(:ncol,:pver)/cpairv(:ncol,:pver,lchnk) @@ -3815,18 +3814,18 @@ subroutine clubb_tend_cam( & call outfld( 'QVTENDICE', qvtend, pcols, lchnk ) call outfld( 'QCTENDICE', qctend, pcols, lchnk ) call outfld( 'NCTENDICE', inctend, pcols, lchnk ) - + where(qctend .ne. 0._r8) fqtend = 1._r8 elsewhere fqtend = 0._r8 end where - + call outfld( 'FQTENDICE', fqtend, pcols, lchnk ) end if - + + ! ------------------------------------------------------------ ! ! ------------------------------------------------------------ ! - ! ------------------------------------------------------------ ! ! ------------------------------------------------------------ ! ! The rest of the code deals with diagnosing variables ! ! for microphysics/radiation computation and macrophysics ! @@ -3834,11 +3833,11 @@ subroutine clubb_tend_cam( & ! ------------------------------------------------------------ ! ! ------------------------------------------------------------ ! - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! COMPUTE THE ICE CLOUD DETRAINMENT ! ! Detrainment of convective condensate into the environment or stratiform cloud ! - ! --------------------------------------------------------------------------------- ! + ! --------------------------------------------------------------------------------- ! ! Initialize the shallow convective detrainment rate, will always be zero dlf2(:,:) = 0.0_r8 @@ -3859,7 +3858,7 @@ subroutine clubb_tend_cam( & call pbuf_get_field(pbuf, dnlfzm_idx, dnlfzm) call pbuf_get_field(pbuf, dnifzm_idx, dnifzm) end if - + do k=1,pver do i=1,ncol if( state1%t(i,k) > meltpt_temp ) then @@ -3867,7 +3866,7 @@ subroutine clubb_tend_cam( & elseif ( state1%t(i,k) < dt_low ) then dum1 = 1.0_r8 else - dum1 = ( meltpt_temp - state1%t(i,k) ) / ( meltpt_temp - dt_low ) + dum1 = ( meltpt_temp - state1%t(i,k) ) / ( meltpt_temp - dt_low ) end if if (zmconv_microp) then @@ -3879,21 +3878,21 @@ subroutine clubb_tend_cam( & ptend_loc%q(i,k,ixnumice) = dnifzm(i,k) + 3._r8 * ( dlf2(i,k) * dum1 ) & / (4._r8*3.14_r8*di_rad**3*500._r8) ! Shallow Convection ptend_loc%s(i,k) = dlf2(i,k) * dum1 * latice - else + else ptend_loc%q(i,k,ixcldliq) = dlf(i,k) * ( 1._r8 - dum1 ) ptend_loc%q(i,k,ixcldice) = dlf(i,k) * dum1 ptend_loc%q(i,k,ixnumliq) = 3._r8 * ( max(0._r8, ( dlf(i,k) - dlf2(i,k) )) * ( 1._r8 - dum1 ) ) & / (4._r8*3.14_r8*dl_rad**3*997._r8) + & ! Deep Convection 3._r8 * ( dlf2(i,k) * ( 1._r8 - dum1 ) ) & - / (4._r8*3.14_r8*10.e-6_r8**3*997._r8) ! Shallow Convection + / (4._r8*3.14_r8*10.e-6_r8**3*997._r8) ! Shallow Convection ptend_loc%q(i,k,ixnumice) = 3._r8 * ( max(0._r8, ( dlf(i,k) - dlf2(i,k) )) * dum1 ) & / (4._r8*3.14_r8*di_rad**3*500._r8) + & ! Deep Convection 3._r8 * ( dlf2(i,k) * dum1 ) & / (4._r8*3.14_r8*50.e-6_r8**3*500._r8) ! Shallow Convection ptend_loc%s(i,k) = dlf(i,k) * dum1 * latice - dlf_liq_out(i,k) = dlf(i,k) * ( 1._r8 - dum1 ) + dlf_liq_out(i,k) = dlf(i,k) * ( 1._r8 - dum1 ) dlf_ice_out(i,k) = dlf(i,k) * dum1 end if @@ -3902,18 +3901,18 @@ subroutine clubb_tend_cam( & ! so that the energy checker doesn't complain. det_s(i) = det_s(i) + ptend_loc%s(i,k)*state1%pdel(i,k)/gravit det_ice(i) = det_ice(i) - ptend_loc%q(i,k,ixcldice)*state1%pdel(i,k)/gravit - + enddo enddo - + det_ice(:ncol) = det_ice(:ncol)/1000._r8 ! divide by density of water call outfld( 'DPDLFLIQ', ptend_loc%q(:,:,ixcldliq), pcols, lchnk) call outfld( 'DPDLFICE', ptend_loc%q(:,:,ixcldice), pcols, lchnk) - + temp2d(:ncol,:pver) = ptend_loc%s(:ncol,:pver)/cpairv(:ncol,:pver, lchnk) call outfld( 'DPDLFT', temp2d, pcols, lchnk) - + call outfld( 'DETNLIQTND', ptend_loc%q(:,:,ixnumliq),pcols, lchnk ) call physics_ptend_sum(ptend_loc,ptend_all,ncol) @@ -3940,20 +3939,20 @@ subroutine clubb_tend_cam( & else relvarmax = 10.0_r8 end if - + relvar(:,:) = relvarmax ! default - if (deep_scheme .ne. 'CLUBB_SGS') then + if (deep_scheme .ne. 'CLUBB_SGS') then where (rcm(:ncol,:pver) /= 0 .and. qclvar(:ncol,:pver) /= 0) & relvar(:ncol,:pver) = min(relvarmax,max(0.001_r8,rcm(:ncol,:pver)**2/qclvar(:ncol,:pver))) end if - + ! ------------------------------------------------- ! ! Optional Accretion enhancement factor ! - ! ------------------------------------------------- ! + ! ------------------------------------------------- ! accre_enhan(:ncol,:pver) = 1._r8 - + ! ------------------------------------------------- ! ! Diagnose some output variables ! ! ------------------------------------------------- ! @@ -3975,54 +3974,54 @@ subroutine clubb_tend_cam( & end if enddo enddo - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! Diagnose some quantities that are computed in macrop_tend here. ! ! These are inputs required for the microphysics calculation. ! ! ! ! FIRST PART COMPUTES THE STRATIFORM CLOUD FRACTION FROM CLUBB CLOUD FRACTION ! - ! --------------------------------------------------------------------------------- ! - - ! initialize variables + ! --------------------------------------------------------------------------------- ! + + ! initialize variables alst(:,:) = 0.0_r8 - qlst(:,:) = 0.0_r8 - + qlst(:,:) = 0.0_r8 + do k=1,pver do i=1,ncol - alst(i,k) = cloud_frac(i,k) + alst(i,k) = cloud_frac(i,k) qlst(i,k) = rcm(i,k)/max(0.01_r8,alst(i,k)) ! Incloud stratus condensate mixing ratio enddo enddo - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! THIS PART COMPUTES CONVECTIVE AND DEEP CONVECTIVE CLOUD FRACTION ! - ! --------------------------------------------------------------------------------- ! - + ! --------------------------------------------------------------------------------- ! + ! Initialize cloud fraction deepcu(:,:) = 0.0_r8 shalcu(:,:) = 0.0_r8 - + do k=1,pver-1 do i=1,ncol - ! diagnose the deep convective cloud fraction, as done in macrophysics based on the - ! deep convective mass flux, read in from pbuf. Since shallow convection is never + ! diagnose the deep convective cloud fraction, as done in macrophysics based on the + ! deep convective mass flux, read in from pbuf. Since shallow convection is never ! called, the shallow convective mass flux will ALWAYS be zero, ensuring that this cloud - ! fraction is purely from deep convection scheme. + ! fraction is purely from deep convection scheme. deepcu(i,k) = max(0.0_r8,min(dp1*log(1.0_r8+dp2*(cmfmc(i,k+1)-cmfmc_sh(i,k+1))),0.6_r8)) shalcu(i,k) = 0._r8 - + if (deepcu(i,k) <= frac_limit .or. dp_icwmr(i,k) < ic_limit) then deepcu(i,k) = 0._r8 end if - - ! using the deep convective cloud fraction, and CLUBB cloud fraction (variable + + ! using the deep convective cloud fraction, and CLUBB cloud fraction (variable ! "cloud_frac"), compute the convective cloud fraction. This follows the formulation - ! found in macrophysics code. Assumes that convective cloud is all nonstratiform cloud + ! found in macrophysics code. Assumes that convective cloud is all nonstratiform cloud ! from CLUBB plus the deep convective cloud fraction concld(i,k) = min(cloud_frac(i,k)-alst(i,k)+deepcu(i,k),0.80_r8) enddo enddo - + if (single_column) then if (trim(scm_clubb_iop_name) == 'ATEX_48hr' .or. & trim(scm_clubb_iop_name) == 'BOMEX_5day' .or. & @@ -4030,20 +4029,20 @@ subroutine clubb_tend_cam( & trim(scm_clubb_iop_name) == 'DYCOMSrf02_06hr' .or. & trim(scm_clubb_iop_name) == 'RICO_3day' .or. & trim(scm_clubb_iop_name) == 'ARM_CC') then - + deepcu(:,:) = 0.0_r8 concld(:,:) = 0.0_r8 - + end if end if - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! COMPUTE THE ICE CLOUD FRACTION PORTION ! ! use the aist_vector function to compute the ice cloud fraction ! - ! --------------------------------------------------------------------------------- ! + ! --------------------------------------------------------------------------------- ! aist(:,:top_lev-1) = 0._r8 - qsatfac(:, :) = 0._r8 ! Zero out entire profile in case qsatfac is left undefined in aist_vector below + qsatfac(:, :) = 0._r8 ! Zero out entire profile in case qsatfac is left undefined in aist_vector below do k = top_lev, pver @@ -4072,41 +4071,41 @@ subroutine clubb_tend_cam( & qsatfac_out=qsatfac(:,k), rhmini_in=rhmini, rhmaxi_in=rhmaxi) end if enddo - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! THIS PART COMPUTES THE LIQUID STRATUS FRACTION ! ! ! ! For now leave the computation of ice stratus fraction from macrop_driver intact ! - ! because CLUBB does nothing with ice. Here I simply overwrite the liquid stratus ! + ! because CLUBB does nothing with ice. Here I simply overwrite the liquid stratus ! ! fraction that was coded in macrop_driver ! - ! --------------------------------------------------------------------------------- ! - + ! --------------------------------------------------------------------------------- ! + ! Recompute net stratus fraction using maximum over-lapping assumption, as done ! in macrophysics code, using alst computed above and aist read in from physics buffer - + do k=1,pver do i=1,ncol ast(i,k) = max(alst(i,k),aist(i,k)) - qist(i,k) = state1%q(i,k,ixcldice)/max(0.01_r8,aist(i,k)) + qist(i,k) = state1%q(i,k,ixcldice)/max(0.01_r8,aist(i,k)) enddo enddo - - ! Probably need to add deepcu cloud fraction to the cloud fraction array, else would just - ! be outputting the shallow convective cloud fraction + + ! Probably need to add deepcu cloud fraction to the cloud fraction array, else would just + ! be outputting the shallow convective cloud fraction do k=1,pver do i=1,ncol cloud_frac(i,k) = min(ast(i,k)+deepcu(i,k),1.0_r8) enddo enddo - - ! --------------------------------------------------------------------------------- ! + + ! --------------------------------------------------------------------------------- ! ! DIAGNOSE THE PBL DEPTH ! ! this is needed for aerosol code ! - ! --------------------------------------------------------------------------------- ! - + ! --------------------------------------------------------------------------------- ! + do i=1,ncol do k=1,pver !use local exner since state%exner is not a proper exner @@ -4115,17 +4114,17 @@ subroutine clubb_tend_cam( & thv(i,k) = th(i,k)*(1.0_r8+zvir*state1%q(i,k,ixq) - state1%q(i,k,ixcldliq)) enddo enddo - + call calc_ustar( ncol, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver), cam_in%wsx(1:ncol), cam_in%wsy(1:ncol), & rrho(1:ncol), ustar2(1:ncol)) ! use correct qflux from coupler call calc_obklen( ncol, th(1:ncol,pver), thv(1:ncol,pver), cam_in%cflx(1:ncol,1), cam_in%shf(1:ncol), & rrho(1:ncol), ustar2(1:ncol), kinheat(1:ncol), kinwat(1:ncol), kbfs(1:ncol), & obklen(1:ncol)) - + dummy2(:) = 0._r8 dummy3(:) = 0._r8 - + where (kbfs(:ncol) == -0.0_r8) kbfs(:ncol) = 0.0_r8 ! Compute PBL depth according to Holtslag-Boville Scheme @@ -4135,14 +4134,14 @@ subroutine clubb_tend_cam( & ! Output the PBL depth call outfld('PBLH', pblh, pcols, lchnk) - + ! Assign the first pver levels of cloud_frac back to cld cld(:,1:pver) = cloud_frac(:,1:pver) - ! --------------------------------------------------------------------------------- ! + ! --------------------------------------------------------------------------------- ! ! END CLOUD FRACTION DIAGNOSIS, begin to store variables back into buffer ! - ! --------------------------------------------------------------------------------- ! - + ! --------------------------------------------------------------------------------- ! + ! Output calls of variables goes here call outfld( 'RELVAR', relvar, pcols, lchnk ) call outfld( 'RHO_CLUBB', rho(:,1:pver), pcols, lchnk ) @@ -4214,60 +4213,60 @@ subroutine clubb_tend_cam( & end if ! Output CLUBB history here - if (l_stats) then - + if (l_stats) then + do j=1,stats_zt(1)%num_output_fields - + temp1 = trim(stats_zt(1)%file%grid_avg_var(j)%name) sub = temp1 if (len(temp1) > 16) sub = temp1(1:16) - + call outfld(trim(sub), out_zt(:,:,j), pcols, lchnk ) enddo - + do j=1,stats_zm(1)%num_output_fields - + temp1 = trim(stats_zm(1)%file%grid_avg_var(j)%name) sub = temp1 if (len(temp1) > 16) sub = temp1(1:16) - + call outfld(trim(sub),out_zm(:,:,j), pcols, lchnk) enddo - if (l_output_rad_files) then + if (l_output_rad_files) then do j=1,stats_rad_zt(1)%num_output_fields call outfld(trim(stats_rad_zt(1)%file%grid_avg_var(j)%name), out_radzt(:,:,j), pcols, lchnk) enddo - + do j=1,stats_rad_zm(1)%num_output_fields call outfld(trim(stats_rad_zm(1)%file%grid_avg_var(j)%name), out_radzm(:,:,j), pcols, lchnk) enddo end if - + do j=1,stats_sfc(1)%num_output_fields call outfld(trim(stats_sfc(1)%file%grid_avg_var(j)%name), out_sfc(:,:,j), pcols, lchnk) enddo - + end if - + call t_stopf("clubb_tend_cam") - + return #endif end subroutine clubb_tend_cam - + ! =============================================================================== ! ! ! ! =============================================================================== ! ! Saturation adjustment for ice ! Add ice mass if supersaturated -subroutine ice_macro_tend(naai,t,p,qv,qi,ni,xxls,deltat,stend,qvtend,qitend,nitend,vlen) +subroutine ice_macro_tend(naai,t,p,qv,qi,ni,xxls,deltat,stend,qvtend,qitend,nitend,vlen) use wv_sat_methods, only: wv_sat_qsat_ice integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: naai !Activated number of ice nuclei + real(r8), dimension(vlen), intent(in) :: naai !Activated number of ice nuclei real(r8), dimension(vlen), intent(in) :: t !temperature (k) real(r8), dimension(vlen), intent(in) :: p !pressure (pa) real(r8), dimension(vlen), intent(in) :: qv !water vapor mixing ratio @@ -4275,11 +4274,11 @@ subroutine ice_macro_tend(naai,t,p,qv,qi,ni,xxls,deltat,stend,qvtend,qitend,nite real(r8), dimension(vlen), intent(in) :: ni !ice number concentration real(r8), intent(in) :: xxls !latent heat of freezing real(r8), intent(in) :: deltat !timestep - real(r8), dimension(vlen), intent(out) :: stend ! 'temperature' tendency + real(r8), dimension(vlen), intent(out) :: stend ! 'temperature' tendency real(r8), dimension(vlen), intent(out) :: qvtend !vapor tendency real(r8), dimension(vlen), intent(out) :: qitend !ice mass tendency - real(r8), dimension(vlen), intent(out) :: nitend !ice number tendency - + real(r8), dimension(vlen), intent(out) :: nitend !ice number tendency + real(r8) :: ESI(vlen) real(r8) :: QSI(vlen) integer :: i @@ -4302,7 +4301,7 @@ subroutine ice_macro_tend(naai,t,p,qv,qi,ni,xxls,deltat,stend,qvtend,qitend,nite qitend(i) = (qv(i)-QSI(i))/deltat qvtend(i) = 0._r8 - qitend(i) stend(i) = qitend(i) * xxls ! moist static energy tend...[J/kg/s] ! - + ! if ice exists (more than 1 L-1) and there is condensation, do not add to number (= growth), else, add 10um ice if (ni(i) < 1.e3_r8 .and. (qi(i)+qitend(i)*deltat) > 1.e-18_r8) then nitend(i) = nitend(i) + 3._r8 * qitend(i)/(4._r8*3.14_r8* 10.e-6_r8**3*997._r8) @@ -4336,7 +4335,7 @@ end subroutine ice_macro_tend ! Code writen March, 1999 by Bjorn Stevens ! -real(r8) function diag_ustar( z, bflx, wnd, z0 ) +real(r8) function diag_ustar( z, bflx, wnd, z0 ) use shr_const_mod, only : shr_const_karman, shr_const_pi, shr_const_g @@ -4405,59 +4404,59 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Description: Initializes the statistics saving functionality of ! the CLUBB model. This is for purpose of CAM-CLUBB interface. Here ! the traditional stats_init of CLUBB is not called, as it is not compatible - ! with CAM output. - + ! with CAM output. + !----------------------------------------------------------------------- use clubb_api_module, only: & - ztscr01, & - ztscr02, & - ztscr03, & - ztscr04, & - ztscr05, & - ztscr06, & - ztscr07, & - ztscr08, & - ztscr09, & - ztscr10, & - ztscr11, & - ztscr12, & - ztscr13, & - ztscr14, & - ztscr15, & - ztscr16, & - ztscr17, & - ztscr18, & - ztscr19, & - ztscr20, & + ztscr01, & + ztscr02, & + ztscr03, & + ztscr04, & + ztscr05, & + ztscr06, & + ztscr07, & + ztscr08, & + ztscr09, & + ztscr10, & + ztscr11, & + ztscr12, & + ztscr13, & + ztscr14, & + ztscr15, & + ztscr16, & + ztscr17, & + ztscr18, & + ztscr19, & + ztscr20, & ztscr21 use clubb_api_module, only: & - zmscr01, & - zmscr02, & - zmscr03, & - zmscr04, & - zmscr05, & - zmscr06, & - zmscr07, & - zmscr08, & - zmscr09, & - zmscr10, & - zmscr11, & - zmscr12, & - zmscr13, & - zmscr14, & + zmscr01, & + zmscr02, & + zmscr03, & + zmscr04, & + zmscr05, & + zmscr06, & + zmscr07, & + zmscr08, & + zmscr09, & + zmscr10, & + zmscr11, & + zmscr12, & + zmscr13, & + zmscr14, & zmscr15, & zmscr16, & zmscr17, & l_stats, & - l_output_rad_files, & - stats_tsamp, & - stats_tout, & - l_stats_samp, & - l_stats_last, & - l_netcdf, & + l_output_rad_files, & + stats_tsamp, & + stats_tout, & + l_stats_samp, & + l_stats_last, & + l_netcdf, & l_grads use clubb_api_module, only: time_precision, & ! @@ -4479,16 +4478,16 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & logical, intent(in) :: l_stats_in ! Stats on? T/F - real(kind=time_precision), intent(in) :: & + real(kind=time_precision), intent(in) :: & stats_tsamp_in, & ! Sampling interval [s] stats_tout_in ! Output interval [s] integer, intent(in) :: nnzp ! Grid points in the vertical [count] - integer, intent(in) :: nnrad_zt ! Grid points in the radiation grid [count] + integer, intent(in) :: nnrad_zt ! Grid points in the radiation grid [count] integer, intent(in) :: nnrad_zm ! Grid points in the radiation grid [count] real(kind=time_precision), intent(in) :: delt ! Timestep (dtmain in CLUBB) [s] - + ! Output Variables type (stats), intent(out) :: stats_zt, & ! stats_zt grid stats_zm, & ! stats_zm grid @@ -4509,11 +4508,11 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & character(len=var_length), dimension(nvarmax_rad_zm) :: clubb_vars_rad_zm ! Variables on the radiation levels character(len=var_length), dimension(nvarmax_sfc) :: clubb_vars_sfc ! Variables at the model surface - namelist /clubb_stats_nl/ & - clubb_vars_zt, & + namelist /clubb_stats_nl/ & + clubb_vars_zt, & clubb_vars_zm, & clubb_vars_rad_zt, & - clubb_vars_rad_zm, & + clubb_vars_rad_zm, & clubb_vars_sfc ! Local Variables @@ -4531,7 +4530,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Set stats_variables variables with inputs from calling subroutine l_stats = l_stats_in - + stats_tsamp = stats_tsamp_in stats_tout = stats_tout_in @@ -4549,7 +4548,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & clubb_vars_rad_zm = '' clubb_vars_sfc = '' - ! Read variables to compute from the namelist + ! Read variables to compute from the namelist if (masterproc) then iunit= getunit() open(unit=iunit,file="atm_in",status='old') @@ -4596,8 +4595,8 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize zt (mass points) i = 1 - do while ( ichar(clubb_vars_zt(i)(1:1)) /= 0 .and. & - len_trim(clubb_vars_zt(i)) /= 0 .and. & + do while ( ichar(clubb_vars_zt(i)(1:1)) /= 0 .and. & + len_trim(clubb_vars_zt(i)) /= 0 .and. & i <= nvarmax_zt ) i = i + 1 enddo @@ -4619,7 +4618,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & allocate( stats_zt%accum_field_values( 1, 1, stats_zt%kk, stats_zt%num_output_fields ) ) allocate( stats_zt%accum_num_samples( 1, 1, stats_zt%kk, stats_zt%num_output_fields ) ) allocate( stats_zt%l_in_update( 1, 1, stats_zt%kk, stats_zt%num_output_fields ) ) - + call stats_zero( stats_zt%kk, stats_zt%num_output_fields, stats_zt%accum_field_values, & stats_zt%accum_num_samples, stats_zt%l_in_update ) @@ -4682,8 +4681,8 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize zm (momentum points) i = 1 - do while ( ichar(clubb_vars_zm(i)(1:1)) /= 0 .and. & - len_trim(clubb_vars_zm(i)) /= 0 .and. & + do while ( ichar(clubb_vars_zm(i)(1:1)) /= 0 .and. & + len_trim(clubb_vars_zm(i)) /= 0 .and. & i <= nvarmax_zm ) i = i + 1 end do @@ -4758,10 +4757,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize rad_zt (radiation points) if (l_output_rad_files) then - + i = 1 - do while ( ichar(clubb_vars_rad_zt(i)(1:1)) /= 0 .and. & - len_trim(clubb_vars_rad_zt(i)) /= 0 .and. & + do while ( ichar(clubb_vars_rad_zt(i)(1:1)) /= 0 .and. & + len_trim(clubb_vars_rad_zt(i)) /= 0 .and. & i <= nvarmax_rad_zt ) i = i + 1 end do @@ -4794,10 +4793,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & stats_rad_zt ) ! Initialize rad_zm (radiation points) - + i = 1 - do while ( ichar(clubb_vars_rad_zm(i)(1:1)) /= 0 .and. & - len_trim(clubb_vars_rad_zm(i)) /= 0 .and. & + do while ( ichar(clubb_vars_rad_zm(i)(1:1)) /= 0 .and. & + len_trim(clubb_vars_rad_zm(i)) /= 0 .and. & i <= nvarmax_rad_zm ) i = i + 1 end do @@ -4825,7 +4824,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & allocate( stats_rad_zm%file%grid_avg_var( stats_rad_zm%num_output_fields ) ) allocate( stats_rad_zm%file%z( stats_rad_zm%kk ) ) - + call stats_init_rad_zm_api( clubb_vars_rad_zm, l_error, & stats_rad_zm ) end if ! l_output_rad_files @@ -4834,8 +4833,8 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize sfc (surface point) i = 1 - do while ( ichar(clubb_vars_sfc(i)(1:1)) /= 0 .and. & - len_trim(clubb_vars_sfc(i)) /= 0 .and. & + do while ( ichar(clubb_vars_sfc(i)(1:1)) /= 0 .and. & + len_trim(clubb_vars_sfc(i)) /= 0 .and. & i <= nvarmax_sfc ) i = i + 1 end do @@ -4877,58 +4876,58 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Now call add fields if (first_call) then - + do i = 1, stats_zt%num_output_fields - + temp1 = trim(stats_zt%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > 16) sub = temp1(1:16) - + call addfld(trim(sub),(/ 'ilev' /),& 'A',trim(stats_zt%file%grid_avg_var(i)%units),trim(stats_zt%file%grid_avg_var(i)%description)) enddo - + do i = 1, stats_zm%num_output_fields - + temp1 = trim(stats_zm%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > 16) sub = temp1(1:16) - + call addfld(trim(sub),(/ 'ilev' /),& 'A',trim(stats_zm%file%grid_avg_var(i)%units),trim(stats_zm%file%grid_avg_var(i)%description)) enddo - if (l_output_rad_files) then + if (l_output_rad_files) then do i = 1, stats_rad_zt%num_output_fields call addfld(trim(stats_rad_zt%file%grid_avg_var(i)%name),(/ 'ilev' /),& 'A',trim(stats_rad_zt%file%grid_avg_var(i)%units),trim(stats_rad_zt%file%grid_avg_var(i)%description)) enddo - + do i = 1, stats_rad_zm%num_output_fields call addfld(trim(stats_rad_zm%file%grid_avg_var(i)%name),(/ 'ilev' /),& 'A',trim(stats_rad_zm%file%grid_avg_var(i)%units),trim(stats_rad_zm%file%grid_avg_var(i)%description)) enddo end if - + do i = 1, stats_sfc%num_output_fields call addfld(trim(stats_sfc%file%grid_avg_var(i)%name),horiz_only,& 'A',trim(stats_sfc%file%grid_avg_var(i)%units),trim(stats_sfc%file%grid_avg_var(i)%description)) enddo - + end if return - end subroutine stats_init_clubb - + end subroutine stats_init_clubb + #endif ! =============================================================================== ! ! ! ! =============================================================================== ! -#ifdef CLUBB_SGS +#ifdef CLUBB_SGS subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, stats_rad_zm, stats_sfc, & out_zt, out_zm, out_radzt, out_radzm, out_sfc) !----------------------------------------------------------------------- @@ -4943,8 +4942,8 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st use clubb_api_module, only: & fstderr, & ! Constant(s) - l_stats_last, & - stats_tsamp, & + l_stats_last, & + stats_tsamp, & stats_tout, & l_output_rad_files, & clubb_at_least_debug_level_api ! Procedure(s) @@ -4954,14 +4953,14 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st implicit none integer :: thecol - + ! Input Variables type (stats), intent(inout) :: stats_zt, & ! stats_zt grid stats_zm, & ! stats_zm grid stats_rad_zt, & ! stats_rad_zt grid stats_rad_zm, & ! stats_rad_zm grid stats_sfc ! stats_sfc - + ! Inout variables real(r8), intent(inout) :: out_zt(:,:,:) ! (pcols,pverp,stats_zt%num_output_fields) real(r8), intent(inout) :: out_zm(:,:,:) ! (pcols,pverp,stats_zt%num_output_fields) @@ -4992,36 +4991,36 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st end if call stats_avg( stats_sfc%kk, stats_sfc%num_output_fields, stats_sfc%accum_field_values, stats_sfc%accum_num_samples ) - ! Here we are not outputting the data, rather reading the stats into + ! Here we are not outputting the data, rather reading the stats into ! arrays which are conformable to CAM output. Also, the data is "flipped" - ! in the vertical level to be the same as CAM output. + ! in the vertical level to be the same as CAM output. do i = 1, stats_zt%num_output_fields - do k = 1, stats_zt%kk + do k = 1, stats_zt%kk out_zt(thecol,pverp-k+1,i) = stats_zt%accum_field_values(1,1,k,i) if(is_nan(out_zt(thecol,k,i))) out_zt(thecol,k,i) = 0.0_r8 - enddo + enddo enddo do i = 1, stats_zm%num_output_fields - do k = 1, stats_zt%kk + do k = 1, stats_zt%kk out_zm(thecol,pverp-k+1,i) = stats_zm%accum_field_values(1,1,k,i) if(is_nan(out_zm(thecol,k,i))) out_zm(thecol,k,i) = 0.0_r8 - enddo + enddo enddo - if (l_output_rad_files) then + if (l_output_rad_files) then do i = 1, stats_rad_zt%num_output_fields - do k = 1, stats_rad_zt%kk + do k = 1, stats_rad_zt%kk out_radzt(thecol,pverp-k+1,i) = stats_rad_zt%accum_field_values(1,1,k,i) if(is_nan(out_radzt(thecol,k,i))) out_radzt(thecol,k,i) = 0.0_r8 - enddo + enddo enddo - + do i = 1, stats_rad_zm%num_output_fields - do k = 1, stats_rad_zm%kk + do k = 1, stats_rad_zm%kk out_radzm(thecol,pverp-k+1,i) = stats_rad_zm%accum_field_values(1,1,k,i) if(is_nan(out_radzm(thecol,k,i))) out_radzm(thecol,k,i) = 0.0_r8 - enddo + enddo enddo ! Fill in values above the CLUBB top. @@ -5031,9 +5030,9 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st out_radzm(thecol,:top_lev-1,:) = 0.0_r8 end if ! l_output_rad_files - + do i = 1, stats_sfc%num_output_fields - out_sfc(thecol,1,i) = stats_sfc%accum_field_values(1,1,1,i) + out_sfc(thecol,1,i) = stats_sfc%accum_field_values(1,1,1,i) if(is_nan(out_sfc(thecol,1,i))) out_sfc(thecol,1,i) = 0.0_r8 enddo @@ -5054,14 +5053,14 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st return end subroutine stats_end_timestep_clubb -#endif - +#endif + ! =============================================================================== ! ! ! ! =============================================================================== ! #ifdef CLUBB_SGS - + !----------------------------------------------------------------------- subroutine stats_zero( kk, num_output_fields, x, n, l_in_update ) @@ -5095,14 +5094,14 @@ subroutine stats_zero( kk, num_output_fields, x, n, l_in_update ) return end subroutine stats_zero - + #endif ! =============================================================================== ! ! ! ! =============================================================================== ! - + #ifdef CLUBB_SGS !----------------------------------------------------------------------- subroutine stats_avg( kk, num_output_fields, x, n ) @@ -5150,7 +5149,7 @@ subroutine grid_size(state, grid_dx, grid_dy) use shr_const_mod, only: shr_const_pi use physics_types, only: physics_state - + type(physics_state), intent(in) :: state real(r8), intent(out) :: grid_dx(pcols), grid_dy(pcols) ! CAM grid [m] @@ -5165,17 +5164,17 @@ subroutine grid_size(state, grid_dx, grid_dy) do i=1,state%ncol column_area = get_area_p(state%lchnk,i) degree = sqrt(column_area)*(180._r8/shr_const_pi) - + ! Now find meters per degree latitude ! Below equation finds distance between two points on an ellipsoid, derived from expansion - ! taking into account ellipsoid using World Geodetic System (WGS84) reference + ! taking into account ellipsoid using World Geodetic System (WGS84) reference mpdeglat = earth_ellipsoid1 - earth_ellipsoid2 * cos(2._r8*state%lat(i)) + earth_ellipsoid3 * cos(4._r8*state%lat(i)) grid_dx(i) = mpdeglat * degree grid_dy(i) = grid_dx(i) ! Assume these are the same - enddo + enddo - end subroutine grid_size + end subroutine grid_size #endif - + end module clubb_intr diff --git a/src/physics/cam/iop_forcing.F90 b/src/physics/cam/iop_forcing.F90 index 55259685b5..7f309caabc 100644 --- a/src/physics/cam/iop_forcing.F90 +++ b/src/physics/cam/iop_forcing.F90 @@ -29,6 +29,8 @@ subroutine scam_use_iop_srf( cam_in ) use physconst, only: stebol, latvap use scamMod use cam_abortutils, only: endrun + use cam_logfile, only: iulog + use spmd_utils, only: masterproc implicit none save @@ -37,6 +39,17 @@ subroutine scam_use_iop_srf( cam_in ) integer :: c ! Chunk index integer :: ncol ! Number of columns + + if (masterproc) write(iulog,*) " Parameters in iop_forcing :" + if (masterproc) write(iulog,*) " scm_iop_lhflxshflxTg =", scm_iop_lhflxshflxTg + if (masterproc) write(iulog,*) " scm_iop_Tg =", scm_iop_Tg + if (masterproc) write(iulog,*) " scm_crm_mode =", scm_crm_mode + if (masterproc) write(iulog,*) " have_lhflx =", have_lhflx + if (masterproc) write(iulog,*) " have_shflx =", have_shflx + if (masterproc) write(iulog,*) " have_Tg =", have_Tg + if (masterproc) write(iulog,*) " Tground =", tground + + if( scm_iop_lhflxshflxTg .and. scm_iop_Tg ) then call endrun( 'scam_use_iop_srf : scm_iop_lhflxshflxTg and scm_iop_Tg must not be specified at the same time.') end if diff --git a/src/physics/cam/ref_pres.F90 b/src/physics/cam/ref_pres.F90 index 742652db11..c80d3cf866 100644 --- a/src/physics/cam/ref_pres.F90 +++ b/src/physics/cam/ref_pres.F90 @@ -1,18 +1,19 @@ module ref_pres !-------------------------------------------------------------------------- -! +! ! Provides access to reference pressures for use by the physics ! parameterizations. The pressures are provided by the dynamical core ! since it determines the grid used by the physics. -! +! ! Note that the init method for this module is called before the init ! method in physpkg; therefore, most physics modules can use these ! reference pressures during their init phases. -! +! !-------------------------------------------------------------------------- use shr_kind_mod, only: r8=>shr_kind_r8 use ppgrid, only: pver, pverp +use scamMod, only: single_column implicit none public @@ -131,7 +132,11 @@ subroutine ref_pres_init(pref_edge_in, pref_mid_in, num_pr_lev_in) top=.true.) ! Find level corresponding to the molecular diffusion bottom. - do_molec_diff = (ptop_ref < do_molec_press) + if (single_column) then + do_molec_diff = .false. + else + do_molec_diff = (ptop_ref < do_molec_press) + end if if (do_molec_diff) then nbot_molec = press_lim_idx(molec_diff_bot_press, & top=.false.)