From 9bfcedc9016a773a0587d4649075c4b28fb4cd32 Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Wed, 3 Jul 2019 15:16:06 -0500 Subject: [PATCH 1/7] Creating new branch based off Balwinder's branch. Merging in Brian Eaton's changes to clubb_intr's handling of pdf_params. larson-group/cam#125 --- components/cam/src/physics/cam/clubb_intr.F90 | 41 +++++++++---------- 1 file changed, 19 insertions(+), 22 deletions(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index ae42ec7dfdc9..85efb04ad8af 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -28,6 +28,9 @@ module clubb_intr use perf_mod, only: t_startf, t_stopf use mpishorthand use cam_history_support, only: fillvalue +#ifdef CLUBB_SGS + use clubb_api_module, only: pdf_parameter +#endif implicit none @@ -203,6 +206,11 @@ module clubb_intr integer, parameter :: ncnst=9 character(len=8) :: cnst_names(ncnst) logical :: do_cnst=.false. + +#ifdef CLUBB_SGS + type(pdf_parameter), target, allocatable :: pdf_params_chnk(:,:,:) ! PDF parameters (thermo. levs.) [units vary] + type(pdf_parameter), target, allocatable :: pdf_params_zm_chnk(:,:,:) ! PDF parameters on momentum levs. [units vary] +#endif logical :: liqcf_fix = .FALSE. ! HW for liquid cloud fraction fix logical :: relvar_fix = .FALSE. !PMA for relvar fix @@ -491,7 +499,7 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) ! From CAM libraries use physics_types, only: physics_state, physics_ptend use cam_history, only: addfld, horiz_only, add_default - use ppgrid, only: pver, pverp, pcols + use ppgrid, only: pver, pverp, pcols, begchunk, endchunk use ref_pres, only: pref_mid use hb_diff, only: init_hb_diff use trb_mtn_stress, only: init_tms @@ -569,6 +577,10 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) !----- Begin Code ----- l_do_expldiff_rtm_thlm = do_expldiff + + allocate( & + pdf_params_chnk(pverp,pcols,begchunk:endchunk), & + pdf_params_zm_chnk(pverp,pcols,begchunk:endchunk) ) ! ----------------------------------------------------------------- ! ! Determine how many constituents CLUBB will transport. Note that @@ -984,9 +996,6 @@ subroutine clubb_tend_cam( & stats_rad_zm, & l_output_rad_files, & pdf_parameter, & - num_pdf_params, & - pack_pdf_params_api, & - unpack_pdf_params_api, & stats_begin_timestep_api, & advance_clubb_core_api, & calculate_thlp2_rad_api, & @@ -1176,8 +1185,6 @@ subroutine clubb_tend_cam( & real(r8) :: rtpthvp(pcols,pverp) ! r_t'th_v' (momentum levels) [kg/kg K] real(r8) :: thlpthvp(pcols,pverp) ! th_l'th_v' (momentum levels) [K^2] real(r8) :: sclrpthvp(pcols,pverp,sclr_dim) ! sclr'th_v' (momentum levels) [{units vary} K] - real(r8) :: pdf_params_ptr(pcols,pverp,num_pdf_params) ! putting the params in the pbuf [variable] - real(r8) :: pdf_params_zm_ptr(pcols,pverp,num_pdf_params) ! putting the zm params in the pbuf [variable] real(r8) :: rcm_in_layer(pcols,pverp) ! CLUBB in-cloud liquid water mixing ratio [kg/kg] real(r8) :: cloud_cover(pcols,pverp) ! CLUBB in-cloud cloud fraction [fraction] real(r8) :: wprcp(pcols,pverp) ! CLUBB liquid water flux [m/s kg/kg] @@ -1221,10 +1228,6 @@ subroutine clubb_tend_cam( & integer :: time_elapsed ! time keep track of stats [s] real(r8), dimension(nparams) :: clubb_params ! These adjustable CLUBB parameters (C1, C2 ...) real(r8), dimension(sclr_dim) :: sclr_tol ! Tolerance on passive scalar [units vary] - type(pdf_parameter), dimension(pverp) :: pdf_params ! PDF parameters (thermo. levs.) [units vary] - type(pdf_parameter), dimension(pverp) :: pdf_params_zm ! PDF parameters on momentum levs. [units vary] - real(r8), dimension(pverp,num_pdf_params) :: pdf_params_packed ! Packed for storage in pbuf - real(r8), dimension(pverp,num_pdf_params) :: pdf_params_zm_packed ! Packed for storage in pbuf character(len=200) :: temp1, sub ! Strings needed for CLUBB output logical :: l_Lscale_plume_centered, l_use_ice_latent @@ -1267,7 +1270,10 @@ subroutine clubb_tend_cam( & real(r8), pointer, dimension(:,:) :: accre_enhan ! accretion enhancement factor [-] real(r8), pointer, dimension(:,:) :: cmeliq real(r8), pointer, dimension(:,:) :: cmfmc_sh ! Shallow convective mass flux--m subc (pcols,pverp) [kg/m2/s/] - + + type(pdf_parameter), pointer :: pdf_params(:) ! PDF parameters (thermo. levs.) [units vary] + type(pdf_parameter), pointer :: pdf_params_zm(:) ! PDF parameters on momentum levs. [units vary] + real(r8), pointer, dimension(:,:) :: naai real(r8), pointer, dimension(:,:) :: prer_evap real(r8), pointer, dimension(:,:) :: qrl @@ -1948,10 +1954,8 @@ subroutine clubb_tend_cam( & ! End cloud-top radiative cooling contribution to CLUBB ! ! --------------------------------------------------------- ! - pdf_params_packed(:,:) = pdf_params_ptr(i,:,:) - pdf_params_zm_packed(:,:) = pdf_params_zm_ptr(i,:,:) - call unpack_pdf_params_api(pdf_params_packed, pverp, pdf_params) - call unpack_pdf_params_api(pdf_params_zm_packed, pverp, pdf_params_zm) + pdf_params => pdf_params_chnk(:,i,lchnk) + pdf_params_zm => pdf_params_zm_chnk(:,i,lchnk) call t_startf('adv_clubb_core_ts_loop') do t=1,nadv ! do needed number of "sub" timesteps for each CAM step @@ -2062,13 +2066,6 @@ subroutine clubb_tend_cam( & enddo endif endif - - - ! Pack up pdf_params and store it in the pbuf - call pack_pdf_params_api(pdf_params, pverp, pdf_params_packed) - call pack_pdf_params_api(pdf_params_zm, pverp, pdf_params_zm_packed) - pdf_params_ptr(i,:,:) = pdf_params_packed(:,:) - pdf_params_zm_ptr(i,:,:) = pdf_params_zm_packed(:,:) ! Arrays need to be "flipped" to CAM grid do k=1,pverp From 8fdf2dfe39c5db9289ba15f903f88f1847cf509c Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Wed, 3 Jul 2019 15:46:24 -0500 Subject: [PATCH 2/7] Moving a number of array zeroings that are done in clubb_tend_cam out of the column loop. They are never changed within clubb_intr and are passed to clubb as intent(in), so they should only need to be set once. larson-group/cam#125 --- components/cam/src/physics/cam/clubb_intr.F90 | 112 +++++++++--------- 1 file changed, 54 insertions(+), 58 deletions(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index 85efb04ad8af..1b38c208f790 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -1336,6 +1336,48 @@ subroutine clubb_tend_cam( & else apply_const = 0._r8 ! Never want this if CLUBB's moments are not advected endif + + ! Define forcings from CAM to CLUBB as zero for momentum and thermo, + ! forcings already applied through CAM + thlm_forcing(1:pverp) = 0._r8 + rtm_forcing(1:pverp) = 0._r8 + um_forcing(1:pverp) = 0._r8 + vm_forcing(1:pverp) = 0._r8 + + wprtp_forcing(1:pverp) = 0._r8 + wpthlp_forcing(1:pverp) = 0._r8 + rtp2_forcing(1:pverp) = 0._r8 + thlp2_forcing(1:pverp) = 0._r8 + rtpthlp_forcing(1:pverp) = 0._r8 + + ! rtp3_in and thlp3_in are not currently used in CLUBB's default code. + rtp3_in(:) = 0.0_r8 + thlp3_in(:) = 0.0_r8 + + ! Define surface sources for transported variables for diffusion, will + ! be zero as these tendencies are done in clubb_surface + do ixind=1,edsclr_dim + wpedsclrp_sfc(ixind) = 0._r8 + enddo + + ice_supersat_frac(1:pverp) = 0._r8 + + ! Higher order scalar inputs, set to zero + wpsclrp_sfc(:) = 0._r8 + hydromet(:,:) = 0._r8 + wphydrometp(:,:) = 0._r8 + wp2hmp(:,:) = 0._r8 + rtphmp_zt(:,:) = 0._r8 + thlphmp_zt(:,:) = 0._r8 + + ! Initialize forcings for transported scalars to zero + sclrm_forcing(:,:) = 0._r8 + edsclrm_forcing(:,:) = 0._r8 + + ! Determine Coriolis force at given latitude. This is never used + ! when CLUBB is implemented in a host model, therefore just set + ! to zero. + fcor = 0._r8 ! Get indicees for cloud and ice mass and cloud and ice number @@ -1513,12 +1555,6 @@ subroutine clubb_tend_cam( & ! determine number of timesteps CLUBB core should be advanced, ! host time step divided by CLUBB time step nadv = max(hdtime/dtime,1._r8) - - ! Initialize forcings for transported scalars to zero - - sclrm_forcing(:,:) = 0._r8 - edsclrm_forcing(:,:) = 0._r8 - sclrm(:,:) = 0._r8 minqn = 0._r8 newfice(:,:) = 0._r8 @@ -1650,11 +1686,6 @@ subroutine clubb_tend_cam( & ! CLUBB's budget stats time_elapsed = hdtime - ! Determine Coriolis force at given latitude. This is never used - ! when CLUBB is implemented in a host model, therefore just set - ! to zero. - fcor = 0._r8 - ! Define the CLUBB momentum grid (in height, units of m) do k=1,pverp zi_g(k) = state1%zi(i,pverp-k+1)-state1%zi(i,pver+1) @@ -1756,27 +1787,6 @@ subroutine clubb_tend_cam( & vpwp_sfc = -vm(i,pver)*ustar**2/ubar endif - - ! Define surface sources for transported variables for diffusion, will - ! be zero as these tendencies are done in clubb_surface - do ixind=1,edsclr_dim - wpedsclrp_sfc(ixind) = 0._r8 - enddo - - ! Define forcings from CAM to CLUBB as zero for momentum and thermo, - ! forcings already applied through CAM - thlm_forcing(1:pverp) = 0._r8 - rtm_forcing(1:pverp) = 0._r8 - um_forcing(1:pverp) = 0._r8 - vm_forcing(1:pverp) = 0._r8 - - wprtp_forcing(1:pverp) = 0._r8 - wpthlp_forcing(1:pverp) = 0._r8 - rtp2_forcing(1:pverp) = 0._r8 - thlp2_forcing(1:pverp) = 0._r8 - rtpthlp_forcing(1:pverp) = 0._r8 - - ice_supersat_frac(1:pverp) = 0._r8 ! Set stats output and increment equal to CLUBB and host dt stats_tsamp = dtime @@ -1847,32 +1857,22 @@ subroutine clubb_tend_cam( & if (k .ne. 1) then pre_in(k) = prer_evap(i,pverp-k+1) endif - - ! Initialize these to prevent crashing behavior - wprcp_out(k) = 0._r8 - rcm_in_layer_out(k) = 0._r8 - cloud_cover_out(k) = 0._r8 - edsclr_in(k,:) = 0._r8 - edsclr_out(k,:) = 0._r8 - khzm_out(k) = 0._r8 - khzt_out(k) = 0._r8 - - ! higher order scalar stuff, put to zero - sclrm(k,:) = 0._r8 - wpsclrp(k,:) = 0._r8 - sclrp2(k,:) = 0._r8 - sclrprtp(k,:) = 0._r8 - sclrpthlp(k,:) = 0._r8 - wpsclrp_sfc(:) = 0._r8 - hydromet(k,:) = 0._r8 - wphydrometp(k,:) = 0._r8 - wp2hmp(k,:) = 0._r8 - rtphmp_zt(k,:) = 0._r8 - thlphmp_zt(k,:) = 0._r8 enddo pre_in(1) = pre_in(2) + + ! Initialize these to prevent crashing behavior + edsclr_in(:,:) = 0._r8 + edsclr_out(:,:) = 0._r8 + + ! Higher order scalar inouts, set to zero + sclrm(:,:) = 0._r8 + wpsclrp(:,:) = 0._r8 + sclrp2(:,:) = 0._r8 + sclrprtp(:,:) = 0._r8 + sclrpthlp(:,:) = 0._r8 + if (clubb_do_adv) then if (macmic_it .eq. 1) then @@ -1894,10 +1894,6 @@ subroutine clubb_tend_cam( & enddo endif endif - - ! rtp3_in and thlp3_in are not currently used in CLUBB's default code. - rtp3_in(:) = 0.0_r8 - thlp3_in(:) = 0.0_r8 ! Do the same for tracers icnt=0 From 548df45fe70fee28b9556bbab99bbb36feea134c Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Wed, 3 Jul 2019 15:50:34 -0500 Subject: [PATCH 3/7] Precalculating many divide operations in clubb_tend_cam to reduce computational expense. These are expected to be bit changing. larson-group/cam#125 --- components/cam/src/physics/cam/clubb_intr.F90 | 65 ++++++++++--------- 1 file changed, 35 insertions(+), 30 deletions(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index 1b38c208f790..c19bfb27fc8a 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -1160,6 +1160,9 @@ subroutine clubb_tend_cam( & real(r8) :: qmin real(r8) :: varmu(pcols) real(r8) :: varmu2 + + real(r8) :: invrs_hdtime ! Preculate 1/hdtime to reduce divide operations + real(r8) :: invrs_gravit ! Preculate 1/gravit to reduce divide operations ! Variables below are needed to compute energy integrals for conservation real(r8) :: ke_a(pcols), ke_b(pcols), te_a(pcols), te_b(pcols) @@ -1328,6 +1331,8 @@ subroutine clubb_tend_cam( & !-----------------------------------------------------------------------------------------------! call t_startf('clubb_tend_cam_init') + invrs_hdtime = 1._r8 / hdtime + invrs_gravit = 1._r8 / gravit frac_limit = 0.01_r8 ic_limit = 1.e-12_r8 @@ -1568,7 +1573,7 @@ subroutine clubb_tend_cam( & do k=1,pver do i=1,ncol - exner_clubb(i,k) = 1._r8/((state1%pmid(i,k)/p0_clubb)**(rair/cpair)) + exner_clubb(i,k) = (p0_clubb/state1%pmid(i,k))**(rair/cpair) enddo enddo @@ -1642,10 +1647,10 @@ subroutine clubb_tend_cam( & wl_b = 0._r8 do k=1,pver do i=1,ncol - se_b(i) = se_b(i) + state1%s(i,k)*state1%pdel(i,k)/gravit - ke_b(i) = ke_b(i) + 0.5_r8*(um(i,k)**2+vm(i,k)**2)*state1%pdel(i,k)/gravit - wv_b(i) = wv_b(i) + state1%q(i,k,ixq)*state1%pdel(i,k)/gravit - wl_b(i) = wl_b(i) + state1%q(i,k,ixcldliq)*state1%pdel(i,k)/gravit + se_b(i) = se_b(i) + state1%s(i,k)*state1%pdel(i,k)*invrs_gravit + ke_b(i) = ke_b(i) + 0.5_r8*(um(i,k)**2+vm(i,k)**2)*state1%pdel(i,k)*invrs_gravit + wv_b(i) = wv_b(i) + state1%q(i,k,ixq)*state1%pdel(i,k)*invrs_gravit + wl_b(i) = wl_b(i) + state1%q(i,k,ixcldliq)*state1%pdel(i,k)*invrs_gravit enddo enddo @@ -1708,7 +1713,7 @@ subroutine clubb_tend_cam( & do k=1,pver p_in_Pa(k+1) = state1%pmid(i,pver-k+1) ! Pressure profile exner(k+1) = 1._r8/exner_clubb(i,pver-k+1) - rho_ds_zt(k+1) = (1._r8/gravit)*(state1%pdel(i,pver-k+1)/dz_g(pver-k+1)) + rho_ds_zt(k+1) = invrs_gravit*state1%pdel(i,pver-k+1)/dz_g(pver-k+1) invrs_rho_ds_zt(k+1) = 1._r8/(rho_ds_zt(k+1)) ! Inverse ds rho at thermo rho(i,k+1) = rho_ds_zt(k+1) ! rho on thermo thv_ds_zt(k+1) = thv(i,pver-k+1) ! thetav on thermo @@ -1733,7 +1738,7 @@ subroutine clubb_tend_cam( & ! Compute mean w wind on thermo grid, convert from omega to w wm_zt(1) = 0._r8 do k=1,pver - wm_zt(k+1) = -1._r8*state1%omega(i,pver-k+1)/(rho(i,k+1)*gravit) + wm_zt(k+1) = -1._r8*state1%omega(i,pver-k+1)*invrs_rho_ds_zt(i,k+1)*invrs_gravit enddo ! ------------------------------------------------- ! @@ -1857,7 +1862,7 @@ subroutine clubb_tend_cam( & if (k .ne. 1) then pre_in(k) = prer_evap(i,pverp-k+1) endif - + enddo pre_in(1) = pre_in(2) @@ -2117,10 +2122,10 @@ subroutine clubb_tend_cam( & do k=1,pver clubb_s(k) = cpair*((thlm(i,k)+(latvap/cpair)*rcm(i,k))/exner_clubb(i,k))+ & gravit*state1%zm(i,k)+state1%phis(i) - se_a(i) = se_a(i) + clubb_s(k)*state1%pdel(i,k)/gravit - ke_a(i) = ke_a(i) + 0.5_r8*(um(i,k)**2+vm(i,k)**2)*state1%pdel(i,k)/gravit - wv_a(i) = wv_a(i) + (rtm(i,k)-rcm(i,k))*state1%pdel(i,k)/gravit - wl_a(i) = wl_a(i) + (rcm(i,k))*state1%pdel(i,k)/gravit + se_a(i) = se_a(i) + clubb_s(k)*state1%pdel(i,k)*invrs_gravit + ke_a(i) = ke_a(i) + 0.5_r8*(um(i,k)**2+vm(i,k)**2)*state1%pdel(i,k)*invrs_gravit + wv_a(i) = wv_a(i) + (rtm(i,k)-rcm(i,k))*state1%pdel(i,k)*invrs_gravit + wl_a(i) = wl_a(i) + (rcm(i,k))*state1%pdel(i,k)*invrs_gravit enddo ! Based on these integrals, compute the total energy before and after CLUBB call @@ -2152,11 +2157,11 @@ subroutine clubb_tend_cam( & ! for all variables and therefore is never called in this loop do k=1,pver - ptend_loc%u(i,k) = (um(i,k)-state1%u(i,k))/hdtime ! east-west wind - ptend_loc%v(i,k) = (vm(i,k)-state1%v(i,k))/hdtime ! north-south wind - ptend_loc%q(i,k,ixq) = (rtm(i,k)-rcm(i,k)-state1%q(i,k,ixq))/hdtime ! water vapor - ptend_loc%q(i,k,ixcldliq) = (rcm(i,k)-state1%q(i,k,ixcldliq))/hdtime ! Tendency of liquid water - ptend_loc%s(i,k) = (clubb_s(k)-state1%s(i,k))/hdtime ! Tendency of static energy + ptend_loc%u(i,k) = (um(i,k)-state1%u(i,k))*invrs_hdtime ! east-west wind + ptend_loc%v(i,k) = (vm(i,k)-state1%v(i,k))*invrs_hdtime ! north-south wind + ptend_loc%q(i,k,ixq) = (rtm(i,k)-rcm(i,k)-state1%q(i,k,ixq))*invrs_hdtime ! water vapor + ptend_loc%q(i,k,ixcldliq) = (rcm(i,k)-state1%q(i,k,ixcldliq))*invrs_hdtime ! Tendency of liquid water + ptend_loc%s(i,k) = (clubb_s(k)-state1%s(i,k))*invrs_hdtime ! Tendency of static energy if (clubb_do_adv) then if (macmic_it .eq. cld_macmic_num_steps) then @@ -2169,15 +2174,15 @@ subroutine clubb_tend_cam( & wpthlp(i,k) = wpthlp(i,k) + wpthlp_const wprtp(i,k) = wprtp(i,k) + wprtp_const - ptend_loc%q(i,k,ixthlp2)=(thlp2(i,k)-state1%q(i,k,ixthlp2))/hdtime ! THLP Variance - ptend_loc%q(i,k,ixrtp2)=(rtp2(i,k)-state1%q(i,k,ixrtp2))/hdtime ! RTP Variance - ptend_loc%q(i,k,ixrtpthlp)=(rtpthlp(i,k)-state1%q(i,k,ixrtpthlp))/hdtime ! RTP THLP covariance - ptend_loc%q(i,k,ixwpthlp)=(wpthlp(i,k)-state1%q(i,k,ixwpthlp))/hdtime ! WPTHLP - ptend_loc%q(i,k,ixwprtp)=(wprtp(i,k)-state1%q(i,k,ixwprtp))/hdtime ! WPRTP - ptend_loc%q(i,k,ixwp2)=(wp2(i,k)-state1%q(i,k,ixwp2))/hdtime ! WP2 - ptend_loc%q(i,k,ixwp3)=(wp3(i,k)-state1%q(i,k,ixwp3))/hdtime ! WP3 - ptend_loc%q(i,k,ixup2)=(up2(i,k)-state1%q(i,k,ixup2))/hdtime ! UP2 - ptend_loc%q(i,k,ixvp2)=(vp2(i,k)-state1%q(i,k,ixvp2))/hdtime ! VP2 + ptend_loc%q(i,k,ixthlp2)=(thlp2(i,k)-state1%q(i,k,ixthlp2))*invrs_hdtime ! THLP Variance + ptend_loc%q(i,k,ixrtp2)=(rtp2(i,k)-state1%q(i,k,ixrtp2))*invrs_hdtime ! RTP Variance + ptend_loc%q(i,k,ixrtpthlp)=(rtpthlp(i,k)-state1%q(i,k,ixrtpthlp))*invrs_hdtime ! RTP THLP covariance + ptend_loc%q(i,k,ixwpthlp)=(wpthlp(i,k)-state1%q(i,k,ixwpthlp))*invrs_hdtime ! WPTHLP + ptend_loc%q(i,k,ixwprtp)=(wprtp(i,k)-state1%q(i,k,ixwprtp))*invrs_hdtime ! WPRTP + ptend_loc%q(i,k,ixwp2)=(wp2(i,k)-state1%q(i,k,ixwp2))*invrs_hdtime ! WP2 + ptend_loc%q(i,k,ixwp3)=(wp3(i,k)-state1%q(i,k,ixwp3))*invrs_hdtime ! WP3 + ptend_loc%q(i,k,ixup2)=(up2(i,k)-state1%q(i,k,ixup2))*invrs_hdtime ! UP2 + ptend_loc%q(i,k,ixvp2)=(vp2(i,k)-state1%q(i,k,ixvp2))*invrs_hdtime ! VP2 else ptend_loc%q(i,k,ixthlp2)=0._r8 ptend_loc%q(i,k,ixrtp2)=0._r8 @@ -2205,7 +2210,7 @@ subroutine clubb_tend_cam( & (ixind /= ixrtpthlp) .and. (ixind /= ixwpthlp) .and.& (ixind /= ixwprtp) .and. (ixind /= ixwp2) .and.& (ixind /= ixwp3) .and. (ixind /= ixup2) .and. (ixind /= ixvp2) ) then - ptend_loc%q(i,k,ixind) = (edsclr_out(k,icnt)-state1%q(i,k,ixind))/hdtime ! transported constituents + ptend_loc%q(i,k,ixind) = (edsclr_out(k,icnt)-state1%q(i,k,ixind))*invrs_hdtime ! transported constituents end if end if enddo @@ -2305,8 +2310,8 @@ subroutine clubb_tend_cam( & ! Only rliq is saved from deep convection, which is the reserved liquid. We need to keep ! track of the integrals of ice and static energy that is effected from conversion to ice ! so that the energy checker doesn't complain. - det_s(i) = det_s(i) + ptend_loc%s(i,k)*state1%pdel(i,k)/gravit - det_ice(i) = det_ice(i) - ptend_loc%q(i,k,ixcldice)*state1%pdel(i,k)/gravit + det_s(i) = det_s(i) + ptend_loc%s(i,k)*state1%pdel(i,k)*invrs_gravit + det_ice(i) = det_ice(i) - ptend_loc%q(i,k,ixcldice)*state1%pdel(i,k)*invrs_gravit enddo enddo @@ -2543,7 +2548,7 @@ subroutine clubb_tend_cam( & ! diagnose surface friction and obukhov length (inputs to diagnose PBL depth) do i=1,ncol - rrho = (1._r8/gravit)*(state1%pdel(i,pver)/dz_g(pver)) + rrho = invrs_gravit*(state1%pdel(i,pver)/dz_g(pver)) call calc_ustar( state1%t(i,pver), state1%pmid(i,pver), cam_in%wsx(i), cam_in%wsy(i), & rrho, ustar2(i) ) call calc_obklen( th(i,pver), thv(i,pver), cam_in%cflx(i,1), cam_in%shf(i), rrho, ustar2(i), & From 4568bdb311c5f13097c36c347db28ee37a371bf6 Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Tue, 9 Jul 2019 17:48:27 -0700 Subject: [PATCH 4/7] Merging in newest clubb changes, requires modifying clubb_intr.f90 to handle new pdf_parameter setup. larson-group/E3SM#7 --- components/cam/src/physics/cam/clubb_intr.F90 | 26 +- .../cam/src/physics/clubb/Skx_module.F90 | 39 +- .../clubb/advance_clubb_core_module.F90 | 499 +++++----- .../physics/clubb/advance_wp2_wp3_module.F90 | 59 +- .../physics/clubb/advance_xm_wpxp_module.F90 | 909 +++++++++++++++--- .../physics/clubb/advance_xp2_xpyp_module.F90 | 2 +- .../src/physics/clubb/clubb_api_module.F90 | 114 ++- .../cam/src/physics/clubb/gmres_cache.F90 | 2 +- .../cam/src/physics/clubb/lapack_wrap.F90 | 8 - .../cam/src/physics/clubb/model_flags.F90 | 2 + components/cam/src/physics/clubb/mt95.f90 | 6 +- .../cam/src/physics/clubb/numerical_check.F90 | 2 +- .../src/physics/clubb/parameter_indices.F90 | 8 +- .../src/physics/clubb/pdf_closure_module.F90 | 101 +- .../physics/clubb/pdf_parameter_module.F90 | 547 +++++------ .../physics/clubb/setup_clubb_pdf_params.F90 | 437 +++++---- .../src/physics/clubb/sigma_sqd_w_module.F90 | 106 +- .../physics/clubb/stats_clubb_utilities.F90 | 3 +- .../cam/src/physics/clubb/stats_variables.F90 | 15 +- .../cam/src/physics/clubb/stats_zm_module.F90 | 38 + .../clubb/variables_diagnostic_module.F90 | 13 +- .../clubb/variables_prognostic_module.F90 | 12 +- 22 files changed, 1905 insertions(+), 1043 deletions(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index c19bfb27fc8a..cb3b3f1799a3 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -208,8 +208,8 @@ module clubb_intr logical :: do_cnst=.false. #ifdef CLUBB_SGS - type(pdf_parameter), target, allocatable :: pdf_params_chnk(:,:,:) ! PDF parameters (thermo. levs.) [units vary] - type(pdf_parameter), target, allocatable :: pdf_params_zm_chnk(:,:,:) ! PDF parameters on momentum levs. [units vary] + type(pdf_parameter), target, allocatable :: pdf_params_chnk(:,:) ! PDF parameters (thermo. levs.) [units vary] + type(pdf_parameter), target, allocatable :: pdf_params_zm_chnk(:,:) ! PDF parameters on momentum levs. [units vary] #endif logical :: liqcf_fix = .FALSE. ! HW for liquid cloud fraction fix @@ -525,7 +525,8 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) stats_rad_zm, & w_tol_sqd, & rt_tol, & - l_do_expldiff_rtm_thlm + l_do_expldiff_rtm_thlm, & + init_pdf_params_api use stats_variables, only: l_output_rad_files use units, only: getunit, freeunit @@ -579,8 +580,15 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) l_do_expldiff_rtm_thlm = do_expldiff allocate( & - pdf_params_chnk(pverp,pcols,begchunk:endchunk), & - pdf_params_zm_chnk(pverp,pcols,begchunk:endchunk) ) + pdf_params_chnk(pcols,begchunk:endchunk), & + pdf_params_zm_chnk(pcols,begchunk:endchunk) ) + + do i = begchunk, endchunk + do j = 1, pcols + call init_pdf_params_api( pverp, pdf_params_chnk(i,j) ) + call init_pdf_params_api( pverp, pdf_params_zm_chnk(i,j) ) + end do + end do ! ----------------------------------------------------------------- ! ! Determine how many constituents CLUBB will transport. Note that @@ -1274,8 +1282,8 @@ subroutine clubb_tend_cam( & real(r8), pointer, dimension(:,:) :: cmeliq real(r8), pointer, dimension(:,:) :: cmfmc_sh ! Shallow convective mass flux--m subc (pcols,pverp) [kg/m2/s/] - type(pdf_parameter), pointer :: pdf_params(:) ! PDF parameters (thermo. levs.) [units vary] - type(pdf_parameter), pointer :: pdf_params_zm(:) ! PDF parameters on momentum levs. [units vary] + type(pdf_parameter), pointer :: pdf_params ! PDF parameters (thermo. levs.) [units vary] + type(pdf_parameter), pointer :: pdf_params_zm ! PDF parameters on momentum levs. [units vary] real(r8), pointer, dimension(:,:) :: naai real(r8), pointer, dimension(:,:) :: prer_evap @@ -1955,8 +1963,8 @@ subroutine clubb_tend_cam( & ! End cloud-top radiative cooling contribution to CLUBB ! ! --------------------------------------------------------- ! - pdf_params => pdf_params_chnk(:,i,lchnk) - pdf_params_zm => pdf_params_zm_chnk(:,i,lchnk) + pdf_params => pdf_params_chnk(i,lchnk) + pdf_params_zm => pdf_params_zm_chnk(i,lchnk) call t_startf('adv_clubb_core_ts_loop') do t=1,nadv ! do needed number of "sub" timesteps for each CAM step diff --git a/components/cam/src/physics/clubb/Skx_module.F90 b/components/cam/src/physics/clubb/Skx_module.F90 index 9ffef15ec6be..ddc8a37a0bf1 100644 --- a/components/cam/src/physics/clubb/Skx_module.F90 +++ b/components/cam/src/physics/clubb/Skx_module.F90 @@ -14,7 +14,7 @@ module Skx_module contains !----------------------------------------------------------------------------- - elemental function Skx_func( xp2, xp3, x_tol ) & + function Skx_func( xp2, xp3, x_tol ) & result( Skx ) ! Description: @@ -34,6 +34,9 @@ elemental function Skx_func( xp2, xp3, x_tol ) & use clubb_precision, only: & core_rknd ! Variable(s) + use grid_class, only: & + gr ! Variable Type + implicit none ! External @@ -41,25 +44,33 @@ elemental function Skx_func( xp2, xp3, x_tol ) & ! Parameter Constants ! Whether to apply clipping to the final result - logical, parameter :: & + logical, parameter :: & l_clipping_kluge = .false. ! Input Variables - real( kind = core_rknd ), intent(in) :: & - xp2, & ! [(x units)^2] - xp3, & ! [(x units)^3] + real( kind = core_rknd ), intent(in) :: & x_tol ! x tolerance value [(x units)] + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + xp2, & ! [(x units)^2] + xp3 ! [(x units)^3] + ! Output Variable - real( kind = core_rknd ) :: & + real( kind = core_rknd ), dimension(gr%nz) :: & Skx ! Skewness of x [-] + ! Local Variable + real( kind = core_rknd ) :: & + Skx_denom_tol + ! ---- Begin Code ---- + Skx_denom_tol = Skw_denom_coef * x_tol**2 + !Skx = xp3 / ( max( xp2, x_tol**two ) )**three_halves ! Calculation of skewness to help reduce the sensitivity of this value to ! small values of xp2. - Skx = xp3 / ( xp2 + Skw_denom_coef * x_tol**2 )**three_halves + Skx = xp3 / ( ( xp2 + Skx_denom_tol ) * sqrt( xp2 + Skx_denom_tol ) ) ! This is no longer needed since clipping is already ! imposed on wp2 and wp3 elsewhere in the code @@ -125,12 +136,10 @@ elemental function LG_2005_ansatz( Skw, wpxp, wp2, & ! Larson and Golaz (2005) eq. 16 nrmlzd_corr_wx & - = ( wpxp / ( sqrt( max( wp2, w_tol_sqd ) ) & - * sqrt( max( xp2, x_tol**2 ) ) ) ) & - / sqrt( one - sigma_sqd_w ) + = wpxp / sqrt( max( wp2, w_tol_sqd ) * max( xp2, x_tol**2 ) * ( one - sigma_sqd_w ) ) ! Larson and Golaz (2005) eq. 11 - nrmlzd_Skw = Skw * ( one - sigma_sqd_w)**(-three_halves) + nrmlzd_Skw = Skw / ( ( one - sigma_sqd_w) * sqrt( one - sigma_sqd_w ) ) ! Larson and Golaz (2005) eq. 33 Skx = nrmlzd_Skw * nrmlzd_corr_wx & @@ -188,8 +197,12 @@ function xp3_LG_2005_ansatz( Skw_zt, wpxp_zt, wp2_zt, & ! Local Variable real( kind = core_rknd ), dimension(gr%nz) :: & - Skx_zt ! Skewness of x on thermodynamic levels [-] + Skx_zt, & ! Skewness of x on thermodynamic levels [-] + Skx_denom_tol + + ! ---- Begin Code ---- + Skx_denom_tol = Skw_denom_coef * x_tol**2 ! Calculate skewness of x using the ansatz of LG05. Skx_zt(1:gr%nz) & @@ -198,7 +211,7 @@ function xp3_LG_2005_ansatz( Skw_zt, wpxp_zt, wp2_zt, & ! Calculate using the reverse of the special sensitivity reduction ! formula in function Skx_func above. - xp3 = Skx_zt * ( xp2_zt + Skw_denom_coef * x_tol**2 )**three_halves + xp3 = Skx_zt * ( xp2_zt + Skx_denom_tol ) * sqrt( xp2_zt + Skx_denom_tol ) return diff --git a/components/cam/src/physics/clubb/advance_clubb_core_module.F90 b/components/cam/src/physics/clubb/advance_clubb_core_module.F90 index 6fc97a0ae331..63503139b483 100644 --- a/components/cam/src/physics/clubb/advance_clubb_core_module.F90 +++ b/components/cam/src/physics/clubb/advance_clubb_core_module.F90 @@ -29,7 +29,7 @@ module advance_clubb_core_module ! ! ! -! Cloud Layers Unified By Binormals (CLUBB) user license +! Cloud Layers Unified By Binormals (CLUBB) user license ! agreement. ! ! Thank you for your interest in CLUBB. We work hard to create a @@ -52,39 +52,39 @@ module advance_clubb_core_module ! contributions of both sets of developers. ! ! 3. You may implement CLUBB as a parameterization in a large- -! scale host model that has 2 or 3 spatial dimensions without -! including "CLUBB" in the combined model name, but please -! acknowledge in presentations and publications that CLUBB has +! scale host model that has 2 or 3 spatial dimensions without +! including "CLUBB" in the combined model name, but please +! acknowledge in presentations and publications that CLUBB has ! been included as a parameterization. ! -! 4. You may not provide all or part of CLUBB to anyone without -! prior permission from Vincent Larson (vlarson@uwm.edu). If -! you wish to share CLUBB with your collaborators without -! seeking permission, please ask your collaborators to register -! as CLUBB users at http://clubb.larson-group.com and to +! 4. You may not provide all or part of CLUBB to anyone without +! prior permission from Vincent Larson (vlarson@uwm.edu). If +! you wish to share CLUBB with your collaborators without +! seeking permission, please ask your collaborators to register +! as CLUBB users at http://clubb.larson-group.com and to ! download CLUBB from there. ! -! 5. You may not use CLUBB for commercial purposes unless you +! 5. You may not use CLUBB for commercial purposes unless you ! receive permission from Vincent Larson. ! ! 6. You may not re-license all or any part of CLUBB. ! ! 7. CLUBB is provided "as is" and without warranty. ! -! We hope that CLUBB will develop into a community resource. We -! encourage users to contribute their CLUBB modifications or -! extensions to the CLUBB development group. We will then -! consider them for inclusion in CLUBB. Such contributions will -! benefit all CLUBB users. We would be pleased to acknowledge -! contributors and list their CLUBB-related papers on our "About -! CLUBB" webpage (http://clubb.larson-group.com/about.php) for +! We hope that CLUBB will develop into a community resource. We +! encourage users to contribute their CLUBB modifications or +! extensions to the CLUBB development group. We will then +! consider them for inclusion in CLUBB. Such contributions will +! benefit all CLUBB users. We would be pleased to acknowledge +! contributors and list their CLUBB-related papers on our "About +! CLUBB" webpage (http://clubb.larson-group.com/about.php) for ! those contributors who so desire. ! ! Thanks so much and best wishes for your research! ! ! The CLUBB Development Group -! (Present and past contributors to the source code include -! Vincent Larson, Chris Golaz, David Schanen, Brian Griffin, +! (Present and past contributors to the source code include +! Vincent Larson, Chris Golaz, David Schanen, Brian Griffin, ! Joshua Fasching, Adam Smith, and Michael Falk). !----------------------------------------------------------------------- @@ -137,7 +137,7 @@ subroutine advance_clubb_core & varmu, & ! intent(in) #endif wphydrometp, wp2hmp, rtphmp_zt, thlphmp_zt, & ! intent(in) - host_dx, host_dy, & ! intent(in) + host_dx, host_dy, & ! intent(in) um, vm, upwp, vpwp, up2, vp2, & ! intent(inout) thlm, rtm, wprtp, wpthlp, & ! intent(inout) wp2, wp3, rtp2, rtp3, thlp2, thlp3, rtpthlp, & ! intent(inout) @@ -177,25 +177,27 @@ subroutine advance_clubb_core & ! Modules to be included - use constants_clubb, only: & - em_min, & - thl_tol, & + use constants_clubb, only: & + em_min, & + thl_tol, & rt_tol, & w_tol, & w_tol_sqd, & - ep2, & - Cp, & - Lv, & - ep1, & + ep2, & + Cp, & + Lv, & + ep1, & fstderr, & zero_threshold, & three_halves, & + one_fourth, & one, & unused_var, & grav, & + vonk, & eps - use parameters_tunable, only: & + use parameters_tunable, only: & taumax, & ! Variable(s) c_K, & mu, & @@ -208,7 +210,11 @@ subroutine advance_clubb_core & c_K10h, & C1, C14, & C5, C4, & - C_wp2_splat + C_wp2_splat,& + C_invrs_tau_bkgnd,& + C_invrs_tau_sfc,& + C_invrs_tau_shear,& + C_invrs_tau_N2 use parameters_model, only: & sclr_dim, & ! Variable(s) @@ -216,7 +222,7 @@ subroutine advance_clubb_core & T0, & sclr_tol - use model_flags, only: & + use model_flags, only: & l_tke_aniso, & ! Variable(s) l_call_pdf_closure_twice, & l_host_applies_sfc_fluxes, & @@ -228,16 +234,19 @@ subroutine advance_clubb_core & l_damp_wp2_using_em, & l_advance_xp3, & l_predict_upwp_vpwp, & - l_diag_Lscale_from_tau + l_diag_Lscale_from_tau, & + l_use_C7_Richardson, & + l_use_C11_Richardson, & + l_use_wp3_pr3 - use grid_class, only: & + use grid_class, only: & gr, & ! Variable(s) zm2zt, & ! Procedure(s) - zt2zm, & + zt2zm, & ddzm, & ddzt - use numerical_check, only: & + use numerical_check, only: & parameterization_check, & ! Procedure(s) calculate_spurious_source @@ -257,10 +266,10 @@ subroutine advance_clubb_core & wp2rcp use variables_diagnostic_module, only: & - thvm, & - em, & + thvm, & + em, & Lscale, & - Lscale_up, & + Lscale_up, & Lscale_down, & tau_zm, & tau_zt, & @@ -272,11 +281,11 @@ subroutine advance_clubb_core & vm_ref use variables_diagnostic_module, only: & - wp2_zt, & - thlp2_zt, & - wpthlp_zt, & - wprtp_zt, & - rtp2_zt, & + wp2_zt, & + thlp2_zt, & + wpthlp_zt, & + wprtp_zt, & + rtp2_zt, & rtpthlp_zt, & up2_zt, & vp2_zt, & @@ -285,8 +294,8 @@ subroutine advance_clubb_core & rtm_ref, & thlm_ref - use variables_diagnostic_module, only: & - wpedsclrp, & + use variables_diagnostic_module, only: & + wpedsclrp, & sclrprcp, & ! sclr'rc' wp2sclrp, & ! w'^2 sclr' wpsclrp2, & ! w'sclr'^2 @@ -297,12 +306,12 @@ subroutine advance_clubb_core & a3_coef, & ! The a3 coefficient [-] a3_coef_zt ! The a3 coefficient interp. to the zt grid [-] - use variables_diagnostic_module, only: & + use variables_diagnostic_module, only: & wp3_on_wp2, & ! Variable(s) wp3_on_wp2_zt use pdf_parameter_module, only: & - pdf_parameter, & ! Variable Type + pdf_parameter, & implicit_coefs_terms #ifdef GFDL @@ -312,32 +321,32 @@ subroutine advance_clubb_core & advance_sclrm_Nd_semi_implicit ! h1g, 2010-06-16 end mod #endif - use advance_xm_wpxp_module, only: & + use advance_xm_wpxp_module, only: & advance_xm_wpxp ! Compute mean/flux terms - use advance_xp2_xpyp_module, only: & + use advance_xp2_xpyp_module, only: & advance_xp2_xpyp ! Computes variance terms - use surface_varnce_module, only: & + use surface_varnce_module, only: & calc_surface_varnce ! Procedure - use mixing_length, only: & + use mixing_length, only: & compute_mixing_length ! Procedure - use advance_windm_edsclrm_module, only: & + use advance_windm_edsclrm_module, only: & advance_windm_edsclrm ! Procedure(s) - use saturation, only: & + use saturation, only: & ! Procedure sat_mixrat_liq ! Saturation mixing ratio - use advance_wp2_wp3_module, only: & + use advance_wp2_wp3_module, only: & advance_wp2_wp3 ! Procedure use advance_xp3_module, only: & advance_xp3 ! Procedure(s) - use clubb_precision, only: & + use clubb_precision, only: & core_rknd ! Variable(s) use error_code, only: & @@ -349,7 +358,7 @@ subroutine advance_clubb_core & Skx_func, & ! Procedure(s) xp3_LG_2005_ansatz - use clip_explicit, only: & + use clip_explicit, only: & clip_covars_denom ! Procedure(s) use T_in_K_module, only: & @@ -359,12 +368,12 @@ subroutine advance_clubb_core & use sigma_sqd_w_module, only: & compute_sigma_sqd_w ! Procedure(s) - use stats_clubb_utilities, only: & + use stats_clubb_utilities, only: & stats_accumulate ! Procedure - use stats_type_utilities, only: & + use stats_type_utilities, only: & stat_update_var_pt, & ! Procedure(s) - stat_update_var, & + stat_update_var, & stat_begin_update, & stat_begin_update_pt, & stat_end_update, & @@ -372,12 +381,12 @@ subroutine advance_clubb_core & use stats_variables, only: & irtp2_bt, & ! Variable(s) - ithlp2_bt, & - irtpthlp_bt, & - iwp2_bt, & - iwp3_bt, & - ivp2_bt, & - iup2_bt, & + ithlp2_bt, & + irtpthlp_bt, & + iwp2_bt, & + iwp3_bt, & + ivp2_bt, & + iup2_bt, & iwprtp_bt, & iwpthlp_bt, & iupwp_bt, & @@ -388,7 +397,8 @@ subroutine advance_clubb_core & ium_bt, & irvm, & irel_humidity, & - iwpthlp_zt + iwpthlp_zt , & + itau_zm_simp use stats_variables, only: & iwprtp_zt, & @@ -425,7 +435,7 @@ subroutine advance_clubb_core & use advance_helper_module, only: & calc_stability_correction, & ! Procedure(s) compute_Cx_fnc_Richardson, & - calc_brunt_vaisala_freq_sqd, & + calc_brunt_vaisala_freq_sqd, & term_wp2_splat, term_wp3_splat use interpolation, only: & @@ -438,7 +448,7 @@ subroutine advance_clubb_core & ! Constant Parameters logical, parameter :: & - l_avg_Lscale = .false. ! Lscale is calculated in subroutine compute_mixing_length + l_avg_Lscale = .false. ! Lscale is calculated in subroutine compute_mixing_length ! if l_avg_Lscale is true, compute_mixing_length is called two additional times with ! perturbed values of rtm and thlm. An average value of Lscale ! from the three calls to compute_mixing_length is then calculated. @@ -448,17 +458,17 @@ subroutine advance_clubb_core & l_iter_xp2_xpyp = .true. ! Set to true when rtp2/thlp2/rtpthlp, et cetera are prognostic real( kind = core_rknd ), parameter :: & - tau_const = 900._core_rknd + tau_const = 1000._core_rknd !!! Input Variables - logical, intent(in) :: & - l_implemented ! True if CLUBB is being run within a large-scale host model, + logical, intent(in) :: & + l_implemented ! True if CLUBB is being run within a large-scale host model, ! rather than a standalone single-column model. - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & dt ! Current timestep duration [s] - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & fcor, & ! Coriolis forcing [s^-1] sfc_elevation ! Elevation of ground level [m above MSL] @@ -466,7 +476,7 @@ subroutine advance_clubb_core & hydromet_dim ! Total number of hydrometeor species [#] ! Input Variables - real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & + real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & thlm_forcing, & ! liquid potential temp. forcing (thermodynamic levels) [K/s] rtm_forcing, & ! total water forcing (thermodynamic levels) [(kg/kg)/s] um_forcing, & ! eastward wind forcing (thermodynamic levels) [m/s/s] @@ -528,7 +538,7 @@ subroutine advance_clubb_core & wpedsclrp_sfc ! Eddy-diffusion passive scalar flux at surface [{units vary} m/s] ! Host model horizontal grid spacing, if part of host model. - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & host_dx, & ! East-west horizontal grid spacing [m] host_dy ! North-south horizontal grid spacing [m] @@ -561,7 +571,7 @@ subroutine advance_clubb_core & sclrprtp, & ! sclr'rt' (momentum levels) [{units vary} (kg/kg)] sclrpthlp ! sclr'thl' (momentum levels) [{units vary} K] - real( kind = core_rknd ), intent(inout), dimension(gr%nz) :: & + real( kind = core_rknd ), intent(inout), dimension(gr%nz) :: & rcm, & ! cloud water mixing ratio, r_c (thermo. levels) [kg/kg] cloud_frac, & ! cloud fraction (thermodynamic levels) [-] wpthvp, & ! < w' th_v' > (momentum levels) [kg/kg K] @@ -572,7 +582,7 @@ subroutine advance_clubb_core & real( kind = core_rknd ), intent(inout), dimension(gr%nz,sclr_dim) :: & sclrpthvp ! < sclr' th_v' > (momentum levels) [units vary] - type(pdf_parameter), dimension(gr%nz), intent(inout) :: & + type(pdf_parameter), intent(inout) :: & pdf_params, & ! Fortran structure of PDF parameters on thermodynamic levels [units vary] pdf_params_zm ! Fortran structure of PDF parameters on momentum levels [units vary] @@ -582,13 +592,13 @@ subroutine advance_clubb_core & #endif ! Eddy passive scalar variable - real( kind = core_rknd ), intent(inout), dimension(gr%nz,edsclr_dim) :: & + real( kind = core_rknd ), intent(inout), dimension(gr%nz,edsclr_dim) :: & edsclrm ! Eddy passive scalar grid-mean (thermo. levels) [units vary] ! Variables that need to be output for use in other parts of the CLUBB ! code, such as microphysics (rcm, pdf_params), forcings (rcm), and/or ! BUGSrad (cloud_cover). - real( kind = core_rknd ), intent(out), dimension(gr%nz) :: & + real( kind = core_rknd ), intent(out), dimension(gr%nz) :: & rcm_in_layer, & ! rcm within cloud layer [kg/kg] cloud_cover ! cloud cover [-] @@ -617,14 +627,14 @@ subroutine advance_clubb_core & #ifdef GFDL ! hlg, 2010-06-16 - real( kind = core_rknd ), intent(inout), dimension(gr%nz, min(1,sclr_dim) , 2) :: & + real( kind = core_rknd ), intent(inout), dimension(gr%nz, min(1,sclr_dim) , 2) :: & RH_crit ! critical relative humidity for droplet and ice nucleation ! ---> h1g, 2012-06-14 logical, intent(in) :: do_liquid_only_in_clubb ! <--- h1g, 2012-06-14 #endif ! Local Variables - integer :: i, k + integer :: i, k #ifdef CLUBB_CAM integer :: ixind @@ -636,7 +646,7 @@ subroutine advance_clubb_core & rc_coef_zm ! Coefficient of X'r_c' in Eq. (34) on m-levs. [K/(kg/kg)] real( kind = core_rknd ), dimension(gr%nz) :: & - Km_zm, & ! Eddy diffusivity for momentum on zm grid levels [m^2/s] + Km_zm, & ! Eddy diffusivity for momentum on zm grid levels [m^2/s] Kmh_zm ! Eddy diffusivity for thermodynamic variables [m^2/s] logical, parameter :: & @@ -651,7 +661,7 @@ subroutine advance_clubb_core & gamma_Skw_fnc, & ! Gamma as a function of skewness [-] sigma_sqd_w, & ! PDF width parameter (momentum levels) [-] sigma_sqd_w_zt, & ! PDF width parameter (thermodynamic levels) [-] - sqrt_em_zt, & ! sqrt( em ) on zt levels; where em is TKE [m/s] + sqrt_em_zt, & ! sqrt( em ) on zt levels; where em is TKE [m/s] Lscale_pert_1, Lscale_pert_2, & ! For avg. calculation of Lscale [m] thlm_pert_1, thlm_pert_2, & ! For avg. calculation of Lscale [K] rtm_pert_1, rtm_pert_2, & ! For avg. calculation of Lscale [kg/kg] @@ -704,7 +714,7 @@ subroutine advance_clubb_core & mu_pert_pos_rt, mu_pert_neg_rt ! For l_Lscale_plume_centered !The following variables are defined for use when l_use_ice_latent = .true. - type(pdf_parameter), dimension(gr%nz) :: & + type(pdf_parameter) :: & pdf_params_frz type(implicit_coefs_terms), dimension(gr%nz) :: & @@ -723,12 +733,23 @@ subroutine advance_clubb_core & rel_humidity ! Relative humidity after PDF closure [-] real( kind = core_rknd ), dimension(gr%nz) :: & - stability_correction, & ! Stability correction factor - tau_N2_zm, & ! Tau with a static stability correction applied to it [s] - tau_C6_zm, & ! Tau values used for the C6 (pr1) term in wpxp [s] - tau_C1_zm, & ! Tau values used for the C1 (dp1) term in wp2 [s] - Cx_fnc_Richardson, & ! Cx_fnc computed from Richardson_num [-] - brunt_vaisala_freq_sqd ! Buoyancy frequency squared, N^2 [s^-2} + stability_correction, & ! Stability correction factor + tau_N2_zm, & ! Tau with a static stability correction applied to it [s] + tau_C6_zm, & ! Tau values used for the C6 (pr1) term in wpxp [s] + tau_C1_zm, & ! Tau values used for the C1 (dp1) term in wp2 [s] + Cx_fnc_Richardson, & ! Cx_fnc computed from Richardson_num [-] + brunt_vaisala_freq_sqd, & ! Buoyancy frequency squared, N^2 [s^-2] + invrs_tau_zm, & ! One divided by tau on zm levels [s^-1] + invrs_tau_N2_zm, & ! One divided by tau, including stability effects [s^-1] + ustar, & ! Friction velocity [m/s] + invrs_tau_zm_simp, & + tau_zm_simp, & + tau_zt_simp + + + real( kind = core_rknd ), parameter :: & + ufmin = 0.01_core_rknd, & ! minimum value of friction velocity [m/s] + z_displace = 20.0_core_rknd ! displacement of log law profile above ground [m] real( kind = core_rknd ) :: Lscale_max @@ -741,10 +762,10 @@ subroutine advance_clubb_core & real( kind = core_rknd ), dimension(gr%nz) :: & wp2_splat, & ! Tendency of due to eddies compressing [m^2/s^3] wp3_splat ! Tendency of due to eddies compressing [m^3/s^4] - + ! Variables associated with upgradient momentum contributions due to cumuli !real( kind = core_rknd ), dimension(gr%nz) :: & - ! Km_Skw_factor ! Factor, with value < 1, that reduces eddy diffusivity, + ! Km_Skw_factor ! Factor, with value < 1, that reduces eddy diffusivity, ! Km_zm, in skewed layers !real( kind = core_rknd ),parameter :: & ! Km_Skw_thresh = zero_threshold, & ! Value of Skw at which Skw correction kicks in @@ -795,7 +816,7 @@ subroutine advance_clubb_core & ! Test input variables !---------------------------------------------------------------- if ( clubb_at_least_debug_level( 2 ) ) then - call parameterization_check & + call parameterization_check & ( thlm_forcing, rtm_forcing, um_forcing, & ! intent(in) vm_forcing, wm_zm, wm_zt, p_in_Pa, & ! intent(in) rho_zm, rho, exner, rho_ds_zm, & ! intent(in) @@ -901,7 +922,7 @@ subroutine advance_clubb_core & newmu = varmu #else newmu = mu -#endif +#endif if ( ipdf_call_placement == ipdf_pre_advance_fields & .or. ipdf_call_placement == ipdf_pre_post_advance_fields ) then @@ -986,7 +1007,8 @@ subroutine advance_clubb_core & ! Compute sigma_sqd_w (dimensionless PDF width parameter) sigma_sqd_w & - = compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, wpthlp, wprtp ) + = compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, & + up2, vp2, wpthlp, wprtp, upwp, vpwp ) ! Smooth in the vertical using interpolation sigma_sqd_w = zt2zm( zm2zt( sigma_sqd_w ) ) @@ -1133,7 +1155,7 @@ subroutine advance_clubb_core & * sqrt( max( pdf_params_frz%varnce_thl_1, thl_tol**2 ) ) ) thlm_pert_neg_rt = pdf_params_frz%thl_2 - ( sign_rtpthlp * Lscale_pert_coef & * sqrt( max( pdf_params_frz%varnce_thl_2, thl_tol**2 ) ) ) - rtm_pert_neg_rt = pdf_params_frz%rt_2 & + rtm_pert_neg_rt = pdf_params_frz%rt_2 & - Lscale_pert_coef * sqrt( max( pdf_params_frz%varnce_rt_2, rt_tol**2 ) ) !Lscale_weight = pdf_params%mixt_frac else where @@ -1143,7 +1165,7 @@ subroutine advance_clubb_core & * sqrt( max( pdf_params_frz%varnce_thl_2, thl_tol**2 ) ) ) thlm_pert_neg_rt = pdf_params_frz%thl_1 - ( sign_rtpthlp * Lscale_pert_coef & * sqrt( max( pdf_params_frz%varnce_thl_1, thl_tol**2 ) ) ) - rtm_pert_neg_rt = pdf_params_frz%rt_1 & + rtm_pert_neg_rt = pdf_params_frz%rt_1 & - Lscale_pert_coef * sqrt( max( pdf_params_frz%varnce_rt_1, rt_tol**2 ) ) !Lscale_weight = 1.0_core_rknd - pdf_params%mixt_frac end where @@ -1155,7 +1177,7 @@ subroutine advance_clubb_core & * sqrt( max( pdf_params%varnce_thl_1, thl_tol**2 ) ) ) thlm_pert_neg_rt = pdf_params%thl_2 - ( sign_rtpthlp * Lscale_pert_coef & * sqrt( max( pdf_params%varnce_thl_2, thl_tol**2 ) ) ) - rtm_pert_neg_rt = pdf_params%rt_2 & + rtm_pert_neg_rt = pdf_params%rt_2 & - Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_2, rt_tol**2 ) ) !Lscale_weight = pdf_params%mixt_frac else where @@ -1165,7 +1187,7 @@ subroutine advance_clubb_core & * sqrt( max( pdf_params%varnce_thl_2, thl_tol**2 ) ) ) thlm_pert_neg_rt = pdf_params%thl_1 - ( sign_rtpthlp * Lscale_pert_coef & * sqrt( max( pdf_params%varnce_thl_1, thl_tol**2 ) ) ) - rtm_pert_neg_rt = pdf_params%rt_1 & + rtm_pert_neg_rt = pdf_params%rt_1 & - Lscale_pert_coef * sqrt( max( pdf_params%varnce_rt_1, rt_tol**2 ) ) !Lscale_weight = 1.0_core_rknd - pdf_params%mixt_frac end where @@ -1200,7 +1222,7 @@ subroutine advance_clubb_core & ! This call to compute_mixing_length must be last. Otherwise, the values of ! Lscale_up and Lscale_down in stats will be based on perturbation length scales ! rather than the mean length scale. - + ! Diagnose CLUBB's turbulent mixing length scale. call compute_mixing_length( thvm, thlm, & !intent(in) rtm, em, Lscale_max, p_in_Pa, & !intent(in) @@ -1227,10 +1249,10 @@ subroutine advance_clubb_core & ! Dissipation time !---------------------------------------------------------------- - ! Calculate CLUBB's turbulent eddy-turnover time scale as + ! Calculate CLUBB's turbulent eddy-turnover time scale as ! CLUBB's length scale divided by a velocity scale. tau_zt = MIN( Lscale / sqrt_em_zt, taumax ) - tau_zm = MIN( ( MAX( zt2zm( Lscale ), zero_threshold ) & + tau_zm = MIN( ( MAX( zt2zm( Lscale ), zero_threshold ) & / SQRT( MAX( em_min, em ) ) ), taumax ) ! End Vince Larson's replacement. @@ -1240,13 +1262,36 @@ subroutine advance_clubb_core & call calc_brunt_vaisala_freq_sqd( thlm, exner, rtm, rcm, p_in_Pa, thvm, & ! intent(in) brunt_vaisala_freq_sqd ) ! intent(out) - tau_zt = tau_const / & - ( one + 0.1_core_rknd * tau_const * & - sqrt( max( zero_threshold, brunt_vaisala_freq_sqd ) ) ) - tau_zm = max( zero_threshold, zt2zm( tau_zt ) ) + ustar = max( ( upwp_sfc**2 + vpwp_sfc**2 )**(one_fourth), ufmin ) + + invrs_tau_zm = C_invrs_tau_bkgnd / tau_const & + + C_invrs_tau_sfc * ( ustar / vonk ) / ( gr%zm - sfc_elevation + z_displace ) & + + C_invrs_tau_shear * zt2zm( zm2zt( sqrt( (ddzt( um ))**2 + (ddzt( vm ))**2 ) ) ) & + + C_invrs_tau_N2 * sqrt( max( zero_threshold, & + zt2zm( zm2zt( brunt_vaisala_freq_sqd ) ) - 1e-4_core_rknd) ) + + invrs_tau_zm_simp = C_invrs_tau_bkgnd / tau_const & + + C_invrs_tau_sfc * ( ustar / vonk ) / ( gr%zm - sfc_elevation + z_displace ) & + + C_invrs_tau_N2 * zt2zm( zm2zt( sqrt( (ddzt( um ))**2 + (ddzt( vm ))**2 ) ) ) + + if ( gr%zm(1) - sfc_elevation + z_displace < eps ) then + stop "Lowest zm grid level is below ground in CLUBB." + end if + + tau_zm = one / invrs_tau_zm + tau_zm_simp = one / invrs_tau_zm_simp + + tau_zt = zm2zt( tau_zm ) + + invrs_tau_N2_zm = invrs_tau_zm & + + C_invrs_tau_N2 * sqrt( max( zero_threshold, brunt_vaisala_freq_sqd ) ) + + tau_N2_zm = tau_zm -! tau_zt = tau_const -! tau_zm = tau_const +! tau_zt = tau_const / & +! ( one + 0.1_core_rknd * tau_const * & +! sqrt( max( zero_threshold, brunt_vaisala_freq_sqd ) ) ) +! tau_zm = max( zero_threshold, zt2zm( tau_zt ) ) Lscale = tau_zt * sqrt_em_zt @@ -1269,10 +1314,10 @@ subroutine advance_clubb_core & ! c_K is 0.548 usually (Duynkerke and Driedonks 1987) ! CLUBB uses a smaller value to better fit empirical data. - ! Calculate CLUBB's eddy diffusivity as + ! Calculate CLUBB's eddy diffusivity as ! CLUBB's length scale times a velocity scale. Kh_zt = c_K * Lscale * sqrt_em_zt - Kh_zm = c_K * max( zt2zm( Lscale ), zero_threshold ) & + Kh_zm = c_K * max( zt2zm( Lscale ), zero_threshold ) & * sqrt( max( em, em_min ) ) #if defined(CLUBB_CAM) || defined(GFDL) @@ -1326,7 +1371,7 @@ subroutine advance_clubb_core & wp2(1), up2(1), vp2(1), & ! intent(out) thlp2(1), rtp2(1), rtpthlp(1), & ! intent(out) sclrp2(1,1:sclr_dim), & ! intent(out) - sclrprtp(1,1:sclr_dim), & ! intent(out) + sclrprtp(1,1:sclr_dim), & ! intent(out) sclrpthlp(1,1:sclr_dim) ) ! intent(out) if ( clubb_at_least_debug_level( 0 ) ) then @@ -1437,20 +1482,26 @@ subroutine advance_clubb_core & tau_C1_zm = tau_N2_zm else - tau_N2_zm = unused_var + tau_N2_zm = unused_var tau_C6_zm = tau_zm tau_C1_zm = tau_zm end if ! l_stability_correction - call compute_Cx_Fnc_Richardson( thlm, um, vm, em, Lscale, exner, rtm, & - rcm, p_in_Pa, thvm, rho_ds_zm, & - Cx_fnc_Richardson ) + ! Cx_fnc_Richardson is only used if one of these flags is true, + ! otherwise its value is irrelevant, set it to 0 to avoid NaN problems + if ( l_use_C7_Richardson .or. l_use_C11_Richardson .or. l_use_wp3_pr3 ) then + call compute_Cx_Fnc_Richardson( thlm, um, vm, em, Lscale, exner, rtm, & + rcm, p_in_Pa, thvm, rho_ds_zm, & + Cx_fnc_Richardson ) + else + Cx_fnc_Richardson = 0.0 + end if - ! Advance the prognostic equations for - ! the scalar grid means (rtm, thlm, sclrm) and - ! scalar turbulent fluxes (wprtp, wpthlp, and wpsclrp) - ! by one time step. + ! Advance the prognostic equations for + ! the scalar grid means (rtm, thlm, sclrm) and + ! scalar turbulent fluxes (wprtp, wpthlp, and wpsclrp) + ! by one time step. ! advance_xm_wpxp_bad_wp2 ! Test error comment, DO NOT modify or move call advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & ! intent(in) Lscale, wp3_on_wp2, wp3_on_wp2_zt, Kh_zt, Kh_zm, & ! intent(in) @@ -1502,7 +1553,7 @@ subroutine advance_clubb_core & ! We found that if we call advance_xp2_xpyp first, we can use a longer timestep. - ! Advance the prognostic equations + ! Advance the prognostic equations ! for scalar variances and covariances, ! plus the horizontal wind variances by one time step, by one time step. call advance_xp2_xpyp( tau_zm, wm_zm, rtm, wprtp, thlm, & ! intent(in) @@ -1551,10 +1602,10 @@ subroutine advance_clubb_core & !---------------------------------------------------------------- - ! Advance the 2nd- and 3rd-order moments + ! Advance the 2nd- and 3rd-order moments ! of vertical velocity (wp2, wp3) by one timestep. !---------------------------------------------------------------- - + ! advance_wp2_wp3_bad_wp2 ! Test error comment, DO NOT modify or move call advance_wp2_wp3 & ( dt, sfc_elevation, sigma_sqd_w, wm_zm, & ! intent(in) @@ -1654,20 +1705,20 @@ subroutine advance_clubb_core & if ( l_use_buoy_mod_Km_zm ) then - tau_factor = ( ( one - C5 ) / C4 ) * tau_zm - Km_zm_denom_term = tau_factor * ( grav / T0 ) * & + tau_factor = ( ( one - C5 ) / C4 ) * tau_zm + Km_zm_denom_term = tau_factor * ( grav / T0 ) * & wpthvp / max( 10._core_rknd*w_tol_sqd, wp2 ) Km_zm_numerator_term = 0.02_core_rknd * 0.5_core_rknd * ( grav / T0 ) & - * tau_zm**2 * ddzt( thlm ) + * tau_zm**2 * ddzt( thlm ) Km_zm = c_K10 * tau_factor * wp2 * & - ( one - min( 0.9_core_rknd, Km_zm_numerator_term ) ) / & + ( one - min( 0.9_core_rknd, Km_zm_numerator_term ) ) / & ( one - min( 0.9_core_rknd, Km_zm_denom_term ) ) ! Old method to account for upgradient contribution due to cumuli !Km_Skw_factor = exp( - (Skw_zm - Km_Skw_thresh) / Km_Skw_factor_efold ) !Km_Skw_factor = max( Km_Skw_factor_min, Km_Skw_factor ) - !Km_Skw_factor = min( one, Km_Skw_factor ) + !Km_Skw_factor = min( one, Km_Skw_factor ) !Km_zm = Km_zm * Km_Skw_factor - else + else Km_zm = Kh_zm * c_K10 ! Coefficient for momentum @@ -1678,7 +1729,7 @@ subroutine advance_clubb_core & if ( l_do_expldiff_rtm_thlm ) then edsclrm(:,edsclr_dim-1)=thlm(:) edsclrm(:,edsclr_dim)=rtm(:) - endif + endif call advance_windm_edsclrm( dt, wm_zt, Km_zm, Kmh_zm, ug, vg, um_ref, vm_ref, & ! intent(in) wp2, up2, vp2, um_forcing, vm_forcing, & ! intent(in) @@ -1724,7 +1775,7 @@ subroutine advance_clubb_core & !######################################################################## ! Given CLUBB's prognosed moments, diagnose CLUBB's PDF parameters ! and quantities integrated over that PDF, including - ! quantities related to clouds, buoyancy, and turbulent advection. + ! quantities related to clouds, buoyancy, and turbulent advection. call pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) thlm, wpthlp, rtp2, rtp3, & ! Intent(in) thlp2, thlp3, rtpthlp, wp2, & ! Intent(in) @@ -1798,7 +1849,7 @@ subroutine advance_clubb_core & endif ! l_predict_upwp_vpwp call stat_end_update( irtp2_bt, rtp2 / dt, & ! intent(in) stats_zm ) ! intent(inout) - call stat_end_update( ithlp2_bt, thlp2 / dt, & ! intent(in) + call stat_end_update( ithlp2_bt, thlp2 / dt, & ! intent(in) stats_zm ) ! intent(inout) call stat_end_update( irtpthlp_bt, rtpthlp / dt, & ! intent(in) stats_zm ) ! intent(inout) @@ -1841,7 +1892,7 @@ subroutine advance_clubb_core & vpwp_zt = zm2zt( vpwp ) end if - call stats_accumulate & + call stats_accumulate & ( um, vm, upwp, vpwp, up2, vp2, & ! intent(in) thlm, rtm, wprtp, wpthlp, & ! intent(in) wp2, wp3, rtp2, rtp3, thlp2, thlp3, rtpthlp, & ! intent(in) @@ -1858,7 +1909,7 @@ subroutine advance_clubb_core & if ( clubb_at_least_debug_level( 2 ) ) then - call parameterization_check & + call parameterization_check & ( thlm_forcing, rtm_forcing, um_forcing, & ! intent(in) vm_forcing, wm_zm, wm_zt, p_in_Pa, & ! intent(in) rho_zm, rho, exner, rho_ds_zm, & ! intent(in) @@ -2038,7 +2089,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) use T_in_K_module, only: & thlm2T_in_K ! Procedure(s) - use saturation, only: & + use saturation, only: & sat_mixrat_liq ! Procedure(s) use array_index, only: & @@ -2113,14 +2164,14 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! respect to ice (zero for liquid) !!! Input Variables - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & dt ! Current timestep duration [s] integer, intent(in) :: & hydromet_dim ! Total number of hydrometeors [#] ! Input Variables - real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & ! rtm, & ! total water mixing ratio, r_t (thermo. levels) [kg/kg] wprtp, & ! w' r_t' (momentum levels) [(kg/kg)m/s] thlm, & ! liq. water pot. temp., th_l (thermo. levels) [K] @@ -2169,12 +2220,12 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) logical, intent(in) :: & l_samp_stats_in_pdf_call ! Sample stats in this call to this subroutine - real( kind = core_rknd ), dimension(gr%nz), intent(inout) :: & + real( kind = core_rknd ), dimension(gr%nz), intent(inout) :: & rtm ! total water mixing ratio, r_t (thermo. levels) [kg/kg] #ifdef GFDL ! hlg, 2010-06-16 - real( kind = core_rknd ), dimension(gr%nz, min(1,sclr_dim) , 2), intent(inout) :: & + real( kind = core_rknd ), dimension(gr%nz, min(1,sclr_dim) , 2), intent(inout) :: & RH_crit ! critical relative humidity for droplet and ice nucleation ! ---> h1g, 2012-06-14 logical, intent(in) :: do_liquid_only_in_clubb @@ -2183,7 +2234,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) !!! Output Variables ! Variables being passed back to and out of advance_clubb_core. - real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & + real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & rcm, & ! mean r_c (thermodynamic levels) [kg/kg] cloud_frac, & ! cloud fraction (thermodynamic levels) [-] ice_supersat_frac, & ! ice supersat. frac. (thermo. levels) [-] @@ -2199,7 +2250,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) rcp2_zt, & ! r_c'^2 (on thermo. grid) [kg^2/kg^2] thlprcp, & ! < th_l' r_c' > (momentum levels) [K kg/kg] rc_coef_zm, & ! Coefficient of X'r_c' on m-levs. [K/(kg/kg)] - rtm_frz, & ! rtm adjusted to include hydrometeors [kg/kg] + rtm_frz, & ! rtm adjusted to include hydrometeors [kg/kg] thlm_frz ! thlm adjusted to include hydrometeors [K] ! Variable being passed back to and out of advance_clubb_core. @@ -2207,7 +2258,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) sclrpthvp ! < sclr' th_v' > (momentum levels) [units vary] ! Variables being passed back to only advance_clubb_core (for statistics). - real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & + real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & wp4, & ! < w'^4 > (momentum levels) [m^4/s^4] wp2rtp, & ! < w'^2 r_t' > (thermodynamic levels) [m^2/s^2 kg/kg] wprtp2, & ! < w' r_t'^2 > (thermodynamic levels) [m/s kg^2/kg^2] @@ -2223,7 +2274,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) vprcp ! < v' r_c' > [(m kg)/(s kg)] ! Variables being passed back to only advance_clubb_core (for statistics). - real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & + real( kind = core_rknd ), dimension(gr%nz), intent(out) :: & Skw_velocity, & ! Skewness velocity [m/s] cloud_frac_zm, & ! Cloud Fraction on momentum levels [-] ice_supersat_frac_zm, & ! Ice supersat. frac. on momentum levels [-] @@ -2240,12 +2291,12 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) wpsclrpthlp ! < w' sclr' th_l' > (thermodynamic levels) [units vary] ! Variable being passed back to and out of advance_clubb_core. - type(pdf_parameter), dimension(gr%nz), intent(out) :: & + type(pdf_parameter), intent(inout) :: & pdf_params, & ! PDF parameters [units vary] pdf_params_frz ! Output for use in pert. Lscale calc. ! Variable being passed back to only advance_clubb_core. - type(pdf_parameter), dimension(gr%nz), intent(out) :: & + type(pdf_parameter), intent(inout) :: & pdf_params_zm ! PDF parameters [units vary] type(implicit_coefs_terms), dimension(gr%nz), intent(out) :: & @@ -2273,13 +2324,13 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) p_in_Pa_zm, & ! Pressure interpolated to momentum levels [Pa] exner_zm ! Exner interpolated to momentum levels [-] - real( kind = core_rknd ), dimension(gr%nz,hydromet_dim) :: & + real( kind = core_rknd ), dimension(gr%nz,hydromet_dim) :: & wphydrometp_zt, & ! Covariance of w and hm (on t-levs.) [(m/s) ] wp2hmp_zm, & ! Moment (on m-levs.) [(m/s)^2 ] rtphmp, & ! Covariance of rt and hm [(kg/kg) ] thlphmp ! Covariance of thl and hm [K ] - real( kind = core_rknd ), dimension(gr%nz,sclr_dim) :: & + real( kind = core_rknd ), dimension(gr%nz,sclr_dim) :: & wpsclrp_zt, & ! w' sclr' on thermo. levels sclrp2_zt, & ! sclr'^2 on thermo. levels sclrprtp_zt, & ! sclr' r_t' on thermo. levels @@ -2289,21 +2340,21 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! momentum grid levels, but pdf_closure outputs them on the thermodynamic ! grid levels. real( kind = core_rknd ), dimension(gr%nz) :: & - wp4_zt, & ! w'^4 (on thermo. grid) [m^4/s^4] + wp4_zt, & ! w'^4 (on thermo. grid) [m^4/s^4] wpthvp_zt, & ! Buoyancy flux (on thermo. grid) [(K m)/s] rtpthvp_zt, & ! r_t' th_v' (on thermo. grid) [(kg K)/kg] thlpthvp_zt, & ! th_l' th_v' (on thermo. grid) [K^2] - wprcp_zt, & ! w' r_c' (on thermo. grid) [(m kg)/(s kg)] - rtprcp_zt, & ! r_t' r_c' (on thermo. grid) [(kg^2)/(kg^2)] + wprcp_zt, & ! w' r_c' (on thermo. grid) [(m kg)/(s kg)] + rtprcp_zt, & ! r_t' r_c' (on thermo. grid) [(kg^2)/(kg^2)] thlprcp_zt, & ! th_l' r_c' (on thermo. grid) [(K kg)/kg] uprcp_zt, & ! u' r_c' (on thermo. grid) [(m kg)/(s kg)] vprcp_zt ! v' r_c' (on thermo. grid) [(m kg)/(s kg)] - real( kind = core_rknd ), dimension(gr%nz, sclr_dim) :: & - sclrpthvp_zt, & ! sclr'th_v' (on thermo. grid) + real( kind = core_rknd ), dimension(gr%nz, sclr_dim) :: & + sclrpthvp_zt, & ! sclr'th_v' (on thermo. grid) sclrprcp_zt ! sclr'rc' (on thermo. grid) - real( kind = core_rknd ), dimension(gr%nz) :: & + real( kind = core_rknd ), dimension(gr%nz) :: & wprtp2_zm, & ! < w' r_t'^2 > on momentum levels [m/s kg^2/kg^2] wp2rtp_zm, & ! < w'^2 r_t' > on momentum levels [m^2/s^2 kg/kg] wpthlp2_zm, & ! < w' th_l'^2 > on momentum levels [m/s K^2] @@ -2312,10 +2363,10 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) wp2thvp_zm, & ! < w'^2 th_v' > on momentum levels [m^2/s^2 K] wp2rcp_zm ! < w'^2 r_c' > on momentum levles [m^2/s^2 kg/kg] - real( kind = core_rknd ), dimension(gr%nz,sclr_dim) :: & - wpsclrprtp_zm, & ! w'sclr'rt' on momentum grid - wpsclrp2_zm, & ! w'sclr'^2 on momentum grid - wpsclrpthlp_zm, & ! w'sclr'thl' on momentum grid + real( kind = core_rknd ), dimension(gr%nz,sclr_dim) :: & + wpsclrprtp_zm, & ! w'sclr'rt' on momentum grid + wpsclrp2_zm, & ! w'sclr'^2 on momentum grid + wpsclrpthlp_zm, & ! w'sclr'thl' on momentum grid wp2sclrp_zm, & ! w'^2 sclr' on momentum grid sclrm_zm ! Passive scalar mean on momentum grid @@ -2350,7 +2401,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) pdf_implicit_coefs_terms_zm_frz ! The following variables are defined for use when l_use_ice_latent = .true. - type(pdf_parameter), dimension(gr%nz) :: & + type(pdf_parameter ):: & pdf_params_zm_frz real( kind = core_rknd ), dimension(gr%nz) :: & @@ -2452,10 +2503,6 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) integer :: i, k - ! Initialize Variables - call init_pdf_params( gr%nz, pdf_params_zm ) - call init_pdf_params( gr%nz, pdf_params_frz ) - !--------------------------------------------------------------------------- ! Interpolate wp3, rtp3, and thlp3 to momentum levels, and wp2, rtp2, and ! thlp2 to thermodynamic levels, and then compute Skw, Skrt, and Skthl for @@ -2512,7 +2559,8 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) end if ! Compute sigma_sqd_w (dimensionless PDF width parameter) - sigma_sqd_w = compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, wpthlp, wprtp ) + sigma_sqd_w = compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, & + up2, vp2, wpthlp, wprtp, upwp, vpwp ) if ( l_stats_samp .and. l_samp_stats_in_pdf_call ) then call stat_update_var( igamma_Skw_fnc, gamma_Skw_fnc, & ! intent(in) @@ -2537,7 +2585,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! Compute skewness velocity for stats output purposes if ( iSkw_velocity > 0 ) then - Skw_velocity = ( 1.0_core_rknd / ( 1.0_core_rknd - sigma_sqd_w(1:gr%nz) ) ) & + Skw_velocity = ( 1.0_core_rknd / ( 1.0_core_rknd - sigma_sqd_w(1:gr%nz) ) ) & * ( wp3_zm(1:gr%nz) / max( wp2(1:gr%nz), w_tol_sqd ) ) end if @@ -2559,7 +2607,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) enddo ! i = 1, hydromet_dim, 1 - call pdf_closure & + call pdf_closure & ( hydromet_dim, p_in_Pa, exner, thv_ds_zt, & ! intent(in) wm_zt, wp2_zt, wp3, sigma_sqd_w_zt, & ! intent(in) Skw_zt, Skthl_zt, Skrt_zt, & ! intent(in) @@ -2618,7 +2666,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! of subgrid clouds do k=1, gr%nz - if ( pdf_params(k)%chi_1/pdf_params(k)%stdev_chi_1 > -1._core_rknd ) then + if ( pdf_params%chi_1(k)/pdf_params%stdev_chi_1(k) > -1._core_rknd ) then ! Recalculate cloud_frac and r_c for each PDF component @@ -2627,29 +2675,29 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) pdf_params%stdev_chi_1, (/(chi_at_liq_sat,i=1,gr%nz)/), & ! Intent(in) cloud_frac_1_refined, rc_1_refined ) ! Intent(out) - call calc_vert_avg_cf_component & + call calc_vert_avg_cf_component & ( gr%nz, k, gr%zt, pdf_params%chi_2, & ! Intent(in) pdf_params%stdev_chi_2, (/(chi_at_liq_sat,i=1,gr%nz)/), & ! Intent(in) cloud_frac_2_refined, rc_2_refined ) ! Intent(out) cloud_frac_refined = compute_mean_binormal & ( cloud_frac_1_refined, cloud_frac_2_refined, & - pdf_params(k)%mixt_frac ) + pdf_params%mixt_frac(k) ) rcm_refined = compute_mean_binormal & - ( rc_1_refined, rc_2_refined, pdf_params(k)%mixt_frac ) + ( rc_1_refined, rc_2_refined, pdf_params%mixt_frac(k) ) if ( l_interactive_refined ) then ! I commented out the lines that modify the values in pdf_params, as it seems that ! these values need to remain consistent with the rest of the PDF. ! Eric Raut Jun 2014 ! Replace pdf_closure estimates with refined estimates - ! pdf_params(k)%rc_1 = rc_1_refined - ! pdf_params(k)%rc_2 = rc_2_refined + ! pdf_params%rc_1(k) = rc_1_refined + ! pdf_params%rc_2(k) = rc_2_refined rcm(k) = rcm_refined - ! pdf_params(k)%cloud_frac_1 = cloud_frac_1_refined - ! pdf_params(k)%cloud_frac_2 = cloud_frac_2_refined + ! pdf_params%cloud_frac_1(k) = cloud_frac_1_refined + ! pdf_params%cloud_frac_2(k) = cloud_frac_2_refined cloud_frac(k) = cloud_frac_refined end if @@ -2658,7 +2706,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! output to stats! cloud_frac_refined = cloud_frac(k) rcm_refined = rcm(k) - end if ! pdf_params(k)%chi_1/pdf_params(k)%stdev_chi_1 > -1._core_rknd + end if ! pdf_params%chi_1(k)/pdf_params%stdev_chi_1(k) > -1._core_rknd ! Stats output if ( l_stats_samp .and. l_samp_stats_in_pdf_call ) then @@ -2721,7 +2769,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) ! Call pdf_closure to output the variables which belong on the momentum grid. - call pdf_closure & + call pdf_closure & ( hydromet_dim, p_in_Pa_zm, exner_zm, thv_ds_zm, & ! intent(in) wm_zm, wp2, wp3_zm, sigma_sqd_w, & ! intent(in) Skw_zm, Skthl_zm, Skrt_zm, & ! intent(in) @@ -2739,7 +2787,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) rtphmp, thlphmp, & ! intent(in) wp4, wprtp2_zm, wp2rtp_zm, & ! intent(out) wpthlp2_zm, wp2thlp_zm, wprtpthlp_zm, & ! intent(out) - cloud_frac_zm, ice_supersat_frac_zm, & ! intent(out) + cloud_frac_zm, ice_supersat_frac_zm, & ! intent(out) rcm_zm, wpthvp, wp2thvp_zm, rtpthvp, & ! intent(out) thlpthvp, wprcp, wp2rcp_zm, rtprcp, & ! intent(out) thlprcp, rcp2, & ! intent(out) @@ -2883,7 +2931,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) !of latent heating due to ice. Thlm and rtm add the effects of ice, and !the terms are all renamed with "_frz" appended. The modified terms will !be fed into the calculations of the turbulence terms. storer-3/14/13 - + !Also added rain for completeness. storer-3/4/14 if ( iirr > 0 ) then @@ -2892,11 +2940,11 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) rrm = zero end if - thlm_frz = thlm - (Lv / (Cp*exner) ) * rrm - (Ls / (Cp*exner) ) * rfrzm + thlm_frz = thlm - (Lv / (Cp*exner) ) * rrm - (Ls / (Cp*exner) ) * rfrzm rtm_frz = rtm + rrm + rfrzm - call pdf_closure & + call pdf_closure & ( hydromet_dim, p_in_Pa, exner, thv_ds_zt, & ! intent(in) wm_zt, wp2_zt, wp3, sigma_sqd_w_zt, & ! intent(in) Skw_zt, Skthl_zt, Skrt_zt, & ! intent(in) @@ -2952,8 +3000,11 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) thlm_zm_frz(gr%nz) = max( thlm_zm_frz(gr%nz), thl_tol ) if ( l_call_pdf_closure_twice ) then + + call init_pdf_params( gr%nz, pdf_params_zm_frz ) + ! Call pdf_closure again to output the variables which belong on the momentum grid. - call pdf_closure & + call pdf_closure & ( hydromet_dim, p_in_Pa_zm, exner_zm, thv_ds_zm, & ! intent(in) wm_zm, wp2, wp3_zm, sigma_sqd_w, & ! intent(in) Skw_zm, Skthl_zm, Skrt_zm, & ! intent(in) @@ -2971,7 +3022,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) rtphmp, thlphmp, & ! intent(in) wp4_frz, wprtp2_zm_frz, wp2rtp_zm_frz, & ! intent(out) wpthlp2_zm_frz, wp2thlp_zm_frz, wprtpthlp_zm_frz, & ! intent(out) - cloud_frac_zm_frz, ice_supersat_frac_zm_frz, & ! intent(out) + cloud_frac_zm_frz, ice_supersat_frac_zm_frz, & ! intent(out) rcm_zm_frz, wpthvp_frz, wp2thvp_zm_frz, rtpthvp_frz, & ! intent(out) thlpthvp_frz, wprcp_frz, wp2rcp_zm_frz, rtprcp_frz, & ! intent(out) thlprcp_frz, rcp2_frz, & ! intent(out) @@ -3035,12 +3086,12 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) end if ! l_use_ice_latent = .true. rsat = sat_mixrat_liq( p_in_Pa, thlm2T_in_K( thlm, exner, rcm ) ) - rel_humidity = (rtm - rcm) / rsat + rel_humidity = (rtm - rcm) / rsat rcm_supersat_adj = zero if ( l_rcm_supersat_adj ) then ! +PAB mods, take remaining supersaturation that may exist - ! after CLUBB PDF call and add it to rcm. Supersaturation + ! after CLUBB PDF call and add it to rcm. Supersaturation ! may exist after PDF call due to issues with calling PDF on the ! thermo grid and momentum grid and the interpolation between the two l_spur_supersat = .false. @@ -3053,7 +3104,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) enddo if ( clubb_at_least_debug_level( 1 ) .and. l_spur_supersat ) then - write(fstderr,*) 'Warning: spurious supersaturation was removed after pdf_closure!' + write(fstderr,*) 'Warning: spurious supersaturation was removed after pdf_closure!' end if end if ! l_rcm_supersat_adj @@ -3064,7 +3115,7 @@ subroutine pdf_closure_driver( dt, hydromet_dim, wprtp, & ! Intent(in) end subroutine pdf_closure_driver !============================================================================= - subroutine setup_clubb_core & + subroutine setup_clubb_core & ( nzmax, T0_in, ts_nudge_in, & ! intent(in) hydromet_dim_in, sclr_dim_in, & ! intent(in) sclr_tol_in, edsclr_dim_in, params, & ! intent(in) @@ -3090,35 +3141,35 @@ subroutine setup_clubb_core & ! None !--------------------------------------------------------------------- - use grid_class, only: & + use grid_class, only: & setup_grid, & ! Procedure gr ! Variable(s) - use parameter_indices, only: & + use parameter_indices, only: & nparams ! Variable(s) - use parameters_tunable, only: & + use parameters_tunable, only: & setup_parameters ! Procedure - use parameters_model, only: & + use parameters_model, only: & setup_parameters_model ! Procedure - use variables_diagnostic_module, only: & + use variables_diagnostic_module, only: & setup_diagnostic_variables ! Procedure - use variables_prognostic_module, only: & + use variables_prognostic_module, only: & setup_prognostic_variables ! Procedure - use constants_clubb, only: & + use constants_clubb, only: & fstderr ! Variable(s) use error_code, only: & clubb_at_least_debug_level, & ! Procedures - initialize_error_headers, & + initialize_error_headers, & err_code, & ! Error Indicator clubb_fatal_error ! Constant - use model_flags, only: & + use model_flags, only: & setup_model_flags, & ! Subroutine l_use_ice_latent, & ! Variable(s) l_predict_upwp_vpwp, & @@ -3182,7 +3233,7 @@ subroutine setup_clubb_core & ! evenly-spaced grid (grid_type = 1), it needs the vertical ! grid spacing, momentum-level starting altitude, and maximum ! altitude as input. - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & deltaz, & ! Change in altitude per level [m] zm_init, & ! Initial grid altitude (momentum level) [m] zm_top ! Maximum grid altitude (momentum level) [m] @@ -3196,23 +3247,23 @@ subroutine setup_clubb_core & ! If the CLUBB model is running by itself, but is using a ! stretched grid entered on momentum levels (grid_type = 3), ! it needs to use the momentum level altitudes as input. - real( kind = core_rknd ), intent(in), dimension(nzmax) :: & + real( kind = core_rknd ), intent(in), dimension(nzmax) :: & momentum_heights, & ! Momentum level altitudes (input) [m] thermodynamic_heights ! Thermodynamic level altitudes (input) [m] ! Model parameters - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & T0_in, ts_nudge_in - integer, intent(in) :: & + integer, intent(in) :: & hydromet_dim_in, & ! Number of hydrometeor species sclr_dim_in, & ! Number of passive scalars edsclr_dim_in ! Number of eddy-diff. passive scalars - real( kind = core_rknd ), intent(in), dimension(sclr_dim_in) :: & + real( kind = core_rknd ), intent(in), dimension(sclr_dim_in) :: & sclr_tol_in ! Thresholds for passive scalars - real( kind = core_rknd ), intent(in), dimension(nparams) :: & + real( kind = core_rknd ), intent(in), dimension(nparams) :: & params ! Including C1, nu1, nu2, etc. ! Flags @@ -3230,7 +3281,7 @@ subroutine setup_clubb_core & logical, intent(in) :: & ! h1g, 2010-06-16 begin mod I_sat_sphum - real( kind = core_rknd ), intent(in) :: & + real( kind = core_rknd ), intent(in) :: & cloud_frac_min ! h1g, 2010-06-16 end mod #endif @@ -3395,7 +3446,7 @@ subroutine setup_clubb_core & // " because uprcp has not been fed through advance_clubb_core." err_code = clubb_fatal_error return - endif + endif endif ! l_predict_upwp_vpwp @@ -3416,7 +3467,7 @@ subroutine setup_clubb_core & write(fstderr,*) "zm_init = ", zm_init write(fstderr,*) "zm_top = ", zm_top write(fstderr,*) "momentum_heights = ", momentum_heights - write(fstderr,*) "thermodynamic_heights = ", & + write(fstderr,*) "thermodynamic_heights = ", & thermodynamic_heights write(fstderr,*) "T0_in = ", T0_in write(fstderr,*) "ts_nudge_in = ", ts_nudge_in @@ -3428,13 +3479,13 @@ subroutine setup_clubb_core & ! Setup flags #ifdef GFDL - call setup_model_flags & + call setup_model_flags & ( l_host_applies_sfc_fluxes, & ! intent(in) - l_uv_nudge, saturation_formula, & ! intent(in) + l_uv_nudge, saturation_formula, & ! intent(in) I_sat_sphum ) ! intent(in) h1g, 2010-06-16 #else - call setup_model_flags & + call setup_model_flags & ( l_host_applies_sfc_fluxes, & ! intent(in) l_uv_nudge, saturation_formula ) ! intent(in) #endif @@ -3453,7 +3504,7 @@ subroutine setup_clubb_core & #endif ! Define tunable constant parameters - call setup_parameters & + call setup_parameters & ( deltaz, params, gr%nz, & ! intent(in) grid_type, momentum_heights(begin_height:end_height), & ! intent(in) thermodynamic_heights(begin_height:end_height) ) ! intent(in) @@ -3469,7 +3520,7 @@ subroutine setup_clubb_core & write(fstderr,*) "zm_init = ", zm_init write(fstderr,*) "zm_top = ", zm_top write(fstderr,*) "momentum_heights = ", momentum_heights - write(fstderr,*) "thermodynamic_heights = ", & + write(fstderr,*) "thermodynamic_heights = ", & thermodynamic_heights write(fstderr,*) "T0_in = ", T0_in write(fstderr,*) "ts_nudge_in = ", ts_nudge_in @@ -3520,10 +3571,10 @@ subroutine cleanup_clubb_core( l_implemented ) !--------------------------------------------------------------------------- use parameters_model, only: sclr_tol ! Variable - use variables_diagnostic_module, only: & + use variables_diagnostic_module, only: & cleanup_diagnostic_variables ! Procedure - use variables_prognostic_module, only: & + use variables_prognostic_module, only: & cleanup_prognostic_variables ! Procedure use grid_class, only: & @@ -3645,12 +3696,12 @@ subroutine trapezoidal_rule_zt & rcm, & ! Liquid water mixing ratio [kg/kg] wp2thvp ! w'^2 th_v' [m^2 K/s^2] - real( kind = core_rknd ), dimension(gr%nz,sclr_dim), intent(inout) :: & - wpsclrprtp, & ! w'sclr'rt' + real( kind = core_rknd ), dimension(gr%nz,sclr_dim), intent(inout) :: & + wpsclrprtp, & ! w'sclr'rt' wpsclrp2, & ! w'sclr'^2 wpsclrpthlp ! w'sclr'thl' - type (pdf_parameter), dimension(gr%nz), intent(inout) :: & + type (pdf_parameter), intent(inout) :: & pdf_params ! PDF parameters [units vary] ! Thermo. level variables brought to momentum levels either by @@ -3665,12 +3716,12 @@ subroutine trapezoidal_rule_zt & rcm_zm, & ! Liquid water mixing ratio on momentum grid [kg/kg] wp2thvp_zm ! w'^2 th_v' on momentum grid [m^2 K/s^2] - real( kind = core_rknd ), dimension(gr%nz,sclr_dim), intent(inout) :: & - wpsclrprtp_zm, & ! w'sclr'rt' on momentum grid - wpsclrp2_zm, & ! w'sclr'^2 on momentum grid + real( kind = core_rknd ), dimension(gr%nz,sclr_dim), intent(inout) :: & + wpsclrprtp_zm, & ! w'sclr'rt' on momentum grid + wpsclrp2_zm, & ! w'sclr'^2 on momentum grid wpsclrpthlp_zm ! w'sclr'thl' on momentum grid - type (pdf_parameter), dimension(gr%nz), intent(inout) :: & + type (pdf_parameter), intent(inout) :: & pdf_params_zm ! PDF parameters on momentum grid [units vary] ! Local variables @@ -4241,7 +4292,7 @@ subroutine compute_cloud_cover & cloud_frac, & ! Cloud fraction [-] rcm ! Liquid water mixing ratio [kg/kg] - type (pdf_parameter), dimension(gr%nz), intent(in) :: & + type (pdf_parameter), intent(in) :: & pdf_params ! PDF Parameters [units vary] ! Output variables @@ -4251,7 +4302,7 @@ subroutine compute_cloud_cover & ! Local variables real( kind = core_rknd ), dimension(gr%nz) :: & - chi_mean, & ! Mean extended cloud water mixing ratio of the + chi_mean, & ! Mean extended cloud water mixing ratio of the ! two Gaussian distributions vert_cloud_frac_upper, & ! Fraction of cloud in top half of grid box vert_cloud_frac_lower, & ! Fraction of cloud in bottom half of grid box @@ -4263,8 +4314,8 @@ subroutine compute_cloud_cover & do k = 1, gr%nz - chi_mean(k) = pdf_params(k)%mixt_frac * pdf_params(k)%chi_1 + & - (1.0_core_rknd-pdf_params(k)%mixt_frac) * pdf_params(k)%chi_2 + chi_mean(k) = pdf_params%mixt_frac(k) * pdf_params%chi_1(k) + & + (1.0_core_rknd-pdf_params%mixt_frac(k)) * pdf_params%chi_2(k) end do @@ -4346,9 +4397,9 @@ subroutine compute_cloud_cover & "Error: Should not arrive here in computation of cloud_cover" write(fstderr,*) "At grid level k = ", k - write(fstderr,*) "pdf_params(k)%mixt_frac = ", pdf_params(k)%mixt_frac - write(fstderr,*) "pdf_params(k)%chi_1 = ", pdf_params(k)%chi_1 - write(fstderr,*) "pdf_params(k)%chi_2 = ", pdf_params(k)%chi_2 + write(fstderr,*) "pdf_params(k)%mixt_frac = ", pdf_params%mixt_frac(k) + write(fstderr,*) "pdf_params(k)%chi_1 = ", pdf_params%chi_1(k) + write(fstderr,*) "pdf_params(k)%chi_2 = ", pdf_params%chi_2(k) write(fstderr,*) "cloud_frac(k) = ", cloud_frac(k) write(fstderr,*) "rcm(k) = ", rcm(k) write(fstderr,*) "rcm(k+1) = ", rcm(k+1) @@ -4391,7 +4442,7 @@ subroutine clip_rcm & use error_code, only: & clubb_at_least_debug_level ! Procedure - use constants_clubb, only: & + use constants_clubb, only: & fstderr, & ! Variable(s) zero_threshold @@ -4539,7 +4590,7 @@ pure subroutine calculate_thlp2_rad & do k = 1, nz if ( rcm_zm(k) > rc_tol ) then - + thlp2_forcing(k) = thlp2_forcing(k) + & thlp2_rad_coef * ( two ) * radht_zm(k) / rcm_zm(k) * thlprcp(k) diff --git a/components/cam/src/physics/clubb/advance_wp2_wp3_module.F90 b/components/cam/src/physics/clubb/advance_wp2_wp3_module.F90 index 160f03b1a755..00370d091b06 100644 --- a/components/cam/src/physics/clubb/advance_wp2_wp3_module.F90 +++ b/components/cam/src/physics/clubb/advance_wp2_wp3_module.F90 @@ -639,7 +639,8 @@ subroutine wp23_solve( dt, sfc_elevation, sigma_sqd_w, wm_zm, & ! Intent(in) lhs ! Implicit contributions to wp2/wp3 (band diag. matrix) real( kind = core_rknd ), dimension(2*gr%nz) :: & - rhs ! RHS of band matrix + rhs, & ! RHS of band matrix + rhs_save ! Saved RHS of band matrix ! real, target, dimension(2*gr%nz) :: real( kind = core_rknd ), dimension(2*gr%nz) :: & @@ -735,6 +736,10 @@ subroutine wp23_solve( dt, sfc_elevation, sigma_sqd_w, wm_zm, & ! Intent(in) l_crank_nich_diff, & ! intent(in) rhs ) ! intent(out) + ! Save the value of rhs, which will be overwritten with the solution as + ! part of the solving routine. + rhs_save = rhs + if (l_gmres) then call wp23_gmres( dt, wp2, wm_zm, wm_zt, a1, a1_zt, a3, a3_zt, & @@ -766,11 +771,26 @@ subroutine wp23_solve( dt, sfc_elevation, sigma_sqd_w, wm_zm, & ! Intent(in) lhs, rhs, solut, rcond ) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "in wp23_solve calling band_solvex for wp2_wp3" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "Error in wp23_solve calling band_solvex for wp2_wp3" + write(fstderr,*) "wp2 & wp3 LU decomp. failed" + write(fstderr,*) "wp2 and wp3 LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", & + gr%zt(k), "LHS = ", lhs(1:5,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", & + gr%zm(k), "LHS = ", lhs(1:5,2*k) + enddo ! k = 1, gr%nz + write(fstderr,*) "wp2 and wp3 RHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", & + gr%zt(k), "RHS = ", rhs_save(2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", & + gr%zm(k), "RHS = ", rhs_save(2*k) + enddo ! k = 1, gr%nz + return + endif + endif ! Est. of the condition number of the w'^2/w^3 LHS matrix call stat_update_var_pt( iwp23_matrix_condt_num, 1, one / rcond, stats_sfc ) @@ -779,14 +799,29 @@ subroutine wp23_solve( dt, sfc_elevation, sigma_sqd_w, wm_zm, & ! Intent(in) ! Perform LU decomp and solve system (LAPACK) call band_solve( "wp2_wp3", 2, 2, 2*gr%nz, nrhs, & - lhs, rhs, solut ) + lhs, rhs, solut ) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "in wp23_solve calling band_solve for wp2_wp3" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "Error in wp23_solve calling band_solve for wp2_wp3" + write(fstderr,*) "wp2 & wp3 LU decomp. failed" + write(fstderr,*) "wp2 and wp3 LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", & + gr%zt(k), "LHS = ", lhs(1:5,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", & + gr%zm(k), "LHS = ", lhs(1:5,2*k) + enddo ! k = 1, gr%nz + write(fstderr,*) "wp2 and wp3 RHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", & + gr%zt(k), "RHS = ", rhs_save(2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", & + gr%zm(k), "RHS = ", rhs_save(2*k) + enddo ! k = 1, gr%nz + return + endif + endif endif diff --git a/components/cam/src/physics/clubb/advance_xm_wpxp_module.F90 b/components/cam/src/physics/clubb/advance_xm_wpxp_module.F90 index 25741fe2d1c5..3b0384dd4862 100644 --- a/components/cam/src/physics/clubb/advance_xm_wpxp_module.F90 +++ b/components/cam/src/physics/clubb/advance_xm_wpxp_module.F90 @@ -32,7 +32,8 @@ module advance_xm_wpxp_module wpxp_terms_bp_pr3_rhs_all, & xm_correction_wpxp_cl, & damp_coefficient, & - diagnose_upxp + diagnose_upxp, & + error_prints_xm_wpxp ! Parameter Constants integer, parameter, private :: & @@ -186,6 +187,10 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & iuprtp, & ivpthlp, & ivprtp, & + iupthvp, & + iuprcp, & + ivpthvp, & + ivprcp, & iC7_Skw_fnc, & iC6rt_Skw_fnc, & iC6thl_Skw_fnc, & @@ -430,8 +435,9 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & zeros_vector ! Array of zeros, of the size of a vertical profile [-] real( kind = core_rknd ), allocatable, dimension(:,:) :: & - rhs, &! Right-hand sides of band diag. matrix. (LAPACK) - solution ! solution vectors of band diag. matrix. (LAPACK) + rhs, & ! Right-hand sides of band diag. matrix. (LAPACK) + rhs_save, & ! Saved Right-hand sides of band diag. matrix. (LAPACK) + solution ! solution vectors of band diag. matrix. (LAPACK) ! Constant parameters as a function of Skw. @@ -440,8 +446,27 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & real( kind = core_rknd ) :: rcond + ! Saved values of predictive fields, prior to being advanced, for use in + ! print statements in case of fatal error. + real( kind = core_rknd ), dimension(gr%nz) :: & + rtm_old, & ! Saved value of r_t [kg/kg] + wprtp_old, & ! Saved value of w'r_t' [(kg/kg) m/s] + thlm_old, & ! Saved value of th_l [K] + wpthlp_old ! Saved value of w'th_l' [K m/s] + + ! Input/Output Variables + real( kind = core_rknd ), dimension(gr%nz,sclr_dim) :: & + sclrm_old, wpsclrp_old ! [Units vary] + + ! Variables used to predict and , as well as and . + real( kind = core_rknd ), dimension(gr%nz) :: & + um_old, & ! Saved value of [m/s] + upwp_old, & ! Saved value of [m^2/s^2] + vm_old, & ! Saved value of [m/s] + vpwp_old ! Saved value of [m^2/s^2] + ! Indices - integer :: i + integer :: i, k !--------------------------------------------------------------------------- @@ -460,12 +485,29 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & ! Allocate rhs and solution vector allocate( rhs(2*gr%nz,nrhs) ) + allocate( rhs_save(2*gr%nz,nrhs) ) allocate( solution(2*gr%nz,nrhs) ) ! This is initialized solely for the purpose of avoiding a compiler ! warning about uninitialized variables. zeros_vector = zero + ! Save values of predictive fields to be printed in case of crash. + rtm_old = rtm + wprtp_old = wprtp + thlm_old = thlm + wpthlp_old = wpthlp + if ( sclr_dim > 0 ) then + sclrm_old = sclrm + wpsclrp_old = wpsclrp + endif ! sclr_dim > 0 + if ( l_predict_upwp_vpwp ) then + um_old = um + upwp_old = upwp + vm_old = vm + vpwp_old = vpwp + endif ! l_predict_upwp_vpwp + ! Compute C6 and C7 as a function of Skw ! The if...then is just here to save compute time if ( abs(C6rt-C6rtb) > abs(C6rt+C6rtb)*eps/2 ) then @@ -851,6 +893,10 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & wpxp_upper_lim, wpxp_lower_lim, & ! In rhs(:,1) ) ! Out + ! Save the value of rhs, which will be overwritten with the solution as + ! part of the solving routine. + rhs_save = rhs + ! Solve r_t / w'r_t' if ( l_stats_samp .and. irtm_matrix_condt_num > 0 ) then call xm_wpxp_solve( nrhs, & ! Intent(in) @@ -863,11 +909,47 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & endif if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,'(a)') "Mean total water & total water flux LU decomp. failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "Mean total water & total water flux LU decomp. failed" + write(fstderr,*) "rtm and wprtp LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k) + enddo ! k = 1, gr%nz + write(fstderr,*) "rtm and wprtp RHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "RHS = ", rhs_save(2*k-1,1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "RHS = ", rhs_save(2*k,1) + enddo ! k = 1, gr%nz + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif call xm_wpxp_clipping_and_stats & ( xm_wpxp_rtm, dt, wp2, rtp2, wm_zt, & ! Intent(in) @@ -879,11 +961,33 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & rtm, rt_tol_mfl, wprtp ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,'(a)') "rtm monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "rtm monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif ! Compute the upper and lower limits of w'th_l' at every level, ! based on the correlation of w and th_l, such that: @@ -918,6 +1022,10 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & wpxp_upper_lim, wpxp_lower_lim, & ! In rhs(:,1) ) ! Out + ! Save the value of rhs, which will be overwritten with the solution as + ! part of the solving routine. + rhs_save = rhs + ! Solve for th_l / w'th_l' if ( l_stats_samp .and. ithlm_matrix_condt_num > 0 ) then call xm_wpxp_solve( nrhs, & ! Intent(in) @@ -930,11 +1038,47 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & endif if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,'(a)') "Liquid pot. temp & thetal flux LU decomp. failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "Liquid pot. temp & thetal flux LU decomp. failed" + write(fstderr,*) "thlm and wpthlp LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k) + enddo ! k = 1, gr%nz + write(fstderr,*) "thlm and wpthlp RHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "RHS = ", rhs_save(2*k-1,1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "RHS = ", rhs_save(2*k,1) + enddo ! k = 1, gr%nz + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif call xm_wpxp_clipping_and_stats & ( xm_wpxp_thlm, dt, wp2, thlp2, wm_zt, & ! Intent(in) @@ -946,11 +1090,33 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & thlm, thl_tol_mfl, wpthlp ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,'(a)') "thlm monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "thlm monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif ! Solve sclrm / wpsclrp ! If sclr_dim is 0, then this loop will execute 0 times. @@ -1006,17 +1172,58 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & wpxp_upper_lim, wpxp_lower_lim, & ! In rhs(:,1) ) ! Out + ! Save the value of rhs, which will be overwritten with the solution as + ! part of the solving routine. + rhs_save = rhs + ! Solve for sclrm / w'sclr' call xm_wpxp_solve( nrhs, & ! Intent(in) lhs, rhs, & ! Intent(inout) solution ) ! Intent(out) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "Passive scalar # ", i, " LU decomp. failed." - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "Passive scalar # ", i, " LU decomp. failed." + write(fstderr,*) "sclrm and wpsclrp LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k) + enddo ! k = 1, gr%nz + write(fstderr,*) "sclrm and wpsclrp RHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "RHS = ", rhs_save(2*k-1,1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "RHS = ", rhs_save(2*k,1) + enddo ! k = 1, gr%nz + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, & + wpthlp_forcing, thlm_ref, rho_ds_zm, & + rho_ds_zt, invrs_rho_ds_zm, & + invrs_rho_ds_zt, thv_ds_zm, rtp2, & + thlp2, w_1_zm, w_2_zm, varnce_w_1_zm, & + varnce_w_2_zm, mixt_frac_zm, & + l_implemented, em, wp2sclrp, & + sclrpthvp, sclrm_forcing, sclrp2, & + exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, & + vprcp, rc_coef, rtm, wprtp, thlm, & + wpthlp, sclrm, wpsclrp, um, upwp, vm, & + vpwp, rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif call xm_wpxp_clipping_and_stats & ( xm_wpxp_scalar, dt, wp2, sclrp2(:,i), & ! Intent(in) @@ -1029,11 +1236,34 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & sclrm(:,i), sclr_tol(i), wpsclrp(:,i) ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "sclrm # ", i, "monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "sclrm # ", i, "monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, & + wpthlp_forcing, thlm_ref, rho_ds_zm, & + rho_ds_zt, invrs_rho_ds_zm, & + invrs_rho_ds_zt, thv_ds_zm, rtp2, & + thlp2, w_1_zm, w_2_zm, varnce_w_1_zm, & + varnce_w_2_zm, mixt_frac_zm, & + l_implemented, em, wp2sclrp, & + sclrpthvp, sclrm_forcing, sclrp2, & + exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, & + vprcp, rc_coef, rtm, wprtp, thlm, & + wpthlp, sclrm, wpsclrp, um, upwp, vm, & + vpwp, rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif enddo ! passive scalars @@ -1201,13 +1431,13 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & endif ! .not. l_implemented ! Add "extra term" and optional Coriolis term for and . - upwp_forcing = zero ! C7_Skw_fnc * wp2 * ddzt( um ) - vpwp_forcing = zero ! C7_Skw_fnc * wp2 * ddzt( vm ) + upwp_forcing = C7_Skw_fnc * wp2 * ddzt( um ) + vpwp_forcing = C7_Skw_fnc * wp2 * ddzt( vm ) if ( l_stats_samp ) then - call stat_update_var( iupwp_pr4, 0.0_core_rknd * C7_Skw_fnc * wp2 * ddzt( um ), & + call stat_update_var( iupwp_pr4, C7_Skw_fnc * wp2 * ddzt( um ), & stats_zm ) - call stat_update_var( ivpwp_pr4, 0.0_core_rknd * C7_Skw_fnc * wp2 * ddzt( vm ), & + call stat_update_var( ivpwp_pr4, C7_Skw_fnc * wp2 * ddzt( vm ), & stats_zm ) endif ! l_stats_samp @@ -1224,24 +1454,27 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & C6rt_Skw_fnc, tau_C6_zm, C7_Skw_fnc, & ! Intent(in) vprtp ) ! Intent(out) + ! Use a crude approximation for buoyancy terms and . + !upthvp = upwp * wpthvp / max( wp2, w_tol_sqd ) + !vpthvp = vpwp * wpthvp / max( wp2, w_tol_sqd ) + !upthvp = 0.3_core_rknd * ( upthlp + 200.0_core_rknd * uprtp ) & + ! + 200._core_rknd * sign( one, upwp) * sqrt( up2 * rcm**2 ) + !vpthvp = 0.3_core_rknd * ( vpthlp + 200.0_core_rknd * vprtp ) & + ! + 200._core_rknd * sign( one, vpwp ) * sqrt( vp2 * rcm**2 ) + upthvp = upthlp + ep1 * thv_ds_zm * uprtp + rc_coef * uprcp + vpthvp = vpthlp + ep1 * thv_ds_zm * vprtp + rc_coef * vprcp + if ( l_stats_samp ) then call stat_update_var( iupthlp, upthlp, stats_zm ) call stat_update_var( iuprtp, uprtp, stats_zm ) call stat_update_var( ivpthlp, vpthlp, stats_zm ) call stat_update_var( ivprtp, vprtp, stats_zm ) + call stat_update_var( iupthvp, upthvp, stats_zm ) + call stat_update_var( iuprcp, uprcp, stats_zm ) + call stat_update_var( ivpthvp, vpthvp, stats_zm ) + call stat_update_var( ivprcp, vprcp, stats_zm ) endif ! l_stats_samp - ! Use a crude approximation for buoyancy terms and . - !upthvp = upwp * wpthvp / max( wp2, w_tol_sqd ) - !vpthvp = vpwp * wpthvp / max( wp2, w_tol_sqd ) - !upthvp = upthlp + ep1 * thv_ds_zm * uprtp + rc_coef * uprcp - !vpthvp = vpthlp + ep1 * thv_ds_zm * vprtp + rc_coef * vprcp - - upthvp = 0.3_core_rknd * ( upthlp + 200.0_core_rknd * uprtp ) & - + 200._core_rknd * sign( one, upwp) * sqrt( up2 * rcm**2 ) - vpthvp = 0.3_core_rknd * ( vpthlp + 200.0_core_rknd * vprtp ) & - + 200._core_rknd * sign( one, vpwp ) * sqrt( vp2 * rcm**2 ) - call xm_wpxp_rhs( xm_wpxp_um, l_iter, dt, um, upwp, & ! In um_tndcy, upwp_forcing, C7_Skw_fnc, & ! In upthvp, C6rt_Skw_fnc, tau_C6_zm, & ! In @@ -1264,6 +1497,10 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & endif ! l_predict_upwp_vpwp + ! Save the value of rhs, which will be overwritten with the solution as + ! part of the solving routine. + rhs_save = rhs + ! Solve for all fields if ( l_stats_samp & .and. ithlm_matrix_condt_num + irtm_matrix_condt_num > 0 ) then @@ -1276,6 +1513,70 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & solution ) ! Intent(out) endif + if ( clubb_at_least_debug_level( 0 ) ) then + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "xm & wpxp LU decomp. failed" + write(fstderr,*) "General xm and wpxp LHS" + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, "height [m] = ", gr%zt(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k-1) + write(fstderr,*) "zm level = ", k, "height [m] = ", gr%zm(k), & + "LHS = ", lhs(1:nsup+nsub+1,2*k) + enddo ! k = 1, gr%nz + do i = 1, nrhs + if ( i == 1 ) then + write(fstderr,*) "rtm and wprtp RHS" + elseif ( i == 2 ) then + write(fstderr,*) "thlm and wpthlp RHS" + else ! i > 2 + if ( sclr_dim > 0 ) then + if ( i <= 2+sclr_dim ) then + write(fstderr,*) "sclrm and wpsclrp RHS for sclr", i-2 + endif ! i <= 2+sclr_dim ) + endif ! sclr_dim > 0 + if ( l_predict_upwp_vpwp ) then + if ( i == 3+sclr_dim ) then + write(fstderr,*) "um and upwp RHS" + elseif ( i == 4+sclr_dim ) then + write(fstderr,*) "vm and vpwp RHS" + endif + endif ! l_predict_upwp_vpwp + endif + do k = 1, gr%nz + write(fstderr,*) "zt level = ", k, & + "height [m] = ", gr%zt(k), & + "RHS = ", rhs_save(2*k-1,i) + write(fstderr,*) "zm level = ", k, & + "height [m] = ", gr%zm(k), & + "RHS = ", rhs_save(2*k,i) + enddo ! k = 1, gr%nz + enddo ! i = 1, nrhs + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif + call xm_wpxp_clipping_and_stats & ( xm_wpxp_rtm, dt, wp2, rtp2, wm_zt, & ! Intent(in) rtm_forcing, rho_ds_zm, rho_ds_zt, & ! Intent(in) @@ -1286,11 +1587,33 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & rtm, rt_tol_mfl, wprtp ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "rtm monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "rtm monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif call xm_wpxp_clipping_and_stats & ( xm_wpxp_thlm, dt, wp2, thlp2, wm_zt, & ! Intent(in) @@ -1302,11 +1625,33 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & thlm, thl_tol_mfl, wpthlp ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,'(a)') "thlm monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "thlm monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif ! ---> h1g, 2010-06-15 ! scalar transport, e.g, droplet and ice number concentration @@ -1329,11 +1674,34 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & sclrm(:,i), sclr_tol(i), wpsclrp(:,i) ) ! Intent(inout) if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - write(fstderr,*) "sclrm # ", i, "monotonic flux limiter: tridag failed" - return - end if - end if + if ( err_code == clubb_fatal_error ) then + write(fstderr,*) "sclrm # ", i, "monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, & + wpthlp_forcing, thlm_ref, rho_ds_zm, & + rho_ds_zt, invrs_rho_ds_zm, & + invrs_rho_ds_zt, thv_ds_zm, rtp2, & + thlp2, w_1_zm, w_2_zm, varnce_w_1_zm, & + varnce_w_2_zm, mixt_frac_zm, & + l_implemented, em, wp2sclrp, & + sclrpthvp, sclrm_forcing, sclrp2, & + exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, & + vprcp, rc_coef, rtm, wprtp, thlm, & + wpthlp, sclrm, wpsclrp, um, upwp, vm, & + vpwp, rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + return + endif + endif end do ! 1..sclr_dim @@ -1353,6 +1721,30 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & if ( clubb_at_least_debug_level( 0 ) ) then if ( err_code == clubb_fatal_error ) then write(fstderr,*) "um monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, & + wpthlp_forcing, thlm_ref, rho_ds_zm, & + rho_ds_zt, invrs_rho_ds_zm, & + invrs_rho_ds_zt, thv_ds_zm, rtp2, & + thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, & + um_forcing, vm_forcing, ug, vg, & + wpthvp, fcor, um_ref, vm_ref, up2, & + vp2, uprcp, vprcp, rc_coef, rtm, & + wprtp, thlm, wpthlp, sclrm, wpsclrp, & + um, upwp, vm, vpwp, rtm_old, & + wprtp_old, thlm_old, wpthlp_old, & + sclrm_old, wpsclrp_old, um_old, & + upwp_old, vm_old, vpwp_old ) return endif endif @@ -1369,6 +1761,30 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & if ( clubb_at_least_debug_level( 0 ) ) then if ( err_code == clubb_fatal_error ) then write(fstderr,*) "vm monotonic flux limiter: tridag failed" + call error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, & + wpthlp_forcing, thlm_ref, rho_ds_zm, & + rho_ds_zt, invrs_rho_ds_zm, & + invrs_rho_ds_zt, thv_ds_zm, rtp2, & + thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, & + um_forcing, vm_forcing, ug, vg, & + wpthvp, fcor, um_ref, vm_ref, up2, & + vp2, uprcp, vprcp, rc_coef, rtm, & + wprtp, thlm, wpthlp, sclrm, wpsclrp, & + um, upwp, vm, vpwp, rtm_old, & + wprtp_old, thlm_old, wpthlp_old, & + sclrm_old, wpsclrp_old, um_old, & + upwp_old, vm_old, vpwp_old ) return endif endif @@ -1380,7 +1796,7 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & ! .and. ( .not. l_explicit_turbulent_adv_wpxp ) ) ! De-allocate memory - deallocate( rhs, solution ) + deallocate( rhs, rhs_save, solution ) if ( rtm_sponge_damp_settings%l_sponge_damping ) then @@ -1463,91 +1879,12 @@ subroutine advance_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & endif ! l_predict_upwp_vpwp - if ( clubb_at_least_debug_level( 0 ) ) then - if ( err_code == clubb_fatal_error ) then - - write(fstderr,*) "xm_wpxp matrix LU decomp. failed" - write(fstderr,*) "Error in advance_xm_wpxp" - - write(fstderr,*) "Intent(in)" - - write(fstderr,*) "gr%zt = ", gr%zt, new_line('c') - write(fstderr,*) "dt = ", dt, new_line('c') - write(fstderr,*) "sigma_sqd_w = ", sigma_sqd_w, new_line('c') - write(fstderr,*) "wm_zm = ", wm_zm, new_line('c') - write(fstderr,*) "wm_zt = ", wm_zt, new_line('c') - write(fstderr,*) "wp2 = ", wp2, new_line('c') - write(fstderr,*) "wp3_on_wp2 = ", wp3_on_wp2, new_line('c') - write(fstderr,*) "wp3_on_wp2_zt = ", wp3_on_wp2_zt, new_line('c') - write(fstderr,*) "Kh_zt = ", Kh_zt, new_line('c') - write(fstderr,*) "tau_C6_zm = ", tau_C6_zm, new_line('c') - write(fstderr,*) "Skw_zm = ", Skw_zm, new_line('c') - write(fstderr,*) "wp2rtp = ", wp2rtp, new_line('c') - write(fstderr,*) "rtpthvp = ", rtpthvp, new_line('c') - write(fstderr,*) "rtm_forcing = ", rtm_forcing, new_line('c') - write(fstderr,*) "wprtp_forcing = ", wprtp_forcing, new_line('c') - write(fstderr,*) "rtm_ref = ", rtm_ref, new_line('c') - write(fstderr,*) "wp2thlp = ", wp2thlp, new_line('c') - write(fstderr,*) "thlpthvp = ", thlpthvp, new_line('c') - write(fstderr,*) "thlm_forcing = ", thlm_forcing, new_line('c') - write(fstderr,*) "wpthlp_forcing = ", wpthlp_forcing, new_line('c') - write(fstderr,*) "thlm_ref = ", thlm_ref, new_line('c') - write(fstderr,*) "rho_ds_zm = ", rho_ds_zm, new_line('c') - write(fstderr,*) "rho_ds_zt = ", rho_ds_zt, new_line('c') - write(fstderr,*) "invrs_rho_ds_zm = ", invrs_rho_ds_zm, new_line('c') - write(fstderr,*) "invrs_rho_ds_zt = ", invrs_rho_ds_zt, new_line('c') - write(fstderr,*) "thv_ds_zm = ", thv_ds_zm, new_line('c') - write(fstderr,*) "rtp2 = ", rtp2, new_line('c') - write(fstderr,*) "thlp2 = ", thlp2, new_line('c') - write(fstderr,*) "w_1_zm = ", w_1_zm, new_line('c') - write(fstderr,*) "w_2_zm = ", w_2_zm, new_line('c') - write(fstderr,*) "varnce_w_1_zm = ", varnce_w_1_zm, new_line('c') - write(fstderr,*) "varnce_w_2_zm = ", varnce_w_2_zm, new_line('c') - write(fstderr,*) "mixt_frac_zm = ", mixt_frac_zm, new_line('c') - write(fstderr,*) "l_implemented = ", l_implemented, new_line('c') - write(fstderr,*) "Lscale = ", Lscale, new_line('c') - write(fstderr,*) "Kh_zm = ", Kh_zm, new_line('c') - write(fstderr,*) "em = ", em, new_line('c') - write(fstderr,*) "exner = ", exner, new_line('c') - write(fstderr,*) "rcm = ", rcm, new_line('c') - write(fstderr,*) "p_in_Pa = ", p_in_Pa, new_line('c') - write(fstderr,*) "thvm = ", thvm, new_line('c') - write(fstderr,*) "Cx_fnc_Richardson = ", Cx_fnc_Richardson, & - new_line('c') - write(fstderr,*) "pdf_implicit_coefs_terms = ", & - pdf_implicit_coefs_terms - write(fstderr,*) new_line('c') - - if ( sclr_dim > 0 ) then - write(fstderr,*) "sclrp2 = ", sclrp2, new_line('c') - write(fstderr,*) "wp2sclrp = ", wp2sclrp, new_line('c') - write(fstderr,*) "sclrpthvp = ", sclrpthvp, new_line('c') - write(fstderr,*) "sclrm_forcing = ", sclrm_forcing, new_line('c') - end if - - write(fstderr,*) "Intent(inout)" - - write(fstderr,*) "rtm = ", rtm, new_line('c') - write(fstderr,*) "wprtp = ", wprtp, new_line('c') - write(fstderr,*) "thlm = ", thlm, new_line('c') - write(fstderr,*) "wpthlp =", wpthlp, new_line('c') - - if ( sclr_dim > 0 ) then - write(fstderr,*) "sclrm = ", sclrm, new_line('c') - write(fstderr,*) "wpsclrp = ", wpsclrp, new_line('c') - end if - - return - - endif ! Fatal error and debug_level >= 1 - endif - return end subroutine advance_xm_wpxp - !================================================================================== + !============================================================================= subroutine xm_wpxp_lhs( l_iter, dt, Kh_zm, wpxp, wm_zm, wm_zt, wp2, & ! In coef_wp2xp_implicit, coef_wp2xp_implicit_zm, & ! In sgn_turbulent_vel, Kw6, tau_C6_zm, C7_Skw_fnc, & ! In @@ -2105,7 +2442,7 @@ subroutine xm_wpxp_lhs( l_iter, dt, Kh_zm, wpxp, wm_zm, wm_zt, wp2, & ! In end subroutine xm_wpxp_lhs - !================================================================================== + !============================================================================= subroutine xm_wpxp_rhs( solve_type, l_iter, dt, xm, wpxp, & ! In xm_forcing, wpxp_forcing, C7_Skw_fnc, & ! In xpthvp, C6x_Skw_fnc, tau_C6_zm, & ! In @@ -2566,7 +2903,7 @@ subroutine xm_wpxp_solve( nrhs, lhs, rhs, solution, rcond ) if ( clubb_at_least_debug_level( 0 ) ) then if ( err_code /= clubb_no_error ) then - write(fstderr,*) "in xm_wpxp_solve" + write(fstderr,*) "Error in xm_wpxp_solve" return end if end if @@ -4057,5 +4394,283 @@ pure subroutine diagnose_upxp( ypwp, xm, wpxp, ym, & end subroutine diagnose_upxp + !============================================================================= + subroutine error_prints_xm_wpxp( dt, sigma_sqd_w, wm_zm, wm_zt, wp2, & + Lscale, wp3_on_wp2, wp3_on_wp2_zt, & + Kh_zt, Kh_zm, tau_C6_zm, Skw_zm, & + wp2rtp, rtpthvp, rtm_forcing, & + wprtp_forcing, rtm_ref, wp2thlp, & + thlpthvp, thlm_forcing, wpthlp_forcing, & + thlm_ref, rho_ds_zm, rho_ds_zt, & + invrs_rho_ds_zm, invrs_rho_ds_zt, & + thv_ds_zm, rtp2, thlp2, w_1_zm, w_2_zm, & + varnce_w_1_zm, varnce_w_2_zm, & + mixt_frac_zm, l_implemented, em, & + wp2sclrp, sclrpthvp, sclrm_forcing, & + sclrp2, exner, rcm, p_in_Pa, thvm, & + Cx_fnc_Richardson, & + pdf_implicit_coefs_terms, um_forcing, & + vm_forcing, ug, vg, wpthvp, fcor, & + um_ref, vm_ref, up2, vp2, uprcp, vprcp, & + rc_coef, rtm, wprtp, thlm, wpthlp, & + sclrm, wpsclrp, um, upwp, vm, vpwp, & + rtm_old, wprtp_old, thlm_old, & + wpthlp_old, sclrm_old, wpsclrp_old, & + um_old, upwp_old, vm_old, vpwp_old ) + + ! Description: + ! Prints values of model fields when fatal errors (LU decomp.) occur. + ! All field that are passed into and out of subroutine advance_xm_wpxp are + ! printed. If additional fields are added to the call to subroutine + ! advance_xm_wpxp, they should also be added here. + + use constants_clubb, only: & + fstderr ! Variable(s) + + use grid_class, only: & + gr ! Variable Type(s) + + use model_flags, only: & + l_predict_upwp_vpwp ! Variable(s) + + use parameters_model, only: & + sclr_dim ! Variable(s) + + use pdf_parameter_module, only: & + implicit_coefs_terms ! Variable Type(s) + + use clubb_precision, only: & + core_rknd ! Variable(s) + + implicit none + + ! Input Variables + real( kind = core_rknd ), intent(in) :: & + dt ! Timestep [s] + + real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & + sigma_sqd_w, & ! sigma_sqd_w on momentum levels [-] + wm_zm, & ! w wind component on momentum levels [m/s] + wm_zt, & ! w wind component on thermodynamic levels [m/s] + wp2, & ! w'^2 (momentum levels) [m^2/s^2] + Lscale, & ! Turbulent mixing length [m] + em, & ! Turbulent Kinetic Energy (TKE) [m^2/s^2] + wp3_on_wp2, & ! Smoothed wp3 / wp2 on momentum levels [m/s] + wp3_on_wp2_zt, & ! Smoothed wp3 / wp2 on thermo. levels [m/s] + Kh_zt, & ! Eddy diffusivity on thermodynamic levels [m^2/s] + Kh_zm, & ! Eddy diffusivity on momentum levels + tau_C6_zm, & ! Time-scale tau on momentum levels applied to C6 term [s] + Skw_zm, & ! Skewness of w on momentum levels [-] + wp2rtp, & ! (thermodynamic levels) [m^2/s^2 kg/kg] + rtpthvp, & ! r_t'th_v' (momentum levels) [(kg/kg) K] + rtm_forcing, & ! r_t forcing (thermodynamic levels) [(kg/kg)/s] + wprtp_forcing, & ! forcing (momentum levels) [(kg/kg)/s^2] + rtm_ref, & ! rtm for nudging [kg/kg] + wp2thlp, & ! (thermodynamic levels) [m^2/s^2 K] + thlpthvp, & ! th_l'th_v' (momentum levels) [K^2] + thlm_forcing, & ! th_l forcing (thermodynamic levels) [K/s] + wpthlp_forcing, & ! forcing (momentum levels) [K/s^2] + thlm_ref, & ! thlm for nudging [K] + rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3] + rho_ds_zt, & ! Dry, static density on thermo. levels [kg/m^3] + invrs_rho_ds_zm, & ! Inv. dry, static density @ moment. levs. [m^3/kg] + invrs_rho_ds_zt, & ! Inv. dry, static density @ thermo. levs. [m^3/kg] + thv_ds_zm, & ! Dry, base-state theta_v on moment. levs. [K] + ! Added for clipping by Vince Larson 29 Sep 2007 + rtp2, & ! r_t'^2 (momentum levels) [(kg/kg)^2] + thlp2, & ! th_l'^2 (momentum levels) [K^2] + ! End of Vince Larson's addition. + w_1_zm, & ! Mean w (1st PDF component) [m/s] + w_2_zm, & ! Mean w (2nd PDF component) [m/s] + varnce_w_1_zm, & ! Variance of w (1st PDF component) [m^2/s^2] + varnce_w_2_zm, & ! Variance of w (2nd PDF component) [m^2/s^2] + mixt_frac_zm ! Weight of 1st PDF component (Sk_w dependent) [-] + + logical, intent(in) :: & + l_implemented ! Flag for CLUBB being implemented in a larger model. + + ! Additional variables for passive scalars + real( kind = core_rknd ), intent(in), dimension(gr%nz,sclr_dim) :: & + wp2sclrp, & ! (thermodynamic levels) [Units vary] + sclrpthvp, & ! (momentum levels) [Units vary] + sclrm_forcing, & ! sclrm forcing (thermodynamic levels) [Units vary] + sclrp2 ! For clipping Vince Larson [Units vary] + + real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & + exner, & ! Exner function [-] + rcm, & ! cloud water mixing ratio, r_c [kg/kg] + p_in_Pa, & ! Air pressure [Pa] + thvm, & ! Virutal potential temperature [K] + Cx_fnc_Richardson ! Cx_fnc computed from Richardson_num [-] + + type(implicit_coefs_terms), dimension(gr%nz), intent(in) :: & + pdf_implicit_coefs_terms ! Implicit coefs / explicit terms [units vary] + + ! Variables used to predict and , as well as and . + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + um_forcing, & ! forcing term (thermodynamic levels) [m/s^2] + vm_forcing, & ! forcing term (thermodynamic levels) [m/s^2] + ug, & ! geostrophic wind (thermodynamic levels) [m/s] + vg, & ! geostrophic wind (thermodynamic levels) [m/s] + wpthvp ! (momentum levels) [m/s K] + + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + uprcp, & ! < u' r_c' > [(m kg)/(s kg)] + vprcp, & ! < v' r_c' > [(m kg)/(s kg)] + rc_coef ! Coefficient on X'r_c' in X'th_v' equation [K/(kg/kg)] + + real( kind = core_rknd ), intent(in) :: & + fcor ! Coriolis parameter [s^-1] + + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + um_ref, & ! Reference u wind component for nudging [m/s] + vm_ref, & ! Reference v wind component for nudging [m/s] + up2, & ! Variance of the u wind component [m^2/s^2] + vp2 ! Variance of the v wind component [m^2/s^2] + + real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & + rtm, & ! r_t (total water mixing ratio) [kg/kg] + wprtp, & ! w'r_t' [(kg/kg) m/s] + thlm, & ! th_l (liquid water potential temperature) [K] + wpthlp ! w'th_l' [K m/s] + + real( kind = core_rknd ), intent(in), dimension(gr%nz,sclr_dim) :: & + sclrm, wpsclrp ! [Units vary] + + ! Variables used to predict and , as well as and . + real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & + um, & ! : mean west-east horiz. velocity (thermo. levs.) [m/s] + upwp, & ! : momentum flux (momentum levels) [m^2/s^2] + vm, & ! : mean south-north horiz. velocity (thermo. levs.) [m/s] + vpwp ! : momentum flux (momentum levels) [m^2/s^2] + + ! Saved values of predictive fields, prior to being advanced, for use in + ! print statements in case of fatal error. + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + rtm_old, & ! Saved value of r_t [kg/kg] + wprtp_old, & ! Saved value of w'r_t' [(kg/kg) m/s] + thlm_old, & ! Saved value of th_l [K] + wpthlp_old ! Saved value of w'th_l' [K m/s] + + ! Input/Output Variables + real( kind = core_rknd ), dimension(gr%nz,sclr_dim), intent(in) :: & + sclrm_old, wpsclrp_old ! [Units vary] + + ! Variables used to predict and , as well as and . + real( kind = core_rknd ), dimension(gr%nz), intent(in) :: & + um_old, & ! Saved value of [m/s] + upwp_old, & ! Saved value of [m^2/s^2] + vm_old, & ! Saved value of [m/s] + vpwp_old ! Saved value of [m^2/s^2] + + + write(fstderr,*) "Error in advance_xm_wpxp", new_line('c') + + write(fstderr,*) "Intent(in)", new_line('c') + + write(fstderr,*) "gr%zt = ", gr%zt, new_line('c') + write(fstderr,*) "gr%zm = ", gr%zm, new_line('c') + write(fstderr,*) "dt = ", dt, new_line('c') + write(fstderr,*) "sigma_sqd_w = ", sigma_sqd_w, new_line('c') + write(fstderr,*) "wm_zm = ", wm_zm, new_line('c') + write(fstderr,*) "wm_zt = ", wm_zt, new_line('c') + write(fstderr,*) "wp2 = ", wp2, new_line('c') + write(fstderr,*) "Lscale = ", Lscale, new_line('c') + write(fstderr,*) "wp3_on_wp2 = ", wp3_on_wp2, new_line('c') + write(fstderr,*) "wp3_on_wp2_zt = ", wp3_on_wp2_zt, new_line('c') + write(fstderr,*) "Kh_zt = ", Kh_zt, new_line('c') + write(fstderr,*) "Kh_zm = ", Kh_zm, new_line('c') + write(fstderr,*) "tau_C6_zm = ", tau_C6_zm, new_line('c') + write(fstderr,*) "Skw_zm = ", Skw_zm, new_line('c') + write(fstderr,*) "wp2rtp = ", wp2rtp, new_line('c') + write(fstderr,*) "rtpthvp = ", rtpthvp, new_line('c') + write(fstderr,*) "rtm_forcing = ", rtm_forcing, new_line('c') + write(fstderr,*) "wprtp_forcing = ", wprtp_forcing, new_line('c') + write(fstderr,*) "rtm_ref = ", rtm_ref, new_line('c') + write(fstderr,*) "wp2thlp = ", wp2thlp, new_line('c') + write(fstderr,*) "thlpthvp = ", thlpthvp, new_line('c') + write(fstderr,*) "thlm_forcing = ", thlm_forcing, new_line('c') + write(fstderr,*) "wpthlp_forcing = ", wpthlp_forcing, new_line('c') + write(fstderr,*) "thlm_ref = ", thlm_ref, new_line('c') + write(fstderr,*) "rho_ds_zm = ", rho_ds_zm, new_line('c') + write(fstderr,*) "rho_ds_zt = ", rho_ds_zt, new_line('c') + write(fstderr,*) "invrs_rho_ds_zm = ", invrs_rho_ds_zm, new_line('c') + write(fstderr,*) "invrs_rho_ds_zt = ", invrs_rho_ds_zt, new_line('c') + write(fstderr,*) "thv_ds_zm = ", thv_ds_zm, new_line('c') + write(fstderr,*) "rtp2 = ", rtp2, new_line('c') + write(fstderr,*) "thlp2 = ", thlp2, new_line('c') + write(fstderr,*) "w_1_zm = ", w_1_zm, new_line('c') + write(fstderr,*) "w_2_zm = ", w_2_zm, new_line('c') + write(fstderr,*) "varnce_w_1_zm = ", varnce_w_1_zm, new_line('c') + write(fstderr,*) "varnce_w_2_zm = ", varnce_w_2_zm, new_line('c') + write(fstderr,*) "mixt_frac_zm = ", mixt_frac_zm, new_line('c') + write(fstderr,*) "l_implemented = ", l_implemented, new_line('c') + write(fstderr,*) "em = ", em, new_line('c') + write(fstderr,*) "exner = ", exner, new_line('c') + write(fstderr,*) "rcm = ", rcm, new_line('c') + write(fstderr,*) "p_in_Pa = ", p_in_Pa, new_line('c') + write(fstderr,*) "thvm = ", thvm, new_line('c') + write(fstderr,*) "Cx_fnc_Richardson = ", Cx_fnc_Richardson, new_line('c') + write(fstderr,*) "pdf_implicit_coefs_terms = ", pdf_implicit_coefs_terms, & + new_line('c') + + if ( sclr_dim > 0 ) then + write(fstderr,*) "sclrp2 = ", sclrp2, new_line('c') + write(fstderr,*) "wp2sclrp = ", wp2sclrp, new_line('c') + write(fstderr,*) "sclrpthvp = ", sclrpthvp, new_line('c') + write(fstderr,*) "sclrm_forcing = ", sclrm_forcing, new_line('c') + endif + + if ( l_predict_upwp_vpwp ) then + write(fstderr,*) "um_forcing = ", um_forcing, new_line('c') + write(fstderr,*) "vm_forcing = ", vm_forcing, new_line('c') + write(fstderr,*) "ug = ", ug, new_line('c') + write(fstderr,*) "vg = ", vg, new_line('c') + write(fstderr,*) "wpthvp = ", wpthvp, new_line('c') + write(fstderr,*) "fcor = ", fcor, new_line('c') + write(fstderr,*) "um_ref = ", um_ref, new_line('c') + write(fstderr,*) "vm_ref = ", vm_ref, new_line('c') + write(fstderr,*) "up2 = ", up2, new_line('c') + write(fstderr,*) "vp2 = ", vp2, new_line('c') + write(fstderr,*) "uprcp = ", uprcp, new_line('c') + write(fstderr,*) "vprcp = ", vprcp, new_line('c') + write(fstderr,*) "rc_coef = ", rc_coef, new_line('c') + endif ! l_predict_upwp_vpwp + + write(fstderr,*) "Intent(inout)", new_line('c') + + write(fstderr,*) "rtm (pre-solve) = ", rtm_old, new_line('c') + write(fstderr,*) "rtm = ", rtm, new_line('c') + write(fstderr,*) "wprtp (pre-solve) = ", wprtp_old, new_line('c') + write(fstderr,*) "wprtp = ", wprtp, new_line('c') + write(fstderr,*) "thlm (pre-solve) = ", thlm_old, new_line('c') + write(fstderr,*) "thlm = ", thlm, new_line('c') + write(fstderr,*) "wpthlp (pre-solve) =", wpthlp_old, new_line('c') + write(fstderr,*) "wpthlp =", wpthlp, new_line('c') + + if ( sclr_dim > 0 ) then + write(fstderr,*) "sclrm (pre-solve) = ", sclrm_old, new_line('c') + write(fstderr,*) "sclrm = ", sclrm, new_line('c') + write(fstderr,*) "wpsclrp (pre-solve) = ", wpsclrp_old, new_line('c') + write(fstderr,*) "wpsclrp = ", wpsclrp, new_line('c') + endif + + if ( l_predict_upwp_vpwp ) then + write(fstderr,*) "um (pre-solve) = ", um_old, new_line('c') + write(fstderr,*) "um = ", um, new_line('c') + write(fstderr,*) "upwp (pre-solve) = ", upwp_old, new_line('c') + write(fstderr,*) "upwp = ", upwp, new_line('c') + write(fstderr,*) "vm (pre-solve) = ", vm_old, new_line('c') + write(fstderr,*) "vm = ", vm, new_line('c') + write(fstderr,*) "vpwp (pre-solve) = ", vpwp_old, new_line('c') + write(fstderr,*) "vpwp = ", vpwp, new_line('c') + endif ! l_predict_upwp_vpwp + + + return + + end subroutine error_prints_xm_wpxp + + !============================================================================= end module advance_xm_wpxp_module diff --git a/components/cam/src/physics/clubb/advance_xp2_xpyp_module.F90 b/components/cam/src/physics/clubb/advance_xp2_xpyp_module.F90 index 97689e45da33..d5daaf3756a3 100644 --- a/components/cam/src/physics/clubb/advance_xp2_xpyp_module.F90 +++ b/components/cam/src/physics/clubb/advance_xp2_xpyp_module.F90 @@ -3592,7 +3592,7 @@ subroutine update_xp2_mc( nz, dt, cloud_frac, rcm, rvm, thlm, & !It is expected that this variable is negative, as !that is the convention in Morrison microphysics - type(pdf_parameter), dimension(nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params ! PDF parameters !input/output variables diff --git a/components/cam/src/physics/clubb/clubb_api_module.F90 b/components/cam/src/physics/clubb/clubb_api_module.F90 index b10aa33db844..fdbfc7362317 100644 --- a/components/cam/src/physics/clubb/clubb_api_module.F90 +++ b/components/cam/src/physics/clubb/clubb_api_module.F90 @@ -148,7 +148,8 @@ module clubb_api_module use parameters_tunable, only : & l_prescribed_avg_deltaz, & ! used in adj_low_res_nu. If .true., avg_deltaz = deltaz - mu + mu, & + params_list use parameter_indices, only: & nparams, & ! Variable(s) @@ -358,7 +359,8 @@ module clubb_api_module stats_tout public :: & - calculate_thlp2_rad_api, mu, update_xp2_mc_api, sat_mixrat_liq_api + calculate_thlp2_rad_api, mu, params_list, & + update_xp2_mc_api, sat_mixrat_liq_api public :: & ! To Convert Between Common CLUBB-related quantities: @@ -377,6 +379,7 @@ module clubb_api_module clubb_no_error, & fill_holes_driver_api, & ! OR fill_holes_vertical_api, & + fill_holes_hydromet_api, & set_clubb_debug_level_api, & vertical_integral_api @@ -430,6 +433,7 @@ module clubb_api_module !#ifdef CLUBB_CAM /* Code for storing pdf_parameter structs in pbuf as array */ pack_pdf_params_api, & unpack_pdf_params_api, & + init_pdf_params_api, & num_pdf_params, & !#endif adj_low_res_nu_api, & @@ -670,7 +674,7 @@ subroutine advance_clubb_core_api( & real( kind = core_rknd ), intent(inout), dimension(gr%nz,sclr_dim) :: & sclrpthvp ! < sclr' th_v' > (momentum levels) [units vary] - type(pdf_parameter), dimension(gr%nz), intent(inout) :: & + type(pdf_parameter), intent(inout) :: & pdf_params, & ! PDF parameters (thermodynamic levels) [units vary] pdf_params_zm ! PDF parameters on momentum levels [units vary] @@ -1158,6 +1162,37 @@ subroutine fill_holes_vertical_api( & field ) end subroutine fill_holes_vertical_api + !============================================================================= + ! fill_holes_hydromet - fills holes in a hydrometeor using mass from another + ! hydrometeor that has the same phase. + !============================================================================= + subroutine fill_holes_hydromet_api( nz, hydromet_dim, hydromet, & ! Intent(in) + hydromet_filled ) ! Intent(out) + + use fill_holes, only: & + fill_holes_hydromet ! Procedure + + implicit none + + ! Input Variables + integer, intent(in) :: & + hydromet_dim, & ! Number of hydrometeor fields + nz ! Number of vertical grid levels + + real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: & + hydromet ! Mean of hydrometeor fields [units vary] + + ! Output Variables + real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: & + hydromet_filled ! Mean of hydrometeor fields after hole filling [un. vary] + + + call fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in) + hydromet_filled ) ! Intent(out) + + + end subroutine fill_holes_hydromet_api + !================================================================================================ ! vertical_integral - Computes the vertical integral. !================================================================================================ @@ -1457,8 +1492,8 @@ end subroutine adj_low_res_nu_api ! pack_pdf_params - Returns a two dimensional real array with all values. !================================================================================================ - subroutine pack_pdf_params_api( & - pdf_params, nz, r_param_array) + subroutine pack_pdf_params_api( pdf_params, nz, r_param_array, & + k_start, k_end ) use pdf_parameter_module, only : pack_pdf_params @@ -1466,16 +1501,23 @@ subroutine pack_pdf_params_api( & implicit none - ! Input a pdf_parameter array with nz instances of pdf_parameter integer, intent(in) :: nz ! Num Vert Model Levs - type (pdf_parameter), dimension(nz), intent(in) :: pdf_params + + ! Input a pdf_parameter array with nz instances of pdf_parameter + type (pdf_parameter), intent(inout) :: pdf_params ! Output a two dimensional real array with all values - real (kind = core_rknd), dimension(nz,num_pdf_params), intent(out) :: & + real (kind = core_rknd), dimension(nz,num_pdf_params), intent(inout) :: & r_param_array - - call pack_pdf_params( & - pdf_params, nz, r_param_array) + + integer, optional, intent(in) :: k_start, k_end + + if( present( k_start ) .and. present( k_end ) ) then + call pack_pdf_params( pdf_params, nz, r_param_array, & + k_start, k_end ) + else + call pack_pdf_params( pdf_params, nz, r_param_array ) + end if end subroutine pack_pdf_params_api @@ -1483,24 +1525,56 @@ end subroutine pack_pdf_params_api ! unpack_pdf_params - Returns a pdf_parameter array with nz instances of pdf_parameter. !================================================================================================ - subroutine unpack_pdf_params_api( & - r_param_array, nz, pdf_params) + subroutine unpack_pdf_params_api( r_param_array, nz, pdf_params, & + k_start, k_end ) use pdf_parameter_module, only : unpack_pdf_params implicit none - - ! Input a two dimensional real array with pdf values + integer, intent(in) :: nz ! Num Vert Model Levs + + ! Input a two dimensional real array with pdf values real (kind = core_rknd), dimension(nz,num_pdf_params), intent(in) :: & r_param_array ! Output a pdf_parameter array with nz instances of pdf_parameter - type (pdf_parameter), dimension(nz), intent(out) :: pdf_params + type (pdf_parameter), intent(inout) :: pdf_params + + integer, optional, intent(in) :: k_start, k_end - call unpack_pdf_params( & - r_param_array, nz, pdf_params) + if( present( k_start ) .and. present( k_end ) ) then + call unpack_pdf_params( r_param_array, nz, pdf_params, & + k_start, k_end ) + else + call unpack_pdf_params( r_param_array, nz, pdf_params ) + end if + + end subroutine unpack_pdf_params_api + + !================================================================================================ + ! init_pdf_params - allocates arrays for pdf_params + !================================================================================================ + subroutine init_pdf_params_api( nz, pdf_params ) + + use pdf_parameter_module, only : init_pdf_params + + implicit none + + ! Input Variable(s) + integer, intent(in) :: & + nz ! Number of vertical grid levels [-] + + ! Output Variable(s) + type(pdf_parameter), intent(out) :: & + pdf_params ! PDF parameters [units vary] + + call init_pdf_params( nz, pdf_params ) + + end subroutine init_pdf_params_api + + !#endif !================================================================================================ @@ -1554,7 +1628,7 @@ subroutine setup_pdf_parameters_api( & corr_array_n_cloud, & ! Prescribed norm. space corr. array in cloud [-] corr_array_n_below ! Prescribed norm. space corr. array below cloud [-] - type(pdf_parameter), dimension(nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params ! PDF parameters [units vary] logical, intent(in) :: & @@ -2102,7 +2176,7 @@ subroutine update_xp2_mc_api( nz, dt, cloud_frac, rcm, rvm, thlm, & !It is expected that this variable is negative, as !that is the convention in Morrison microphysics - type(pdf_parameter), dimension(nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params ! PDF parameters !input/output variables diff --git a/components/cam/src/physics/clubb/gmres_cache.F90 b/components/cam/src/physics/clubb/gmres_cache.F90 index c1563d587e70..571afde2bc16 100644 --- a/components/cam/src/physics/clubb/gmres_cache.F90 +++ b/components/cam/src/physics/clubb/gmres_cache.F90 @@ -45,7 +45,7 @@ module gmres_cache ! for the non-interlaced matrices (gr%nz grid ! levels) -!$omp threadprivate( gmres_tmp_intlc, gmres_temp_norm ) +!$omp threadprivate( gmres_temp_intlc, gmres_temp_norm ) integer, public :: & gmres_tempsize_norm, & ! Size of the temporary array for diff --git a/components/cam/src/physics/clubb/lapack_wrap.F90 b/components/cam/src/physics/clubb/lapack_wrap.F90 index 3161c2142553..c38e6066a084 100644 --- a/components/cam/src/physics/clubb/lapack_wrap.F90 +++ b/components/cam/src/physics/clubb/lapack_wrap.F90 @@ -498,8 +498,6 @@ subroutine band_solvex( solve_type, nsup, nsub, ndim, nrhs, & case( :-1 ) write(fstderr,*) "in band_solvex for ", trim( solve_type ), & ": illegal value for argument", -info - write(fstderr,*) "lhs = ", lhs, new_line('c') - write(fstderr,*) "rhs = ", rhs, new_line('c') err_code = clubb_fatal_error case( 0 ) @@ -519,8 +517,6 @@ subroutine band_solvex( solve_type, nsup, nsub, ndim, nrhs, & else write(fstderr,*) "in band_solvex for", trim( solve_type ), & ": singular matrix, solution not computed" - write(fstderr,*) "lhs = ", lhs, new_line('c') - write(fstderr,*) "rhs = ", rhs, new_line('c') err_code = clubb_fatal_error end if @@ -754,8 +750,6 @@ subroutine band_solve( solve_type, nsup, nsub, ndim, nrhs, & case( :-1 ) write(fstderr,*) "in band_solve for ", trim( solve_type ), & ": illegal value for argument", -info - write(fstderr,*) "lhs = ", lhs, new_line('c') - write(fstderr,*) "rhs = ", rhs, new_line('c') err_code = clubb_fatal_error case( 0 ) ! Success! @@ -770,8 +764,6 @@ subroutine band_solve( solve_type, nsup, nsub, ndim, nrhs, & case( 1: ) write(fstderr,*) "in band_solve for ", trim( solve_type ), & ": singular matrix, solution not computed" - write(fstderr,*) "lhs = ", lhs, new_line('c') - write(fstderr,*) "rhs = ", rhs, new_line('c') err_code = clubb_fatal_error end select diff --git a/components/cam/src/physics/clubb/model_flags.F90 b/components/cam/src/physics/clubb/model_flags.F90 index f4dde319621f..679aa6acc11d 100644 --- a/components/cam/src/physics/clubb/model_flags.F90 +++ b/components/cam/src/physics/clubb/model_flags.F90 @@ -142,6 +142,8 @@ module model_flags logical, public :: & l_high_accuracy_parab_cyl_fnc = .false. +!$omp threadprivate(l_high_accuracy_parab_cyl_fnc) + ! These flags determine whether we want to use an upwind differencing approximation ! rather than a centered differencing for turbulent or mean advection terms. ! wpxp_ta affects wprtp, wpthlp, & wpsclrp diff --git a/components/cam/src/physics/clubb/mt95.f90 b/components/cam/src/physics/clubb/mt95.f90 index cdbb5270ae87..eea189f172f2 100644 --- a/components/cam/src/physics/clubb/mt95.f90 +++ b/components/cam/src/physics/clubb/mt95.f90 @@ -1,10 +1,10 @@ ! A C-program for MT19937, with initialization improved 2002/1/26. ! Coded by Takuji Nishimura and Makoto Matsumoto. -! Code converted to Fortran 95 by José Rui Faustino de Sousa +! Code converted to Fortran 95 by Jose Rui Faustino de Sousa ! Date: 2002-02-01 -! Enhanced version by José Rui Faustino de Sousa +! Enhanced version by Jose Rui Faustino de Sousa ! Date: 2003-04-30 ! Interface: @@ -1312,7 +1312,7 @@ subroutine genrand_res53_7d( r ) end subroutine genrand_res53_7d ! These real versions are due to Isaku Wada, 2002/01/09 added - ! Altered by José Sousa genrand_real[1-3] will not return exactely + ! Altered by Jose Sousa genrand_real[1-3] will not return exactely ! the same values but should have the same properties and are faster end module mt95 diff --git a/components/cam/src/physics/clubb/numerical_check.F90 b/components/cam/src/physics/clubb/numerical_check.F90 index 3d5f2e2ff149..ffbd5ce119a8 100644 --- a/components/cam/src/physics/clubb/numerical_check.F90 +++ b/components/cam/src/physics/clubb/numerical_check.F90 @@ -137,7 +137,7 @@ subroutine pdf_closure_check( wp4, wprtp2, wp2rtp, wpthlp2, & crt_1, crt_2, & cthl_1, cthl_2 - type(pdf_parameter), dimension(gr%nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params ! PDF parameters [units vary] ! Input (Optional passive scalar variables) diff --git a/components/cam/src/physics/clubb/parameter_indices.F90 b/components/cam/src/physics/clubb/parameter_indices.F90 index d3e3c37da296..a32693e64c31 100644 --- a/components/cam/src/physics/clubb/parameter_indices.F90 +++ b/components/cam/src/physics/clubb/parameter_indices.F90 @@ -26,7 +26,7 @@ module parameter_indices private ! Default Scope integer, parameter, public :: & - nparams = 78 ! Total tunable parameters + nparams = 82 ! Total tunable parameters !*************************************************************** ! ***** IMPORTANT ***** @@ -120,7 +120,11 @@ module parameter_indices ithlp2_rad_coef = 75, & ithlp2_rad_cloud_frac_thresh = 76, & iup2_vp2_factor = 77, & - iSkw_max_mag = 78 + iSkw_max_mag = 78, & + iC_invrs_tau_bkgnd = 79, & + iC_invrs_tau_sfc = 80, & + iC_invrs_tau_shear = 81, & + iC_invrs_tau_N2 = 82 end module parameter_indices diff --git a/components/cam/src/physics/clubb/pdf_closure_module.F90 b/components/cam/src/physics/clubb/pdf_closure_module.F90 index f3338c40da00..d8405715ce04 100644 --- a/components/cam/src/physics/clubb/pdf_closure_module.F90 +++ b/components/cam/src/physics/clubb/pdf_closure_module.F90 @@ -261,7 +261,7 @@ subroutine pdf_closure( hydromet_dim, p_in_Pa, exner, thv_ds, & uprcp, & ! u' r_c' [(m kg)/(s kg)] vprcp ! v' r_c' [(m kg)/(s kg)] - type(pdf_parameter), dimension(gr%nz), intent(out) :: & + type(pdf_parameter), intent(inout) :: & pdf_params ! pdf paramters [units vary] type(implicit_coefs_terms), dimension(gr%nz), intent(out) :: & @@ -908,6 +908,7 @@ subroutine pdf_closure( hydromet_dim, p_in_Pa, exner, thv_ds, & ! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:anl_int_cloud_terms cloud_frac = mixt_frac * cloud_frac_1 + ( one - mixt_frac ) * cloud_frac_2 rcm = mixt_frac * rc_1 + ( one - mixt_frac ) * rc_2 + rcm = max( zero_threshold, rcm ) if ( iiPDF_type == iiPDF_ADG1 .or. iiPDF_type == iiPDF_ADG2 ) then @@ -1088,58 +1089,54 @@ subroutine pdf_closure( hydromet_dim, p_in_Pa, exner, thv_ds, & ! Save PDF parameters - do k = 1, gr%nz - - pdf_params(k)%w_1 = w_1(k) - pdf_params(k)%w_2 = w_2(k) - pdf_params(k)%varnce_w_1 = varnce_w_1(k) - pdf_params(k)%varnce_w_2 = varnce_w_2(k) - pdf_params(k)%rt_1 = rt_1(k) - pdf_params(k)%rt_2 = rt_2(k) - pdf_params(k)%varnce_rt_1 = varnce_rt_1(k) - pdf_params(k)%varnce_rt_2 = varnce_rt_2(k) - pdf_params(k)%thl_1 = thl_1(k) - pdf_params(k)%thl_2 = thl_2(k) - pdf_params(k)%varnce_thl_1 = varnce_thl_1(k) - pdf_params(k)%varnce_thl_2 = varnce_thl_2(k) - pdf_params(k)%corr_w_rt_1 = corr_w_rt_1(k) - pdf_params(k)%corr_w_rt_2 = corr_w_rt_2(k) - pdf_params(k)%corr_w_thl_1 = corr_w_thl_1(k) - pdf_params(k)%corr_w_thl_2 = corr_w_thl_2(k) - pdf_params(k)%corr_rt_thl_1 = corr_rt_thl_1(k) - pdf_params(k)%corr_rt_thl_2 = corr_rt_thl_2(k) - pdf_params(k)%alpha_thl = alpha_thl(k) - pdf_params(k)%alpha_rt = alpha_rt(k) - pdf_params(k)%crt_1 = crt_1(k) - pdf_params(k)%crt_2 = crt_2(k) - pdf_params(k)%cthl_1 = cthl_1(k) - pdf_params(k)%cthl_2 = cthl_2(k) - pdf_params(k)%chi_1 = chi_1(k) - pdf_params(k)%chi_2 = chi_2(k) - pdf_params(k)%stdev_chi_1 = stdev_chi_1(k) - pdf_params(k)%stdev_chi_2 = stdev_chi_2(k) - pdf_params(k)%stdev_eta_1 = stdev_eta_1(k) - pdf_params(k)%stdev_eta_2 = stdev_eta_2(k) - pdf_params(k)%covar_chi_eta_1 = covar_chi_eta_1(k) - pdf_params(k)%covar_chi_eta_2 = covar_chi_eta_2(k) - pdf_params(k)%corr_w_chi_1 = corr_w_chi_1(k) - pdf_params(k)%corr_w_chi_2 = corr_w_chi_2(k) - pdf_params(k)%corr_w_eta_1 = corr_w_eta_1(k) - pdf_params(k)%corr_w_eta_2 = corr_w_eta_2(k) - pdf_params(k)%corr_chi_eta_1 = corr_chi_eta_1(k) - pdf_params(k)%corr_chi_eta_2 = corr_chi_eta_2(k) - pdf_params(k)%rsatl_1 = rsatl_1(k) - pdf_params(k)%rsatl_2 = rsatl_2(k) - pdf_params(k)%rc_1 = rc_1(k) - pdf_params(k)%rc_2 = rc_2(k) - pdf_params(k)%cloud_frac_1 = cloud_frac_1(k) - pdf_params(k)%cloud_frac_2 = cloud_frac_2(k) - pdf_params(k)%mixt_frac = mixt_frac(k) - - pdf_params(k)%ice_supersat_frac_1 = ice_supersat_frac_1(k) - pdf_params(k)%ice_supersat_frac_2 = ice_supersat_frac_2(k) + pdf_params%w_1 = w_1 + pdf_params%w_2 = w_2 + pdf_params%varnce_w_1 = varnce_w_1 + pdf_params%varnce_w_2 = varnce_w_2 + pdf_params%rt_1 = rt_1 + pdf_params%rt_2 = rt_2 + pdf_params%varnce_rt_1 = varnce_rt_1 + pdf_params%varnce_rt_2 = varnce_rt_2 + pdf_params%thl_1 = thl_1 + pdf_params%thl_2 = thl_2 + pdf_params%varnce_thl_1 = varnce_thl_1 + pdf_params%varnce_thl_2 = varnce_thl_2 + pdf_params%corr_w_rt_1 = corr_w_rt_1 + pdf_params%corr_w_rt_2 = corr_w_rt_2 + pdf_params%corr_w_thl_1 = corr_w_thl_1 + pdf_params%corr_w_thl_2 = corr_w_thl_2 + pdf_params%corr_rt_thl_1 = corr_rt_thl_1 + pdf_params%corr_rt_thl_2 = corr_rt_thl_2 + pdf_params%alpha_thl = alpha_thl + pdf_params%alpha_rt = alpha_rt + pdf_params%crt_1 = crt_1 + pdf_params%crt_2 = crt_2 + pdf_params%cthl_1 = cthl_1 + pdf_params%cthl_2 = cthl_2 + pdf_params%chi_1 = chi_1 + pdf_params%chi_2 = chi_2 + pdf_params%stdev_chi_1 = stdev_chi_1 + pdf_params%stdev_chi_2 = stdev_chi_2 + pdf_params%stdev_eta_1 = stdev_eta_1 + pdf_params%stdev_eta_2 = stdev_eta_2 + pdf_params%covar_chi_eta_1 = covar_chi_eta_1 + pdf_params%covar_chi_eta_2 = covar_chi_eta_2 + pdf_params%corr_w_chi_1 = corr_w_chi_1 + pdf_params%corr_w_chi_2 = corr_w_chi_2 + pdf_params%corr_w_eta_1 = corr_w_eta_1 + pdf_params%corr_w_eta_2 = corr_w_eta_2 + pdf_params%corr_chi_eta_1 = corr_chi_eta_1 + pdf_params%corr_chi_eta_2 = corr_chi_eta_2 + pdf_params%rsatl_1 = rsatl_1 + pdf_params%rsatl_2 = rsatl_2 + pdf_params%rc_1 = rc_1 + pdf_params%rc_2 = rc_2 + pdf_params%cloud_frac_1 = cloud_frac_1 + pdf_params%cloud_frac_2 = cloud_frac_2 + pdf_params%mixt_frac = mixt_frac + pdf_params%ice_supersat_frac_1 = ice_supersat_frac_1 + pdf_params%ice_supersat_frac_2 = ice_supersat_frac_2 - end do if ( clubb_at_least_debug_level( 2 ) ) then diff --git a/components/cam/src/physics/clubb/pdf_parameter_module.F90 b/components/cam/src/physics/clubb/pdf_parameter_module.F90 index 9efdc6547245..cd5daa467115 100644 --- a/components/cam/src/physics/clubb/pdf_parameter_module.F90 +++ b/components/cam/src/physics/clubb/pdf_parameter_module.F90 @@ -24,7 +24,7 @@ module pdf_parameter_module ! CLUBB's PDF parameters. type pdf_parameter - real( kind = core_rknd ) :: & + real( kind = core_rknd ), dimension(:), allocatable :: & w_1, & ! Mean of w (1st PDF component) [m/s] w_2, & ! Mean of w (2nd PDF component) [m/s] varnce_w_1, & ! Variance of w (1st PDF component) [m^2/s^2] @@ -69,14 +69,12 @@ module pdf_parameter_module rc_2, & ! Mean of r_c (2nd PDF component) [kg/kg] cloud_frac_1, & ! Cloud fraction (1st PDF component) [-] cloud_frac_2, & ! Cloud fraction (2nd PDF component) [-] - mixt_frac ! Weight of 1st PDF component (Sk_w dependent) [-] - - real( kind = core_rknd ) :: & + mixt_frac, & ! Weight of 1st PDF component (Sk_w dependent) [-] ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.) [-] ice_supersat_frac_2 ! Ice supersaturation fraction (2nd PDF comp.) [-] end type pdf_parameter - + ! The implicit coefficients, semi-implicit coefficients and terms, and ! explicit terms for turbulent advection of turbulent fields are calculated ! from the PDF and the resulting PDF parameters. @@ -115,7 +113,7 @@ module pdf_parameter_module !#endif /* CLUBB_CAM */ contains - + !============================================================================= subroutine init_pdf_params( nz, pdf_params ) @@ -135,349 +133,266 @@ subroutine init_pdf_params( nz, pdf_params ) nz ! Number of vertical grid levels [-] ! Output Variable(s) - type(pdf_parameter), dimension(nz), intent(out) :: & + type(pdf_parameter), intent(out) :: & pdf_params ! PDF parameters [units vary] - integer :: k - - - ! Initialize all PDF parameters in variable type pdf_parameter. - do k = 1, nz - - pdf_params(k)%w_1 = zero - pdf_params(k)%w_2 = zero - pdf_params(k)%varnce_w_1 = zero - pdf_params(k)%varnce_w_2 = zero - pdf_params(k)%rt_1 = zero - pdf_params(k)%rt_2 = zero - pdf_params(k)%varnce_rt_1 = zero - pdf_params(k)%varnce_rt_2 = zero - pdf_params(k)%thl_1 = zero - pdf_params(k)%thl_2 = zero - pdf_params(k)%varnce_thl_1 = zero - pdf_params(k)%varnce_thl_2 = zero - pdf_params(k)%corr_w_rt_1 = zero - pdf_params(k)%corr_w_rt_2 = zero - pdf_params(k)%corr_w_thl_1 = zero - pdf_params(k)%corr_w_thl_2 = zero - pdf_params(k)%corr_rt_thl_1 = zero - pdf_params(k)%corr_rt_thl_2 = zero - pdf_params(k)%alpha_thl = zero - pdf_params(k)%alpha_rt = zero - pdf_params(k)%crt_1 = zero - pdf_params(k)%crt_2 = zero - pdf_params(k)%cthl_1 = zero - pdf_params(k)%cthl_2 = zero - pdf_params(k)%chi_1 = zero - pdf_params(k)%chi_2 = zero - pdf_params(k)%stdev_chi_1 = zero - pdf_params(k)%stdev_chi_2 = zero - pdf_params(k)%stdev_eta_1 = zero - pdf_params(k)%stdev_eta_2 = zero - pdf_params(k)%covar_chi_eta_1 = zero - pdf_params(k)%covar_chi_eta_2 = zero - pdf_params(k)%corr_w_chi_1 = zero - pdf_params(k)%corr_w_chi_2 = zero - pdf_params(k)%corr_w_eta_1 = zero - pdf_params(k)%corr_w_eta_2 = zero - pdf_params(k)%corr_chi_eta_1 = zero - pdf_params(k)%corr_chi_eta_2 = zero - pdf_params(k)%rsatl_1 = zero - pdf_params(k)%rsatl_2 = zero - pdf_params(k)%rc_1 = zero - pdf_params(k)%rc_2 = zero - pdf_params(k)%cloud_frac_1 = zero - pdf_params(k)%cloud_frac_2 = zero - pdf_params(k)%mixt_frac = zero - pdf_params(k)%ice_supersat_frac_1 = zero - pdf_params(k)%ice_supersat_frac_2 = zero - - end do + allocate( pdf_params%w_1(nz), & + pdf_params%w_2(nz), & + pdf_params%varnce_w_1(nz), & + pdf_params%varnce_w_2(nz), & + pdf_params%rt_1(nz), & + pdf_params%rt_2(nz), & + pdf_params%varnce_rt_1(nz), & + pdf_params%varnce_rt_2(nz), & + pdf_params%thl_1(nz), & + pdf_params%thl_2(nz), & + pdf_params%varnce_thl_1(nz), & + pdf_params%varnce_thl_2(nz), & + pdf_params%corr_w_rt_1(nz), & + pdf_params%corr_w_rt_2(nz), & + pdf_params%corr_w_thl_1(nz), & + pdf_params%corr_w_thl_2(nz), & + pdf_params%corr_rt_thl_1(nz), & + pdf_params%corr_rt_thl_2(nz), & + pdf_params%alpha_thl(nz), & + pdf_params%alpha_rt(nz), & + pdf_params%crt_1(nz), & + pdf_params%crt_2(nz), & + pdf_params%cthl_1(nz), & + pdf_params%cthl_2(nz), & + pdf_params%chi_1(nz), & + pdf_params%chi_2(nz), & + pdf_params%stdev_chi_1(nz), & + pdf_params%stdev_chi_2(nz), & + pdf_params%stdev_eta_1(nz), & + pdf_params%stdev_eta_2(nz), & + pdf_params%covar_chi_eta_1(nz), & + pdf_params%covar_chi_eta_2(nz), & + pdf_params%corr_w_chi_1(nz), & + pdf_params%corr_w_chi_2(nz), & + pdf_params%corr_w_eta_1(nz), & + pdf_params%corr_w_eta_2(nz), & + pdf_params%corr_chi_eta_1(nz), & + pdf_params%corr_chi_eta_2(nz), & + pdf_params%rsatl_1(nz), & + pdf_params%rsatl_2(nz), & + pdf_params%rc_1(nz), & + pdf_params%rc_2(nz), & + pdf_params%cloud_frac_1(nz), & + pdf_params%cloud_frac_2(nz), & + pdf_params%mixt_frac(nz), & + pdf_params%ice_supersat_frac_1(nz), & + pdf_params%ice_supersat_frac_2(nz) ) + + pdf_params%w_1 = zero + pdf_params%w_2 = zero + pdf_params%varnce_w_1 = zero + pdf_params%varnce_w_2 = zero + pdf_params%rt_1 = zero + pdf_params%rt_2 = zero + pdf_params%varnce_rt_1 = zero + pdf_params%varnce_rt_2 = zero + pdf_params%thl_1 = zero + pdf_params%thl_2 = zero + pdf_params%varnce_thl_1 = zero + pdf_params%varnce_thl_2 = zero + pdf_params%corr_w_rt_1 = zero + pdf_params%corr_w_rt_2 = zero + pdf_params%corr_w_thl_1 = zero + pdf_params%corr_w_thl_2 = zero + pdf_params%corr_rt_thl_1 = zero + pdf_params%corr_rt_thl_2 = zero + pdf_params%alpha_thl = zero + pdf_params%alpha_rt = zero + pdf_params%crt_1 = zero + pdf_params%crt_2 = zero + pdf_params%cthl_1 = zero + pdf_params%cthl_2 = zero + pdf_params%chi_1 = zero + pdf_params%chi_2 = zero + pdf_params%stdev_chi_1 = zero + pdf_params%stdev_chi_2 = zero + pdf_params%stdev_eta_1 = zero + pdf_params%stdev_eta_2 = zero + pdf_params%covar_chi_eta_1 = zero + pdf_params%covar_chi_eta_2 = zero + pdf_params%corr_w_chi_1 = zero + pdf_params%corr_w_chi_2 = zero + pdf_params%corr_w_eta_1 = zero + pdf_params%corr_w_eta_2 = zero + pdf_params%corr_chi_eta_1 = zero + pdf_params%corr_chi_eta_2 = zero + pdf_params%rsatl_1 = zero + pdf_params%rsatl_2 = zero + pdf_params%rc_1 = zero + pdf_params%rc_2 = zero + pdf_params%cloud_frac_1 = zero + pdf_params%cloud_frac_2 = zero + pdf_params%mixt_frac = zero + pdf_params%ice_supersat_frac_1 = zero + pdf_params%ice_supersat_frac_2 = zero return end subroutine init_pdf_params - + !============================================================================= ! The CLUBB_CAM preprocessor directives are being commented out because this ! code is now also used for WRF-CLUBB. !#ifdef CLUBB_CAM /* Code for storing pdf_parameter structs in pbuf as array */ - subroutine pack_pdf_params(pdf_params, nz, r_param_array) + subroutine pack_pdf_params(pdf_params, nz, r_param_array, & + k_start_in, k_end_in ) implicit none - ! Input a pdf_parameter array with nz instances of pdf_parameter + integer, intent(in) :: nz ! Num Vert Model Levs - type (pdf_parameter), dimension(nz), intent(in) :: pdf_params + + ! Input a pdf_parameter array with nz instances of pdf_parameter + type (pdf_parameter), intent(in) :: pdf_params ! Output a two dimensional real array with all values real (kind = core_rknd), dimension(nz,num_pdf_params), intent(out) :: & r_param_array - ! Local Loop vars - integer :: k, p - - do k = 1,nz - do p = 1,num_pdf_params - - r_param_array(k,p) = get_param_at_ind(pdf_params(k), p) - - end do ! p - end do ! k + integer, optional, intent(in) :: k_start_in, k_end_in + + integer :: k_start, k_end + + if( present( k_start_in ) .and. present( k_end_in ) ) then + k_start = k_start_in + k_end = k_end_in + else + k_start = 1 + k_end = nz + end if + + r_param_array(:,1) = pdf_params%w_1(k_start:k_end) + r_param_array(:,2) = pdf_params%w_2(k_start:k_end) + r_param_array(:,3) = pdf_params%varnce_w_1(k_start:k_end) + r_param_array(:,4) = pdf_params%varnce_w_2(k_start:k_end) + r_param_array(:,5) = pdf_params%rt_1(k_start:k_end) + r_param_array(:,6) = pdf_params%rt_2(k_start:k_end) + r_param_array(:,7) = pdf_params%varnce_rt_1(k_start:k_end) + r_param_array(:,8) = pdf_params%varnce_rt_2(k_start:k_end) + r_param_array(:,9) = pdf_params%thl_1(k_start:k_end) + r_param_array(:,10) = pdf_params%thl_2(k_start:k_end) + r_param_array(:,11) = pdf_params%varnce_thl_1(k_start:k_end) + r_param_array(:,12) = pdf_params%varnce_thl_2(k_start:k_end) + r_param_array(:,13) = pdf_params%corr_w_rt_1(k_start:k_end) + r_param_array(:,14) = pdf_params%corr_w_rt_2(k_start:k_end) + r_param_array(:,15) = pdf_params%corr_w_thl_1(k_start:k_end) + r_param_array(:,16) = pdf_params%corr_w_thl_2(k_start:k_end) + r_param_array(:,17) = pdf_params%corr_rt_thl_1(k_start:k_end) + r_param_array(:,18) = pdf_params%corr_rt_thl_2(k_start:k_end) + r_param_array(:,19) = pdf_params%alpha_thl(k_start:k_end) + r_param_array(:,20) = pdf_params%alpha_rt(k_start:k_end) + r_param_array(:,21) = pdf_params%crt_1(k_start:k_end) + r_param_array(:,22) = pdf_params%crt_2(k_start:k_end) + r_param_array(:,23) = pdf_params%cthl_1(k_start:k_end) + r_param_array(:,24) = pdf_params%cthl_2(k_start:k_end) + r_param_array(:,25) = pdf_params%chi_1(k_start:k_end) + r_param_array(:,26) = pdf_params%chi_2(k_start:k_end) + r_param_array(:,27) = pdf_params%stdev_chi_1(k_start:k_end) + r_param_array(:,28) = pdf_params%stdev_chi_2(k_start:k_end) + r_param_array(:,29) = pdf_params%stdev_eta_1(k_start:k_end) + r_param_array(:,30) = pdf_params%stdev_eta_2(k_start:k_end) + r_param_array(:,31) = pdf_params%covar_chi_eta_1(k_start:k_end) + r_param_array(:,32) = pdf_params%covar_chi_eta_2(k_start:k_end) + r_param_array(:,33) = pdf_params%corr_w_chi_1(k_start:k_end) + r_param_array(:,34) = pdf_params%corr_w_chi_2(k_start:k_end) + r_param_array(:,35) = pdf_params%corr_w_eta_1(k_start:k_end) + r_param_array(:,36) = pdf_params%corr_w_eta_2(k_start:k_end) + r_param_array(:,37) = pdf_params%corr_chi_eta_1(k_start:k_end) + r_param_array(:,38) = pdf_params%corr_chi_eta_2(k_start:k_end) + r_param_array(:,39) = pdf_params%rsatl_1(k_start:k_end) + r_param_array(:,40) = pdf_params%rsatl_2(k_start:k_end) + r_param_array(:,41) = pdf_params%rc_1(k_start:k_end) + r_param_array(:,42) = pdf_params%rc_2(k_start:k_end) + r_param_array(:,43) = pdf_params%cloud_frac_1(k_start:k_end) + r_param_array(:,44) = pdf_params%cloud_frac_2(k_start:k_end) + r_param_array(:,45) = pdf_params%mixt_frac(k_start:k_end) + r_param_array(:,46) = pdf_params%ice_supersat_frac_1(k_start:k_end) + r_param_array(:,47) = pdf_params%ice_supersat_frac_2(k_start:k_end) end subroutine pack_pdf_params - subroutine unpack_pdf_params(r_param_array, nz, pdf_params) + subroutine unpack_pdf_params(r_param_array, nz, pdf_params, & + k_start_in, k_end_in ) implicit none - ! Input a two dimensional real array with pdf values + integer, intent(in) :: nz ! Num Vert Model Levs + + ! Input a two dimensional real array with pdf values real (kind = core_rknd), dimension(nz,num_pdf_params), intent(in) :: & r_param_array ! Output a pdf_parameter array with nz instances of pdf_parameter - type (pdf_parameter), dimension(nz), intent(out) :: pdf_params - - ! Local Loop vars - integer :: k, p - ! temp var - real (kind = core_rknd) :: value + type (pdf_parameter), intent(inout) :: pdf_params - do k = 1,nz - do p = 1,num_pdf_params - - value = r_param_array(k,p) - call set_param_at_ind(pdf_params(k), p, value) - - end do ! p - end do ! k + integer, optional, intent(in) :: k_start_in, k_end_in + + integer :: k_start, k_end + + if( present( k_start_in ) .and. present( k_end_in ) ) then + k_start = k_start_in + k_end = k_end_in + else + k_start = 1 + k_end = nz + end if + + pdf_params%w_1(k_start:k_end) = r_param_array(:,1) + pdf_params%w_2(k_start:k_end) = r_param_array(:,2) + pdf_params%varnce_w_1(k_start:k_end) = r_param_array(:,3) + pdf_params%varnce_w_2(k_start:k_end) = r_param_array(:,4) + pdf_params%rt_1(k_start:k_end) = r_param_array(:,5) + pdf_params%rt_2(k_start:k_end) = r_param_array(:,6) + pdf_params%varnce_rt_1(k_start:k_end) = r_param_array(:,7) + pdf_params%varnce_rt_2(k_start:k_end)= r_param_array(:,8) + pdf_params%thl_1(k_start:k_end) = r_param_array(:,9) + pdf_params%thl_2(k_start:k_end) = r_param_array(:,10) + pdf_params%varnce_thl_1(k_start:k_end) = r_param_array(:,11) + pdf_params%varnce_thl_2(k_start:k_end) = r_param_array(:,12) + pdf_params%corr_w_rt_1(k_start:k_end) = r_param_array(:,13) + pdf_params%corr_w_rt_2(k_start:k_end) = r_param_array(:,14) + pdf_params%corr_w_thl_1(k_start:k_end) = r_param_array(:,15) + pdf_params%corr_w_thl_2(k_start:k_end) = r_param_array(:,16) + pdf_params%corr_rt_thl_1(k_start:k_end) = r_param_array(:,17) + pdf_params%corr_rt_thl_2(k_start:k_end) = r_param_array(:,18) + pdf_params%alpha_thl(k_start:k_end) = r_param_array(:,19) + pdf_params%alpha_rt(k_start:k_end) = r_param_array(:,20) + pdf_params%crt_1(k_start:k_end) = r_param_array(:,21) + pdf_params%crt_2(k_start:k_end) = r_param_array(:,22) + pdf_params%cthl_1(k_start:k_end) = r_param_array(:,23) + pdf_params%cthl_2(k_start:k_end) = r_param_array(:,24) + pdf_params%chi_1(k_start:k_end) = r_param_array(:,25) + pdf_params%chi_2(k_start:k_end) = r_param_array(:,26) + pdf_params%stdev_chi_1(k_start:k_end) = r_param_array(:,27) + pdf_params%stdev_chi_2(k_start:k_end) = r_param_array(:,28) + pdf_params%stdev_eta_1(k_start:k_end) = r_param_array(:,29) + pdf_params%stdev_eta_2(k_start:k_end) = r_param_array(:,30) + pdf_params%covar_chi_eta_1(k_start:k_end) = r_param_array(:,31) + pdf_params%covar_chi_eta_2(k_start:k_end) = r_param_array(:,32) + pdf_params%corr_w_chi_1(k_start:k_end) = r_param_array(:,33) + pdf_params%corr_w_chi_2(k_start:k_end) = r_param_array(:,34) + pdf_params%corr_w_eta_1(k_start:k_end) = r_param_array(:,35) + pdf_params%corr_w_eta_2(k_start:k_end) = r_param_array(:,36) + pdf_params%corr_chi_eta_1(k_start:k_end) = r_param_array(:,37) + pdf_params%corr_chi_eta_2(k_start:k_end) = r_param_array(:,38) + pdf_params%rsatl_1(k_start:k_end) = r_param_array(:,39) + pdf_params%rsatl_2(k_start:k_end) = r_param_array(:,40) + pdf_params%rc_1(k_start:k_end) = r_param_array(:,41) + pdf_params%rc_2(k_start:k_end) = r_param_array(:,42) + pdf_params%cloud_frac_1(k_start:k_end) = r_param_array(:,43) + pdf_params%cloud_frac_2(k_start:k_end) = r_param_array(:,44) + pdf_params%mixt_frac(k_start:k_end) = r_param_array(:,45) + pdf_params%ice_supersat_frac_1(k_start:k_end) = r_param_array(:,46) + pdf_params%ice_supersat_frac_2(k_start:k_end) = r_param_array(:,47) end subroutine unpack_pdf_params - real( kind = core_rknd ) function get_param_at_ind(pp_struct, ind) - - use constants_clubb, only: & - fstderr ! File I/O Constant - - use error_code, only : & - err_code, & ! Error Indicator - clubb_fatal_error ! Constant - - implicit none - type (pdf_parameter), intent(in) :: pp_struct - integer, intent(in) :: ind - - SELECT CASE (ind) - CASE (1) - get_param_at_ind = pp_struct%w_1 - CASE (2) - get_param_at_ind = pp_struct%w_2 - CASE (3) - get_param_at_ind = pp_struct%varnce_w_1 - CASE (4) - get_param_at_ind = pp_struct%varnce_w_2 - CASE (5) - get_param_at_ind = pp_struct%rt_1 - CASE (6) - get_param_at_ind = pp_struct%rt_2 - CASE (7) - get_param_at_ind = pp_struct%varnce_rt_1 - CASE (8) - get_param_at_ind = pp_struct%varnce_rt_2 - CASE (9) - get_param_at_ind = pp_struct%thl_1 - CASE (10) - get_param_at_ind = pp_struct%thl_2 - CASE (11) - get_param_at_ind = pp_struct%varnce_thl_1 - CASE (12) - get_param_at_ind = pp_struct%varnce_thl_2 - CASE (13) - get_param_at_ind = pp_struct%corr_w_rt_1 - CASE (14) - get_param_at_ind = pp_struct%corr_w_rt_2 - CASE (15) - get_param_at_ind = pp_struct%corr_w_thl_1 - CASE (16) - get_param_at_ind = pp_struct%corr_w_thl_2 - CASE (17) - get_param_at_ind = pp_struct%corr_rt_thl_1 - CASE (18) - get_param_at_ind = pp_struct%corr_rt_thl_2 - CASE (19) - get_param_at_ind = pp_struct%alpha_thl - CASE (20) - get_param_at_ind = pp_struct%alpha_rt - CASE (21) - get_param_at_ind = pp_struct%crt_1 - CASE (22) - get_param_at_ind = pp_struct%crt_2 - CASE (23) - get_param_at_ind = pp_struct%cthl_1 - CASE (24) - get_param_at_ind = pp_struct%cthl_2 - CASE (25) - get_param_at_ind = pp_struct%chi_1 - CASE (26) - get_param_at_ind = pp_struct%chi_2 - CASE (27) - get_param_at_ind = pp_struct%stdev_chi_1 - CASE (28) - get_param_at_ind = pp_struct%stdev_chi_2 - CASE (29) - get_param_at_ind = pp_struct%stdev_eta_1 - CASE (30) - get_param_at_ind = pp_struct%stdev_eta_2 - CASE (31) - get_param_at_ind = pp_struct%covar_chi_eta_1 - CASE (32) - get_param_at_ind = pp_struct%covar_chi_eta_2 - CASE (33) - get_param_at_ind = pp_struct%corr_w_chi_1 - CASE (34) - get_param_at_ind = pp_struct%corr_w_chi_2 - CASE (35) - get_param_at_ind = pp_struct%corr_w_eta_1 - CASE (36) - get_param_at_ind = pp_struct%corr_w_eta_2 - CASE (37) - get_param_at_ind = pp_struct%corr_chi_eta_1 - CASE (38) - get_param_at_ind = pp_struct%corr_chi_eta_2 - CASE (39) - get_param_at_ind = pp_struct%rsatl_1 - CASE (40) - get_param_at_ind = pp_struct%rsatl_2 - CASE (41) - get_param_at_ind = pp_struct%rc_1 - CASE (42) - get_param_at_ind = pp_struct%rc_2 - CASE (43) - get_param_at_ind = pp_struct%cloud_frac_1 - CASE (44) - get_param_at_ind = pp_struct%cloud_frac_2 - CASE (45) - get_param_at_ind = pp_struct%mixt_frac - CASE (46) - get_param_at_ind = pp_struct%ice_supersat_frac_1 - CASE (47) - get_param_at_ind = pp_struct%ice_supersat_frac_2 - CASE DEFAULT - write(fstderr,*) "Invalid index in get_param_at_ind" - err_code = clubb_fatal_error - return - END SELECT - - RETURN - end function get_param_at_ind - - subroutine set_param_at_ind(pp_struct, ind, val) - implicit none - type (pdf_parameter), intent(inout) :: pp_struct - integer, intent(in) :: ind - real (kind = core_rknd), intent(in) :: val - - SELECT CASE (ind) - CASE (1) - pp_struct%w_1 = val - CASE (2) - pp_struct%w_2 = val - CASE (3) - pp_struct%varnce_w_1 = val - CASE (4) - pp_struct%varnce_w_2 = val - CASE (5) - pp_struct%rt_1 = val - CASE (6) - pp_struct%rt_2 = val - CASE (7) - pp_struct%varnce_rt_1 = val - CASE (8) - pp_struct%varnce_rt_2 = val - CASE (9) - pp_struct%thl_1 = val - CASE (10) - pp_struct%thl_2 = val - CASE (11) - pp_struct%varnce_thl_1 = val - CASE (12) - pp_struct%varnce_thl_2 = val - CASE (13) - pp_struct%corr_w_rt_1 = val - CASE (14) - pp_struct%corr_w_rt_2 = val - CASE (15) - pp_struct%corr_w_thl_1 = val - CASE (16) - pp_struct%corr_w_thl_2 = val - CASE (17) - pp_struct%corr_rt_thl_1 = val - CASE (18) - pp_struct%corr_rt_thl_2 = val - CASE (19) - pp_struct%alpha_thl = val - CASE (20) - pp_struct%alpha_rt = val - CASE (21) - pp_struct%crt_1 = val - CASE (22) - pp_struct%crt_2 = val - CASE (23) - pp_struct%cthl_1 = val - CASE (24) - pp_struct%cthl_2 = val - CASE (25) - pp_struct%chi_1 = val - CASE (26) - pp_struct%chi_2 = val - CASE (27) - pp_struct%stdev_chi_1 = val - CASE (28) - pp_struct%stdev_chi_2 = val - CASE (29) - pp_struct%stdev_eta_1 = val - CASE (30) - pp_struct%stdev_eta_2 = val - CASE (31) - pp_struct%covar_chi_eta_1 = val - CASE (32) - pp_struct%covar_chi_eta_2 = val - CASE (33) - pp_struct%corr_w_chi_1 = val - CASE (34) - pp_struct%corr_w_chi_2 = val - CASE (35) - pp_struct%corr_w_eta_1 = val - CASE (36) - pp_struct%corr_w_eta_2 = val - CASE (37) - pp_struct%corr_chi_eta_1 = val - CASE (38) - pp_struct%corr_chi_eta_2 = val - CASE (39) - pp_struct%rsatl_1 = val - CASE (40) - pp_struct%rsatl_2 = val - CASE (41) - pp_struct%rc_1 = val - CASE (42) - pp_struct%rc_2 = val - CASE (43) - pp_struct%cloud_frac_1 = val - CASE (44) - pp_struct%cloud_frac_2 = val - CASE (45) - pp_struct%mixt_frac = val - CASE (46) - pp_struct%ice_supersat_frac_1 = val - CASE (47) - pp_struct%ice_supersat_frac_2 = val - CASE DEFAULT - ! do nothing ! - END SELECT - - end subroutine set_param_at_ind - !#endif /* CLUBB_CAM */ end module pdf_parameter_module diff --git a/components/cam/src/physics/clubb/setup_clubb_pdf_params.F90 b/components/cam/src/physics/clubb/setup_clubb_pdf_params.F90 index ac46f07b64a1..c86ab29a6d76 100644 --- a/components/cam/src/physics/clubb/setup_clubb_pdf_params.F90 +++ b/components/cam/src/physics/clubb/setup_clubb_pdf_params.F90 @@ -183,7 +183,7 @@ subroutine setup_pdf_parameters( nz, pdf_dim, dt, & ! Inten corr_array_n_cloud, & ! Prescribed normal space corr. array in cloud [-] corr_array_n_below ! Prescribed normal space corr. array below cl. [-] - type(pdf_parameter), dimension(nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params ! PDF parameters [units vary] logical, intent(in) :: & @@ -528,7 +528,19 @@ subroutine setup_pdf_parameters( nz, pdf_dim, dt, & ! Inten Ncnm, mixt_frac, precip_frac, & ! Intent(in) precip_frac_1, precip_frac_2, & ! Intent(in) precip_frac_tol, & ! Intent(in) - pdf_params, pdf_dim, & ! Intent(in) + pdf_params%w_1, & ! Intent(in) + pdf_params%w_2, & ! Intent(in) + pdf_params%varnce_w_1, & ! Intent(in) + pdf_params%varnce_w_2, & ! Intent(in) + pdf_params%chi_1, & ! Intent(in) + pdf_params%chi_2, & ! Intent(in) + pdf_params%stdev_chi_1, & ! Intent(in) + pdf_params%stdev_chi_2, & ! Intent(in) + pdf_params%stdev_eta_1, & ! Intent(in) + pdf_params%stdev_eta_2, & ! Intent(in) + pdf_params%thl_1, & ! Intent(in) + pdf_params%thl_2, & ! Intent(in) + pdf_dim, & ! Intent(in) mu_x_1, mu_x_2, & ! Intent(out) sigma_x_1, sigma_x_2, & ! Intent(out) hm_1, hm_2, & ! Intent(out) @@ -587,7 +599,13 @@ subroutine setup_pdf_parameters( nz, pdf_dim, dt, & ! Inten sigma_x_1, sigma_x_2, & sigma_x_1_n, sigma_x_2_n, & corr_array_n_cloud, corr_array_n_below, & - pdf_params, pdf_dim, & + pdf_params%corr_chi_eta_1, & + pdf_params%corr_chi_eta_2, & + pdf_params%corr_w_chi_1, & + pdf_params%corr_w_chi_2, & + pdf_params%corr_w_eta_1, & + pdf_params%corr_w_eta_2, & + pdf_dim, & corr_array_1_n, corr_array_2_n ) endif ! l_diagnose_correlations @@ -671,10 +689,19 @@ subroutine setup_pdf_parameters( nz, pdf_dim, dt, & ! Inten if ( l_stats_samp ) then if ( irtp2_from_chi > 0 ) then + rtp2_zt_from_chi & - = compute_rtp2_from_chi( pdf_params(:), & - corr_array_1_n(iiPDF_chi,iiPDF_eta,:), & - corr_array_2_n(iiPDF_chi,iiPDF_eta,:) ) + = compute_rtp2_from_chi( pdf_params%stdev_chi_1(:), pdf_params%stdev_chi_2(:), & + pdf_params%stdev_eta_1(:), pdf_params%stdev_eta_2(:), & + pdf_params%rt_1(:), pdf_params%rt_2(:), & + pdf_params%crt_1(:), pdf_params%crt_2(:), & + pdf_params%mixt_frac(:), & + corr_array_1_n(iiPDF_chi,iiPDF_eta,:), & + corr_array_2_n(iiPDF_chi,iiPDF_eta,:) ) + + + + call stat_update_var( irtp2_from_chi, zt2zm( rtp2_zt_from_chi ), & stats_zm ) endif @@ -709,7 +736,13 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) Ncnm, mixt_frac, precip_frac, & ! Intent(in) precip_frac_1, precip_frac_2, & ! Intent(in) precip_frac_tol, & ! Intent(in) - pdf_params, pdf_dim, & ! Intent(in) + w_1, w_2, & ! Intent(in) + varnce_w_1, varnce_w_2, & ! Intent(in) + chi_1, chi_2, & ! Intent(in) + stdev_chi_1, stdev_chi_2, & ! Intent(in) + stdev_eta_1, stdev_eta_2, & ! Intent(in) + thl_1, thl_2, & ! Intent(in) + pdf_dim, & ! Intent(in) mu_x_1, mu_x_2, & ! Intent(out) sigma_x_1, sigma_x_2, & ! Intent(out) hm_1, hm_2, & ! Intent(out) @@ -744,9 +777,6 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) use index_mapping, only: & pdf2hydromet_idx ! Procedure(s) - use pdf_parameter_module, only: & - pdf_parameter ! Variable(s) type - use parameters_tunable, only: & omicron, & ! Variable(s) zeta_vrnce_rat @@ -778,11 +808,16 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) precip_frac_1, & ! Precipitation fraction (1st PDF component) [-] precip_frac_2 ! Precipitation fraction (2nd PDF component) [-] - real( kind = core_rknd ), intent(in) :: & - precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-] + real( kind = core_rknd ), dimension(nz), intent(in) :: & + w_1, w_2, & ! Mean of w [m/s] + varnce_w_1, varnce_w_2, & ! Variance of w [m^2/s^2] + chi_1, chi_2, & ! Mean of chi [kg/kg] + stdev_chi_1, stdev_chi_2, & ! Standard dev. of chi [kg/kg] + stdev_eta_1, stdev_eta_2, & ! Standard dev. of eta [kg/kg] + thl_1, thl_2 ! Mean of th_l [K] - type(pdf_parameter), dimension(nz), intent(in) :: & - pdf_params ! PDF parameters [units vary] + real( kind = core_rknd ), intent(in) :: & + precip_frac_tol ! Minimum precip. frac. when hydromet. are present [-] integer, intent(in) :: & pdf_dim ! Number of PDF variables @@ -828,35 +863,35 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) !!! Vertical velocity, w. ! Mean of vertical velocity, w, in PDF component 1. - mu_x_1(iiPDF_w,:) = pdf_params%w_1 + mu_x_1(iiPDF_w,:) = w_1 ! Mean of vertical velocity, w, in PDF component 2. - mu_x_2(iiPDF_w,:) = pdf_params%w_2 + mu_x_2(iiPDF_w,:) = w_2 ! Standard deviation of vertical velocity, w, in PDF component 1. - sigma_x_1(iiPDF_w,:) = sqrt( pdf_params%varnce_w_1 ) + sigma_x_1(iiPDF_w,:) = sqrt( varnce_w_1 ) ! Standard deviation of vertical velocity, w, in PDF component 2. - sigma_x_2(iiPDF_w,:) = sqrt( pdf_params%varnce_w_2 ) + sigma_x_2(iiPDF_w,:) = sqrt( varnce_w_2 ) !!! Extended liquid water mixing ratio, chi. ! Mean of extended liquid water mixing ratio, chi (old s), ! in PDF component 1. - mu_x_1(iiPDF_chi,:) = pdf_params%chi_1 + mu_x_1(iiPDF_chi,:) = chi_1 ! Mean of extended liquid water mixing ratio, chi (old s), ! in PDF component 2. - mu_x_2(iiPDF_chi,:) = pdf_params%chi_2 + mu_x_2(iiPDF_chi,:) = chi_2 ! Standard deviation of extended liquid water mixing ratio, chi (old s), ! in PDF component 1. - sigma_x_1(iiPDF_chi,:) = pdf_params%stdev_chi_1 + sigma_x_1(iiPDF_chi,:) = stdev_chi_1 ! Standard deviation of extended liquid water mixing ratio, chi (old s), ! in PDF component 2. - sigma_x_2(iiPDF_chi,:) = pdf_params%stdev_chi_2 + sigma_x_2(iiPDF_chi,:) = stdev_chi_2 !!! Coordinate orthogonal to chi, eta. @@ -876,10 +911,10 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) mu_x_2(iiPDF_eta,:) = zero ! Standard deviation of eta (old t) in PDF component 1. - sigma_x_1(iiPDF_eta,:) = pdf_params%stdev_eta_1 + sigma_x_1(iiPDF_eta,:) = stdev_eta_1 ! Standard deviation of eta (old t) in PDF component 2. - sigma_x_2(iiPDF_eta,:) = pdf_params%stdev_eta_2 + sigma_x_2(iiPDF_eta,:) = stdev_eta_2 !!! Simplified cloud nuclei concentration, Ncn. @@ -934,7 +969,7 @@ subroutine compute_mean_stdev( nz, hydromet, hydrometp2_zt, & ! Intent(in) mixt_frac, precip_frac, & precip_frac_1, precip_frac_2, & hydromet_tol(hm_idx), precip_frac_tol, & - pdf_params%thl_1, pdf_params%thl_2, & + thl_1, thl_2, & omicron, zeta_vrnce_rat, & mu_x_1(ivar,:), mu_x_2(ivar,:), & sigma_x_1(ivar,:), sigma_x_2(ivar,:), & @@ -957,7 +992,13 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & mu_x_1, mu_x_2, sigma_x_1, sigma_x_2, & sigma_x_1_n, sigma_x_2_n, & corr_array_n_cloud, corr_array_n_below, & - pdf_params, pdf_dim, & + corr_chi_eta_1, & + corr_chi_eta_2, & + corr_w_chi_1, & + corr_w_chi_2, & + corr_w_eta_1, & + corr_w_eta_2, & + pdf_dim, & corr_array_1_n, corr_array_2_n ) ! Description: @@ -986,9 +1027,6 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & use clubb_precision, only: & core_rknd ! Variable(s) - use pdf_parameter_module, only: & - pdf_parameter ! Variable(s) type - use array_index, only: & iiPDF_chi, & ! Variable(s) iiPDF_eta, & @@ -1030,8 +1068,13 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & corr_array_n_cloud, & ! Prescribed correlation array in cloud [-] corr_array_n_below ! Prescribed correlation array below cloud [-] - type(pdf_parameter), dimension(nz), intent(in) :: & - pdf_params ! PDF parameters [units vary] + real( kind = core_rknd ), dimension(nz), intent(in) :: & + corr_chi_eta_1, & + corr_chi_eta_2, & + corr_w_chi_1, & + corr_w_chi_2, & + corr_w_eta_1, & + corr_w_eta_2 ! Output Variables real( kind = core_rknd ), dimension(pdf_dim,pdf_dim,nz), & @@ -1054,7 +1097,7 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & logical :: & l_limit_corr_chi_eta ! Flag to limit the correlation of chi and eta [-] - integer :: ivar, jvar ! Loop iterators + integer :: ivar, jvar, hm_idx ! Indices ! ---- Begin Code ---- @@ -1090,8 +1133,10 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & ! hydrometeor for each PDF component and each hydrometeor type. do jvar = iiPDF_Ncn+1, pdf_dim + hm_idx = pdf2hydromet_idx(jvar) + call calc_corr_w_hm_n( nz, wm_zt, & - wphydrometp_zt(:,pdf2hydromet_idx(jvar)), & + wphydrometp_zt(:,hm_idx), & mu_x_1(iiPDF_w,:), mu_x_2(iiPDF_w,:), & mu_x_1(jvar,:), mu_x_2(jvar,:), & sigma_x_1(iiPDF_w,:), sigma_x_2(iiPDF_w,:), & @@ -1099,7 +1144,7 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & sigma_x_1_n(jvar,:), sigma_x_2_n(jvar,:), & mixt_frac, precip_frac_1, precip_frac_2, & corr_w_hm_1_n(jvar,:), corr_w_hm_2_n(jvar,:), & - hydromet_tol(pdf2hydromet_idx(jvar)) ) + hydromet_tol(hm_idx) ) enddo ! jvar = iiPDF_Ncn+1, pdf_dim @@ -1129,14 +1174,14 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & ! Correlation of chi (old s) and eta (old t) corr_array_1_n(iiPDF_eta,iiPDF_chi,:) & - = component_corr_chi_eta( nz, pdf_params%corr_chi_eta_1, & + = component_corr_chi_eta( nz, corr_chi_eta_1, & rc_1, cloud_frac_1, & corr_array_n_cloud(iiPDF_eta,iiPDF_chi), & corr_array_n_below(iiPDF_eta,iiPDF_chi), & l_limit_corr_chi_eta ) corr_array_2_n(iiPDF_eta,iiPDF_chi,:) & - = component_corr_chi_eta( nz, pdf_params%corr_chi_eta_2, & + = component_corr_chi_eta( nz, corr_chi_eta_2, & rc_2, cloud_frac_2, & corr_array_n_cloud(iiPDF_eta,iiPDF_chi), & corr_array_n_below(iiPDF_eta,iiPDF_chi), & @@ -1144,12 +1189,12 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & ! Correlation of chi (old s) and w corr_array_1_n(iiPDF_w,iiPDF_chi,:) & - = component_corr_w_x( nz, pdf_params%corr_w_chi_1, rc_1, cloud_frac_1, & + = component_corr_w_x( nz, corr_w_chi_1, rc_1, cloud_frac_1, & corr_array_n_cloud(iiPDF_w,iiPDF_chi), & corr_array_n_below(iiPDF_w,iiPDF_chi) ) corr_array_2_n(iiPDF_w,iiPDF_chi,:) & - = component_corr_w_x( nz, pdf_params%corr_w_chi_2, rc_2, cloud_frac_2, & + = component_corr_w_x( nz, corr_w_chi_2, rc_2, cloud_frac_2, & corr_array_n_cloud(iiPDF_w,iiPDF_chi), & corr_array_n_below(iiPDF_w,iiPDF_chi) ) @@ -1181,12 +1226,12 @@ subroutine comp_corr_norm( nz, wm_zt, rc_1, rc_2, cloud_frac_1, & ! Correlation of eta (old t) and w corr_array_1_n(iiPDF_w,iiPDF_eta,:) & - = component_corr_w_x( nz, pdf_params%corr_w_eta_1, rc_1, cloud_frac_1, & + = component_corr_w_x( nz, corr_w_eta_1, rc_1, cloud_frac_1, & corr_array_n_cloud(iiPDF_w,iiPDF_eta), & corr_array_n_below(iiPDF_w,iiPDF_eta) ) corr_array_2_n(iiPDF_w,iiPDF_eta,:) & - = component_corr_w_x( nz, pdf_params%corr_w_eta_2, rc_2, cloud_frac_2, & + = component_corr_w_x( nz, corr_w_eta_2, rc_2, cloud_frac_2, & corr_array_n_cloud(iiPDF_w,iiPDF_eta), & corr_array_n_below(iiPDF_w,iiPDF_eta) ) @@ -2638,7 +2683,7 @@ subroutine norm_transform_mean_stdev( nz, hm_1, hm_2, & sigma_x_2_n ! Std. dev. array (normal space): PDF vars (comp. 2) [u.v.] ! Local Variable - integer :: ivar ! Loop index + integer :: ivar, hm_idx ! Indices ! The means and standard deviations in each PDF component of w, chi (old s), @@ -2714,10 +2759,11 @@ subroutine norm_transform_mean_stdev( nz, hm_1, hm_2, & ! Normal space precipitating hydrometeor means and standard deviations. do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Normal space mean of a precipitating hydrometeor, hm, in PDF ! component 1. - where ( hm_1(:,pdf2hydromet_idx(ivar)) & - >= hydromet_tol(pdf2hydromet_idx(ivar)) ) + where ( hm_1(:,hm_idx) >= hydromet_tol(hm_idx) ) mu_x_1_n(ivar,:) = mean_L2N( nz, mu_x_1(ivar,:), & sigma2_on_mu2_ip_1(ivar,:) ) @@ -2741,8 +2787,7 @@ subroutine norm_transform_mean_stdev( nz, hm_1, hm_2, & ! Normal space mean of a precipitating hydrometeor, hm, in PDF ! component 2. - where ( hm_2(:,pdf2hydromet_idx(ivar)) & - >= hydromet_tol(pdf2hydromet_idx(ivar)) ) + where ( hm_2(:,hm_idx) >= hydromet_tol(hm_idx) ) mu_x_2_n(ivar,:) = mean_L2N( nz, mu_x_2(ivar,:), & sigma2_on_mu2_ip_2(ivar,:) ) @@ -3264,7 +3309,7 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & l_stats_samp ! Flag to record statistical output. ! Local Variable - integer :: ivar, jvar ! Loop indices + integer :: ivar, jvar, hm_idx, hm_idx_ivar, hm_idx_jvar ! Indices !!! Output the statistics for hydrometeor PDF parameters. @@ -3288,18 +3333,18 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Mean of the precipitating hydrometeor (in-precip) ! in PDF component 1. - if ( imu_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( imu_hm_1(pdf2hydromet_idx(ivar)), & - mu_x_1(ivar,:), stats_zt ) + if ( imu_hm_1(hm_idx) > 0 ) then + call stat_update_var( imu_hm_1(hm_idx), mu_x_1(ivar,:), stats_zt ) endif ! Mean of the precipitating hydrometeor (in-precip) ! in PDF component 2. - if ( imu_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( imu_hm_2(pdf2hydromet_idx(ivar)), & - mu_x_2(ivar,:), stats_zt ) + if ( imu_hm_2(hm_idx) > 0 ) then + call stat_update_var( imu_hm_2(hm_idx), mu_x_2(ivar,:), stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 @@ -3316,18 +3361,20 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Standard deviation of the precipitating hydrometeor (in-precip) ! in PDF component 1. - if ( isigma_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( isigma_hm_1(pdf2hydromet_idx(ivar)), & - sigma_x_1(ivar,:), stats_zt ) + if ( isigma_hm_1(hm_idx) > 0 ) then + call stat_update_var( isigma_hm_1(hm_idx), sigma_x_1(ivar,:), & + stats_zt ) endif ! Standard deviation of the precipitating hydrometeor (in-precip) ! in PDF component 2. - if ( isigma_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( isigma_hm_2(pdf2hydromet_idx(ivar)), & - sigma_x_2(ivar,:), stats_zt ) + if ( isigma_hm_2(hm_idx) > 0 ) then + call stat_update_var( isigma_hm_2(hm_idx), sigma_x_2(ivar,:), & + stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 @@ -3404,17 +3451,19 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of w and the precipitating hydrometeor ! in PDF component 1. - if ( icorr_w_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_w_hm_1(pdf2hydromet_idx(ivar)), & + if ( icorr_w_hm_1(hm_idx) > 0 ) then + call stat_update_var( icorr_w_hm_1(hm_idx), & corr_array_1(ivar,iiPDF_w,:), stats_zt ) endif ! Correlation (in-precip) of w and the precipitating hydrometeor ! in PDF component 2. - if ( icorr_w_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_w_hm_2(pdf2hydromet_idx(ivar)), & + if ( icorr_w_hm_2(hm_idx) > 0 ) then + call stat_update_var( icorr_w_hm_2(hm_idx), & corr_array_2(ivar,iiPDF_w,:), stats_zt ) endif @@ -3466,17 +3515,19 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of chi (old s) and the precipitating ! hydrometeor in PDF component 1. - if ( icorr_chi_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_chi_hm_1(pdf2hydromet_idx(ivar)), & + if ( icorr_chi_hm_1(hm_idx) > 0 ) then + call stat_update_var( icorr_chi_hm_1(hm_idx), & corr_array_1(ivar,iiPDF_chi,:), stats_zt ) endif ! Correlation (in-precip) of chi (old s) and the precipitating ! hydrometeor in PDF component 2. - if ( icorr_chi_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_chi_hm_2(pdf2hydromet_idx(ivar)), & + if ( icorr_chi_hm_2(hm_idx) > 0 ) then + call stat_update_var( icorr_chi_hm_2(hm_idx), & corr_array_2(ivar,iiPDF_chi,:), stats_zt ) endif @@ -3496,17 +3547,19 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of eta (old t) and the precipitating ! hydrometeor in PDF component 1. - if ( icorr_eta_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_eta_hm_1(pdf2hydromet_idx(ivar)), & + if ( icorr_eta_hm_1(hm_idx) > 0 ) then + call stat_update_var( icorr_eta_hm_1(hm_idx), & corr_array_1(ivar,iiPDF_eta,:), stats_zt ) endif ! Correlation (in-precip) of eta (old t) and the precipitating ! hydrometeor in PDF component 2. - if ( icorr_eta_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_eta_hm_2(pdf2hydromet_idx(ivar)), & + if ( icorr_eta_hm_2(hm_idx) > 0 ) then + call stat_update_var( icorr_eta_hm_2(hm_idx), & corr_array_2(ivar,iiPDF_eta,:), stats_zt ) endif @@ -3526,44 +3579,47 @@ subroutine pdf_param_hm_stats( nz, pdf_dim, hm_1, hm_2, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of N_cn and the precipitating ! hydrometeor in PDF component 1. - if ( icorr_Ncn_hm_1(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_Ncn_hm_1(pdf2hydromet_idx(ivar)), & + if ( icorr_Ncn_hm_1(hm_idx) > 0 ) then + call stat_update_var( icorr_Ncn_hm_1(hm_idx), & corr_array_1(ivar,iiPDF_Ncn,:), stats_zt ) endif ! Correlation (in-precip) of N_cn and the precipitating ! hydrometeor in PDF component 2. - if ( icorr_Ncn_hm_2(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_Ncn_hm_2(pdf2hydromet_idx(ivar)), & + if ( icorr_Ncn_hm_2(hm_idx) > 0 ) then + call stat_update_var( icorr_Ncn_hm_2(hm_idx), & corr_array_2(ivar,iiPDF_Ncn,:), stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 do ivar = iiPDF_Ncn+1, pdf_dim, 1 - do jvar = ivar+1, pdf_dim, 1 - - ! Correlation (in-precip) of two different hydrometeors (hmx and - ! hmy) in PDF component 1. - if ( icorr_hmx_hmy_1(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)) & - > 0 ) then - call stat_update_var( & - icorr_hmx_hmy_1(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)), & - corr_array_1(jvar,ivar,:), stats_zt ) - endif - - ! Correlation (in-precip) of two different hydrometeors (hmx and - ! hmy) in PDF component 2. - if ( icorr_hmx_hmy_2(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)) & - > 0 ) then - call stat_update_var( & - icorr_hmx_hmy_2(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)), & - corr_array_2(jvar,ivar,:), stats_zt ) - endif - - enddo ! jvar = ivar+1, pdf_dim, 1 + + hm_idx_ivar = pdf2hydromet_idx(ivar) + + do jvar = ivar+1, pdf_dim, 1 + + hm_idx_jvar = pdf2hydromet_idx(jvar) + + ! Correlation (in-precip) of two different hydrometeors (hmx and + ! hmy) in PDF component 1. + if ( icorr_hmx_hmy_1(hm_idx_jvar,hm_idx_ivar) > 0 ) then + call stat_update_var( icorr_hmx_hmy_1(hm_idx_jvar,hm_idx_ivar), & + corr_array_1(jvar,ivar,:), stats_zt ) + endif + + ! Correlation (in-precip) of two different hydrometeors (hmx and + ! hmy) in PDF component 2. + if ( icorr_hmx_hmy_2(hm_idx_jvar,hm_idx_ivar) > 0 ) then + call stat_update_var( icorr_hmx_hmy_2(hm_idx_jvar,hm_idx_ivar), & + corr_array_2(jvar,ivar,:), stats_zt ) + endif + + enddo ! jvar = ivar+1, pdf_dim, 1 enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 endif ! l_stats_samp @@ -3657,7 +3713,7 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & mu_Ncn_1_n, & ! Mean of ln Ncn (1st PDF component) [ln(num/kg)] mu_Ncn_2_n ! Mean of ln Ncn (2nd PDF component) [ln(num/kg)] - integer :: ivar, jvar ! Loop indices + integer :: ivar, jvar, hm_idx, hm_idx_ivar, hm_idx_jvar ! Indices !!! Output the statistics for normal space hydrometeor PDF parameters. @@ -3667,8 +3723,10 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Mean (in-precip) of ln hm in PDF component 1. - if ( imu_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then + if ( imu_hm_1_n(hm_idx) > 0 ) then where ( mu_x_1_n(ivar,:) > real( -huge( 0.0 ), kind = core_rknd ) ) mu_hm_1_n = mu_x_1_n(ivar,:) elsewhere @@ -3679,12 +3737,11 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & ! Set to -huge for single precision. mu_hm_1_n = real( -huge( 0.0 ), kind = core_rknd ) endwhere - call stat_update_var( imu_hm_1_n(pdf2hydromet_idx(ivar)), & - mu_hm_1_n, stats_zt ) + call stat_update_var( imu_hm_1_n(hm_idx), mu_hm_1_n, stats_zt ) endif ! Mean (in-precip) of ln hm in PDF component 2. - if ( imu_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then + if ( imu_hm_2_n(hm_idx) > 0 ) then where ( mu_x_2_n(ivar,:) > real( -huge( 0.0 ), kind = core_rknd ) ) mu_hm_2_n = mu_x_2_n(ivar,:) elsewhere @@ -3695,8 +3752,7 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & ! Set to -huge for single precision. mu_hm_2_n = real( -huge( 0.0 ), kind = core_rknd ) endwhere - call stat_update_var( imu_hm_2_n(pdf2hydromet_idx(ivar)), & - mu_hm_2_n, stats_zt ) + call stat_update_var( imu_hm_2_n(hm_idx), mu_hm_2_n, stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 @@ -3735,16 +3791,18 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Standard deviation (in-precip) of ln hm in PDF component 1. - if ( isigma_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( isigma_hm_1_n(pdf2hydromet_idx(ivar)), & - sigma_x_1_n(ivar,:), stats_zt ) + if ( isigma_hm_1_n(hm_idx) > 0 ) then + call stat_update_var( isigma_hm_1_n(hm_idx), sigma_x_1_n(ivar,:), & + stats_zt ) endif ! Standard deviation (in-precip) of ln hm in PDF component 2. - if ( isigma_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( isigma_hm_2_n(pdf2hydromet_idx(ivar)), & - sigma_x_2_n(ivar,:), stats_zt ) + if ( isigma_hm_2_n(hm_idx) > 0 ) then + call stat_update_var( isigma_hm_2_n(hm_idx), sigma_x_2_n(ivar,:), & + stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 @@ -3763,15 +3821,17 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of w and ln hm in PDF component 1. - if ( icorr_w_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_w_hm_1_n(pdf2hydromet_idx(ivar)), & + if ( icorr_w_hm_1_n(hm_idx) > 0 ) then + call stat_update_var( icorr_w_hm_1_n(hm_idx), & corr_array_1_n(ivar,iiPDF_w,:), stats_zt ) endif ! Correlation (in-precip) of w and ln hm in PDF component 2. - if ( icorr_w_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_w_hm_2_n(pdf2hydromet_idx(ivar)), & + if ( icorr_w_hm_2_n(hm_idx) > 0 ) then + call stat_update_var( icorr_w_hm_2_n(hm_idx), & corr_array_2_n(ivar,iiPDF_w,:), stats_zt ) endif @@ -3791,15 +3851,17 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of chi (old s) and ln hm in PDF component 1. - if ( icorr_chi_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_chi_hm_1_n(pdf2hydromet_idx(ivar)), & + if ( icorr_chi_hm_1_n(hm_idx) > 0 ) then + call stat_update_var( icorr_chi_hm_1_n(hm_idx), & corr_array_1_n(ivar,iiPDF_chi,:), stats_zt ) endif ! Correlation (in-precip) of chi( old s) and ln hm in PDF component 2. - if ( icorr_chi_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_chi_hm_2_n(pdf2hydromet_idx(ivar)), & + if ( icorr_chi_hm_2_n(hm_idx) > 0 ) then + call stat_update_var( icorr_chi_hm_2_n(hm_idx), & corr_array_2_n(ivar,iiPDF_chi,:), stats_zt ) endif @@ -3821,15 +3883,17 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of eta (old t) and ln hm in PDF component 1. - if ( icorr_eta_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_eta_hm_1_n(pdf2hydromet_idx(ivar)), & + if ( icorr_eta_hm_1_n(hm_idx) > 0 ) then + call stat_update_var( icorr_eta_hm_1_n(hm_idx), & corr_array_1_n(ivar,iiPDF_eta,:), stats_zt ) endif ! Correlation (in-precip) of eta (old t) and ln hm in PDF component 2. - if ( icorr_eta_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_eta_hm_2_n(pdf2hydromet_idx(ivar)), & + if ( icorr_eta_hm_2_n(hm_idx) > 0 ) then + call stat_update_var( icorr_eta_hm_2_n(hm_idx), & corr_array_2_n(ivar,iiPDF_eta,:), stats_zt ) endif @@ -3851,44 +3915,49 @@ subroutine pdf_param_ln_hm_stats( nz, pdf_dim, mu_x_1_n, & do ivar = iiPDF_Ncn+1, pdf_dim, 1 + hm_idx = pdf2hydromet_idx(ivar) + ! Correlation (in-precip) of ln N_cn and ln hm in PDF ! component 1. - if ( icorr_Ncn_hm_1_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_Ncn_hm_1_n(pdf2hydromet_idx(ivar)), & + if ( icorr_Ncn_hm_1_n(hm_idx) > 0 ) then + call stat_update_var( icorr_Ncn_hm_1_n(hm_idx), & corr_array_1_n(ivar,iiPDF_Ncn,:), stats_zt ) endif ! Correlation (in-precip) of ln N_cn and ln hm in PDF ! component 2. - if ( icorr_Ncn_hm_2_n(pdf2hydromet_idx(ivar)) > 0 ) then - call stat_update_var( icorr_Ncn_hm_2_n(pdf2hydromet_idx(ivar)), & + if ( icorr_Ncn_hm_2_n(hm_idx) > 0 ) then + call stat_update_var( icorr_Ncn_hm_2_n(hm_idx), & corr_array_2_n(ivar,iiPDF_Ncn,:), stats_zt ) endif enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 do ivar = iiPDF_Ncn+1, pdf_dim, 1 - do jvar = ivar+1, pdf_dim, 1 - - ! Correlation (in-precip) of ln hmx and ln hmy (two different - ! hydrometeors) in PDF component 1. - if (icorr_hmx_hmy_1_n(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar))& - > 0 ) then - call stat_update_var( & - icorr_hmx_hmy_1_n(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)), & - corr_array_1_n(jvar,ivar,:), stats_zt ) - endif - - ! Correlation (in-precip) of ln hmx and ln hmy (two different - ! hydrometeors) in PDF component 2. - if (icorr_hmx_hmy_2_n(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar))& - > 0 ) then - call stat_update_var( & - icorr_hmx_hmy_2_n(pdf2hydromet_idx(jvar),pdf2hydromet_idx(ivar)), & - corr_array_2_n(jvar,ivar,:), stats_zt ) - endif - - enddo ! jvar = ivar+1, pdf_dim, 1 + + hm_idx_ivar = pdf2hydromet_idx(ivar) + + do jvar = ivar+1, pdf_dim, 1 + + hm_idx_jvar= pdf2hydromet_idx(jvar) + + ! Correlation (in-precip) of ln hmx and ln hmy (two different + ! hydrometeors) in PDF component 1. + if ( icorr_hmx_hmy_1_n(hm_idx_jvar,hm_idx_ivar) > 0 ) then + call stat_update_var( & + icorr_hmx_hmy_1_n(hm_idx_jvar,hm_idx_ivar), & + corr_array_1_n(jvar,ivar,:), stats_zt ) + endif + + ! Correlation (in-precip) of ln hmx and ln hmy (two different + ! hydrometeors) in PDF component 2. + if ( icorr_hmx_hmy_2_n(hm_idx_jvar,hm_idx_ivar) > 0 ) then + call stat_update_var( & + icorr_hmx_hmy_2_n(hm_idx_jvar,hm_idx_ivar), & + corr_array_2_n(jvar,ivar,:), stats_zt ) + endif + + enddo ! jvar = ivar+1, pdf_dim, 1 enddo ! ivar = iiPDF_Ncn+1, pdf_dim, 1 endif ! l_stats_samp @@ -3967,58 +4036,58 @@ subroutine pack_hydromet_pdf_params( nz, hm_1, hm_2, pdf_dim, mu_x_1, & ! In hydromet_pdf_params ! Hydrometeor PDF parameters [units vary] ! Local Variables - integer :: ivar, jvar ! Loop indices + integer :: ivar, jvar, pdf_idx ! Indices ! Pack remaining means and standard deviations into hydromet_pdf_params. do ivar = 1, hydromet_dim, 1 + pdf_idx = hydromet2pdf_idx(ivar) + ! Mean of a hydrometeor (overall) in the 1st PDF component. hydromet_pdf_params%hm_1(ivar) = hm_1(:,ivar) ! Mean of a hydrometeor (overall) in the 2nd PDF component. hydromet_pdf_params%hm_2(ivar) = hm_2(:,ivar) ! Mean of a hydrometeor (in-precip) in the 1st PDF component. - hydromet_pdf_params%mu_hm_1(ivar) = mu_x_1(hydromet2pdf_idx(ivar),:) + hydromet_pdf_params%mu_hm_1(ivar) = mu_x_1(pdf_idx,:) ! Mean of a hydrometeor (in-precip) in the 2nd PDF component. - hydromet_pdf_params%mu_hm_2(ivar) = mu_x_2(hydromet2pdf_idx(ivar),:) + hydromet_pdf_params%mu_hm_2(ivar) = mu_x_2(pdf_idx,:) ! Standard deviation of a hydrometeor (in-precip) in the ! 1st PDF component. - hydromet_pdf_params%sigma_hm_1(ivar) = sigma_x_1(hydromet2pdf_idx(ivar),:) + hydromet_pdf_params%sigma_hm_1(ivar) = sigma_x_1(pdf_idx,:) ! Standard deviation of a hydrometeor (in-precip) in the ! 2nd PDF component. - hydromet_pdf_params%sigma_hm_2(ivar) = sigma_x_2(hydromet2pdf_idx(ivar),:) + hydromet_pdf_params%sigma_hm_2(ivar) = sigma_x_2(pdf_idx,:) ! Correlation (in-precip) of w and a hydrometeor in the 1st PDF ! component. - hydromet_pdf_params%corr_w_hm_1(ivar) & - = corr_array_1( hydromet2pdf_idx(ivar), iiPDF_w, : ) + hydromet_pdf_params%corr_w_hm_1(ivar) = corr_array_1(pdf_idx,iiPDF_w,:) ! Correlation (in-precip) of w and a hydrometeor in the 2nd PDF ! component. - hydromet_pdf_params%corr_w_hm_2(ivar) & - = corr_array_2( hydromet2pdf_idx(ivar), iiPDF_w, : ) + hydromet_pdf_params%corr_w_hm_2(ivar) = corr_array_2(pdf_idx,iiPDF_w,:) ! Correlation (in-precip) of chi and a hydrometeor in the 1st PDF ! component. hydromet_pdf_params%corr_chi_hm_1(ivar) & - = corr_array_1( hydromet2pdf_idx(ivar), iiPDF_chi, : ) + = corr_array_1(pdf_idx,iiPDF_chi,:) ! Correlation (in-precip) of chi and a hydrometeor in the 2nd PDF ! component. hydromet_pdf_params%corr_chi_hm_2(ivar) & - = corr_array_2( hydromet2pdf_idx(ivar), iiPDF_chi, : ) + = corr_array_2(pdf_idx,iiPDF_chi,:) ! Correlation (in-precip) of eta and a hydrometeor in the 1st PDF ! component. hydromet_pdf_params%corr_eta_hm_1(ivar) & - = corr_array_1( hydromet2pdf_idx(ivar), iiPDF_eta, : ) + = corr_array_1(pdf_idx,iiPDF_eta,:) ! Correlation (in-precip) of eta and a hydrometeor in the 2nd PDF ! component. hydromet_pdf_params%corr_eta_hm_2(ivar) & - = corr_array_2( hydromet2pdf_idx(ivar), iiPDF_eta, : ) + = corr_array_2(pdf_idx,iiPDF_eta,:) ! Correlation (in-precip) of two hydrometeors, hmx and hmy, in the 1st ! PDF component. @@ -4027,7 +4096,7 @@ subroutine pack_hydromet_pdf_params( nz, hm_1, hm_2, pdf_dim, mu_x_1, & ! In do jvar = ivar+1, hydromet_dim, 1 hydromet_pdf_params%corr_hmx_hmy_1(jvar,ivar) & - = corr_array_1( hydromet2pdf_idx(jvar), hydromet2pdf_idx(ivar), : ) + = corr_array_1(hydromet2pdf_idx(jvar),pdf_idx,:) hydromet_pdf_params%corr_hmx_hmy_1(ivar,jvar) & = hydromet_pdf_params%corr_hmx_hmy_1(jvar,ivar) @@ -4041,7 +4110,7 @@ subroutine pack_hydromet_pdf_params( nz, hm_1, hm_2, pdf_dim, mu_x_1, & ! In do jvar = ivar+1, hydromet_dim, 1 hydromet_pdf_params%corr_hmx_hmy_2(jvar,ivar) & - = corr_array_2( hydromet2pdf_idx(jvar), hydromet2pdf_idx(ivar), : ) + = corr_array_2(hydromet2pdf_idx(jvar),pdf_idx,:) hydromet_pdf_params%corr_hmx_hmy_2(ivar,jvar) & = hydromet_pdf_params%corr_hmx_hmy_2(jvar,ivar) @@ -4073,8 +4142,12 @@ subroutine pack_hydromet_pdf_params( nz, hm_1, hm_2, pdf_dim, mu_x_1, & ! In end subroutine pack_hydromet_pdf_params !============================================================================= - elemental function compute_rtp2_from_chi( pdf_params, corr_chi_eta_1, & - corr_chi_eta_2 ) & + elemental function compute_rtp2_from_chi( sigma_chi_1, sigma_chi_2, & + sigma_eta_1, sigma_eta_2, & + rt_1, rt_2, & + crt_1, crt_2, & + mixt_frac, & + corr_chi_eta_1, corr_chi_eta_2 ) & result( rtp2_zt_from_chi ) ! Description: @@ -4097,14 +4170,19 @@ elemental function compute_rtp2_from_chi( pdf_params, corr_chi_eta_1, & one, & two - use pdf_parameter_module, only: & - pdf_parameter ! Type - implicit none ! Input Variables - type(pdf_parameter), intent(in) :: & - pdf_params + real( kind = core_rknd ), intent(in) :: & + sigma_chi_1, & ! Standard deviation of chi (1st PDF comp.) [kg/kg] + sigma_chi_2, & ! Standard deviation of chi (2nd PDF comp.) [kg/kg] + sigma_eta_1, & ! Standard deviation of eta (1st PDF comp.) [kg/kg] + sigma_eta_2, & ! Standard deviation of eta (2nd PDF comp.) [kg/kg] + crt_1, & ! Coef. of r_t in chi/eta eqns. (1st comp.) [-] + crt_2, & ! Coef. of r_t in chi/eta eqns. (2nd comp.) [-] + rt_1, & ! Mean of rt (1st PDF component) [kg/kg] + rt_2, & ! Mean of rt (2nd PDF component) [kg/kg] + mixt_frac ! Weight of 1st gaussian PDF component [-] real( kind = core_rknd ), intent(in) :: & corr_chi_eta_1, & ! Correlation of chi and eta in 1st PDF component [-] @@ -4119,34 +4197,15 @@ elemental function compute_rtp2_from_chi( pdf_params, corr_chi_eta_1, & varnce_rt_1_zt_from_chi, varnce_rt_2_zt_from_chi real( kind = core_rknd ) :: & - sigma_chi_1, & ! Standard deviation of chi (1st PDF comp.) [kg/kg] - sigma_chi_2, & ! Standard deviation of chi (2nd PDF comp.) [kg/kg] - sigma_eta_1, & ! Standard deviation of eta (1st PDF comp.) [kg/kg] - sigma_eta_2, & ! Standard deviation of eta (2nd PDF comp.) [kg/kg] - crt_1, & ! Coef. of r_t in chi/eta eqns. (1st comp.) [-] - crt_2, & ! Coef. of r_t in chi/eta eqns. (2nd comp.) [-] - rt_1, & ! Mean of rt (1st PDF component) [kg/kg] - rt_2, & ! Mean of rt (2nd PDF component) [kg/kg] - rtm, & ! Mean of rt (overall) [kg/kg] + rtm, & ! Mean of rt (overall) [kg/kg] sigma_rt_1_from_chi, & ! Standard deviation of rt (1st PDF comp.) [kg/kg] - sigma_rt_2_from_chi, & ! Standard deviation of rt (2nd PDF comp.) [kg/kg] - mixt_frac ! Weight of 1st gaussian PDF component [-] + sigma_rt_2_from_chi ! Standard deviation of rt (2nd PDF comp.) [kg/kg] + !----------------------------------------------------------------------- !----- Begin Code ----- - ! Enter some PDF parameters - sigma_chi_1 = pdf_params%stdev_chi_1 - sigma_chi_2 = pdf_params%stdev_chi_2 - sigma_eta_1 = pdf_params%stdev_eta_1 - sigma_eta_2 = pdf_params%stdev_eta_2 - rt_1 = pdf_params%rt_1 - rt_2 = pdf_params%rt_2 - crt_1 = pdf_params%crt_1 - crt_2 = pdf_params%crt_2 - mixt_frac = pdf_params%mixt_frac - varnce_rt_1_zt_from_chi & = ( corr_chi_eta_1 * sigma_chi_1 * sigma_eta_1 & + one_half * sigma_chi_1**2 + one_half * sigma_eta_1**2 ) & diff --git a/components/cam/src/physics/clubb/sigma_sqd_w_module.F90 b/components/cam/src/physics/clubb/sigma_sqd_w_module.F90 index 8335e5fd81ad..14effe7f0584 100644 --- a/components/cam/src/physics/clubb/sigma_sqd_w_module.F90 +++ b/components/cam/src/physics/clubb/sigma_sqd_w_module.F90 @@ -10,22 +10,51 @@ module sigma_sqd_w_module private ! Default scope contains -!--------------------------------------------------------------------------------------------------- - elemental function compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, wpthlp, wprtp ) & + + !============================================================================= + elemental function compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, & + up2, vp2, wpthlp, wprtp, upwp, vpwp ) & result( sigma_sqd_w ) -! Description: -! Compute the variable sigma_sqd_w (PDF width parameter) -! -! References: -! Eqn 22 in ``Equations for CLUBB'' -!--------------------------------------------------------------------------------------------------- + + ! Description: + ! Compute the variable sigma_sqd_w (PDF width parameter). + ! + ! The value of sigma_sqd_w is restricted in the ADG1 PDF in order to keep + ! the marginal PDFs of all responder variables (variables that do not set + ! the mixture fraction) valid. The limits on sigma_sqd_w in order to keep + ! the PDF of a responder variable, x, valid are: + ! + ! 0 <= sigma_sqd_w <= 1 - ^2 / ( * ). + ! + ! The overall limits on sigma_sqd_w must be applied based on the most + ! restrictive case so that all Double Gaussian PDF responder variables, x, + ! have realizable PDFs. The overall limits on sigma_sqd_w are: + ! + ! 0 <= sigma_sqd_w <= 1 - max( ^2 / ( * ), for all x). + ! + ! The equation used for sigma_sqd_w is: + ! + ! sigma_sqd_w = gamma_Skw_fnc + ! * ( 1 - max( ^2 / ( * ), for all x) ); + ! + ! where 0 <= gamma_Skw_fnc <= 1. + + ! References: + ! Eqn 22 in ``Equations for CLUBB'' + !----------------------------------------------------------------------- + use constants_clubb, only: & - w_tol, & ! Constant(s) - rt_tol, & - thl_tol + one, & ! Constant(s) + w_tol, & + rt_tol, & + thl_tol, & + w_tol_sqd + + use model_flags, only: & + l_predict_upwp_vpwp ! Variable(s) use clubb_precision, only: & - core_rknd ! Variable(s) + core_rknd ! Variable(s) implicit none @@ -34,33 +63,56 @@ elemental function compute_sigma_sqd_w( gamma_Skw_fnc, wp2, thlp2, rtp2, wpthlp, ! Input Variables real( kind = core_rknd ), intent(in) :: & - gamma_Skw_fnc, & ! Gamma as a function of skewness [-] - wp2, & ! Variance of vertical velocity [m^2/s^2] - thlp2, & ! Variance of liquid pot. temp. [K^2] - rtp2, & ! Variance of total water [kg^2/kg^2] - wpthlp, & ! Flux of liquid pot. temp. [m/s K] - wprtp ! Flux of total water [m/s kg/kg] + gamma_Skw_fnc, & ! Gamma as a function of skewness [-] + wp2, & ! Variance of vertical velocity [m^2/s^2] + thlp2, & ! Variance of liquid water potential temp. [K^2] + rtp2, & ! Variance of total water mixing ratio [kg^2/kg^2] + up2, & ! Variance of west-east horizontal velocity [m^2/s^2] + vp2, & ! Variance of south-north horizontal velocity [m^2/s^2] + wpthlp, & ! Flux of liquid water potential temp. [m/s K] + wprtp, & ! Flux of total water mixing ratio [m/s kg/kg] + upwp, & ! Flux of west-east horizontal velocity [m^2/s^2] + vpwp ! Flux of south-north horizontal velocity [m^2/s^2] ! Output Variable real( kind = core_rknd ) :: sigma_sqd_w ! PDF width parameter [-] + ! Local Variable + real( kind = core_rknd ) :: & + max_corr_w_x_sqd ! Max. val. of wpxp^2/(wp2*xp2) for all vars. x [-] + ! ---- Begin Code ---- !---------------------------------------------------------------- ! Compute sigma_sqd_w with new formula from Vince !---------------------------------------------------------------- - sigma_sqd_w = gamma_Skw_fnc * & - ( 1.0_core_rknd - min( & - max( ( wpthlp / ( sqrt( wp2 * thlp2 ) & - + 0.01_core_rknd * w_tol * thl_tol ) )**2, & - ( wprtp / ( sqrt( wp2 * rtp2 ) & - + 0.01_core_rknd * w_tol * rt_tol ) )**2 & - ), & ! max - 1.0_core_rknd ) & ! min - Known magic number (eq. 22 from "Equations for CLUBB") - ) + ! Find the maximum value of ^2 / ( * ) for all + ! variables x that are Double Gaussian PDF responder variables. This + ! includes rt and theta-l. When l_predict_upwp_vpwp is enabled, u and v are + ! also calculated as part of the PDF, and they are included as well. + ! Additionally, when sclr_dim > 0, passive scalars (sclr) are also included. + max_corr_w_x_sqd = max( ( wpthlp / ( sqrt( wp2 * thlp2 ) & + + 0.01_core_rknd * w_tol * thl_tol ) )**2, & + ( wprtp / ( sqrt( wp2 * rtp2 ) & + + 0.01_core_rknd * w_tol * rt_tol ) )**2 ) + + if ( l_predict_upwp_vpwp ) then + max_corr_w_x_sqd = max( max_corr_w_x_sqd, & + ( upwp / ( sqrt( up2 * wp2 ) & + + 0.01_core_rknd * w_tol_sqd ) )**2, & + ( vpwp / ( sqrt( vp2 * wp2 ) & + + 0.01_core_rknd * w_tol_sqd ) )**2 ) + endif ! l_predict_upwp_vpwp + + ! Calculate the value of sigma_sqd_w . + sigma_sqd_w = gamma_Skw_fnc * ( one - min( max_corr_w_x_sqd, one ) ) + return + end function compute_sigma_sqd_w + !============================================================================= + end module sigma_sqd_w_module diff --git a/components/cam/src/physics/clubb/stats_clubb_utilities.F90 b/components/cam/src/physics/clubb/stats_clubb_utilities.F90 index fec5200a8ea6..a72382ab1774 100644 --- a/components/cam/src/physics/clubb/stats_clubb_utilities.F90 +++ b/components/cam/src/physics/clubb/stats_clubb_utilities.F90 @@ -1698,6 +1698,7 @@ subroutine stats_end_timestep( ) end if call write_grads( stats_sfc%file ) else ! l_netcdf + #ifdef NETCDF call write_netcdf( stats_zt%file ) call write_netcdf( stats_zm%file ) @@ -2123,7 +2124,7 @@ subroutine stats_accumulate & real( kind = core_rknd ), intent(in), dimension(gr%nz) :: & sigma_sqd_w ! PDF width parameter (momentum levels) [-] - type(pdf_parameter), dimension(gr%nz), intent(in) :: & + type(pdf_parameter), intent(in) :: & pdf_params, & ! PDF parameters (thermodynamic levels) [units vary] pdf_params_zm ! PDF parameters on momentum levels [units vary] diff --git a/components/cam/src/physics/clubb/stats_variables.F90 b/components/cam/src/physics/clubb/stats_variables.F90 index f4e9aac02c2b..67ef19efcbd0 100755 --- a/components/cam/src/physics/clubb/stats_variables.F90 +++ b/components/cam/src/physics/clubb/stats_variables.F90 @@ -118,6 +118,11 @@ module stats_variables !$omp itau_zt, iKh_zt, iwp2thvp, iwp2rcp, iwprtpthlp, irc_coef, & !$omp isigma_sqd_w_zt, irho ) + integer, public :: & + itau_zm_simp = 0 +!$omp threadprivate( itau_zm_simp ) + + integer, dimension(:), allocatable, public :: & ihm_1, & ihm_2 @@ -214,7 +219,7 @@ module stats_variables isilhs_variance_category, & ilh_samp_frac_category -!$omp threadprivate( isilhs_variance_category ) +!$omp threadprivate( isilhs_variance_category, ilh_samp_frac_category ) integer, public :: & iNcm = 0, & ! Brian @@ -435,7 +440,7 @@ module stats_variables !$omp threadprivate(irrm_bt, irrm_ma, irrm_ta, irrm_sd) !$omp threadprivate(irrm_ts, irrm_sd_morr) !$omp threadprivate(irrm_cond, irrm_auto, irrm_accr) -!$omp threadprivate(irrm_cond_adj, irrm_src_adj ) +!$omp threadprivate(irrm_cond_adj, irrm_src_adj, irrm_mc_nonadj) !$omp threadprivate(irrm_mc, irrm_hf, irrm_wvhf, irrm_cl) integer, public :: & @@ -995,6 +1000,10 @@ module stats_variables iuprtp = 0, & ivpthlp = 0, & ivprtp = 0, & + iupthvp = 0, & + iuprcp = 0, & + ivpthvp = 0, & + ivprcp = 0, & iSkw_zm = 0, & iSkthl_zm = 0, & iSkrt_zm = 0 @@ -1028,6 +1037,7 @@ module stats_variables !$omp threadprivate(iwprcp, irc_coef_zm, ithlprcp, irtprcp, ircp2) !$omp threadprivate(iupwp, ivpwp) !$omp threadprivate(iupthlp, iuprtp, ivpthlp, ivprtp) +!$omp threadprivate(iupthvp, iuprcp, ivpthvp, ivprcp) !$omp threadprivate(iSkw_zm, iSkthl_zm, iSkrt_zm) !$omp threadprivate(irho_zm, isigma_sqd_w, irho_ds_zm, ithv_ds_zm, iem, ishear) !$omp threadprivate(imean_w_up, imean_w_down) @@ -1051,6 +1061,7 @@ module stats_variables !$omp threadprivate(igamma_Skw_fnc, iC6rt_Skw_fnc, iC6thl_Skw_fnc) !$omp threadprivate(iC7_Skw_fnc, iC1_Skw_fnc) +!$omp threadprivate(ibrunt_vaisala_freq_sqd, iRichardson_num, ishear_sqd) integer, public :: & icoef_wp4_implicit = 0 diff --git a/components/cam/src/physics/clubb/stats_zm_module.F90 b/components/cam/src/physics/clubb/stats_zm_module.F90 index a45a505d4b61..7ea84a203189 100644 --- a/components/cam/src/physics/clubb/stats_zm_module.F90 +++ b/components/cam/src/physics/clubb/stats_zm_module.F90 @@ -46,6 +46,7 @@ subroutine stats_init_zm( vars_zm, l_error ) irtpthvp, & ithlpthvp, & itau_zm, & + itau_zm_simp, & iKh_zm, & iK_hm, & iwprcp, & @@ -64,6 +65,10 @@ subroutine stats_init_zm( vars_zm, l_error ) iuprtp, & ivpthlp, & ivprtp, & + iupthvp, & + iuprcp, & + ivpthvp, & + ivprcp, & irho_zm, & isigma_sqd_w, & irho_ds_zm, & @@ -645,6 +650,15 @@ subroutine stats_init_zm( vars_zm, l_error ) l_silhs=.false., grid_kind=stats_zm ) k = k + 1 + case ('tau_zm_simp') + itau_zm_simp = k + + call stat_assign( var_index=itau_zm_simp, var_name="tau_zm_simp", & + var_description="simple tau on momentum levels [s]", var_units="s", & + l_silhs=.false., grid_kind=stats_zm ) + k = k + 1 + + case ('Kh_zm') iKh_zm = k @@ -746,6 +760,30 @@ subroutine stats_init_zm( vars_zm, l_error ) var_description="v'rt', Northward total water flux [(m/s)(kg/kg)]", & var_units="(m/s)(kg/kg)", l_silhs=.false., grid_kind=stats_zm ) k = k + 1 + case ('upthvp') + iupthvp = k + call stat_assign( var_index=iupthvp, var_name="upthvp", & + var_description="u'thv', Eastward theta_v flux [(m/s)K]", & + var_units="(m/s)K", l_silhs=.false., grid_kind=stats_zm ) + k = k + 1 + case ('uprcp') + iuprcp = k + call stat_assign( var_index=iuprcp, var_name="uprcp", & + var_description="u'rc', Eastward liquid water flux [(m/s)(kg/kg)]", & + var_units="(m/s)(kg/kg)", l_silhs=.false., grid_kind=stats_zm ) + k = k + 1 + case ('vpthvp') + ivpthvp = k + call stat_assign( var_index=ivpthvp, var_name="vpthvp", & + var_description="v'thv', Northward theta_v flux [(m/s)K]", & + var_units="(m/s)K", l_silhs=.false., grid_kind=stats_zm ) + k = k + 1 + case ('vprcp') + ivprcp = k + call stat_assign( var_index=ivprcp, var_name="vprcp", & + var_description="v'rc', Northward liquid water flux [(m/s)(kg/kg)]", & + var_units="(m/s)(kg/kg)", l_silhs=.false., grid_kind=stats_zm ) + k = k + 1 case ('rho_zm') irho_zm = k call stat_assign( var_index=irho_zm, var_name="rho_zm", & diff --git a/components/cam/src/physics/clubb/variables_diagnostic_module.F90 b/components/cam/src/physics/clubb/variables_diagnostic_module.F90 index 60e77b60545f..dbf92c6d72eb 100644 --- a/components/cam/src/physics/clubb/variables_diagnostic_module.F90 +++ b/components/cam/src/physics/clubb/variables_diagnostic_module.F90 @@ -55,7 +55,7 @@ module variables_diagnostic_module !$omp threadprivate(rsat) - type(pdf_parameter), allocatable, dimension(:), target, public :: & + type(pdf_parameter), allocatable, public, save :: & pdf_params_zm, & ! pdf_params on momentum levels [units vary] pdf_params_zm_frz !used when l_use_ice_latent = .true. @@ -282,8 +282,10 @@ subroutine setup_diagnostic_variables( nz ) allocate( radht(1:nz) ) ! SW + LW heating rate ! pdf_params on momentum levels - allocate( pdf_params_zm(1:nz) ) - allocate( pdf_params_zm_frz(1:nz) ) + allocate( pdf_params_zm ) + allocate( pdf_params_zm_frz ) + call init_pdf_params( nz, pdf_params_zm ) + call init_pdf_params( nz, pdf_params_zm_frz ) ! Second order moments @@ -411,11 +413,6 @@ subroutine setup_diagnostic_variables( nz ) Frad_SW_down = 0.0_core_rknd Frad_LW_down = 0.0_core_rknd - - ! pdf_params on momentum levels - call init_pdf_params( nz, pdf_params_zm ) - call init_pdf_params( nz, pdf_params_zm_frz ) - ! Second order moments thlprcp = 0.0_core_rknd rtprcp = 0.0_core_rknd diff --git a/components/cam/src/physics/clubb/variables_prognostic_module.F90 b/components/cam/src/physics/clubb/variables_prognostic_module.F90 index 0e55c7e63e6c..8cb946d04627 100644 --- a/components/cam/src/physics/clubb/variables_prognostic_module.F90 +++ b/components/cam/src/physics/clubb/variables_prognostic_module.F90 @@ -169,7 +169,7 @@ module variables_prognostic_module !$omp threadprivate(sigma_sqd_w) - type(pdf_parameter), target, allocatable, dimension(:), public :: & + type(pdf_parameter), allocatable, public, save :: & pdf_params, & pdf_params_frz !for use when l_use_ice_latent = .true. @@ -298,8 +298,10 @@ subroutine setup_prognostic_variables( nz ) allocate( sigma_sqd_w(1:nz) ) ! PDF width parameter (momentum levels) ! Variables for pdf closure scheme - allocate( pdf_params(1:nz) ) - allocate( pdf_params_frz(1:nz) ) + allocate( pdf_params ) + allocate( pdf_params_frz ) + call init_pdf_params( nz, pdf_params ) + call init_pdf_params( nz, pdf_params_frz ) !--------- Set initial values for array variables --------- @@ -361,10 +363,6 @@ subroutine setup_prognostic_variables( nz ) sigma_sqd_w = 0.0_core_rknd ! PDF width parameter (momentum levels) - ! Variables for PDF closure scheme - call init_pdf_params( nz, pdf_params ) - call init_pdf_params( nz, pdf_params_frz ) - ! Surface fluxes wpthlp_sfc = 0.0_core_rknd wprtp_sfc = 0.0_core_rknd From ee40d7d3fe7f8ca99a481f229808955ca4e19d93 Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Wed, 10 Jul 2019 15:43:13 -0500 Subject: [PATCH 5/7] Modifying parameters_tunable to update it but still include Balwinders changes. larson-group/E3SM#7 --- .../src/physics/clubb/parameters_tunable.F90 | 69 ++++++++++++++++--- 1 file changed, 58 insertions(+), 11 deletions(-) diff --git a/components/cam/src/physics/clubb/parameters_tunable.F90 b/components/cam/src/physics/clubb/parameters_tunable.F90 index 642d2583e996..0a1cee709ca5 100644 --- a/components/cam/src/physics/clubb/parameters_tunable.F90 +++ b/components/cam/src/physics/clubb/parameters_tunable.F90 @@ -212,6 +212,14 @@ module parameters_tunable !$omp threadprivate(Skw_max_mag) + real( kind = core_rknd ), public :: & + C_invrs_tau_bkgnd = 1.0_core_rknd ,& ! + C_invrs_tau_sfc = 0.1_core_rknd ,& ! + C_invrs_tau_shear = 0.02_core_rknd,& ! + C_invrs_tau_N2 = 0.1_core_rknd ! +!$omp threadprivate(C_invrs_tau_bkgnd,C_invrs_tau_sfc) +!$omp threadprivate(C_invrs_tau_shear,C_invrs_tau_N2) + ! Parameters for the new PDF (w, rt, and theta-l). ! ! Brian Griffin added a tunable parameter for the PDF of w, @@ -350,7 +358,8 @@ module parameters_tunable omicron, zeta_vrnce_rat, upsilon_precip_frac_rat, & lambda0_stability_coef, mult_coef, taumin, taumax, mu, Lscale_mu_coef, & Lscale_pert_coef, alpha_corr, Skw_denom_coef, c_K10, c_K10h, & - thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag + thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, C_invrs_tau_shear, C_invrs_tau_N2 ! These are referenced together often enough that it made sense to ! make a list of them. Note that lmin_coef is the input parameter, @@ -403,7 +412,9 @@ module parameters_tunable "Skw_denom_coef ", "c_K10 ", & "c_K10h ", "thlp2_rad_coef ", & "thlp2_rad_cloud_frac_thresh ", "up2_vp2_factor ", & - "Skw_max_mag " /) + "Skw_max_mag ", "C_invrs_tau_bkgnd ", & + "C_invrs_tau_sfc ", "C_invrs_tau_shear ", & + "C_invrs_tau_N2 "/) real( kind = core_rknd ), parameter, private :: & init_value = -999._core_rknd ! Initial value for the parameters, used to detect missing values @@ -650,7 +661,9 @@ subroutine setup_parameters & lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, & + C_invrs_tau_shear, C_invrs_tau_N2) ! It was decided after some experimentation, that the best @@ -1166,7 +1179,9 @@ subroutine read_parameters( iunit, filename, params ) lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, params ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, & + C_invrs_tau_shear, C_invrs_tau_N2, params ) l_error = .false. @@ -1243,7 +1258,9 @@ subroutine read_param_max & omicron, zeta_vrnce_rat, upsilon_precip_frac_rat, & lambda0_stability_coef, mult_coef, taumin, taumax, mu, Lscale_mu_coef, & Lscale_pert_coef, alpha_corr, Skw_denom_coef, c_K10, c_K10h, & - thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag + thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc,& + C_invrs_tau_shear, C_invrs_tau_N2 ! Initialize values to -999. call init_parameters_999( ) @@ -1271,7 +1288,9 @@ subroutine read_param_max & lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, param_max ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc,& + C_invrs_tau_shear, C_invrs_tau_N2, param_max ) l_error = .false. @@ -1319,7 +1338,9 @@ subroutine pack_parameters & lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, params ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag,& + C_invrs_tau_bkgnd, C_invrs_tau_sfc, & + C_invrs_tau_shear, C_invrs_tau_N2, params ) ! Description: ! Takes the list of scalar variables and puts them into a 1D vector. @@ -1413,6 +1434,10 @@ subroutine pack_parameters & ithlp2_rad_cloud_frac_thresh, & iup2_vp2_factor, & iSkw_max_mag, & + iC_invrs_tau_bkgnd, & + iC_invrs_tau_sfc, & + iC_invrs_tau_shear, & + iC_invrs_tau_N2, & nparams implicit none @@ -1432,7 +1457,8 @@ subroutine pack_parameters & omicron, zeta_vrnce_rat, upsilon_precip_frac_rat, & lambda0_stability_coef, mult_coef, taumin, taumax, Lscale_mu_coef, & Lscale_pert_coef, alpha_corr, Skw_denom_coef, c_K10, c_K10h, & - thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag + thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, C_invrs_tau_shear, C_invrs_tau_N2 ! Output variables real( kind = core_rknd ), intent(out), dimension(nparams) :: params @@ -1527,6 +1553,10 @@ subroutine pack_parameters & params(ithlp2_rad_cloud_frac_thresh) = thlp2_rad_cloud_frac_thresh params(iup2_vp2_factor) = up2_vp2_factor params(iSkw_max_mag) = Skw_max_mag + params(iC_invrs_tau_bkgnd) = C_invrs_tau_bkgnd + params(iC_invrs_tau_sfc) = C_invrs_tau_sfc + params(iC_invrs_tau_shear) = C_invrs_tau_shear + params(iC_invrs_tau_N2) = C_invrs_tau_N2 return end subroutine pack_parameters @@ -1548,7 +1578,9 @@ subroutine unpack_parameters & lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, & + C_invrs_tau_shear, C_invrs_tau_N2 ) ! Description: ! Takes the 1D vector and returns the list of scalar variables. @@ -1642,6 +1674,10 @@ subroutine unpack_parameters & ithlp2_rad_cloud_frac_thresh, & iup2_vp2_factor, & iSkw_max_mag, & + iC_invrs_tau_bkgnd,& + iC_invrs_tau_sfc,& + iC_invrs_tau_shear,& + iC_invrs_tau_N2,& nparams implicit none @@ -1664,7 +1700,8 @@ subroutine unpack_parameters & omicron, zeta_vrnce_rat, upsilon_precip_frac_rat, & lambda0_stability_coef, mult_coef, taumin, taumax, Lscale_mu_coef, & Lscale_pert_coef, alpha_corr, Skw_denom_coef, c_K10, c_K10h, & - thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag + thlp2_rad_coef, thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, C_invrs_tau_shear, C_invrs_tau_N2 C1 = params(iC1) C1b = params(iC1b) @@ -1756,6 +1793,10 @@ subroutine unpack_parameters & thlp2_rad_cloud_frac_thresh = params(ithlp2_rad_cloud_frac_thresh) up2_vp2_factor = params(iup2_vp2_factor) Skw_max_mag = params(iSkw_max_mag) + C_invrs_tau_bkgnd = params(iC_invrs_tau_bkgnd) + C_invrs_tau_sfc = params(iC_invrs_tau_sfc ) + C_invrs_tau_shear = params(iC_invrs_tau_shear) + C_invrs_tau_N2 = params(iC_invrs_tau_N2) return end subroutine unpack_parameters @@ -1790,7 +1831,9 @@ subroutine get_parameters( params ) lambda0_stability_coef, mult_coef, taumin, taumax, & Lscale_mu_coef, Lscale_pert_coef, alpha_corr, & Skw_denom_coef, c_K10, c_K10h, thlp2_rad_coef, & - thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, params ) + thlp2_rad_cloud_frac_thresh, up2_vp2_factor, Skw_max_mag, & + C_invrs_tau_bkgnd, C_invrs_tau_sfc, & + C_invrs_tau_shear, C_invrs_tau_N2, params ) return @@ -1888,6 +1931,10 @@ subroutine init_parameters_999( ) thlp2_rad_cloud_frac_thresh = init_value up2_vp2_factor = init_value Skw_max_mag = init_value + C_invrs_tau_bkgnd = init_value + C_invrs_tau_sfc = init_value + C_invrs_tau_shear = init_value + C_invrs_tau_N2 = init_value return From 25ca78c0b054437362641aaa4d03b8164fe827e9 Mon Sep 17 00:00:00 2001 From: zhunguo <48636672+zhunguo@users.noreply.github.com> Date: Wed, 10 Jul 2019 16:13:44 -0500 Subject: [PATCH 6/7] Update clubb_intr.F90 --- components/cam/src/physics/cam/clubb_intr.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index cb3b3f1799a3..39987ad0948b 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -1746,7 +1746,7 @@ subroutine clubb_tend_cam( & ! Compute mean w wind on thermo grid, convert from omega to w wm_zt(1) = 0._r8 do k=1,pver - wm_zt(k+1) = -1._r8*state1%omega(i,pver-k+1)*invrs_rho_ds_zt(i,k+1)*invrs_gravit + wm_zt(k+1) = -1._r8*state1%omega(i,pver-k+1)*invrs_rho_ds_zt(k+1)*invrs_gravit enddo ! ------------------------------------------------- ! From 869072581cb175ccae72d2ec1c7928868dcfa386 Mon Sep 17 00:00:00 2001 From: Gunther Huebler Date: Wed, 10 Jul 2019 18:15:58 -0500 Subject: [PATCH 7/7] Fixing the init_pdf_params_api calls, the indices were incorrect. Also renaming them for clarity. larson-group/E3SM#7 --- components/cam/src/physics/cam/clubb_intr.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/components/cam/src/physics/cam/clubb_intr.F90 b/components/cam/src/physics/cam/clubb_intr.F90 index cb3b3f1799a3..d348ca2554cd 100644 --- a/components/cam/src/physics/cam/clubb_intr.F90 +++ b/components/cam/src/physics/cam/clubb_intr.F90 @@ -564,7 +564,7 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) character(len=128) :: errstring ! error status for CLUBB init integer :: err_code, iunit ! Code for when CLUBB fails - integer :: i, j, k, l ! Indices + integer :: i, j, k, l, idx_chunk, idx_pcols ! Indices integer :: read_status ! Length of a string integer :: ntop_eddy ! Top interface level to which eddy vertical diffusion is applied ( = 1 ) integer :: nbot_eddy ! Bottom interface level to which eddy vertical diffusion is applied ( = pver ) @@ -583,10 +583,10 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in) pdf_params_chnk(pcols,begchunk:endchunk), & pdf_params_zm_chnk(pcols,begchunk:endchunk) ) - do i = begchunk, endchunk - do j = 1, pcols - call init_pdf_params_api( pverp, pdf_params_chnk(i,j) ) - call init_pdf_params_api( pverp, pdf_params_zm_chnk(i,j) ) + do idx_chunk = begchunk, endchunk + do idx_pcols = 1, pcols + call init_pdf_params_api( pverp, pdf_params_chnk(idx_pcols,idx_chunk) ) + call init_pdf_params_api( pverp, pdf_params_zm_chnk(idx_pcols,idx_chunk) ) end do end do