From 812d47682c20ba07a484fd883cdbd0894cf0b8df Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 17 May 2024 14:49:08 -0600 Subject: [PATCH 01/18] Remove unused variables --- src/physics/cam/zm_conv_intr.F90 | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index febf576443..337c722b43 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -58,8 +58,6 @@ module zm_conv_intr dp_cldice_idx, & dlfzm_idx, & ! detrained convective cloud water mixing ratio. difzm_idx, & ! detrained convective cloud ice mixing ratio. - dnlfzm_idx, & ! detrained convective cloud water num concen. - dnifzm_idx, & ! detrained convective cloud ice num concen. prec_dp_idx, & snow_dp_idx, & mconzm_idx ! convective mass flux @@ -445,8 +443,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8), pointer, dimension(:,:) :: dp_cldice real(r8), pointer :: dlf(:,:) ! detrained convective cloud water mixing ratio. real(r8), pointer :: dif(:,:) ! detrained convective cloud ice mixing ratio. - real(r8), pointer :: dnlf(:,:) ! detrained convective cloud water num concen. - real(r8), pointer :: dnif(:,:) ! detrained convective cloud ice num concen. real(r8), pointer :: lambdadpcu(:,:) ! slope of cloud liquid size distr real(r8), pointer :: mudpcu(:,:) ! width parameter of droplet size distr real(r8), pointer :: mconzm(:,:) !convective mass fluxes @@ -543,8 +539,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call pbuf_get_field(pbuf, difzm_idx, dif) call pbuf_get_field(pbuf, mconzm_idx, mconzm) - allocate(dnlf(pcols,pver), dnif(pcols,pver)) - ! ! Begin with Zhang-McFarlane (1996) convection parameterization ! @@ -567,8 +561,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & zdu(:,:) = 0._r8 rprd(:,:) = 0._r8 dif(:,:) = 0._r8 - dnlf(:,:) = 0._r8 - dnif(:,:) = 0._r8 mu(:,:) = 0._r8 eu(:,:) = 0._r8 du(:,:) = 0._r8 @@ -596,8 +588,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & dp(:ncol,:), dsubcld(:ncol), jt(:ncol), maxg(:ncol), ideep(:ncol), & ql(:ncol,:), rliq(:ncol), landfrac(:ncol), & org_ncol(:ncol,:), orgt_ncol(:ncol,:), zm_org2d_ncol(:ncol,:), & - dif(:ncol,:), dnlf(:ncol,:), dnif(:ncol,:), & - rice(:ncol), errmsg, errflg) + dif(:ncol,:), rice(:ncol), errmsg, errflg) if (zmconv_org) then @@ -765,7 +756,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call t_startf ('zm_conv_momtran_run') call zm_conv_momtran_run (ncol, pver, pverp, & - l_windt,state1%u(:ncol,:), state1%v(:ncol,:), 2, mu(:ncol,:), md(:ncol,:), & + l_windt,state1%u(:ncol,:), state1%v(:ncol,:), mu(:ncol,:), md(:ncol,:), & zmconv_momcu, zmconv_momcd, & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & @@ -845,7 +836,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & deallocate(zm_org2d) end if - deallocate(dnlf, dnif) end subroutine zm_conv_tend !========================================================================================= From 600b8b69d0098e2fc990de06c3ffea9fbf0f3c97 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 22 May 2024 16:10:38 -0600 Subject: [PATCH 02/18] Remove DIFZM as it is not used --- src/physics/cam/clubb_intr.F90 | 86 +++++++++++++++---------------- src/physics/cam/macrop_driver.F90 | 2 - src/physics/cam/zm_conv_intr.F90 | 7 +-- 3 files changed, 43 insertions(+), 52 deletions(-) diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index c5bdcd71ce..272f8caf02 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -186,7 +186,7 @@ module clubb_intr real(r8) :: clubb_bv_efold = unset_r8 real(r8) :: clubb_wpxp_Ri_exp = unset_r8 real(r8) :: clubb_z_displace = unset_r8 - + integer :: & clubb_iiPDF_type, & ! Selected option for the two-component normal ! (double Gaussian) PDF type to use for the w, rt, @@ -314,7 +314,7 @@ module clubb_intr clubb_l_mono_flux_lim_um, & ! Flag to turn on monotonic flux limiter for um clubb_l_mono_flux_lim_vm, & ! Flag to turn on monotonic flux limiter for vm clubb_l_mono_flux_lim_spikefix, & ! Flag to implement monotonic flux limiter code that - ! eliminates spurious drying tendencies at model top + ! eliminates spurious drying tendencies at model top clubb_l_intr_sfc_flux_smooth = .false. ! Add a locally calculated roughness to upwp and vpwp sfc fluxes ! Constant parameters @@ -433,7 +433,6 @@ module clubb_intr integer :: & dlfzm_idx = -1, & ! ZM detrained convective cloud water mixing ratio. - difzm_idx = -1, & ! ZM detrained convective cloud ice mixing ratio. dnlfzm_idx = -1, & ! ZM detrained convective cloud water num concen. dnifzm_idx = -1 ! ZM detrained convective cloud ice num concen. @@ -477,7 +476,7 @@ 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, dtype_i4, dyn_time_lvls use subcol_utils, only: subcol_get_scheme @@ -844,7 +843,7 @@ subroutine clubb_readnl(nlfile) !----- 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 stats_metadata%l_stats = .false. ! Initialize to false stats_metadata%l_output_rad_files = .false. ! Initialize to false @@ -1201,7 +1200,7 @@ subroutine clubb_readnl(nlfile) ! Overwrite defaults if they are true if (clubb_history) stats_metadata%l_stats = .true. - if (clubb_rad_history) stats_metadata%l_output_rad_files = .true. + if (clubb_rad_history) stats_metadata%l_output_rad_files = .true. if (clubb_cloudtop_cooling) do_cldcool = .true. if (clubb_rainevap_turb) do_rainturb = .true. @@ -1529,7 +1528,7 @@ subroutine clubb_ini_cam(pbuf2d) stats_metadata%l_stats_samp = .false. stats_metadata%l_grads = .false. - ! Overwrite defaults if needbe + ! Overwrite defaults if needbe if (stats_metadata%l_stats) stats_metadata%l_stats_samp = .true. ! Define physics buffers indexes @@ -1679,7 +1678,7 @@ subroutine clubb_ini_cam(pbuf2d) clubb_params(ibv_efold) = clubb_bv_efold clubb_params(iwpxp_Ri_exp) = clubb_wpxp_Ri_exp clubb_params(iz_displace) = clubb_z_displace - + ! 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 @@ -2426,7 +2425,6 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & ! ZM microphysics real(r8), pointer :: dlfzm(:,:) ! ZM detrained convective cloud water mixing ratio. - real(r8), pointer :: difzm(:,:) ! ZM detrained convective cloud ice mixing ratio. real(r8), pointer :: dnlfzm(:,:) ! ZM detrained convective cloud water num concen. real(r8), pointer :: dnifzm(:,:) ! ZM detrained convective cloud ice num concen. @@ -2489,9 +2487,9 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & character(len=*), parameter :: subr='clubb_tend_cam' real(r8), parameter :: rad2deg=180.0_r8/pi real(r8) :: tmp_lon1, tmp_lonN - + type(grid) :: gr - + type(nu_vertical_res_dep) :: nu_vert_res_dep ! Vertical resolution dependent nu values real(r8) :: lmin @@ -3033,10 +3031,10 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & stats_nsamp = nint(stats_metadata%stats_tsamp/dtime) stats_nout = nint(stats_metadata%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 @@ -3333,7 +3331,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & do t=1,nadv ! do needed number of "sub" timesteps for each CAM step - + ! Increment the statistics then begin stats timestep if (stats_metadata%l_stats) then call stats_begin_timestep_api( t, stats_nsamp, stats_nout, & @@ -3808,7 +3806,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & rtm_integral_ltend(:) = rtm_integral_ltend(:)/gravit rtm_integral_vtend(:) = rtm_integral_vtend(:)/gravit - + if (clubb_do_adv) then if (macmic_it == cld_macmic_num_steps) then @@ -4370,8 +4368,8 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end if ! Output CLUBB history here - if (stats_metadata%l_stats) then - + if (stats_metadata%l_stats) then + do j=1,stats_zt(1)%num_output_fields temp1 = trim(stats_zt(1)%file%grid_avg_var(j)%name) @@ -4390,7 +4388,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & call outfld(trim(sub),out_zm(:,:,j), pcols, lchnk) enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%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 @@ -4758,8 +4756,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 @@ -4802,8 +4800,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 @@ -4839,10 +4837,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize rad_zt (radiation points) if (stats_metadata%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 @@ -4876,10 +4874,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & stats_metadata, stats_rad_zt(j) ) ! 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 @@ -4908,7 +4906,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & allocate( stats_rad_zm(j)%file%grid_avg_var( stats_rad_zm(j)%num_output_fields ) ) allocate( stats_rad_zm(j)%file%z( stats_rad_zm(j)%kk ) ) - + call stats_init_rad_zm_api( clubb_vars_rad_zm, & l_error, & stats_metadata, stats_rad_zm(j) ) @@ -4918,8 +4916,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 @@ -4961,30 +4959,30 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & endif ! Now call add fields - + do i = 1, stats_zt(1)%num_output_fields - + temp1 = trim(stats_zt(1)%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > max_fieldname_len) sub = temp1(1:max_fieldname_len) - + call addfld( trim(sub), (/ 'ilev' /), 'A', & trim(stats_zt(1)%file%grid_avg_var(i)%units), & trim(stats_zt(1)%file%grid_avg_var(i)%description) ) enddo - + do i = 1, stats_zm(1)%num_output_fields - + temp1 = trim(stats_zm(1)%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > max_fieldname_len) sub = temp1(1:max_fieldname_len) - + call addfld( trim(sub), (/ 'ilev' /), 'A', & trim(stats_zm(1)%file%grid_avg_var(i)%units), & trim(stats_zm(1)%file%grid_avg_var(i)%description) ) enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%l_output_rad_files) then do i = 1, stats_rad_zt(1)%num_output_fields temp1 = trim(stats_rad_zt(1)%file%grid_avg_var(i)%name) @@ -4994,7 +4992,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_rad_zt(1)%file%grid_avg_var(i)%units), & trim(stats_rad_zt(1)%file%grid_avg_var(i)%description) ) enddo - + do i = 1, stats_rad_zm(1)%num_output_fields temp1 = trim(stats_rad_zm(1)%file%grid_avg_var(i)%name) sub = temp1 @@ -5004,7 +5002,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_rad_zm(1)%file%grid_avg_var(i)%description) ) enddo endif - + do i = 1, stats_sfc(1)%num_output_fields temp1 = trim(stats_sfc(1)%file%grid_avg_var(i)%name) sub = temp1 @@ -5013,7 +5011,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_sfc(1)%file%grid_avg_var(i)%units), & trim(stats_sfc(1)%file%grid_avg_var(i)%description) ) enddo - + return @@ -5102,7 +5100,7 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st enddo enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%l_output_rad_files) then do i = 1, stats_rad_zt%num_output_fields do k = 1, stats_rad_zt%kk out_radzt(thecol,pverp-k+1,i) = stats_rad_zt%accum_field_values(1,1,k,i) diff --git a/src/physics/cam/macrop_driver.F90 b/src/physics/cam/macrop_driver.F90 index 92d52fff8c..66638508f4 100644 --- a/src/physics/cam/macrop_driver.F90 +++ b/src/physics/cam/macrop_driver.F90 @@ -87,7 +87,6 @@ module macrop_driver integer :: & dlfzm_idx = -1, & ! ZM detrained convective cloud water mixing ratio. - difzm_idx = -1, & ! ZM detrained convective cloud ice mixing ratio. dnlfzm_idx = -1, & ! ZM detrained convective cloud water num concen. dnifzm_idx = -1 ! ZM detrained convective cloud ice num concen. @@ -486,7 +485,6 @@ subroutine macrop_driver_tend( & ! ZM microphysics real(r8), pointer :: dlfzm(:,:) ! ZM detrained convective cloud water mixing ratio. - real(r8), pointer :: difzm(:,:) ! ZM detrained convective cloud ice mixing ratio. real(r8), pointer :: dnlfzm(:,:) ! ZM detrained convective cloud water num concen. real(r8), pointer :: dnifzm(:,:) ! ZM detrained convective cloud ice num concen. diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 337c722b43..4c9a6a3550 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -57,7 +57,6 @@ module zm_conv_intr ixorg, & dp_cldice_idx, & dlfzm_idx, & ! detrained convective cloud water mixing ratio. - difzm_idx, & ! detrained convective cloud ice mixing ratio. prec_dp_idx, & snow_dp_idx, & mconzm_idx ! convective mass flux @@ -149,8 +148,6 @@ subroutine zm_conv_register ! detrained convective cloud water mixing ratio. call pbuf_add_field('DLFZM', 'physpkg', dtype_r8, (/pcols,pver/), dlfzm_idx) ! detrained convective cloud ice mixing ratio. - call pbuf_add_field('DIFZM', 'physpkg', dtype_r8, (/pcols,pver/), difzm_idx) - ! convective mass fluxes call pbuf_add_field('CMFMC_DP', 'physpkg', dtype_r8, (/pcols,pverp/), mconzm_idx) !CACNOTE - Is zm_org really a constituent or was it just a handy structure to use for an allocatable which persists in the run? @@ -442,7 +439,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8), pointer, dimension(:,:) :: dp_cldliq real(r8), pointer, dimension(:,:) :: dp_cldice real(r8), pointer :: dlf(:,:) ! detrained convective cloud water mixing ratio. - real(r8), pointer :: dif(:,:) ! detrained convective cloud ice mixing ratio. + real(r8) :: dif(pcols,pver) ! detrained convective cloud ice mixing ratio. real(r8), pointer :: lambdadpcu(:,:) ! slope of cloud liquid size distr real(r8), pointer :: mudpcu(:,:) ! width parameter of droplet size distr real(r8), pointer :: mconzm(:,:) !convective mass fluxes @@ -536,10 +533,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call pbuf_get_field(pbuf, zm_ideep_idx, ideep) call pbuf_get_field(pbuf, dlfzm_idx, dlf) - call pbuf_get_field(pbuf, difzm_idx, dif) call pbuf_get_field(pbuf, mconzm_idx, mconzm) -! ! Begin with Zhang-McFarlane (1996) convection parameterization ! call t_startf ('zm_convr_run') From 61fa56c40d8f33f83b3c7477b75dfc2c96e86ea4 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Thu, 23 May 2024 16:01:59 -0600 Subject: [PATCH 03/18] Remove DIFZM output --- src/physics/cam/zm_conv_intr.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 4c9a6a3550..c2db86a407 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -304,7 +304,6 @@ subroutine zm_conv_init(pref_edge) call addfld ('ZMICVU', (/ 'lev' /), 'A', 'm/s', 'ZM in-cloud V updrafts') call addfld ('ZMICVD', (/ 'lev' /), 'A', 'm/s', 'ZM in-cloud V downdrafts') - call addfld ('DIFZM' ,(/ 'lev' /), 'A','kg/kg/s ','Detrained ice water from ZM convection') call addfld ('DLFZM' ,(/ 'lev' /), 'A','kg/kg/s ','Detrained liquid water from ZM convection') call phys_getopts( history_budget_out = history_budget, & @@ -439,7 +438,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8), pointer, dimension(:,:) :: dp_cldliq real(r8), pointer, dimension(:,:) :: dp_cldice real(r8), pointer :: dlf(:,:) ! detrained convective cloud water mixing ratio. - real(r8) :: dif(pcols,pver) ! detrained convective cloud ice mixing ratio. real(r8), pointer :: lambdadpcu(:,:) ! slope of cloud liquid size distr real(r8), pointer :: mudpcu(:,:) ! width parameter of droplet size distr real(r8), pointer :: mconzm(:,:) !convective mass fluxes @@ -555,7 +553,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & cape(:) = 0._r8 zdu(:,:) = 0._r8 rprd(:,:) = 0._r8 - dif(:,:) = 0._r8 mu(:,:) = 0._r8 eu(:,:) = 0._r8 du(:,:) = 0._r8 @@ -583,7 +580,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & dp(:ncol,:), dsubcld(:ncol), jt(:ncol), maxg(:ncol), ideep(:ncol), & ql(:ncol,:), rliq(:ncol), landfrac(:ncol), & org_ncol(:ncol,:), orgt_ncol(:ncol,:), zm_org2d_ncol(:ncol,:), & - dif(:ncol,:), rice(:ncol), errmsg, errflg) + rice(:ncol), errmsg, errflg) if (zmconv_org) then @@ -636,7 +633,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call outfld('ZMDQ ',ptend_loc%q(1,1,1) ,pcols ,lchnk ) call t_stopf ('zm_convr_run') - call outfld('DIFZM' ,dif ,pcols, lchnk) call outfld('DLFZM' ,dlf ,pcols, lchnk) pcont(:ncol) = state%ps(:ncol) From 4a5cdb578de93bcf050739543381e08aca253c2b Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 29 May 2024 16:28:25 -0600 Subject: [PATCH 04/18] Move limcnv calculation into ZM layer --- src/physics/cam/zm_conv_intr.F90 | 26 ++------------------------ 1 file changed, 2 insertions(+), 24 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index c2db86a407..98c8fad5e9 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -323,32 +323,10 @@ subroutine zm_conv_init(pref_edge) call add_default('ZMMTT ', history_budget_histfile_num, ' ') end if -! -! Limit deep convection to regions below 40 mb -! Note this calculation is repeated in the shallow convection interface -! - limcnv = 0 ! null value to check against below - if (pref_edge(1) >= 4.e3_r8) then - limcnv = 1 - else - do k=1,plev - if (pref_edge(k) < 4.e3_r8 .and. pref_edge(k+1) >= 4.e3_r8) then - limcnv = k - exit - end if - end do - if ( limcnv == 0 ) limcnv = plevp - end if - - if (masterproc) then - write(iulog,*)'ZM_CONV_INIT: Deep convection will be capped at intfc ',limcnv, & - ' which is ',pref_edge(limcnv),' pascals' - end if - no_deep_pbl = phys_deepconv_pbl() !CACNOTE - Need to check errflg and report errors - call zm_convr_init(cpair, epsilo, gravit, latvap, tmelt, rair, & - limcnv,zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_ke_lnd, & + call zm_convr_init(plev, plevp, cpair, epsilo, gravit, latvap, tmelt, rair, & + pref_edge,zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_ke_lnd, & zmconv_momcu, zmconv_momcd, zmconv_num_cin, zmconv_org, & no_deep_pbl, zmconv_tiedke_add, & zmconv_capelmt, zmconv_dmpdz,zmconv_parcel_pbl, zmconv_tau, & From 2713176f7fe488c48294ce42feca25d67353b299 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Tue, 18 Jun 2024 12:54:25 -0600 Subject: [PATCH 05/18] Add constituent hooks for ccpp compatability --- src/physics/cam/physpkg.F90 | 6 +- src/physics/cam/zm_conv_intr.F90 | 8 +- src/physics/cam_dev/physpkg.F90 | 2 +- src/physics/simple/physpkg.F90 | 7 +- .../cam_ccpp/ccpp_constituent_prop_mod.F90 | 98 ++++++++++++++++++- 5 files changed, 108 insertions(+), 13 deletions(-) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index cb7322254f..c9ce88d436 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -740,6 +740,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use co2_cycle, only: co2_init, co2_transport use convect_deep, only: convect_deep_init use convect_shallow, only: convect_shallow_init + use constituents, only: cnst_get_ind use cam_diagnostics, only: diag_init use gw_drag, only: gw_init use cam3_aero_data, only: cam3_aero_data_on, cam3_aero_data_init @@ -798,7 +799,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! local variables integer :: lchnk - integer :: ierr + integer :: ierr, ixq logical :: history_budget ! output tendencies and state variables for ! temperature, water vapor, cloud @@ -983,7 +984,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Initialize CAM CCPP constituent properties array ! for use in CCPP-ized physics schemes: - call ccpp_const_props_init() + call cnst_get_ind('Q', ixq) + call ccpp_const_props_init(ixq) ! Initialize qneg3 and qneg4 call qneg_init() diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 98c8fad5e9..a353684f89 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -452,7 +452,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8) :: icwdu(pcols,pver) real(r8) :: icwdv(pcols,pver) real(r8) :: seten(pcols, pver) - logical :: l_windt(2) + logical :: l_windt real(r8) :: tfinal1, tfinal2 integer :: ii @@ -714,8 +714,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call physics_ptend_init(ptend_loc, state1%psetcols, 'zm_conv_momtran_run', ls=.true., lu=.true., lv=.true.) - l_windt(1) = .true. - l_windt(2) = .true. + l_windt = .true. !REMOVECAM - no longer need these when CAM is retired and pcols no longer exists ptend_loc%s(:,:) = 0._r8 ptend_loc%u(:,:) = 0._r8 @@ -816,6 +815,8 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) use time_manager, only: get_nstep use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc use constituents, only: pcnst, cnst_is_convtran2 + use ccpp_constituent_prop_mod, only: ccpp_const_props + ! Arguments type(physics_state), intent(in ) :: state ! Physics state variables @@ -845,6 +846,7 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) integer, pointer :: jt(:) ! (pcols) integer, pointer :: maxg(:) ! (pcols) integer, pointer :: ideep(:) ! (pcols) + !----------------------------------------------------------------------------------- diff --git a/src/physics/cam_dev/physpkg.F90 b/src/physics/cam_dev/physpkg.F90 index aef997716f..e99d6d9dc3 100644 --- a/src/physics/cam_dev/physpkg.F90 +++ b/src/physics/cam_dev/physpkg.F90 @@ -952,7 +952,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Initialize CAM CCPP constituent properties array ! for use in CCPP-ized physics schemes: - call ccpp_const_props_init() + call ccpp_const_props_init(ixq) ! Initialize qneg3 and qneg4 call qneg_init() diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index a296fd2fdb..f2be0114cc 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -208,6 +208,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use nudging, only: Nudge_Model, nudging_init use cam_snapshot, only: cam_snapshot_init use cam_budget, only: cam_budget_init + use constituents, only: cnst_get_ind + use ccpp_constituent_prop_mod, only: ccpp_const_props_init @@ -220,7 +222,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) type(cam_out_t),intent(inout) :: cam_out(begchunk:endchunk) ! local variables - integer :: lchnk + integer :: lchnk, ixq !----------------------------------------------------------------------- call physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, pcols) @@ -281,7 +283,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Initialize CAM CCPP constituent properties array ! for use in CCPP-ized physics schemes: - call ccpp_const_props_init() + call cnst_get_ind('Q', ixq) + call ccpp_const_props_init(ixq) ! Initialize qneg3 and qneg4 call qneg_init() diff --git a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 index edbc4d59e9..ab795f7c69 100644 --- a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 +++ b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 @@ -8,12 +8,19 @@ module ccpp_constituent_prop_mod type, public :: ccpp_constituent_prop_ptr_t logical, private :: thermo_active = .false. logical, private :: water_species = .false. - contains + logical, private :: species_is_dry + character(len=256) :: std_name = '' + + contains procedure :: standard_name => ccp_get_standard_name + procedure :: set_standard_name => ccp_set_standard_name procedure :: is_thermo_active => ccp_is_thermo_active procedure :: is_water_species => ccp_is_water_species procedure :: set_thermo_active => ccp_set_thermo_active procedure :: set_water_species => ccp_set_water_species + procedure :: is_dry => ccp_is_dry + procedure :: set_dry => ccp_set_dry + end type ccpp_constituent_prop_ptr_t ! CCPP properties init routine @@ -37,20 +44,43 @@ subroutine ccp_get_standard_name(this, std_name, errcode, errmsg) integer, optional, intent(out) :: errcode character(len=*), optional, intent(out) :: errmsg - std_name = 'Not Used!' + std_name = this%std_name ! Provide err values if requested: if(present(errcode)) then errcode = 0 end if if(present(errmsg)) then - errmsg = 'Still Not Used!' + errmsg = '' end if end subroutine ccp_get_standard_name !------ + subroutine ccp_set_standard_name(this, std_name, errcode, errmsg) + ! Return this constituent's standard name + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(inout) :: this + character(len=*), intent(in) :: std_name + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + this%std_name = std_name + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_set_standard_name + + !------ + subroutine ccp_is_thermo_active(this, val_out, errcode, errmsg) ! Dummy arguments @@ -97,6 +127,29 @@ end subroutine ccp_is_water_species !------ + subroutine ccp_is_dry(this, val_out, errcode, errmsg) + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(in) :: this + logical, intent(out) :: val_out + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Pass back water species property: + val_out = this%species_is_dry + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_is_dry + + !------ + subroutine ccp_set_thermo_active(this, thermo_flag, errcode, errmsg) ! Set whether this constituent is thermodynamically active, which ! means that certain physics schemes will use this constitutent @@ -147,18 +200,41 @@ subroutine ccp_set_water_species(this, water_flag, errcode, errmsg) end subroutine ccp_set_water_species + subroutine ccp_set_dry(this, dry_flag, errcode, errmsg) + ! Set whether this constituent is a dry species or not using the dry_flag which is passed in + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(inout) :: this + logical, intent(in) :: dry_flag + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Set dry_flag for this constituent: + this%species_is_dry = dry_flag + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_set_dry + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++ !CAM-equivalent CCPP constituents initialization routine !+++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine ccpp_const_props_init() +subroutine ccpp_const_props_init(ix_qv) ! Use statements: - use constituents, only: pcnst + use constituents, only: pcnst, cnst_get_type_byind use cam_abortutils, only: handle_allocate_error use air_composition, only: dry_air_species_num use air_composition, only: thermodynamic_active_species_idx + integer, intent(in) :: ix_qv ! Local variables: integer :: ierr integer :: m @@ -185,6 +261,18 @@ subroutine ccpp_const_props_init() end if end do + ! Set "set_dry" property: + do m=1,pcnst + if (cnst_get_type_byind(m).eq.'dry') then + call ccpp_const_props(m)%set_dry(.true.) + else + call ccpp_const_props(m)%set_dry(.false.) + end if + end do + + ! Set "std_name" property: + call ccpp_const_props(ix_qv)%set_standard_name('water_vapor_wrt_moist_air_and_condensed_water') + end subroutine ccpp_const_props_init end module ccpp_constituent_prop_mod From 3504c518b03146c51c0d75216c7719766c2416a8 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 19 Jul 2024 15:20:12 -0600 Subject: [PATCH 06/18] Remove the 'use' statements --- bld/configure | 1 + src/physics/cam/wv_sat_methods.F90 | 759 -------------- src/physics/cam/wv_saturation.F90 | 1484 ---------------------------- src/physics/cam/zm_conv_intr.F90 | 9 +- src/utils/error_messages.F90 | 151 --- src/utils/namelist_utils.F90 | 6 - src/utils/units.F90 | 36 - 7 files changed, 8 insertions(+), 2438 deletions(-) delete mode 100644 src/physics/cam/wv_sat_methods.F90 delete mode 100644 src/physics/cam/wv_saturation.F90 delete mode 100644 src/utils/error_messages.F90 delete mode 100644 src/utils/namelist_utils.F90 delete mode 100644 src/utils/units.F90 diff --git a/bld/configure b/bld/configure index e2fb784495..5806d97171 100755 --- a/bld/configure +++ b/bld/configure @@ -2304,6 +2304,7 @@ sub write_filepath #Add the CCPP'ized subdirectories print $fh "$camsrcdir/src/atmos_phys/zhang_mcfarlane\n"; + print $fh "$camsrcdir/src/atmos_phys/to_be_ccppized\n"; # Dynamics package and test utilities print $fh "$camsrcdir/src/dynamics/$dyn\n"; diff --git a/src/physics/cam/wv_sat_methods.F90 b/src/physics/cam/wv_sat_methods.F90 deleted file mode 100644 index bb2ffeb45b..0000000000 --- a/src/physics/cam/wv_sat_methods.F90 +++ /dev/null @@ -1,759 +0,0 @@ -module wv_sat_methods - -! This portable module contains all CAM methods for estimating -! the saturation vapor pressure of water. -! -! wv_saturation provides CAM-specific interfaces and utilities -! based on these formulae. -! -! Typical usage of this module: -! -! Init: -! call wv_sat_methods_init(r8, , errstring) -! -! Get scheme index from a name string: -! scheme_idx = wv_sat_get_scheme_idx(scheme_name) -! if (.not. wv_sat_valid_idx(scheme_idx)) -! -! Get pressures: -! es = wv_sat_svp_water(t, scheme_idx) -! es = wv_sat_svp_ice(t, scheme_idx) -! -! Use ice/water transition range: -! es = wv_sat_svp_trice(t, ttrice, scheme_idx) -! -! Note that elemental functions cannot be pointed to, nor passed -! as arguments. If you need to do either, it is recommended to -! wrap the function so that it can be given an explicit (non- -! elemental) interface. - -implicit none -private -save - -integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real - -integer, parameter :: VLENS = 128 ! vector length for a GPU kernel - -real(r8) :: tmelt ! Melting point of water at 1 atm (K) -real(r8) :: h2otrip ! Triple point temperature of water (K) -real(r8) :: tboil ! Boiling point of water at 1 atm (K) - -real(r8) :: ttrice ! Ice-water transition range - -real(r8) :: epsilo ! Ice-water transition range -real(r8) :: omeps ! 1._r8 - epsilo - -! Indices representing individual schemes -integer, parameter :: Invalid_idx = -1 -integer, parameter :: GoffGratch_idx = 1 -integer, parameter :: MurphyKoop_idx = 2 -integer, parameter :: Bolton_idx = 3 - -! Index representing the current default scheme. -integer, parameter :: initial_default_idx = GoffGratch_idx -integer :: default_idx = initial_default_idx - -!$acc declare create (epsilo, tmelt, tboil, omeps, h2otrip, ttrice) - -public wv_sat_methods_init -public wv_sat_get_scheme_idx -public wv_sat_valid_idx - -public wv_sat_set_default -public wv_sat_reset_default - -public wv_sat_qsat_water, wv_sat_qsat_water_vect -public wv_sat_qsat_ice, wv_sat_qsat_ice_vect - -public wv_sat_svp_trans, wv_sat_svp_trans_vect - -! pressure -> humidity conversion -public wv_sat_svp_to_qsat, wv_sat_svp_to_qsat_vect - -! Combined qsat operations -public wv_sat_qsat_trans - -public wv_sat_svp_water, wv_sat_svp_water_vect -public wv_sat_svp_ice, wv_sat_svp_ice_vect - -contains - -!--------------------------------------------------------------------- -! ADMINISTRATIVE FUNCTIONS -!--------------------------------------------------------------------- - -! Get physical constants -subroutine wv_sat_methods_init(kind, tmelt_in, h2otrip_in, tboil_in, & - ttrice_in, epsilo_in, errstring) - integer, intent(in) :: kind - real(r8), intent(in) :: tmelt_in - real(r8), intent(in) :: h2otrip_in - real(r8), intent(in) :: tboil_in - real(r8), intent(in) :: ttrice_in - real(r8), intent(in) :: epsilo_in - character(len=*), intent(out) :: errstring - - errstring = ' ' - - if (kind /= r8) then - write(errstring,*) 'wv_sat_methods_init: ERROR: ', & - kind,' was input kind but ',r8,' is internal kind.' - return - end if - - if (ttrice_in < 0._r8) then - write(errstring,*) 'wv_sat_methods_init: ERROR: ', & - ttrice_in,' was input for ttrice, but negative range is invalid.' - return - end if - - tmelt = tmelt_in - h2otrip = h2otrip_in - tboil = tboil_in - ttrice = ttrice_in - epsilo = epsilo_in - - omeps = 1._r8 - epsilo - - !$acc update device (tmelt,h2otrip,tboil,ttrice,epsilo,omeps) - -end subroutine wv_sat_methods_init - -! Look up index by name. -pure function wv_sat_get_scheme_idx(name) result(idx) - character(len=*), intent(in) :: name - integer :: idx - - select case (name) - case("GoffGratch") - idx = GoffGratch_idx - case("MurphyKoop") - idx = MurphyKoop_idx - case("Bolton") - idx = Bolton_idx - case default - idx = Invalid_idx - end select - -end function wv_sat_get_scheme_idx - -! Check validity of an index from the above routine. -pure function wv_sat_valid_idx(idx) result(status) - integer, intent(in) :: idx - logical :: status - - status = (idx /= Invalid_idx) - -end function wv_sat_valid_idx - -! Set default scheme (otherwise, Goff & Gratch is default) -! Returns a logical representing success (.true.) or -! failure (.false.). -function wv_sat_set_default(name) result(status) - character(len=*), intent(in) :: name - logical :: status - - ! Don't want to overwrite valid default with invalid, - ! so assign to temporary and check it first. - integer :: tmp_idx - - tmp_idx = wv_sat_get_scheme_idx(name) - - status = wv_sat_valid_idx(tmp_idx) - - if (status) default_idx = tmp_idx - -end function wv_sat_set_default - -! Reset default scheme to initial value. -! The same thing can be accomplished with wv_sat_set_default; -! the real reason to provide this routine is to reset the -! module for testing purposes. -subroutine wv_sat_reset_default() - - default_idx = initial_default_idx - -end subroutine wv_sat_reset_default - -!--------------------------------------------------------------------- -! UTILITIES -!--------------------------------------------------------------------- - -! Get saturation specific humidity given pressure and SVP. -! Specific humidity is limited to range 0-1. -function wv_sat_svp_to_qsat(es, p) result(qs) - real(r8), intent(in) :: es ! SVP - real(r8), intent(in) :: p ! Current pressure. - real(r8) :: qs - - ! If pressure is less than SVP, set qs to maximum of 1. - if ( (p - es) <= 0._r8 ) then - qs = 1.0_r8 - else - qs = epsilo*es / (p - omeps*es) - end if - -end function wv_sat_svp_to_qsat - -! Get saturation specific humidity given pressure and SVP. -! Specific humidity is limited to range 0-1. -subroutine wv_sat_svp_to_qsat_vect(es, p, qs, vlen) - - integer, intent(in) :: vlen - real(r8), intent(in) :: es(vlen) ! SVP - real(r8), intent(in) :: p(vlen) ! Current pressure. - real(r8), intent(out) :: qs(vlen) - integer :: i - - ! If pressure is less than SVP, set qs to maximum of 1. - - !$acc data present (es,p,qs) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i=1,vlen - if ( (p(i) - es(i)) <= 0._r8 ) then - qs(i) = 1.0_r8 - else - qs(i) = epsilo*es(i) / (p(i) - omeps*es(i)) - end if - end do - !$acc end parallel - - !$acc end data -end subroutine wv_sat_svp_to_qsat_vect - -subroutine wv_sat_qsat_water(t, p, es, qs, idx) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over water at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - integer, intent(in), optional :: idx ! Scheme index - - es = wv_sat_svp_water(t, idx) - - qs = wv_sat_svp_to_qsat(es, p) - - ! Ensures returned es is consistent with limiters on qs. - es = min(es, p) - -end subroutine wv_sat_qsat_water - -subroutine wv_sat_qsat_water_vect(t, p, es, qs, vlen, idx) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over water at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - !------------------------------------------------------------------! - ! Inputs - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature - real(r8), intent(in) :: p(vlen) ! Pressure - ! Outputs - real(r8), intent(out) :: es(vlen) ! Saturation vapor pressure - real(r8), intent(out) :: qs(vlen) ! Saturation specific humidity - - integer, intent(in), optional :: idx ! Scheme index - integer :: i - - !$acc data present (t,p,es,qs) - - call wv_sat_svp_water_vect(t, es, vlen, idx) - call wv_sat_svp_to_qsat_vect(es, p, qs, vlen) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i=1,vlen - ! Ensures returned es is consistent with limiters on qs. - es(i) = min(es(i), p(i)) - enddo - !$acc end parallel - - !$acc end data -end subroutine wv_sat_qsat_water_vect - -subroutine wv_sat_qsat_ice(t, p, es, qs, idx) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - integer, intent(in), optional :: idx ! Scheme index - - es = wv_sat_svp_ice(t, idx) - - qs = wv_sat_svp_to_qsat(es, p) - - ! Ensures returned es is consistent with limiters on qs. - es = min(es, p) - -end subroutine wv_sat_qsat_ice - -subroutine wv_sat_qsat_ice_vect(t, p, es, qs, vlen, idx) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - !------------------------------------------------------------------! - ! Inputs - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature - real(r8), intent(in) :: p(vlen) ! Pressure - ! Outputs - real(r8), intent(out) :: es(vlen) ! Saturation vapor pressure - real(r8), intent(out) :: qs(vlen) ! Saturation specific humidity - - integer, intent(in), optional :: idx ! Scheme index - integer :: i - - !$acc data present (t,p,es,qs) - - call wv_sat_svp_ice_vect(t, es, vlen, idx) - call wv_sat_svp_to_qsat_vect(es, p, qs, vlen) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i=1,vlen - ! Ensures returned es is consistent with limiters on qs. - es(i) = min(es(i), p(i)) - enddo - !$acc end parallel - - !$acc end data -end subroutine wv_sat_qsat_ice_vect - -subroutine wv_sat_qsat_trans(t, p, es, qs, idx) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - integer, intent(in), optional :: idx ! Scheme index - - es = wv_sat_svp_trans(t, idx) - - qs = wv_sat_svp_to_qsat(es, p) - - ! Ensures returned es is consistent with limiters on qs. - es = min(es, p) - -end subroutine wv_sat_qsat_trans - -!--------------------------------------------------------------------- -! SVP INTERFACE FUNCTIONS -!--------------------------------------------------------------------- - -function wv_sat_svp_water(t, idx) result(es) - real(r8), intent(in) :: t - integer, intent(in), optional :: idx - real(r8) :: es - - integer :: use_idx - - if (present(idx)) then - use_idx = idx - else - use_idx = default_idx - end if - - select case (use_idx) - case(GoffGratch_idx) - es = GoffGratch_svp_water(t) - case(MurphyKoop_idx) - es = MurphyKoop_svp_water(t) - case(Bolton_idx) - es = Bolton_svp_water(t) - end select - -end function wv_sat_svp_water - -subroutine wv_sat_svp_water_vect(t, es, vlen, idx) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) - integer, intent(in), optional :: idx - real(r8), intent(out) :: es(vlen) - integer :: i - integer :: use_idx - - !$acc data present (t,es) - - if (present(idx)) then - use_idx = idx - else - use_idx = default_idx - end if - - select case (use_idx) - case(GoffGratch_idx) - call GoffGratch_svp_water_vect(t,es,vlen) - case(MurphyKoop_idx) - call MurphyKoop_svp_water_vect(t,es,vlen) - case(Bolton_idx) - call Bolton_svp_water_vect(t,es,vlen) - end select - - !$acc end data -end subroutine wv_sat_svp_water_vect - -function wv_sat_svp_ice(t, idx) result(es) - real(r8), intent(in) :: t - integer, intent(in), optional :: idx - real(r8) :: es - - integer :: use_idx - - if (present(idx)) then - use_idx = idx - else - use_idx = default_idx - end if - - select case (use_idx) - case(GoffGratch_idx) - es = GoffGratch_svp_ice(t) - case(MurphyKoop_idx) - es = MurphyKoop_svp_ice(t) - case(Bolton_idx) - es = Bolton_svp_water(t) - end select - -end function wv_sat_svp_ice - -subroutine wv_sat_svp_ice_vect(t, es, vlen, idx) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) - integer, intent(in), optional :: idx - real(r8), intent(out) :: es(vlen) - integer :: i - - integer :: use_idx - - !$acc data present (t,es) - - if (present(idx)) then - use_idx = idx - else - use_idx = default_idx - end if - - select case (use_idx) - case(GoffGratch_idx) - call GoffGratch_svp_ice_vect(t,es,vlen) - case(MurphyKoop_idx) - call MurphyKoop_svp_ice_vect(t,es,vlen) - case(Bolton_idx) - call Bolton_svp_water_vect(t,es,vlen) - end select - - !$acc end data -end subroutine wv_sat_svp_ice_vect - -function wv_sat_svp_trans(t, idx) result(es) - - real(r8), intent(in) :: t - integer, intent(in), optional :: idx - real(r8) :: es - - real(r8) :: esice ! Saturation vapor pressure over ice - real(r8) :: weight ! Intermediate scratch variable for es transition - -! -! Water -! - if (t >= (tmelt - ttrice)) then - es = wv_sat_svp_water(t,idx) - else - es = 0.0_r8 - end if - -! -! Ice -! - if (t < tmelt) then - - esice = wv_sat_svp_ice(t,idx) - - if ( (tmelt - t) > ttrice ) then - weight = 1.0_r8 - else - weight = (tmelt - t)/ttrice - end if - - es = weight*esice + (1.0_r8 - weight)*es - end if - -end function wv_sat_svp_trans - -subroutine wv_sat_svp_trans_vect(t, es, vlen, idx) - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) - integer, intent(in), optional :: idx - real(r8), intent(out) :: es(vlen) - - real(r8) :: esice(vlen) ! Saturation vapor pressure over ice - real(r8) :: weight ! Intermediate scratch variable for es transition - integer :: i - - !$acc data present (t,es) & - !$acc create (esice) - -! -! Water -! - call wv_sat_svp_water_vect(t,es,vlen,idx) - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - if (t(i) < (tmelt - ttrice)) then - es(i) = 0.0_r8 - end if - end do - !$acc end parallel -! -! Ice -! - call wv_sat_svp_ice_vect(t,esice,vlen,idx) - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - if (t(i) < tmelt) then - if ( (tmelt - t(i)) > ttrice ) then - weight = 1.0_r8 - else - weight = (tmelt - t(i))/ttrice - end if - - es(i) = weight*esice(i) + (1.0_r8 - weight)*es(i) - end if - end do - !$acc end parallel - - !$acc end data -end subroutine wv_sat_svp_trans_vect - -!--------------------------------------------------------------------- -! SVP METHODS -!--------------------------------------------------------------------- - -! Goff & Gratch (1946) - -function GoffGratch_svp_water(t) result(es) - real(r8), intent(in) :: t ! Temperature in Kelvin - real(r8) :: es ! SVP in Pa - - ! uncertain below -70 C - es = 10._r8**(-7.90298_r8*(tboil/t-1._r8)+ & - 5.02808_r8*log10(tboil/t)- & - 1.3816e-7_r8*(10._r8**(11.344_r8*(1._r8-t/tboil))-1._r8)+ & - 8.1328e-3_r8*(10._r8**(-3.49149_r8*(tboil/t-1._r8))-1._r8)+ & - log10(1013.246_r8))*100._r8 - -end function GoffGratch_svp_water - -subroutine GoffGratch_svp_water_vect(t, es, vlen) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature in Kelvin - real(r8), intent(out) :: es(vlen) ! SVP in Pa - real(r8) :: log_tboil - integer :: i - - !$acc data present (t,es) - - ! Goff, J. A., and S. Gratch. “Low-Pressure Properties of Water from -160F - ! to 212F.” Trans. Am. Soc. Heat. Vent. Eng. 52 (1946): 95–121. - ! uncertain below -70 C - - log_tboil = log10(tboil) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i=1,vlen - es(i) = 10._r8**(-7.90298_r8*(tboil/t(i)-1._r8)+ & - 5.02808_r8*(log_tboil-log10(t(i)))- & - 1.3816e-7_r8*(10._r8**(11.344_r8*(1._r8-t(i)/tboil))-1._r8)+ & - 8.1328e-3_r8*(10._r8**(-3.49149_r8*(tboil/t(i)-1._r8))-1._r8)+ & - log10(1013.246_r8))*100._r8 - enddo - !$acc end parallel - - !$acc end data -end subroutine GoffGratch_svp_water_vect - -function GoffGratch_svp_ice(t) result(es) - real(r8), intent(in) :: t ! Temperature in Kelvin - real(r8) :: es ! SVP in Pa - - ! good down to -100 C - es = 10._r8**(-9.09718_r8*(h2otrip/t-1._r8)-3.56654_r8* & - log10(h2otrip/t)+0.876793_r8*(1._r8-t/h2otrip)+ & - log10(6.1071_r8))*100._r8 - -end function GoffGratch_svp_ice - -subroutine GoffGratch_svp_ice_vect(t, es, vlen) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature in Kelvin - real(r8), intent(out) :: es(vlen) ! SVP in Pa - real(r8), parameter :: log_param = log10(6.1071_r8) - integer :: i - ! good down to -100 C - - !$acc data present (t,es) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i=1,vlen - es(i) = 10._r8**(-9.09718_r8*(h2otrip/t(i)-1._r8)-3.56654_r8* & - log10(h2otrip/t(i))+0.876793_r8*(1._r8-t(i)/h2otrip)+ & - log_param)*100._r8 - enddo - !$acc end parallel - - !$acc end data -end subroutine GoffGratch_svp_ice_vect - -! Murphy & Koop (2005) - -function MurphyKoop_svp_water(t) result(es) - real(r8), intent(in) :: t ! Temperature in Kelvin - real(r8) :: es ! SVP in Pa - - ! (good for 123 < T < 332 K) - es = exp(54.842763_r8 - (6763.22_r8 / t) - (4.210_r8 * log(t)) + & - (0.000367_r8 * t) + (tanh(0.0415_r8 * (t - 218.8_r8)) * & - (53.878_r8 - (1331.22_r8 / t) - (9.44523_r8 * log(t)) + & - 0.014025_r8 * t))) - -end function MurphyKoop_svp_water - -subroutine MurphyKoop_svp_water_vect(t, es, vlen) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature in Kelvin - real(r8), intent(out) :: es(vlen) ! SVP in Pa - - integer :: i - ! Murphy, D. M., and T. Koop. “Review of the Vapour Pressure of Ice and - ! Supercooled Water for Atmospheric Applications.” Q. J. R. Meteorol. - ! Soc. 131, no. 608 (2005): 1539–65. 10.1256/qj.04.94 - ! (good for 123 < T < 332 K) - - !$acc data present (t,es) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - es(i) = exp(54.842763_r8 - (6763.22_r8 / t(i)) - (4.210_r8 * log(t(i))) + & - (0.000367_r8 * t(i)) + (tanh(0.0415_r8 * (t(i) - 218.8_r8)) * & - (53.878_r8 - (1331.22_r8 / t(i)) - (9.44523_r8 * log(t(i))) + & - 0.014025_r8 * t(i)))) - end do - !$acc end parallel - - !$acc end data -end subroutine MurphyKoop_svp_water_vect - -function MurphyKoop_svp_ice(t) result(es) - real(r8), intent(in) :: t ! Temperature in Kelvin - real(r8) :: es ! SVP in Pa - - ! (good down to 110 K) - es = exp(9.550426_r8 - (5723.265_r8 / t) + (3.53068_r8 * log(t)) & - - (0.00728332_r8 * t)) - -end function MurphyKoop_svp_ice - -subroutine MurphyKoop_svp_ice_vect(t, es, vlen) - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature in Kelvin - real(r8), intent(out) :: es(vlen) ! SVP in Pa - - integer :: i - ! (good down to 110 K) - - !$acc data present (t,es) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - es(i) = exp(9.550426_r8 - (5723.265_r8 / t(i)) + (3.53068_r8 * log(t(i))) & - - (0.00728332_r8 * t(i))) - end do - !$acc end parallel - - !$acc end data -end subroutine MurphyKoop_svp_ice_vect - -! Bolton (1980) -! zm_conv deep convection scheme contained this SVP calculation. -! It appears to be from D. Bolton, 1980, Monthly Weather Review. -! Unlike the other schemes, no distinct ice formula is associated -! with it. (However, a Bolton ice formula exists in CLUBB.) - -! The original formula used degrees C, but this function -! takes Kelvin and internally converts. - -function Bolton_svp_water(t) result(es) - real(r8),parameter :: c1 = 611.2_r8 - real(r8),parameter :: c2 = 17.67_r8 - real(r8),parameter :: c3 = 243.5_r8 - - real(r8), intent(in) :: t ! Temperature in Kelvin - real(r8) :: es ! SVP in Pa - - es = c1*exp( (c2*(t - tmelt))/((t - tmelt)+c3) ) - -end function Bolton_svp_water - -subroutine Bolton_svp_water_vect(t, es,vlen) - real(r8),parameter :: c1 = 611.2_r8 - real(r8),parameter :: c2 = 17.67_r8 - real(r8),parameter :: c3 = 243.5_r8 - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature in Kelvin - real(r8), intent(out) :: es(vlen) ! SVP in Pa - - integer :: i - - !$acc data present (t,es) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - es(i) = c1*exp( (c2*(t(i) - tmelt))/((t(i) - tmelt)+c3) ) - end do - !$acc end parallel - - !$acc end data -end subroutine Bolton_svp_water_vect - -end module wv_sat_methods diff --git a/src/physics/cam/wv_saturation.F90 b/src/physics/cam/wv_saturation.F90 deleted file mode 100644 index ac94482e20..0000000000 --- a/src/physics/cam/wv_saturation.F90 +++ /dev/null @@ -1,1484 +0,0 @@ -module wv_saturation - -!--------------------------------------------------------------------! -! Module Overview: ! -! ! -! This module provides an interface to wv_sat_methods, providing ! -! saturation vapor pressure and related calculations to CAM. ! -! ! -! The original wv_saturation codes were introduced by J. J. Hack, ! -! February 1990. The code has been extensively rewritten since then, ! -! including a total refactoring in Summer 2012. ! -! ! -!--------------------------------------------------------------------! -! Methods: ! -! ! -! Pure water/ice saturation vapor pressures are calculated on the ! -! fly, with the specific method determined by a runtime option. ! -! Mixed phase SVP is interpolated from the internal table, estbl, ! -! which is created during initialization. ! -! ! -! The default method for calculating SVP is determined by a namelist ! -! option, and used whenever svp_water/ice or qsat are called. ! -! ! -!--------------------------------------------------------------------! - -use shr_kind_mod, only: r8 => shr_kind_r8 -use physconst, only: epsilo, & - latvap, & - latice, & - rh2o, & - cpair, & - tmelt, & - h2otrip - -use wv_sat_methods, only: & - svp_to_qsat => wv_sat_svp_to_qsat, & - svp_to_qsat_vect => wv_sat_svp_to_qsat_vect - -implicit none -private -save - -! Public interfaces -! Namelist, initialization, finalization -public wv_sat_readnl -public wv_sat_init -public wv_sat_final - -! Saturation vapor pressure calculations -public svp_water, svp_water_vect -public svp_ice, svp_ice_vect - -! Mixed phase (water + ice) saturation vapor pressure table lookup -public estblf - -public svp_to_qsat - -! Subroutines that return both SVP and humidity -! Optional arguments do temperature derivatives -interface qsat - module procedure qsat_line - module procedure qsat_vect - module procedure qsat_2D -end interface -public qsat ! Mixed phase -interface qsat_water - module procedure qsat_water_line - module procedure qsat_water_vect - module procedure qsat_water_2D -end interface -public qsat_water ! SVP over water only -interface qsat_ice - module procedure qsat_ice_line - module procedure qsat_ice_vect - module procedure qsat_ice_2D -end interface -public qsat_ice ! SVP over ice only - -! Wet bulb temperature solver -public :: findsp_vc, findsp - -! Data - -! This value is slightly high, but it seems to be the value for the -! steam point of water originally (and most frequently) used in the -! Goff & Gratch scheme. -real(r8), parameter :: tboil = 373.16_r8 - -! Table of saturation vapor pressure values (estbl) from tmin to -! tmax+1 Kelvin, in one degree increments. ttrice defines the -! transition region, estbl contains a combination of ice & water -! values. -! Make these public parameters in case another module wants to see the -! extent of the table. - real(r8), public, parameter :: tmin = 127.16_r8 - real(r8), public, parameter :: tmax = 375.16_r8 - - real(r8), parameter :: ttrice = 20.00_r8 ! transition range from es over H2O to es over ice - - integer :: plenest ! length of estbl - real(r8), allocatable :: estbl(:) ! table values of saturation vapor pressure - - real(r8) :: omeps ! 1.0_r8 - epsilo - - real(r8) :: c3 ! parameter used by findsp - - ! Set coefficients for polynomial approximation of difference - ! between saturation vapor press over water and saturation pressure - ! over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial is - ! valid in the range -40 < t < 0 (degrees C). - real(r8) :: pcf(5) = (/ & - 5.04469588506e-01_r8, & - -5.47288442819e+00_r8, & - -3.67471858735e-01_r8, & - -8.95963532403e-03_r8, & - -7.78053686625e-05_r8 /) - -! --- Degree 6 approximation --- -! real(r8) :: pcf(6) = (/ & -! 7.63285250063e-02, & -! 5.86048427932e+00, & -! 4.38660831780e-01, & -! 1.37898276415e-02, & -! 2.14444472424e-04, & -! 1.36639103771e-06 /) - - integer, parameter :: VLENS = 128 ! vector length for a GPU kernel - - !$acc declare create (plenest,estbl,omeps,c3,pcf) - -contains - -!--------------------------------------------------------------------- -! ADMINISTRATIVE FUNCTIONS -!--------------------------------------------------------------------- - -subroutine wv_sat_readnl(nlfile) - !------------------------------------------------------------------! - ! Purpose: ! - ! Get runtime options for wv_saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_get_scheme_idx, & - wv_sat_valid_idx, & - wv_sat_set_default - - use spmd_utils, only: masterproc - use namelist_utils, only: find_group_name - use units, only: getunit, freeunit - use mpishorthand - use cam_abortutils, only: endrun - - character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - - ! Local variables - integer :: unitn, ierr - - character(len=32) :: wv_sat_scheme = "GoffGratch" - - character(len=*), parameter :: subname = 'wv_sat_readnl' - - namelist /wv_sat_nl/ wv_sat_scheme - !----------------------------------------------------------------------------- - - if (masterproc) then - unitn = getunit() - open( unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'wv_sat_nl', status=ierr) - if (ierr == 0) then - read(unitn, wv_sat_nl, iostat=ierr) - if (ierr /= 0) then - call endrun(subname // ':: ERROR reading namelist') - return - end if - end if - close(unitn) - call freeunit(unitn) - - end if - -#ifdef SPMD - call mpibcast(wv_sat_scheme, len(wv_sat_scheme) , mpichar, 0, mpicom) -#endif - - if (.not. wv_sat_set_default(wv_sat_scheme)) then - call endrun('wv_sat_readnl :: Invalid wv_sat_scheme.') - return - end if - -end subroutine wv_sat_readnl - -subroutine wv_sat_init - !------------------------------------------------------------------! - ! Purpose: ! - ! Initialize module (e.g. setting parameters, initializing the ! - ! SVP lookup table). ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_methods_init, & - wv_sat_get_scheme_idx, & - wv_sat_valid_idx - use spmd_utils, only: masterproc - use cam_logfile, only: iulog - use cam_abortutils, only: endrun - use shr_assert_mod, only: shr_assert_in_domain - use error_messages, only: handle_errmsg - - integer :: status - - ! For wv_sat_methods error reporting. - character(len=256) :: errstring - - ! For generating internal SVP table. - real(r8) :: t ! Temperature - integer :: i ! Increment counter - - ! Precalculated because so frequently used. - omeps = 1.0_r8 - epsilo - - ! Transition range method is only valid for transition temperatures at: - ! -40 deg C < T < 0 deg C - call shr_assert_in_domain(ttrice, ge=0._r8, le=40._r8, varname="ttrice",& - msg="wv_sat_init: Invalid transition temperature range.") - -! This parameter uses a hardcoded 287.04_r8? - c3 = 287.04_r8*(7.5_r8*log(10._r8))/cpair - -! Init "methods" module containing actual SVP formulae. - - call wv_sat_methods_init(r8, tmelt, h2otrip, tboil, ttrice, & - epsilo, errstring) - - call handle_errmsg(errstring, subname="wv_sat_methods_init") - - ! Add two to make the table slightly too big, just in case. - plenest = ceiling(tmax-tmin) + 2 - - ! Allocate SVP table. - allocate(estbl(plenest), stat=status) - if (status /= 0) then - call endrun('wv_sat_init :: ERROR allocating saturation vapor pressure table') - return - end if - - do i = 1, plenest - estbl(i) = svp_trans(tmin + real(i-1,r8)) - end do - - !$acc update device (plenest,estbl,omeps,c3,pcf) - - if (masterproc) then - write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***' - end if - -end subroutine wv_sat_init - -subroutine wv_sat_final - !------------------------------------------------------------------! - ! Purpose: ! - ! Deallocate global variables in module. ! - !------------------------------------------------------------------! - use cam_abortutils, only: endrun - - integer :: status - - if (allocated(estbl)) then - - deallocate(estbl, stat=status) - - if (status /= 0) then - call endrun('wv_sat_final :: ERROR deallocating table') - return - end if - - end if - -end subroutine wv_sat_final - -!--------------------------------------------------------------------- -! DEFAULT SVP FUNCTIONS -!--------------------------------------------------------------------- - -! Compute saturation vapor pressure over water -function svp_water(t) result(es) - - use wv_sat_methods, only: & - wv_sat_svp_water - - real(r8), intent(in) :: t ! Temperature (K) - real(r8) :: es ! SVP (Pa) - - es = wv_sat_svp_water(t) - -end function svp_water - -! Compute saturation vapor pressure over water -subroutine svp_water_vect(t, es, vlen) - - use wv_sat_methods, only: & - wv_sat_svp_water_vect - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature (K) - real(r8), intent(out) :: es(vlen) ! SVP (Pa) - - !$acc data copyin (t) copyout (es) - - call wv_sat_svp_water_vect(t, es, vlen) - - !$acc end data -end subroutine svp_water_vect - -! Compute saturation vapor pressure over ice -function svp_ice(t) result(es) - - use wv_sat_methods, only: & - wv_sat_svp_ice - - real(r8), intent(in) :: t ! Temperature (K) - real(r8) :: es ! SVP (Pa) - - es = wv_sat_svp_ice(t) - -end function svp_ice - -! Compute saturation vapor pressure over ice -subroutine svp_ice_vect(t, es, vlen) - - use wv_sat_methods, only: & - wv_sat_svp_ice_vect - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature (K) - real(r8), intent(out) :: es(vlen) ! SVP (Pa) - - !$acc data copyin(t) copyout(es) - - call wv_sat_svp_ice_vect(t, es, vlen) - - !$acc end data -end subroutine svp_ice_vect - -! Compute saturation vapor pressure with an ice-water transition -function svp_trans(t) result(es) - - use wv_sat_methods, only: & - wv_sat_svp_trans - - real(r8), intent(in) :: t ! Temperature (K) - real(r8) :: es ! SVP (Pa) - - es = wv_sat_svp_trans(t) - -end function svp_trans - -! Compute saturation vapor pressure with an ice-water transition -subroutine svp_trans_vect(t, es, vlen) - - use wv_sat_methods, only: & - wv_sat_svp_trans_vect - - integer, intent(in) :: vlen - real(r8), intent(in) :: t(vlen) ! Temperature (K) - real(r8), intent(out) :: es(vlen) ! SVP (Pa) - - !$acc data copyin(t) copyout(es) - - call wv_sat_svp_trans_vect(t, es, vlen) - - !$acc end data -end subroutine svp_trans_vect - -!--------------------------------------------------------------------- -! UTILITIES -!--------------------------------------------------------------------- - -! Does linear interpolation from nearest values found -! in the table (estbl). -elemental function estblf(t) result(es) - - real(r8), intent(in) :: t ! Temperature - real(r8) :: es ! SVP (Pa) - - integer :: i ! Index for t in the table - real(r8) :: t_tmp ! intermediate temperature for es look-up - - real(r8) :: weight ! Weight for interpolation - - t_tmp = max(min(t,tmax)-tmin, 0._r8) ! Number of table entries above tmin - i = int(t_tmp) + 1 ! Corresponding index. - weight = t_tmp - aint(t_tmp, r8) ! Fractional part of t_tmp (for interpolation). - es = (1._r8 - weight)*estbl(i) + weight*estbl(i+1) - -end function estblf - -! Does linear interpolation from nearest values found -! in the table (estbl). -subroutine estblf_vect(t, es, vlen) - - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(out) :: es ! SVP (Pa) - - integer :: i ! Index for t in the table - integer :: j - real(r8) :: t_tmp ! intermediate temperature for es look-up - - real(r8) :: weight ! Weight for interpolation - - !$acc data present (t,es) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector private(t_tmp,weight,i) - do j = 1, vlen - t_tmp = max(min(t(j),tmax)-tmin, 0._r8) ! Number of table entries above tmin - i = int(t_tmp) + 1 ! Corresponding index. - weight = t_tmp - aint(t_tmp, r8) ! Fractional part of t_tmp (for interpolation). - es(j) = (1._r8 - weight)*estbl(i) + weight*estbl(i+1) - end do - !$acc end parallel - - !$acc end data -end subroutine estblf_vect - -! Get enthalpy based only on temperature -! and specific humidity. -elemental function tq_enthalpy(t, q, hltalt) result(enthalpy) - - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: q ! Specific humidity - real(r8), intent(in) :: hltalt ! Modified hlat for T derivatives - - real(r8) :: enthalpy - - enthalpy = cpair * t + hltalt * q - -end function tq_enthalpy - -! Get enthalpy based only on temperature -! and specific humidity. -subroutine tq_enthalpy_vect(t, q, hltalt, enthalpy, vlen) - - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(in) :: q ! Specific humidity - real(r8), dimension(vlen), intent(in) :: hltalt ! Modified hlat for T derivatives - - real(r8), dimension(vlen), intent(out) :: enthalpy - - integer :: i - - !$acc data present(t,q,hltalt,enthalpy) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - enthalpy(i) = cpair * t(i) + hltalt(i) * q(i) - end do - !$acc end parallel - - !$acc end data -end subroutine tq_enthalpy_vect - -!--------------------------------------------------------------------- -! LATENT HEAT OF VAPORIZATION CORRECTIONS -!--------------------------------------------------------------------- - -elemental subroutine no_ip_hltalt(t, hltalt) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate latent heat of vaporization of pure liquid water at ! - ! a given temperature. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - ! Outputs - real(r8), intent(out) :: hltalt ! Appropriately modified hlat - - hltalt = latvap - - ! Account for change of latvap with t above freezing where - ! constant slope is given by -2369 j/(kg c) = cpv - cw - if (t >= tmelt) then - hltalt = hltalt - 2369.0_r8*(t-tmelt) - end if - -end subroutine no_ip_hltalt - -subroutine no_ip_hltalt_vect(t, hltalt, vlen) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate latent heat of vaporization of pure liquid water at ! - ! a given temperature. ! - !------------------------------------------------------------------! - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - ! Outputs - real(r8), dimension(vlen), intent(out) :: hltalt ! Appropriately modified hlat - - integer :: i - - !$acc data present(t,hltalt) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - hltalt(i) = latvap - ! Account for change of latvap with t above freezing where - ! constant slope is given by -2369 j/(kg c) = cpv - cw - if (t(i) >= tmelt) then - hltalt(i) = hltalt(i) - 2369.0_r8*(t(i)-tmelt) - end if - end do - !$acc end parallel - - !$acc end data -end subroutine no_ip_hltalt_vect - -elemental subroutine calc_hltalt(t, hltalt, tterm) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate latent heat of vaporization of water at a given ! - ! temperature, taking into account the ice phase if temperature ! - ! is below freezing. ! - ! Optional argument also calculates a term used to calculate ! - ! d(es)/dT within the water-ice transition range. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - ! Outputs - real(r8), intent(out) :: hltalt ! Appropriately modified hlat - ! Term to account for d(es)/dT in transition region. - real(r8), intent(out), optional :: tterm - - ! Local variables - real(r8) :: tc ! Temperature in degrees C - real(r8) :: weight ! Weight for es transition from water to ice - ! Loop iterator - integer :: i - - if (present(tterm)) tterm = 0.0_r8 - - call no_ip_hltalt(t,hltalt) - if (t < tmelt) then - ! Weighting of hlat accounts for transition from water to ice. - tc = t - tmelt - - if (tc >= -ttrice) then - weight = -tc/ttrice - - ! polynomial expression approximates difference between es - ! over water and es over ice from 0 to -ttrice (C) (max of - ! ttrice is 40): required for accurate estimate of es - ! derivative in transition range from ice to water - if (present(tterm)) then - do i = size(pcf), 1, -1 - tterm = pcf(i) + tc*tterm - end do - tterm = tterm/ttrice - end if - - else - weight = 1.0_r8 - end if - - hltalt = hltalt + weight*latice - - end if - -end subroutine calc_hltalt - -subroutine calc_hltalt_vect(t, hltalt, vlen, tterm) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate latent heat of vaporization of water at a given ! - ! temperature, taking into account the ice phase if temperature ! - ! is below freezing. ! - ! Optional argument also calculates a term used to calculate ! - ! d(es)/dT within the water-ice transition range. ! - !------------------------------------------------------------------! - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - ! Outputs - real(r8), dimension(vlen), intent(out) :: hltalt ! Appropriately modified hlat - ! Term to account for d(es)/dT in transition region. - real(r8), dimension(vlen), intent(out), optional :: tterm - - ! Local variables - real(r8) :: tc ! Temperature in degrees C - real(r8) :: weight ! Weight for es transition from water to ice - logical :: present_tterm - ! Loop iterator - integer :: i, j, size_pcf - - present_tterm = present(tterm) - size_pcf = size(pcf) - - !$acc data present(t,hltalt,tterm) - - if (present_tterm) then - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - tterm(i) = 0.0_r8 - end do - !$acc end parallel - end if - - call no_ip_hltalt_vect(t,hltalt,vlen) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector private(tc,weight) - do j = 1, vlen - if (t(j) < tmelt) then - ! Weighting of hlat accounts for transition from water to ice. - tc = t(j) - tmelt - - if (tc >= -ttrice) then - weight = -tc/ttrice - - ! polynomial expression approximates difference between es - ! over water and es over ice from 0 to -ttrice (C) (max of - ! ttrice is 40): required for accurate estimate of es - ! derivative in transition range from ice to water - if (present_tterm) then - !$acc loop seq - do i = size_pcf, 1, -1 - tterm(j) = pcf(i) + tc*tterm(j) - end do - tterm(j) = tterm(j)/ttrice - end if - - else - weight = 1.0_r8 - end if - - hltalt(j) = hltalt(j) + weight*latice - - end if - end do - !$acc end parallel - - !$acc end data -end subroutine calc_hltalt_vect - -!--------------------------------------------------------------------- -! OPTIONAL OUTPUTS -!--------------------------------------------------------------------- - -! Temperature derivative outputs, for qsat_* -subroutine deriv_outputs_line(t, p, es, qs, hltalt, tterm, & - gam, dqsdt) - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - real(r8), intent(in) :: es ! Saturation vapor pressure - real(r8), intent(in) :: qs ! Saturation specific humidity - real(r8), intent(in) :: hltalt ! Modified latent heat - real(r8), intent(in) :: tterm ! Extra term for d(es)/dT in - ! transition region. - - ! Outputs - real(r8), intent(out), optional :: gam ! (hltalt/cpair)*(d(qs)/dt) - real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) - - ! Local variables - real(r8) :: desdt ! d(es)/dt - real(r8) :: dqsdt_loc ! local copy of dqsdt - - if (qs == 1.0_r8) then - dqsdt_loc = 0._r8 - else - desdt = hltalt*es/(rh2o*t*t) + tterm - dqsdt_loc = qs*p*desdt/(es*(p-omeps*es)) - end if - - if (present(dqsdt)) dqsdt = dqsdt_loc - if (present(gam)) gam = dqsdt_loc * (hltalt/cpair) - -end subroutine deriv_outputs_line - -! Temperature derivative outputs, for qsat_* -subroutine deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam, dqsdt) - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(in) :: p ! Pressure - real(r8), dimension(vlen), intent(in) :: es ! Saturation vapor pressure - real(r8), dimension(vlen), intent(in) :: qs ! Saturation specific humidity - real(r8), dimension(vlen), intent(in) :: hltalt ! Modified latent heat - real(r8), dimension(vlen), intent(in) :: tterm ! Extra term for d(es)/dT in - ! transition region. - - ! Outputs - real(r8), dimension(vlen), intent(out), optional :: gam ! (hltalt/cpair)*(d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt) - - ! Local variables - real(r8) :: desdt ! d(es)/dt - real(r8) :: dqsdt_loc ! local copy of dqsdt - logical :: present_dqsdt, present_gam - integer :: i - - present_dqsdt = present(dqsdt) - present_gam = present(gam) - - !$acc data present(t,p,es,qs,hltalt,tterm,gam,dqsdt) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector private(dqsdt_loc,desdt) - do i = 1, vlen - if (qs(i) == 1.0_r8) then - dqsdt_loc = 0._r8 - else - desdt = hltalt(i)*es(i)/(rh2o*t(i)*t(i)) + tterm(i) - dqsdt_loc = qs(i)*p(i)*desdt/(es(i)*(p(i)-omeps*es(i))) - end if - - if (present_dqsdt) dqsdt(i) = dqsdt_loc - if (present_gam) gam(i) = dqsdt_loc * (hltalt(i)/cpair) - end do - !$acc end parallel - - !$acc end data -end subroutine deriv_outputs_vect - -!--------------------------------------------------------------------- -! QSAT (SPECIFIC HUMIDITY) PROCEDURES -!--------------------------------------------------------------------- - -subroutine qsat_line(t, p, es, qs, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Look up and return saturation vapor pressure from precomputed ! - ! table, then calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8) :: hltalt ! Modified latent heat for T derivatives - real(r8) :: tterm ! Account for d(es)/dT in transition region - - es = estblf(t) - - qs = svp_to_qsat(es, p) - - ! Ensures returned es is consistent with limiters on qs. - es = min(es, p) - - ! Calculate optional arguments. - if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call calc_hltalt(t, hltalt, tterm) - - if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) - - call deriv_outputs_line(t, p, es, qs, hltalt, tterm, & - gam=gam, dqsdt=dqsdt) - - end if - -end subroutine qsat_line - -subroutine qsat_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Look up and return saturation vapor pressure from precomputed ! - ! table, then calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(vlen) :: tterm ! Account for d(es)/dT in transition region - integer :: i - logical :: present_gam, present_dqsdt, present_enthalpy - - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (hltalt,tterm) - - call estblf_vect(t, es, vlen) - - call svp_to_qsat_vect(es, p, qs, vlen) - - ! Ensures returned es is consistent with limiters on qs. - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - es(i) = min(es(i), p(i)) - end do - !$acc end parallel - - ! Calculate optional arguments. - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call calc_hltalt_vect(t, hltalt, vlen, tterm) - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_vect - -subroutine qsat_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Look up and return saturation vapor pressure from precomputed ! - ! table, then calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - ! Inputs - integer, intent(in) :: dim1, dim2 - real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature - real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(dim1,dim2) :: tterm ! Account for d(es)/dT in transition region - integer :: i, k, vlen - logical :: present_gam, present_dqsdt, present_enthalpy - - vlen = dim1 * dim2 - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (hltalt,tterm) - - call estblf_vect(t, es, vlen) - - call svp_to_qsat_vect(es, p, qs, vlen) - - ! Ensures returned es is consistent with limiters on qs. - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector collapse(2) - do k = 1, dim2 - do i = 1, dim1 - es(i,k) = min(es(i,k), p(i,k)) - end do - end do - !$acc end parallel - - ! Calculate optional arguments. - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call calc_hltalt_vect(t, hltalt, vlen, tterm) - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_2D - -subroutine qsat_water_line(t, p, es, qs, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over water at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_water - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8) :: hltalt ! Modified latent heat for T derivatives - - call wv_sat_qsat_water(t, p, es, qs) - - if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call no_ip_hltalt(t, hltalt) - - if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) - - ! For pure water/ice transition term is 0. - call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, & - gam=gam, dqsdt=dqsdt) - - end if - -end subroutine qsat_water_line - -subroutine qsat_water_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over water at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_water_vect - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(vlen) :: tterm - integer :: i - logical :: present_gam, present_dqsdt, present_enthalpy - - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (tterm,hltalt) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - tterm(i) = 0._r8 - end do - !$acc end parallel - - call wv_sat_qsat_water_vect(t, p, es, qs, vlen) - - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call no_ip_hltalt_vect(t, hltalt, vlen) - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - ! For pure water/ice transition term is 0. - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_water_vect - -subroutine qsat_water_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over water at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_water_vect - - ! Inputs - integer, intent(in) :: dim1, dim2 - real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature - real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(dim1,dim2) :: tterm - integer :: i, k, vlen - logical :: present_gam, present_dqsdt, present_enthalpy - - vlen = dim1 * dim2 - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (hltalt,tterm) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector collapse(2) - do k = 1, dim2 - do i = 1, dim1 - tterm(i,k) = 0._r8 - end do - end do - !$acc end parallel - - call wv_sat_qsat_water_vect(t, p, es, qs, vlen) - - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - ! "generalized" analytic expression for t derivative of es - ! accurate to within 1 percent for 173.16 < t < 373.16 - call no_ip_hltalt_vect(t, hltalt, vlen) - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - ! For pure water/ice transition term is 0. - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_water_2D - -subroutine qsat_ice_line(t, p, es, qs, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_ice - - ! Inputs - real(r8), intent(in) :: t ! Temperature - real(r8), intent(in) :: p ! Pressure - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure - real(r8), intent(out) :: qs ! Saturation specific humidity - - real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8) :: hltalt ! Modified latent heat for T derivatives - - call wv_sat_qsat_ice(t, p, es, qs) - - if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then - - ! For pure ice, just add latent heats. - hltalt = latvap + latice - - if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) - - ! For pure water/ice transition term is 0. - call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, & - gam=gam, dqsdt=dqsdt) - - end if - -end subroutine qsat_ice_line - -subroutine qsat_ice_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_ice_vect - - ! Inputs - integer, intent(in) :: vlen - real(r8), dimension(vlen), intent(in) :: t ! Temperature - real(r8), dimension(vlen), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(vlen) :: tterm - integer :: i - logical :: present_gam, present_dqsdt, present_enthalpy - - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (hltalt,tterm) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - tterm(i) = 0._r8 - end do - !$acc end parallel - - call wv_sat_qsat_ice_vect(t, p, es, qs, vlen) - - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector - do i = 1, vlen - ! For pure ice, just add latent heats. - hltalt(i) = latvap + latice - end do - !$acc end parallel - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - ! For pure water/ice transition term is 0. - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_ice_vect - -subroutine qsat_ice_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy) - !------------------------------------------------------------------! - ! Purpose: ! - ! Calculate SVP over ice at a given temperature, and then ! - ! calculate and return saturation specific humidity. ! - ! Optionally return various temperature derivatives or enthalpy ! - ! at saturation. ! - !------------------------------------------------------------------! - - use wv_sat_methods, only: wv_sat_qsat_ice_vect - - ! Inputs - integer, intent(in) :: dim1, dim2 - real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature - real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure - ! Outputs - real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure - real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity - - real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt) - real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q - - ! Local variables - real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives - real(r8), dimension(dim1,dim2) :: tterm - integer :: i, k, vlen - logical :: present_gam, present_dqsdt, present_enthalpy - - vlen = dim1 * dim2 - present_gam = present(gam) - present_dqsdt = present(dqsdt) - present_enthalpy = present(enthalpy) - - !$acc data copyin (t,p) & - !$acc copyout (es,qs,gam,dqsdt,enthalpy) & - !$acc create (hltalt,tterm) - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector collapse(2) - do k = 1, dim2 - do i = 1, dim1 - tterm(i,k) = 0._r8 - end do - end do - !$acc end parallel - - call wv_sat_qsat_ice_vect(t, p, es, qs, vlen) - - if (present_gam .or. present_dqsdt .or. present_enthalpy) then - - !$acc parallel vector_length(VLENS) default(present) - !$acc loop gang vector collapse(2) - do k = 1, dim2 - do i = 1, dim1 - ! For pure ice, just add latent heats. - hltalt(i,k) = latvap + latice - end do - end do - !$acc end parallel - - if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen) - - ! For pure water/ice transition term is 0. - call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, & - gam=gam, dqsdt=dqsdt) - - end if - - !$acc end data -end subroutine qsat_ice_2D - -!--------------------------------------------------------------------- -! FINDSP (WET BULB TEMPERATURE) PROCEDURES -!--------------------------------------------------------------------- - -subroutine findsp_vc(q, t, p, use_ice, tsp, qsp) - - use cam_logfile, only: iulog - use cam_abortutils, only: endrun - - ! Wrapper for findsp which is 1D and handles the output status. - ! Changing findsp to elemental restricted debugging output. - ! If that output is needed again, it's preferable *not* to copy findsp, - ! but to change the existing version. - - ! input arguments - real(r8), intent(in) :: q(:) ! water vapor (kg/kg) - real(r8), intent(in) :: t(:) ! temperature (K) - real(r8), intent(in) :: p(:) ! pressure (Pa) - logical, intent(in) :: use_ice ! flag to include ice phase in calculations - - ! output arguments - real(r8), intent(out) :: tsp(:) ! saturation temp (K) - real(r8), intent(out) :: qsp(:) ! saturation mixing ratio (kg/kg) - - integer :: status(size(q)) ! flag representing state of output - ! 0 => Successful convergence - ! 1 => No calculation done: pressure or specific - ! humidity not within usable range - ! 2 => Run failed to converge - ! 4 => Temperature fell below minimum - ! 8 => Enthalpy not conserved - - integer :: n, i - - n = size(q) - - ! Currently, only 2 and 8 seem to be treated as fatal errors. - do i = 1,n - call findsp(q(i), t(i), p(i), use_ice, tsp(i), qsp(i), status(i)) - if (status(i) == 2) then - write(iulog,*) ' findsp not converging at i = ', i - write(iulog,*) ' t, q, p ', t(i), q(i), p(i) - write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i) - call endrun ('wv_saturation::FINDSP -- not converging') - else if (status(i) == 8) then - write(iulog,*) ' the enthalpy is not conserved at i = ', i - write(iulog,*) ' t, q, p ', t(i), q(i), p(i) - write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i) - call endrun ('wv_saturation::FINDSP -- enthalpy is not conserved') - endif - end do - -end subroutine findsp_vc - -subroutine findsp (q, t, p, use_ice, tsp, qsp, status) -!----------------------------------------------------------------------- -! -! Purpose: -! find the wet bulb temperature for a given t and q -! in a longitude height section -! wet bulb temp is the temperature and spec humidity that is -! just saturated and has the same enthalpy -! if q > qs(t) then tsp > t and qsp = qs(tsp) < q -! if q < qs(t) then tsp < t and qsp = qs(tsp) > q -! -! Method: -! a Newton method is used -! first guess uses an algorithm provided by John Petch from the UKMO -! we exclude points where the physical situation is unrealistic -! e.g. where the temperature is outside the range of validity for the -! saturation vapor pressure, or where the water vapor pressure -! exceeds the ambient pressure, or the saturation specific humidity is -! unrealistic -! -! Author: P. Rasch -! -!----------------------------------------------------------------------- -! -! input arguments -! - - real(r8), intent(in) :: q ! water vapor (kg/kg) - real(r8), intent(in) :: t ! temperature (K) - real(r8), intent(in) :: p ! pressure (Pa) - logical, intent(in) :: use_ice ! flag to include ice phase in calculations -! -! output arguments -! - real(r8), intent(out) :: tsp ! saturation temp (K) - real(r8), intent(out) :: qsp ! saturation mixing ratio (kg/kg) - integer, intent(out) :: status ! flag representing state of output - ! 0 => Successful convergence - ! 1 => No calculation done: pressure or specific - ! humidity not within usable range - ! 2 => Run failed to converge - ! 4 => Temperature fell below minimum - ! 8 => Enthalpy not conserved -! -! local variables -! - integer, parameter :: iter = 8 ! max number of times to iterate the calculation - integer :: l ! iterator - - real(r8) es ! sat. vapor pressure - real(r8) gam ! change in sat spec. hum. wrt temperature (times hltalt/cpair) - real(r8) dgdt ! work variable - real(r8) g ! work variable - real(r8) hltalt ! lat. heat. of vap. - real(r8) qs ! spec. hum. of water vapor - -! work variables - real(r8) t1, q1, dt, dq - real(r8) qvd - real(r8) r1b, c1, c2 - real(r8), parameter :: dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration - real(r8), parameter :: dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration - real(r8) enin, enout - - ! Saturation specific humidity at this temperature - if (use_ice) then - call qsat(t, p, es, qs) - else - call qsat_water(t, p, es, qs) - end if - - ! make sure a meaningful calculation is possible - if (p <= 5._r8*es .or. qs <= 0._r8 .or. qs >= 0.5_r8 & - .or. t < tmin .or. t > tmax) then - status = 1 - ! Keep initial parameters when conditions aren't suitable - tsp = t - qsp = q - enin = 1._r8 - enout = 1._r8 - - return - end if - - ! Prepare to iterate - status = 2 - - ! Get initial enthalpy - if (use_ice) then - call calc_hltalt(t,hltalt) - else - call no_ip_hltalt(t,hltalt) - end if - enin = tq_enthalpy(t, q, hltalt) - - ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch) - c1 = hltalt*c3 - c2 = (t + 36._r8)**2 - r1b = c2/(c2 + c1*qs) - qvd = r1b * (q - qs) - tsp = t + ((hltalt/cpair)*qvd) - - ! Generate qsp, gam, and enout from tsp. - if (use_ice) then - call qsat(tsp, p, es, qsp, gam=gam, enthalpy=enout) - else - call qsat_water(tsp, p, es, qsp, gam=gam, enthalpy=enout) - end if - - ! iterate on first guess - do l = 1, iter - - g = enin - enout - dgdt = -cpair * (1 + gam) - - ! New tsp - t1 = tsp - g/dgdt - dt = abs(t1 - tsp)/t1 - tsp = t1 - - ! bail out if past end of temperature range - if ( tsp < tmin ) then - tsp = tmin - ! Get latent heat and set qsp to a value - ! that preserves enthalpy. - if (use_ice) then - call calc_hltalt(tsp,hltalt) - else - call no_ip_hltalt(tsp,hltalt) - end if - qsp = (enin - cpair*tsp)/hltalt - enout = tq_enthalpy(tsp, qsp, hltalt) - status = 4 - exit - end if - - ! Re-generate qsp, gam, and enout from new tsp. - if (use_ice) then - call qsat(tsp, p, es, q1, gam=gam, enthalpy=enout) - else - call qsat_water(tsp, p, es, q1, gam=gam, enthalpy=enout) - end if - dq = abs(q1 - qsp)/max(q1,1.e-12_r8) - qsp = q1 - - ! if converged at this point, exclude it from more iterations - if (dt < dttol .and. dq < dqtol) then - status = 0 - exit - endif - end do - - ! Test for enthalpy conservation - if (abs((enin-enout)/(enin+enout)) > 1.e-4_r8) status = 8 - -end subroutine findsp - -end module wv_saturation diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index a353684f89..3aa4c1e30d 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -359,6 +359,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & use physconst, only: gravit, latice, latvap, tmelt, cpwv, cpliq, rh2o use phys_control, only: cam_physpkg_is + use ccpp_constituent_prop_mod, only: ccpp_const_props ! Arguments @@ -788,7 +789,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ptend_loc%lq,state1%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ztodt) + nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ztodt, ccpp_const_props, errflg, errmsg) call t_stopf ('convtran1') call outfld('ZMDICE ',ptend_loc%q(1,1,ixcldice) ,pcols ,lchnk ) @@ -847,6 +848,9 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) integer, pointer :: maxg(:) ! (pcols) integer, pointer :: ideep(:) ! (pcols) + character(len=512) :: errmsg + integer :: errflg + !----------------------------------------------------------------------------------- @@ -886,11 +890,12 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) ptend%q(:,:,:) = 0._r8 !REMOVECAM_END +!CACNOTE - Need to check errflg and report errors call zm_conv_convtran_run (ncol, pver, & ptend%lq,state%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ztodt) + nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ztodt, ccpp_const_props, errflg, errmsg) call t_stopf ('convtran2') end if diff --git a/src/utils/error_messages.F90 b/src/utils/error_messages.F90 deleted file mode 100644 index a2a64bca91..0000000000 --- a/src/utils/error_messages.F90 +++ /dev/null @@ -1,151 +0,0 @@ -module error_messages - - !----------------------------------------------------------------------- - ! - ! Purpose: - ! General purpose routines for issuing error messages. - ! - ! Author: B. Eaton - ! - !----------------------------------------------------------------------- - use cam_abortutils, only: endrun - use cam_logfile, only: iulog - - implicit none - save - private - public :: & - alloc_err, &! Issue error message after non-zero return from an allocate statement. - handle_err, &! Issue error message after non-zero return from anything - handle_ncerr ! Handle error returns from netCDF library procedures. - - ! If an error message string is not empty, abort with that string as the - ! error message. - public :: handle_errmsg - -!############################################################################## -contains -!############################################################################## - - subroutine alloc_err( istat, routine, name, nelem ) - - !----------------------------------------------------------------------- - ! Purpose: - ! Issue error message after non-zero return from an allocate statement. - ! - ! Author: B. Eaton - !----------------------------------------------------------------------- - - integer, intent(in) ::& - istat ! status from allocate statement - character(len=*), intent(in) ::& - routine, &! routine that called allocate - name ! name of array - integer, intent(in) ::& - nelem ! number of elements attempted to allocate - !----------------------------------------------------------------------- - - if ( istat .ne. 0 ) then - write(iulog,*)'ERROR trying to allocate memory in routine: ' & - //trim(routine) - write(iulog,*)' Variable name: '//trim(name) - write(iulog,*)' Number of elements: ',nelem - call endrun ('ALLOC_ERR') - end if - - return - - end subroutine alloc_err - -!############################################################################## - - subroutine handle_err( istat, msg ) - - !----------------------------------------------------------------------- - ! Purpose: - ! Issue error message after non-zero return from anything. - ! - ! Author: T. Henderson - !----------------------------------------------------------------------- - - integer, intent(in) :: istat ! status, zero = "no error" - character(len=*), intent(in) :: msg ! error message to print - !----------------------------------------------------------------------- - - if ( istat .ne. 0 ) then - call endrun (trim(msg)) - end if - - return - - end subroutine handle_err - -!############################################################################## - - subroutine handle_ncerr( ret, mes, line ) - - !----------------------------------------------------------------------- - ! Purpose: - ! Check netCDF library function return code. If error detected - ! issue error message then abort. - ! - ! Author: B. Eaton - !----------------------------------------------------------------------- - -!----------------------------------------------------------------------- - use netcdf -!----------------------------------------------------------------------- - - integer, intent(in) ::& - ret ! return code from netCDF library routine - character(len=*), intent(in) ::& - mes ! message to be printed if error detected - integer, intent(in), optional :: line - !----------------------------------------------------------------------- - - if ( ret .ne. NF90_NOERR ) then - if(present(line)) then - write(iulog,*) mes, line - else - write(iulog,*) mes - end if - write(iulog,*) nf90_strerror( ret ) - call endrun ('HANDLE_NCERR') - endif - - return - - end subroutine handle_ncerr - -!############################################################################## - - subroutine handle_errmsg(errmsg, subname, extra_msg) - - ! String that is asserted to be null. - character(len=*), intent(in) :: errmsg - ! Name of procedure generating the message. - character(len=*), intent(in), optional :: subname - ! Additional message from the procedure calling this one. - character(len=*), intent(in), optional :: extra_msg - - if (trim(errmsg) /= "") then - - if (present(extra_msg)) & - write(iulog,*) "handle_errmsg: & - &Message from caller: ",trim(extra_msg) - - if (present(subname)) then - call endrun("ERROR: handle_errmsg: "// & - trim(subname)//": "//trim(errmsg)) - else - call endrun("ERROR: handle_errmsg: "// & - "Error message received from routine: "//trim(errmsg)) - end if - - end if - - end subroutine handle_errmsg - -!############################################################################## - -end module error_messages diff --git a/src/utils/namelist_utils.F90 b/src/utils/namelist_utils.F90 deleted file mode 100644 index c12dfad2d6..0000000000 --- a/src/utils/namelist_utils.F90 +++ /dev/null @@ -1,6 +0,0 @@ -module namelist_utils - -use shr_nl_mod, only: & - find_group_name => shr_nl_find_group_name - -end module namelist_utils diff --git a/src/utils/units.F90 b/src/utils/units.F90 deleted file mode 100644 index 8497cbf687..0000000000 --- a/src/utils/units.F90 +++ /dev/null @@ -1,36 +0,0 @@ -module units - -use shr_file_mod, only: shr_file_getUnit, shr_file_freeUnit - -implicit none -private - -public :: getunit, freeunit - -!------------------------------------------------------------------------------- -contains -!------------------------------------------------------------------------------- - -integer function getunit(iu) - - ! return an available unit number for i/o - - integer, intent(in), optional :: iu ! desired unit number - - getunit = shr_file_getUnit(iu) - -end function getunit - -!------------------------------------------------------------------------------- - -subroutine freeunit(iu) - - ! release the unit - - integer, intent(in) :: iu ! unit number to be freed - - call shr_file_freeUnit(iu) - -end subroutine freeunit - -end module units From a6b5db5c923955be7ade6971ebd1f743192fbd0f Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 2 Aug 2024 17:09:59 -0600 Subject: [PATCH 07/18] Bug fixes to get cleanup3 to run in CAM --- bld/configure | 1 + src/physics/cam/convect_shallow.F90 | 22 +++++++++++++++++++++- src/physics/cam/zm_conv_intr.F90 | 22 +++++++++++++++++++--- 3 files changed, 41 insertions(+), 4 deletions(-) diff --git a/bld/configure b/bld/configure index 5806d97171..9707679a48 100755 --- a/bld/configure +++ b/bld/configure @@ -2305,6 +2305,7 @@ sub write_filepath #Add the CCPP'ized subdirectories print $fh "$camsrcdir/src/atmos_phys/zhang_mcfarlane\n"; print $fh "$camsrcdir/src/atmos_phys/to_be_ccppized\n"; + print $fh "$camsrcdir/src/atmos_phys/cloud_fraction\n"; # Dynamics package and test utilities print $fh "$camsrcdir/src/dynamics/$dyn\n"; diff --git a/src/physics/cam/convect_shallow.F90 b/src/physics/cam/convect_shallow.F90 index ffd1db8f5f..a99cb3ad17 100644 --- a/src/physics/cam/convect_shallow.F90 +++ b/src/physics/cam/convect_shallow.F90 @@ -18,6 +18,9 @@ module convect_shallow use cam_history, only : outfld, addfld, horiz_only use cam_logfile, only : iulog use phys_control, only : phys_getopts + use cloud_fraction_fice, only: cloud_fraction_fice_run + use ref_pres, only: trop_cloud_top_lev + use phys_control, only: phys_getopts implicit none private @@ -472,10 +475,18 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & real(r8), pointer, dimension(:,:) :: cmfmc2 ! (pcols,pverp) Updraft mass flux by shallow convection [ kg/s/m2 ] real(r8), pointer, dimension(:,:) :: sh_e_ed_ratio ! (pcols,pver) fer/(fer+fdr) from uwschu + real(r8), dimension(pcols,pver) :: fsnow_conv + real(r8), dimension(pcols,pver) :: fice + logical :: lq(pcnst) type(unicon_out_t) :: unicon_out + character(len=16) :: macrop_scheme + integer :: top_lev + + + ! ----------------------- ! ! Main Computation Begins ! ! ----------------------- ! @@ -871,8 +882,16 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & tend_s_snwprd(:,:) = 0._r8 tend_s_snwevmlt(:,:) = 0._r8 snow(:) = 0._r8 + fice(:,:) = 0._r8 + fsnow_conv(:,:) = 0._r8 !REMOVECAM_END + top_lev = 1 + call phys_getopts (macrop_scheme_out = macrop_scheme) + if ( .not. (macrop_scheme == "rk" .or. macrop_scheme == "SPCAM_sam1mom")) top_lev = trop_cloud_top_lev + + call cloud_fraction_fice_run(ncol, state1%t(1:ncol,:), tmelt, top_lev, pver, fice(1:ncol,:), fsnow_conv(1:ncol,:)) + call zm_conv_evap_run(state1%ncol, pver, pverp, & gravit, latice, latvap, tmelt, & cpair, zmconv_ke, zmconv_ke_lnd, zmconv_org, & @@ -880,7 +899,8 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & landfracdum(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprdsh(:ncol,:), cld(:ncol,:), ztodt, & - precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:) ) + precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:), & + fsnow_conv(:ncol,:) ) ! ---------------------------------------------- ! ! record history variables from zm_conv_evap_run ! diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 3aa4c1e30d..41acb25dcf 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -14,6 +14,7 @@ module zm_conv_intr use zm_convr, only: zm_convr_init, zm_convr_run use zm_conv_convtran, only: zm_conv_convtran_run use zm_conv_momtran, only: zm_conv_momtran_run + use cloud_fraction_fice, only: cloud_fraction_fice_run use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_mode_num, rad_cnst_get_aer_mmr, & rad_cnst_get_aer_props, rad_cnst_get_mode_props !, & @@ -24,6 +25,8 @@ module zm_conv_intr use perf_mod use cam_logfile, only: iulog use constituents, only: cnst_add + use ref_pres, only: trop_cloud_top_lev + use phys_control, only: phys_getopts implicit none private @@ -457,6 +460,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8) :: tfinal1, tfinal2 integer :: ii + real(r8) :: fice(pcols,pver) + real(r8) :: fsnow_conv(pcols,pver) real(r8),pointer :: zm_org2d(:,:) real(r8),allocatable :: orgt_alloc(:,:), org_alloc(:,:) @@ -464,6 +469,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8) :: orgt_ncol(state%ncol,pver), org_ncol(state%ncol,pver) logical :: lq(pcnst) + character(len=16) :: macrop_scheme + integer :: top_lev !---------------------------------------------------------------------- @@ -660,8 +667,16 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & flxprec(:,:) = 0._r8 flxsnow(:,:) = 0._r8 snow(:) = 0._r8 + fice(:,:) = 0._r8 + fsnow_conv(:,:) = 0._r8 !REMOVECAM_END + top_lev = 1 + call phys_getopts (macrop_scheme_out = macrop_scheme) + if ( .not. (macrop_scheme == "rk" .or. macrop_scheme == "SPCAM_sam1mom")) top_lev = trop_cloud_top_lev + + call cloud_fraction_fice_run(ncol, state1%t(:ncol,:), tmelt, top_lev, pver, fice(:ncol,:), fsnow_conv(:ncol,:)) + call zm_conv_evap_run(state1%ncol, pver, pverp, & gravit, latice, latvap, tmelt, & cpair, zmconv_ke, zmconv_ke_lnd, zmconv_org, & @@ -669,7 +684,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & landfrac(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprd(:ncol,:), cld(:ncol,:), ztodt, & - prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:) ) + prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:),& + fsnow_conv(:ncol,:) ) evapcdp(:ncol,:pver) = ptend_loc%q(:ncol,:pver,1) @@ -789,7 +805,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ptend_loc%lq,state1%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ztodt, ccpp_const_props, errflg, errmsg) + nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ztodt, ccpp_const_props, errmsg, errflg) call t_stopf ('convtran1') call outfld('ZMDICE ',ptend_loc%q(1,1,ixcldice) ,pcols ,lchnk ) @@ -895,7 +911,7 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) ptend%lq,state%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ztodt, ccpp_const_props, errflg, errmsg) + nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ztodt, ccpp_const_props, errmsg, errflg) call t_stopf ('convtran2') end if From 1e78173e18ed477a46e2a9370498a6a2cb599e9b Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 5 Aug 2024 12:01:09 -0600 Subject: [PATCH 08/18] Fixes after the atmos_phys merge --- src/physics/cam/convect_shallow.F90 | 3 +-- src/physics/cam/zm_conv_intr.F90 | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/physics/cam/convect_shallow.F90 b/src/physics/cam/convect_shallow.F90 index a99cb3ad17..3252b4cf03 100644 --- a/src/physics/cam/convect_shallow.F90 +++ b/src/physics/cam/convect_shallow.F90 @@ -899,8 +899,7 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & landfracdum(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprdsh(:ncol,:), cld(:ncol,:), ztodt, & - precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:), & - fsnow_conv(:ncol,:) ) + precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:)) ! ---------------------------------------------- ! ! record history variables from zm_conv_evap_run ! diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 41acb25dcf..932566f79b 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -684,8 +684,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & landfrac(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprd(:ncol,:), cld(:ncol,:), ztodt, & - prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:),& - fsnow_conv(:ncol,:) ) + prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:)) evapcdp(:ncol,:pver) = ptend_loc%q(:ncol,:pver,1) From 49bfd2a0ef16a5b26a5b7bbccf3a6332e86fe158 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Thu, 19 Sep 2024 17:29:30 -0600 Subject: [PATCH 09/18] Updated calling lists for zm routines --- src/physics/cam/convect_shallow.F90 | 6 +++++- src/physics/cam/zm_conv_intr.F90 | 20 +++++++++++--------- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/physics/cam/convect_shallow.F90 b/src/physics/cam/convect_shallow.F90 index eff7b700eb..8874242a5b 100644 --- a/src/physics/cam/convect_shallow.F90 +++ b/src/physics/cam/convect_shallow.F90 @@ -482,7 +482,10 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & type(unicon_out_t) :: unicon_out + character(len=40) :: scheme_name character(len=16) :: macrop_scheme + character(len=512):: errmsg + integer :: errflg integer :: top_lev @@ -899,7 +902,8 @@ subroutine convect_shallow_tend( ztodt , cmfmc , & landfracdum(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprdsh(:ncol,:), cld(:ncol,:), ztodt, & - precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:)) + precc(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:),& + scheme_name, errmsg, errflg) ! ---------------------------------------------- ! ! record history variables from zm_conv_evap_run ! diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index b3c8afc6ac..d58b0b197d 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -488,6 +488,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & logical :: lq(pcnst) character(len=16) :: macrop_scheme + character(len=40) :: scheme_name integer :: top_lev !---------------------------------------------------------------------- @@ -573,11 +574,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & mu(:ncol,:), md(:ncol,:), du(:ncol,:), eu(:ncol,:), ed(:ncol,:), & dp(:ncol,:), dsubcld(:ncol), jt(:ncol), maxg(:ncol), ideep(:ncol), & ql(:ncol,:), rliq(:ncol), landfrac(:ncol), & - rice(:ncol), errmsg, errflg) - - - lengath = count(ideep > 0) - if (lengath > ncol) lengath = ncol ! should not happen, but force it to not be larger than ncol for safety sake + rice(:ncol), lengath, scheme_name, errmsg, errflg) jctop(:) = real(pver,r8) jcbot(:) = 1._r8 @@ -683,7 +680,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & landfrac(:ncol), & ptend_loc%s(:ncol,:), tend_s_snwprd(:ncol,:), tend_s_snwevmlt(:ncol,:), ptend_loc%q(:ncol,:pver,1), & rprd(:ncol,:), cld(:ncol,:), ztodt, & - prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:)) + prec(:ncol), snow(:ncol), ntprprd(:ncol,:), ntsnprd(:ncol,:), fsnow_conv(:ncol,:), flxprec(:ncol,:), flxsnow(:ncol,:),& + scheme_name, errmsg, errflg) evapcdp(:ncol,:pver) = ptend_loc%q(:ncol,:pver,1) @@ -739,7 +737,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & nstep, ptend_loc%u(:ncol,:), ptend_loc%v(:ncol,:),& pguallu(:ncol,:), pguallv(:ncol,:), pgdallu(:ncol,:), pgdallv(:ncol,:), & - icwuu(:ncol,:), icwuv(:ncol,:), icwdu(:ncol,:), icwdv(:ncol,:), ztodt, seten(:ncol,:) ) + icwuu(:ncol,:), icwuv(:ncol,:), icwdu(:ncol,:), icwdv(:ncol,:), ztodt, seten(:ncol,:) ,& + scheme_name, errmsg, errflg) call t_stopf ('zm_conv_momtran_run') ptend_loc%s(:ncol,:pver) = seten(:ncol,:pver) @@ -793,7 +792,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ptend_loc%lq,state1%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ccpp_const_props, errmsg, errflg) + nstep, fracis(:ncol,:,:), ptend_loc%q(:ncol,:,:), fake_dpdry(:ncol,:), ccpp_const_props, & + scheme_name, errmsg, errflg) call t_stopf ('convtran1') call outfld('ZMDICE ',ptend_loc%q(1,1,ixcldice) ,pcols ,lchnk ) @@ -849,6 +849,7 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) integer, pointer :: maxg(:) ! (pcols) integer, pointer :: ideep(:) ! (pcols) + character(len=40) :: scheme_name character(len=512) :: errmsg integer :: errflg @@ -896,7 +897,8 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) ptend%lq,state%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & - nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ccpp_const_props, errmsg, errflg) + nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ccpp_const_props, & + scheme_name, errmsg, errflg) call t_stopf ('convtran2') end if From 0fdf1cc5055eea4356c40c8f3555642d2946642e Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 18 Dec 2024 16:21:31 -0700 Subject: [PATCH 10/18] Add diagnostic messages --- src/physics/cam/zm_conv_intr.F90 | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index bd168efeb6..4df3a7e3ad 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -344,7 +344,6 @@ subroutine zm_conv_init(pref_edge) end if no_deep_pbl = phys_deepconv_pbl() -!CACNOTE - Need to check errflg and report errors call zm_convr_init(plev, plevp, cpair, epsilo, gravit, latvap, tmelt, rair, & pref_edge,zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_ke_lnd, & zmconv_momcu, zmconv_momcd, zmconv_num_cin, & @@ -352,6 +351,10 @@ subroutine zm_conv_init(pref_edge) zmconv_capelmt, zmconv_dmpdz,zmconv_parcel_pbl, zmconv_tau, & masterproc, iulog, errmsg, errflg) + if (errflg /= 0) then + call endrun('From zm_convr_init:' // errmsg) + end if + cld_idx = pbuf_get_index('CLD') fracis_idx = pbuf_get_index('FRACIS') @@ -377,6 +380,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & use constituents, only: pcnst, cnst_get_ind, cnst_is_convtran1 use check_energy, only: check_energy_chng use physconst, only: gravit, latice, latvap, tmelt, cpwv, cpliq, rh2o + use phys_grid, only: get_rlat_all_p, get_rlon_all_p use phys_control, only: cam_physpkg_is use ccpp_constituent_prop_mod, only: ccpp_const_props @@ -456,6 +460,8 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & real(r8) :: pcont(pcols), pconb(pcols), freqzm(pcols) + real(r8) :: lat_all(pcols), long_all(pcols) + ! history output fields real(r8) :: cape(pcols) ! w convective available potential energy. real(r8) :: mu_out(pcols,pver) @@ -482,6 +488,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & logical :: lq(pcnst) character(len=16) :: macrop_scheme character(len=40) :: scheme_name + character(len=40) :: str integer :: top_lev !---------------------------------------------------------------------- @@ -557,9 +564,13 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ideep(:) = 0._r8 !REMOVECAM_END -!CACNOTE - Need to check errflg and report errors + + call get_rlat_all_p(lchnk, ncol, lat_all) + call get_rlon_all_p(lchnk, ncol, long_all) + call zm_convr_run(ncol, pver, & pverp, gravit, latice, cpwv, cpliq, rh2o, & + lat_all, long_all, & state%t(:ncol,:), state%q(:ncol,:,1), prec(:ncol), & pblh(:ncol), state%zm(:ncol,:), state%phis(:ncol), state%zi(:ncol,:), ptend_loc%q(:ncol,:,1), & ptend_loc%s(:ncol,:), state%pmid(:ncol,:), state%pint(:ncol,:), state%pdel(:ncol,:), & @@ -570,6 +581,11 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ql(:ncol,:), rliq(:ncol), landfrac(:ncol), & rice(:ncol), lengath, scheme_name, errmsg, errflg) + if (errflg /= 0) then + write(str,*) 'From zm_convr_run: at chunk',lchnk + call endrun(str // errmsg) + end if + jctop(:) = real(pver,r8) jcbot(:) = 1._r8 do i = 1,lengath @@ -878,13 +894,17 @@ subroutine zm_conv_tend_2( state, ptend, ztodt, pbuf) ptend%q(:,:,:) = 0._r8 !REMOVECAM_END -!CACNOTE - Need to check errflg and report errors call zm_conv_convtran_run (ncol, pver, & ptend%lq,state%q(:ncol,:,:), pcnst, mu(:ncol,:), md(:ncol,:), & du(:ncol,:), eu(:ncol,:), ed(:ncol,:), dp(:ncol,:), dsubcld(:ncol), & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & nstep, fracis(:ncol,:,:), ptend%q(:ncol,:,:), dpdry(:ncol,:), ccpp_const_props, & scheme_name, errmsg, errflg) + + if (errflg /= 0) then + call endrun('From zm_conv_convtran_run:' // errmsg) + end if + call t_stopf ('convtran2') end if From 062373633a97186ad0ce3a38582d370b78fa7601 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 3 Jan 2025 10:48:51 -0700 Subject: [PATCH 11/18] Eliminate final whitespace which snuck in file --- src/physics/cam/clubb_intr.F90 | 92 +++++++++++++++++----------------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index f238dbd07f..20964ab845 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -31,7 +31,7 @@ module clubb_intr #ifdef CLUBB_SGS use clubb_api_module, only: pdf_parameter, implicit_coefs_terms - use clubb_api_module, only: clubb_config_flags_type, grid, stats, & + use clubb_api_module, only: clubb_config_flags_type, grid, stats, & nu_vertical_res_dep, stats_metadata_type, & hm_metadata_type, sclr_idx_type @@ -52,7 +52,7 @@ module clubb_intr stats_sfc(pcols) ! stats_sfc type (hm_metadata_type) :: & hm_metadata - + type (stats_metadata_type) :: & stats_metadata @@ -95,7 +95,7 @@ module clubb_intr ! These are zero by default, but will be set by SILHS before they are used by subcolumns integer :: & - hydromet_dim = 0, & + hydromet_dim = 0, & pdf_dim = 0 @@ -212,7 +212,7 @@ module clubb_intr real(r8) :: clubb_bv_efold = unset_r8 real(r8) :: clubb_wpxp_Ri_exp = unset_r8 real(r8) :: clubb_z_displace = unset_r8 - + integer :: & clubb_iiPDF_type, & ! Selected option for the two-component normal ! (double Gaussian) PDF type to use for the w, rt, @@ -340,7 +340,7 @@ module clubb_intr clubb_l_mono_flux_lim_um, & ! Flag to turn on monotonic flux limiter for um clubb_l_mono_flux_lim_vm, & ! Flag to turn on monotonic flux limiter for vm clubb_l_mono_flux_lim_spikefix, & ! Flag to implement monotonic flux limiter code that - ! eliminates spurious drying tendencies at model top + ! eliminates spurious drying tendencies at model top clubb_l_host_applies_sfc_fluxes ! Whether the host model applies the surface fluxes logical :: & @@ -349,7 +349,7 @@ module clubb_intr ! Constant parameters logical, parameter, private :: & l_implemented = .true. ! Implemented in a host model (always true) - + logical, parameter, private :: & apply_to_heat = .false. ! Apply WACCM energy fixer to heat or not (.true. = yes (duh)) @@ -517,7 +517,7 @@ 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, dtype_i4, dyn_time_lvls use subcol_utils, only: subcol_get_scheme @@ -891,7 +891,7 @@ subroutine clubb_readnl(nlfile) !----- 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 stats_metadata%l_stats = .false. ! Initialize to false stats_metadata%l_output_rad_files = .false. ! Initialize to false @@ -1254,7 +1254,7 @@ subroutine clubb_readnl(nlfile) ! Overwrite defaults if they are true if (clubb_history) stats_metadata%l_stats = .true. - if (clubb_rad_history) stats_metadata%l_output_rad_files = .true. + if (clubb_rad_history) stats_metadata%l_output_rad_files = .true. if (clubb_cloudtop_cooling) do_cldcool = .true. if (clubb_rainevap_turb) do_rainturb = .true. @@ -1731,7 +1731,7 @@ subroutine clubb_ini_cam(pbuf2d) clubb_params_single_col(iwpxp_Ri_exp) = clubb_wpxp_Ri_exp clubb_params_single_col(iz_displace) = clubb_z_displace - ! Override clubb default + ! Override clubb default if ( trim(subcol_scheme) == 'SILHS' ) then clubb_config_flags%saturation_formula = saturation_flatau else @@ -1740,7 +1740,7 @@ subroutine clubb_ini_cam(pbuf2d) ! Define model constant parameters call setup_parameters_model_api( theta0, ts_nudge, clubb_params_single_col(iSkw_max_mag) ) - + ! 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 @@ -2581,15 +2581,15 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & character(len=*), parameter :: subr='clubb_tend_cam' real(r8), parameter :: rad2deg=180.0_r8/pi real(r8) :: tmp_lon1, tmp_lonN - + type(grid) :: gr - + type(nu_vertical_res_dep) :: nu_vert_res_dep ! Vertical resolution dependent nu values real(r8) :: lmin real(r8), dimension(state%ncol,nparams) :: & clubb_params ! Adjustable CLUBB parameters (C1, C2 ...) - + #endif det_s(:) = 0.0_r8 det_ice(:) = 0.0_r8 @@ -3138,10 +3138,10 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & stats_nsamp = nint(stats_metadata%stats_tsamp/dtime) stats_nout = nint(stats_metadata%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 @@ -3448,7 +3448,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & if (macmic_it==1) wpthlp_clubb_gw_mc(:ncol,:) = 0._r8 do t=1,nadv ! do needed number of "sub" timesteps for each CAM step - + ! Increment the statistics then begin stats timestep if (stats_metadata%l_stats) then call stats_begin_timestep_api( t, stats_nsamp, stats_nout, & @@ -4510,8 +4510,8 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & end if ! Output CLUBB history here - if (stats_metadata%l_stats) then - + if (stats_metadata%l_stats) then + do j=1,stats_zt(1)%num_output_fields temp1 = trim(stats_zt(1)%file%grid_avg_var(j)%name) @@ -4530,7 +4530,7 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & call outfld(trim(sub),out_zm(:,:,j), pcols, lchnk) enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%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 @@ -4898,8 +4898,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 @@ -4944,8 +4944,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 @@ -4983,10 +4983,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & ! Initialize rad_zt (radiation points) if (stats_metadata%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 @@ -5020,10 +5020,10 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & stats_metadata, stats_rad_zt(j) ) ! 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 @@ -5052,7 +5052,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & allocate( stats_rad_zm(j)%file%grid_avg_var( stats_rad_zm(j)%num_output_fields ) ) allocate( stats_rad_zm(j)%file%z( stats_rad_zm(j)%kk ) ) - + call stats_init_rad_zm_api( clubb_vars_rad_zm, & l_error, & stats_metadata, stats_rad_zm(j) ) @@ -5062,8 +5062,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 @@ -5105,30 +5105,30 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & endif ! Now call add fields - + do i = 1, stats_zt(1)%num_output_fields - + temp1 = trim(stats_zt(1)%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > max_fieldname_len) sub = temp1(1:max_fieldname_len) - + call addfld( trim(sub), (/ 'ilev' /), 'A', & trim(stats_zt(1)%file%grid_avg_var(i)%units), & trim(stats_zt(1)%file%grid_avg_var(i)%description) ) enddo - + do i = 1, stats_zm(1)%num_output_fields - + temp1 = trim(stats_zm(1)%file%grid_avg_var(i)%name) sub = temp1 if (len(temp1) > max_fieldname_len) sub = temp1(1:max_fieldname_len) - + call addfld( trim(sub), (/ 'ilev' /), 'A', & trim(stats_zm(1)%file%grid_avg_var(i)%units), & trim(stats_zm(1)%file%grid_avg_var(i)%description) ) enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%l_output_rad_files) then do i = 1, stats_rad_zt(1)%num_output_fields temp1 = trim(stats_rad_zt(1)%file%grid_avg_var(i)%name) @@ -5138,7 +5138,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_rad_zt(1)%file%grid_avg_var(i)%units), & trim(stats_rad_zt(1)%file%grid_avg_var(i)%description) ) enddo - + do i = 1, stats_rad_zm(1)%num_output_fields temp1 = trim(stats_rad_zm(1)%file%grid_avg_var(i)%name) sub = temp1 @@ -5148,7 +5148,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_rad_zm(1)%file%grid_avg_var(i)%description) ) enddo endif - + do i = 1, stats_sfc(1)%num_output_fields temp1 = trim(stats_sfc(1)%file%grid_avg_var(i)%name) sub = temp1 @@ -5157,7 +5157,7 @@ subroutine stats_init_clubb( l_stats_in, stats_tsamp_in, stats_tout_in, & trim(stats_sfc(1)%file%grid_avg_var(i)%units), & trim(stats_sfc(1)%file%grid_avg_var(i)%description) ) enddo - + return @@ -5246,7 +5246,7 @@ subroutine stats_end_timestep_clubb(thecol, stats_zt, stats_zm, stats_rad_zt, st enddo enddo - if (stats_metadata%l_output_rad_files) then + if (stats_metadata%l_output_rad_files) then do i = 1, stats_rad_zt%num_output_fields do k = 1, stats_rad_zt%kk out_radzt(thecol,pverp-k+1,i) = stats_rad_zt%accum_field_values(1,1,k,i) From 6f63133e9f23e9d869eb0488373f88675a78cec5 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 6 Jan 2025 09:37:02 -0700 Subject: [PATCH 12/18] Address reviewer comments --- src/physics/cam/cloud_fraction.F90 | 152 +++++------------- src/physics/cam/macrop_driver.F90 | 5 +- src/physics/cam/rk_stratiform.F90 | 12 +- .../cam_ccpp/ccpp_constituent_prop_mod.F90 | 2 +- 4 files changed, 52 insertions(+), 119 deletions(-) diff --git a/src/physics/cam/cloud_fraction.F90 b/src/physics/cam/cloud_fraction.F90 index 3285862fae..821b84695d 100644 --- a/src/physics/cam/cloud_fraction.F90 +++ b/src/physics/cam/cloud_fraction.F90 @@ -5,7 +5,7 @@ module cloud_fraction use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver, pverp - use ref_pres, only: pref_mid + use ref_pres, only: pref_mid use spmd_utils, only: masterproc use cam_logfile, only: iulog use cam_abortutils, only: endrun @@ -22,7 +22,6 @@ module cloud_fraction cldfrc_init, &! Inititialization of cloud_fraction run-time parameters cldfrc_getparams, &! public access of tuning parameters cldfrc, &! Computation of cloud fraction - cldfrc_fice, &! Calculate fraction of condensate in ice phase (radiation partitioning) dp1, &! parameter for deep convection cloud fraction needed in clubb_intr dp2 ! parameter for deep convection cloud fraction needed in clubb_intr @@ -32,9 +31,9 @@ module cloud_fraction ! Top level integer :: top_lev = 1 - ! Physics buffer indices - integer :: sh_frac_idx = 0 - integer :: dp_frac_idx = 0 + ! Physics buffer indices + integer :: sh_frac_idx = 0 + integer :: dp_frac_idx = 0 ! Namelist variables logical :: cldfrc_freeze_dry ! switch for Vavrus correction @@ -154,8 +153,8 @@ subroutine cldfrc_register !----------------------------------------------------------------------- - call pbuf_add_field('SH_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), sh_frac_idx) - call pbuf_add_field('DP_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), dp_frac_idx) + call pbuf_add_field('SH_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), sh_frac_idx) + call pbuf_add_field('DP_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), dp_frac_idx) end subroutine cldfrc_register @@ -215,7 +214,7 @@ subroutine cldfrc_init inversion_cld_off = .false. endif - if ( masterproc ) then + if ( masterproc ) then write(iulog,*)'tuning parameters cldfrc_init: inversion_cld_off',inversion_cld_off write(iulog,*)'tuning parameters cldfrc_init: dp1',dp1,'dp2',dp2,'sh1',sh1,'sh2',sh2 if (shallow_scheme .ne. 'UW' .or. shallow_scheme .eq. 'SPCAM_m2005' ) then @@ -249,38 +248,38 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & cmfmc ,cmfmc2 ,landfrac,snowh ,concld ,cldst , & ts ,sst ,ps ,zdu ,ocnfrac ,& rhu00 ,cldice ,icecldf ,liqcldf ,relhum ,dindex ) - !----------------------------------------------------------------------- - ! - ! Purpose: - ! Compute cloud fraction - ! - ! - ! Method: + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute cloud fraction + ! + ! + ! Method: ! This calculate cloud fraction using a relative humidity threshold - ! The threshold depends upon pressure, and upon the presence or absence - ! of convection as defined by a reasonably large vertical mass flux + ! The threshold depends upon pressure, and upon the presence or absence + ! of convection as defined by a reasonably large vertical mass flux ! entering that layer from below. - ! + ! ! Author: Many. Last modified by Jim McCaa - ! + ! !----------------------------------------------------------------------- use cam_history, only: outfld use physconst, only: cappa, gravit, rair, tmelt use wv_saturation, only: qsat, qsat_water, svp_ice_vect use phys_grid, only: get_rlat_all_p, get_rlon_all_p - + !RBN - Need this to write shallow,deep fraction to phys buffer. !PJR - we should probably make seperate modules for determining convective ! clouds and make this one just responsible for relative humidity clouds - + use physics_buffer, only: physics_buffer_desc, pbuf_get_field ! Arguments integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: dindex ! 0 or 1 to perturb rh - + type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures real(r8), intent(in) :: temp(pcols,pver) ! temperature @@ -307,7 +306,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & real(r8), intent(out) :: clc(pcols) ! column convective cloud amount real(r8), intent(out) :: cldst(pcols,pver) ! cloud fraction real(r8), intent(out) :: rhu00(pcols,pver) ! RH threshold for cloud - real(r8), intent(out) :: relhum(pcols,pver) ! RH + real(r8), intent(out) :: relhum(pcols,pver) ! RH real(r8), intent(out) :: icecldf(pcols,pver) ! ice cloud fraction real(r8), intent(out) :: liqcldf(pcols,pver) ! liquid cloud fraction (combined into cloud) @@ -376,7 +375,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! The idea is that the RH limits for condensation are strict only for ! water saturation ! - ! Ice clouds are formed by explicit parameterization of ice nucleation. + ! Ice clouds are formed by explicit parameterization of ice nucleation. ! Closure for ice cloud fraction is done on available cloud ice, such that ! the in-cloud ice content matches an empirical fit ! thus, icecldf = min(cldice/icicval,1) where icicval = f(temp,cldice,numice) @@ -385,17 +384,17 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! No dA/dt term for ice? ! ! There are three co-existing cloud types: convective, inversion related low-level - ! stratocumulus, and layered cloud (based on relative humidity). Layered and - ! stratocumulus clouds do not compete with convective cloud for which one creates - ! the most cloud. They contribute collectively to the total grid-box average cloud - ! amount. This is reflected in the way in which the total cloud amount is evaluated + ! stratocumulus, and layered cloud (based on relative humidity). Layered and + ! stratocumulus clouds do not compete with convective cloud for which one creates + ! the most cloud. They contribute collectively to the total grid-box average cloud + ! amount. This is reflected in the way in which the total cloud amount is evaluated ! (a sum as opposed to a logical "or" operation) ! !================================================================================== ! set defaults for rhu00 rhu00(:,:) = 2.0_r8 ! define rh perturbation in order to estimate rhdfda - rhpert = 0.01_r8 + rhpert = 0.01_r8 !set Wang and Sassen IWC paramters a=26.87_r8 @@ -460,7 +459,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! ! Estimate of local convective cloud cover based on convective mass flux - ! Modify local large-scale relative humidity to account for presence of + ! Modify local large-scale relative humidity to account for presence of ! convective cloud when evaluating relative humidity based layered cloud amount ! concld(:ncol,top_lev:pver) = 0.0_r8 @@ -468,7 +467,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! cloud mass flux in SI units of kg/m2/s; should produce typical numbers of 20% ! shallow and deep convective cloudiness are evaluated separately (since processes ! are evaluated separately) and summed - ! + ! #ifndef PERGRO do k=top_lev,pver do i=1,ncol @@ -488,7 +487,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! ****** Compute layer cloudiness ****** ! !==================================================================== - ! Begin the evaluation of layered cloud amount based on (modified) RH + ! Begin the evaluation of layered cloud amount based on (modified) RH !==================================================================== ! numkcld = pver @@ -517,7 +516,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! SJV: decrease cloud amount if very low water vapor content ! (thus very cold): "freeze dry" if (cldfrc_freeze_dry) then - rhcloud(i,k) = rhcloud(i,k)*max(0.15_r8,min(1.0_r8,q(i,k)/0.0030_r8)) + rhcloud(i,k) = rhcloud(i,k)*max(0.15_r8,min(1.0_r8,q(i,k)/0.0030_r8)) endif else if ( pmid(i,k).lt.premit ) then @@ -537,7 +536,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & ! linear rh threshold transition between thresholds for low & high cloud ! rhwght = (premib-(max(pmid(i,k),premit)))/(premib-premit) - + if (land(i) .and. (snowh(i) <= 0.000001_r8)) then rhlim = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght) else @@ -591,7 +590,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & !--------ICE CLOUD OPTION 3--------Wood & Field 2000 (JAS) ! eq 6: cloud fraction = 1 - exp (-K * qc/qsati) - + icecldf(i,k)=1._r8 - exp(-Kc*cldice(i,k)/(qs(i,k)*(esi(i,k)/esl(i,k)))) icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8)) else @@ -634,7 +633,7 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & cloud(i,k) = rhcloud(i,k) end if end do - end do + end do ! ! Add in the marine strat ! MARINE STRATUS SHOULD BE A SPECIAL CASE OF LAYERED CLOUD @@ -644,20 +643,20 @@ subroutine cldfrc(lchnk ,ncol , pbuf, & !=================================================================================== ! ! SOME OBSERVATIONS ABOUT THE FOLLOWING SECTION OF CODE (missed in earlier look) - ! K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON - ! LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND + ! K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON + ! LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND ! DISCONTINUOUS IN SPACE (I.E., STRATUS WILL END SUDDENLY WITH NO TRANSITION) ! ! IT APPEARS THAT STRAT IS EVALUATED ACCORDING TO KLEIN AND HARTMANN; HOWEVER, ! THE ACTUAL STRATUS AMOUNT (CLDST) APPEARS TO DEPEND DIRECTLY ON THE RH BELOW - ! THE STRONGEST PART OF THE LOW LEVEL INVERSION. + ! THE STRONGEST PART OF THE LOW LEVEL INVERSION. !PJR answers: 1) the rh limitation is a physical/mathematical limitation ! cant have more cloud than there is RH ! allowed the cloud to exist two layers below the inversion ! because the numerics frequently make 50% relative humidity ! in level below the inversion which would allow no cloud ! 2) since the cloud is only allowed over ocean, it should - ! be very insensitive to surface pressure (except due to + ! be very insensitive to surface pressure (except due to ! spectral ringing, which also causes so many other problems ! I didnt worry about it. ! @@ -738,77 +737,4 @@ end subroutine cldfrc !================================================================================================ - subroutine cldfrc_fice(ncol, t, fice, fsnow) -! -! Compute the fraction of the total cloud water which is in ice phase. -! The fraction depends on temperature only. -! This is the form that was used for radiation, the code came from cldefr originally -! -! Author: B. A. Boville Sept 10, 2002 -! modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection ) -!----------------------------------------------------------------------- - use physconst, only: tmelt - -! Arguments - integer, intent(in) :: ncol ! number of active columns - real(r8), intent(in) :: t(:,:) ! temperature - - real(r8), intent(out) :: fice(:,:) ! Fractional ice content within cloud - real(r8), intent(out) :: fsnow(:,:) ! Fractional snow content for convection - -! Local variables - real(r8) :: tmax_fice ! max temperature for cloud ice formation - real(r8) :: tmin_fice ! min temperature for cloud ice formation - real(r8) :: tmax_fsnow ! max temperature for transition to convective snow - real(r8) :: tmin_fsnow ! min temperature for transition to convective snow - - integer :: i,k ! loop indexes - -!----------------------------------------------------------------------- - - tmax_fice = tmelt - 10._r8 ! max temperature for cloud ice formation - tmin_fice = tmax_fice - 30._r8 ! min temperature for cloud ice formation - tmax_fsnow = tmelt ! max temperature for transition to convective snow - tmin_fsnow = tmelt - 5._r8 ! min temperature for transition to convective snow - - fice(:,:top_lev-1) = 0._r8 - fsnow(:,:top_lev-1) = 0._r8 - -! Define fractional amount of cloud that is ice - do k=top_lev,pver - do i=1,ncol - -! If warmer than tmax then water phase - if (t(i,k) > tmax_fice) then - fice(i,k) = 0.0_r8 - -! If colder than tmin then ice phase - else if (t(i,k) < tmin_fice) then - fice(i,k) = 1.0_r8 - -! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax - else - fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice) - end if - -! snow fraction partitioning - -! If warmer than tmax then water phase - if (t(i,k) > tmax_fsnow) then - fsnow(i,k) = 0.0_r8 - -! If colder than tmin then ice phase - else if (t(i,k) < tmin_fsnow) then - fsnow(i,k) = 1.0_r8 - -! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax - else - fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow) - end if - - end do - end do - - end subroutine cldfrc_fice - end module cloud_fraction diff --git a/src/physics/cam/macrop_driver.F90 b/src/physics/cam/macrop_driver.F90 index 66638508f4..3441c5059e 100644 --- a/src/physics/cam/macrop_driver.F90 +++ b/src/physics/cam/macrop_driver.F90 @@ -394,7 +394,8 @@ subroutine macrop_driver_tend( & ! ! !-------------------------------------------------------- ! - use cloud_fraction, only: cldfrc, cldfrc_fice + use cloud_fraction, only: cldfrc + use cloud_fraction_fice, only: cloud_fraction_fice_run use physics_types, only: physics_state, physics_ptend use physics_types, only: physics_ptend_init, physics_update use physics_types, only: physics_ptend_sum, physics_state_copy @@ -870,8 +871,8 @@ subroutine macrop_driver_tend( & fice(:,:) = 0._r8 fsnow(:,:) = 0._r8 !REMOVECAM_END - call cldfrc_fice( ncol, state_loc%t(:ncol,:), fice(:ncol,:), fsnow(:ncol,:) ) + call cloud_fraction_fice_run(ncol, state_loc%t(:ncol,:), tmelt, top_lev, pver, fice(:ncol,:), fsnow(:ncol,:)) lq(:) = .FALSE. diff --git a/src/physics/cam/rk_stratiform.F90 b/src/physics/cam/rk_stratiform.F90 index 002300dbfd..148ea3fd28 100644 --- a/src/physics/cam/rk_stratiform.F90 +++ b/src/physics/cam/rk_stratiform.F90 @@ -427,7 +427,8 @@ subroutine rk_stratiform_tend( & ! ! !-------------------------------------------------------- ! - use cloud_fraction, only: cldfrc, cldfrc_fice + use cloud_fraction, only: cldfrc + use cloud_fraction_fice, only: cloud_fraction_fice_run use physics_types, only: physics_state, physics_ptend use physics_types, only: physics_ptend_init, physics_update use physics_types, only: physics_ptend_sum, physics_state_copy @@ -440,7 +441,7 @@ subroutine rk_stratiform_tend( & use phys_control, only: cam_physpkg_is use tropopause, only: tropopause_find_cam use phys_grid, only: get_rlat_all_p - use physconst, only: pi + use physconst, only: pi, tmelt ! Arguments type(physics_state), intent(in) :: state ! State variables @@ -577,6 +578,9 @@ subroutine rk_stratiform_tend( & real(r8) :: dlat(pcols) real(r8), parameter :: rad2deg = 180._r8/pi + integer :: top_lev + + ! ====================================================================== lchnk = state%lchnk @@ -812,7 +816,9 @@ subroutine rk_stratiform_tend( & fice(:,:) = 0._r8 fsnow(:,:) = 0._r8 !REMOVECAM_END - call cldfrc_fice(ncol, state1%t(1:ncol,:), fice(1:ncol,:), fsnow(1:ncol,:)) + top_lev = 1 + call cloud_fraction_fice_run(ncol, state1%t(:ncol,:), tmelt, top_lev, pver, fice(:ncol,:), fsnow(:ncol,:)) + ! Perform repartitioning of stratiform condensate. ! Corresponding heating tendency will be added later. diff --git a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 index ab795f7c69..fd08cacfb7 100644 --- a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 +++ b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 @@ -271,7 +271,7 @@ subroutine ccpp_const_props_init(ix_qv) end do ! Set "std_name" property: - call ccpp_const_props(ix_qv)%set_standard_name('water_vapor_wrt_moist_air_and_condensed_water') + call ccpp_const_props(ix_qv)%set_standard_name('water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water') end subroutine ccpp_const_props_init From e900c7112c61f8d892db972d965dd573a3ef85ae Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 22 Jan 2025 15:32:13 -0700 Subject: [PATCH 13/18] Use Fotran unit handling --- src/physics/cam/zm_conv_intr.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 4df3a7e3ad..1d0e732994 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -152,7 +152,6 @@ subroutine zm_conv_readnl(nlfile) use spmd_utils, only: mpicom, masterproc, masterprocid, mpi_real8, mpi_integer, mpi_logical use namelist_utils, only: find_group_name - use units, only: getunit, freeunit character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input @@ -168,8 +167,7 @@ subroutine zm_conv_readnl(nlfile) !----------------------------------------------------------------------------- if (masterproc) then - unitn = getunit() - open( unitn, file=trim(nlfile), status='old' ) + open( newunit=unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'zmconv_nl', status=ierr) if (ierr == 0) then read(unitn, zmconv_nl, iostat=ierr) @@ -178,7 +176,6 @@ subroutine zm_conv_readnl(nlfile) end if end if close(unitn) - call freeunit(unitn) end if From 1180831a90e99525c3dc7f364a6e7498c958320b Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 29 Jan 2025 13:19:26 -0700 Subject: [PATCH 14/18] Update atmos_phys tag --- .gitmodules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index 752c3bf4f7..1d38fe3287 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_07_001 + fxtag = atmos_phys0_08_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics From 85b14f35270001ce74c472cf4e1ede43844730e5 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 10 Feb 2025 09:25:55 -0700 Subject: [PATCH 15/18] remove unnecessary directory for configure to use --- bld/configure | 1 - 1 file changed, 1 deletion(-) diff --git a/bld/configure b/bld/configure index 21d4c5e0d1..a64d6e4605 100755 --- a/bld/configure +++ b/bld/configure @@ -2344,7 +2344,6 @@ sub write_filepath print $fh "$camsrcdir/src/atmos_phys/schemes/check_energy\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/utilities\n"; - print $fh "$camsrcdir/src/atmos_phys/schemes/to_be_ccppized\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/cloud_fraction\n"; # Dynamics package and test utilities From 3ab5a5dae3d2deb0997d3aa531f54eabe930c83d Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 10 Feb 2025 10:46:54 -0700 Subject: [PATCH 16/18] Address reviewer comments --- bld/config_files/definition.xml | 4 ---- bld/configure | 18 ------------------ src/physics/cam/convect_shallow.F90 | 1 - src/physics/cam/zm_conv_intr.F90 | 6 +++--- .../cam_ccpp/ccpp_constituent_prop_mod.F90 | 2 +- 5 files changed, 4 insertions(+), 27 deletions(-) diff --git a/bld/config_files/definition.xml b/bld/config_files/definition.xml index 0b7b6bca45..9feec64a03 100644 --- a/bld/config_files/definition.xml +++ b/bld/config_files/definition.xml @@ -83,10 +83,6 @@ Switch to turn on UNICON package: 0 => off, 1 => on Switch to turn on/off advecting CLUBB moments: 0 => no, 1 => yes - -Switch to turn on/off parameterization for sub-grid scale convective organization for the ZM deep convective scheme based -on Mapes and Neale (2011): 0 => no, 1 => yes - PBL package: uw (University of Washington), hb (Holtslag and Boville), hbr (Holtslag, Boville, and Rasch), clubb_sgs, spcam_sam1om, spcam_m2005, none. diff --git a/bld/configure b/bld/configure index a64d6e4605..82bae67421 100755 --- a/bld/configure +++ b/bld/configure @@ -116,8 +116,6 @@ OPTIONS The user does not need to set this if one of the waccm chemistry options is chosen. -waccmx Build CAM/WACCM with WACCM upper Thermosphere/Ionosphere extended package - -zmconv_org Include parameterization for sub-grid scale convective organization for the ZM deep convective scheme based - on Mapes and Neale (2011) Options relevent to SCAM mode: @@ -317,7 +315,6 @@ GetOptions( "version" => \$opts{'version'}, "waccm_phys" => \$opts{'waccm_phys'}, "waccmx" => \$opts{'waccmx'}, - "zmconv_org" => \$opts{'zmconv_org'}, ) or usage(); # Give usage message. @@ -939,16 +936,6 @@ if (defined $opts{'clubb_opts'}) { my $clubb_do_adv = $cfg_ref->get('clubb_do_adv'); if ($print>=2) { print "clubb_do_adv: $clubb_do_adv$eol"; } -#----------------------------------------------------------------------------------------------- -# ZM convective organization - -if (defined $opts{'zmconv_org'}) { - $cfg_ref->set('zmconv_org', $opts{'zmconv_org'}); -} - -my $zmconv_org = $cfg_ref->get('zmconv_org'); -if ($print>=2) { print "zmconv_org: $zmconv_org$eol"; } - #----------------------------------------------------------------------------------------------- # Macro-physics package @@ -1569,11 +1556,6 @@ else { if ($print>=2) { print "Advected constituents added by $microphys_pkg microphysics: 10$eol"; } } - if ($zmconv_org == 1 ) { - $nadv += 1; - if ($print>=2) { print "Advected constituents added by $microphys_pkg microphysics: 8$eol"; } - } - if ($clubb_do_adv) { $nadv += 9; if ($print>=2) { print "Advected constituents added by $microphys_pkg microphysics: 8$eol"; } diff --git a/src/physics/cam/convect_shallow.F90 b/src/physics/cam/convect_shallow.F90 index 0a3ba8f4e3..e757e1eef4 100644 --- a/src/physics/cam/convect_shallow.F90 +++ b/src/physics/cam/convect_shallow.F90 @@ -20,7 +20,6 @@ module convect_shallow use phys_control, only : phys_getopts use cloud_fraction_fice, only: cloud_fraction_fice_run use ref_pres, only: trop_cloud_top_lev - use phys_control, only: phys_getopts implicit none private diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index f26f3a5bd1..cf44e2d8fc 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -138,7 +138,7 @@ subroutine zm_conv_register call pbuf_add_field('PREC_DP', 'physpkg',dtype_r8,(/pcols/), prec_dp_idx) call pbuf_add_field('SNOW_DP', 'physpkg',dtype_r8,(/pcols/), snow_dp_idx) - ! detrained convective cloud water mixing ratio. + ! convective mass fluxes call pbuf_add_field('DLFZM', 'physpkg', dtype_r8, (/pcols,pver/), dlfzm_idx) ! detrained convective cloud ice mixing ratio. call pbuf_add_field('CMFMC_DP', 'physpkg', dtype_r8, (/pcols,pverp/), mconzm_idx) @@ -578,7 +578,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & rice(:ncol), lengath, scheme_name, errmsg, errflg) if (errflg /= 0) then - write(str,*) 'From zm_convr_run: at chunk',lchnk + write(str,*) 'From zm_convr_run: at chunk ',lchnk, ' : ' call endrun(str // errmsg) end if @@ -737,7 +737,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & jt(:ncol), maxg(:ncol), ideep(:ncol), 1, lengath, & nstep, ptend_loc%u(:ncol,:), ptend_loc%v(:ncol,:),& pguallu(:ncol,:), pguallv(:ncol,:), pgdallu(:ncol,:), pgdallv(:ncol,:), & - icwuu(:ncol,:), icwuv(:ncol,:), icwdu(:ncol,:), icwdv(:ncol,:), ztodt, seten(:ncol,:) ,& + icwuu(:ncol,:), icwuv(:ncol,:), icwdu(:ncol,:), icwdv(:ncol,:), ztodt, seten(:ncol,:), & scheme_name, errmsg, errflg) call t_stopf ('zm_conv_momtran_run') diff --git a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 index fd08cacfb7..29ae61fc53 100644 --- a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 +++ b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 @@ -59,7 +59,7 @@ end subroutine ccp_get_standard_name !------ subroutine ccp_set_standard_name(this, std_name, errcode, errmsg) - ! Return this constituent's standard name + ! Set this constituent's standard name ! Dummy arguments class(ccpp_constituent_prop_ptr_t), intent(inout) :: this From 94072d1b3284432dcaf0962462d3ff1a304e741c Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 10 Feb 2025 12:22:21 -0700 Subject: [PATCH 17/18] Update comments --- src/physics/cam/zm_conv_intr.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index cf44e2d8fc..44488ac737 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -138,9 +138,9 @@ subroutine zm_conv_register call pbuf_add_field('PREC_DP', 'physpkg',dtype_r8,(/pcols/), prec_dp_idx) call pbuf_add_field('SNOW_DP', 'physpkg',dtype_r8,(/pcols/), snow_dp_idx) - ! convective mass fluxes + ! detrained convective cloud water mixing ratio. call pbuf_add_field('DLFZM', 'physpkg', dtype_r8, (/pcols,pver/), dlfzm_idx) - ! detrained convective cloud ice mixing ratio. + ! convective mass fluxes call pbuf_add_field('CMFMC_DP', 'physpkg', dtype_r8, (/pcols,pverp/), mconzm_idx) From 79aee2f79ffaa56b20817e993a92bfce8b54847c Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Tue, 11 Feb 2025 13:42:55 -0700 Subject: [PATCH 18/18] Updated ChangeLog for cam6_4_065 --- doc/ChangeLog | 88 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index 0aaa217c2a..1ad2135740 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,91 @@ + +=============================================================== + +Tag name: cam6_4_065 +Originator(s): cacraig +Date: Feb 11, 2025 +One-line Summary: ZM CCPP'ization round 4 (completes CCPP conversion of ZM) + +Github PR URL: https://github.com/ESCOMP/CAM/pull/1218 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + - Convert ZM to CCPP and move into atmospheric_physics github repo: https://github.com/ESCOMP/CAM/issues/873 + +Describe any changes made to build system: N/A + +Describe any changes made to the namelist: + - Removed zmconv_org namelist as that partially implemented capability has been removed + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: nusbaume, jimmielin + +List all files eliminated: N/A +D src/physics/cam/wv_sat_methods.F90 +D src/physics/cam/wv_saturation.F90 +D src/utils/error_messages.F90 +D src/utils/namelist_utils.F90 + - Moved to atmospheric_physics (and currently reside in the to_be_ccppized directory) + +List all files added and what they do: N/A + +List all existing files that have been modified, and describe the changes: +M .gitmodules + - update atmospheric_physics to bring in ZM changes + +M bld/build-namelist +M bld/config_files/definition.xml +M bld/configure +M bld/namelist_files/namelist_defaults_cam.xml +M bld/namelist_files/namelist_definition.xml + - remove zmconv_org namelist + +M src/physics/cam/cloud_fraction.F90 + - moved cldfrc_fice to atmospheric_physics and ccppized + +M src/physics/cam/clubb_intr.F90 + - removed difzm declarations as no longer needed + +M src/physics/cam/convect_shallow.F90 +M src/physics/cam/macrop_driver.F90 +M src/physics/cam/physpkg.F90 +M src/physics/cam/rk_stratiform.F90 +M src/physics/cam/zm_conv_intr.F90 +M src/physics/cam7/physpkg.F90 +M src/physics/simple/physpkg.F90 + - various mods to get this to work with the routines that are ccppized + +M src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 + - Add routine: + ccp_set_standard_name to set constituent's standard name + ccp_is_dry to return whether species is dry + ccp_set_dry to set constituent's dry property based on what is passed in + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +derecho/intel/aux_cam: all BFB, except: + ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL) + SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF) + - pre-existing failure due to HEMCO not having reproducible results (issues #1018 and #856) + + SMS_D_Ln9.f19_f19_mg17.FXHIST.derecho_intel.cam-outfrq9s_amie (Overall: FAIL) + SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL) + - pre-existing failures due to build-namelist error requiring CLM/CTSM external update + +derecho/nvhpc/aux_cam: BFB + +izumi/nag/aux_cam: BFB + +izumi/gnu/aux_cam: + ERC_D_Ln9.f10_f10_mg37.QPSPCAMM.izumi_gnu.cam-outfrq3s (Overall: FAIL) details: + - New failure, but since SPCAM is being removed(PR 1217) will create this new pre-existing failure + +=============================================================== =============================================================== Tag name: cam6_4_064