diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index a6c2d8f771..a24ba9b7db 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -1551,7 +1551,7 @@ state integer STEPCU - misc 1 - r "S state real RTHRATEN ikj misc 1 - rd "RTHRATEN" "THETA TENDENCY DUE TO RADIATION" "K s-1" state real RTHRATENLW ikj misc 1 - r "RTHRATLW" "UNCOUPLED THETA TENDENCY DUE TO LONG WAVE RADIATION" "K s-1" state real RTHRATENSW ikj misc 1 - r "RTHRATSW" "UNCOUPLED THETA TENDENCY DUE TO SHORT WAVE RADIATION" "K s-1" -state real CLDFRA ikj misc 1 - rh "CLDFRA" "CLOUD FRACTION" "" +state real CLDFRA ikj misc 1 - i0rh "CLDFRA" "CLOUD FRACTION" "" state real CONVCLD ij misc 1 - r "CONVCLD" "BMJ CONVECTIVE CLOUD" "kg m-2" state real CCLDFRA ikj misc 1 - r "CCLDFRA" "CONVECTIVE CLOUD FRACTION" "" state real CLDFRA_OLD ikj misc 1 - r "CLDFRA_OLD" "previous time level cldfra" "" @@ -2458,6 +2458,7 @@ rconfig integer tracer_pblmix namelist,physics max_domains 1 rconfig logical use_aero_icbc namelist,physics 1 .false. rh "use_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero" "logical flag" rconfig logical use_rap_aero_icbc namelist,physics 1 .false. r "use_rap_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero from RAP" "logical flag" rconfig integer use_mp_re namelist,physics 1 1 h "use_mp_re" "use effective radii computed in some mp schemes in RRTMG" "flag" +rconfig logical insert_init_cloud namelist,physics 1 .false. irh "insert_init_cloud" "Insert initial cloud using icloud=3 method (Thompson)" "logical flag" # The following two options are hooked into various microphysics schemes to allow for ensemble perturbations of CCN and hail/graupel PSDs - GAC (AFWA) rconfig real ccn_conc namelist,physics 1 1.0E8 h "ccn_conc" "CCN concentration" "# m-3" diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index 00a559912d..09eec7232b 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -69,6 +69,8 @@ SUBROUTINE init_domain_rk ( grid & ) USE module_optional_input + USE module_radiation_driver + USE module_dm, ONLY : wrf_dm_max_real IMPLICIT NONE ! Input space and data. No gridded meteorological data has been stored, though. @@ -121,6 +123,13 @@ SUBROUTINE init_domain_rk ( grid & INTEGER :: auto_levels_opt REAL :: max_dz, dzbot, dzstretch_s, dzstretch_u, z1, airmass + INTEGER:: i_end, j_end + REAL, ALLOCATABLE, DIMENSION(:):: temp_P, temp_T, temp_R, temp_Dz + REAL, ALLOCATABLE, DIMENSION(:):: temp_Qv, temp_Qc, temp_Qi, temp_Qs + REAL, ALLOCATABLE, DIMENSION(:):: temp_CF, temp_Nc, temp_Ni + REAL:: max_relh, temp_xland, gridkm + LOGICAL :: debug_flag = .FALSE. + REAL :: dclat ! INTEGER , PARAMETER :: nl_max = 1000 @@ -4370,6 +4379,125 @@ SUBROUTINE init_domain_rk ( grid & enddo ENDIF +!+---+-----------------------------------------------------------------+ +!..We can consider that in circumstance of a 'cold start' we can make +!.. an attempt to insert some initial clouds to get a better starting +!.. radiation representation due to clouds using the icloud3 cloud fraction +!.. scheme. +!+---+-----------------------------------------------------------------+ + + if (config_flags%insert_init_cloud .AND. & + (P_QC .gt. 1 .AND. P_QI .gt. 1)) then + + ALLOCATE(temp_P(kts:kte-1)) + ALLOCATE(temp_Dz(kts:kte-1)) + ALLOCATE(temp_T(kts:kte-1)) + ALLOCATE(temp_R(kts:kte-1)) + ALLOCATE(temp_Qv(kts:kte-1)) + ALLOCATE(temp_Qc(kts:kte-1)) + ALLOCATE(temp_Nc(kts:kte-1)) + ALLOCATE(temp_Qi(kts:kte-1)) + ALLOCATE(temp_Ni(kts:kte-1)) + ALLOCATE(temp_Qs(kts:kte-1)) + ALLOCATE(temp_CF(kts:kte-1)) + + i_end = MIN(ite,ide-1) + j_end = MIN(jte,jde-1) + max_relh = 1.5 +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + max_relh = wrf_dm_max_real ( MAXVAL(grid%u_1(its:i_end,:,jts:j_end)) ) +#else + max_relh = MAXVAL ( grid%u_1(its:i_end,:,jts:j_end) ) +#endif + max_relh = max_relh*0.01 + + gridkm = SQRT(config_flags%dx*config_flags%dy)*0.001 + + !..As it occurs up above, temporarily utilizing the v_1 variable, + !.. to hold temperature, which it does when time_loop=0. + + IF ( internal_time_loop .GT. 1 ) THEN + grid%v_1 = grid%t_2+t0 + CALL theta_to_t ( grid%v_1 , grid%p_hyd , p00 , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + ENDIF + + do j = jts, j_end + do i = its, i_end + IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE + debug_flag = .false. +! if (i.eq.9 .and. j.eq.9) debug_flag = .true. + + temp_xland = grid%xland(i,j) + if (grid%lakemask(i,j) .eq. 1) temp_xland = 1 + do k = kts, kte-1 + temp_Dz(k) = (grid%ph_2(i,k+1,j)+grid%phb(i,k+1,j) - (grid%ph_2(i,k,j)+grid%phb(i,k,j)))/g + temp_P(k) = grid%p_hyd(i,k,j) + temp_T(k) = grid%v_1(i,k,j) ! Around line num 1800 v_1 used to hold temperature. + temp_R(k) = 1./grid%alt(i,k,j) + temp_Qv(k) = moist(i,k,j,P_QV) + temp_Qc(k) = MAX(0., moist(i,k,j,P_QC)) + temp_Qi(k) = MAX(0., moist(i,k,j,P_QI)) + if (P_QS .gt. 1) then + temp_Qs(k) = MAX(0., moist(i,k,j,P_QS)) + else + temp_Qs(k) = 0. + endif + if (P_QNI .gt. 1) then + temp_Ni(k) = MAX(0., scalar(i,k,j,P_QNI)) + else + temp_Ni(k) = 0. + endif + if (P_QNC .gt. 1) then + temp_Nc(k) = MAX(0., scalar(i,k,j,P_QNC)) + else + temp_Nc(k) = 0. + endif + temp_CF(k) = 0. + enddo + + call cal_cldfra3(temp_CF,temp_Qv,temp_Qc,temp_Qi,temp_Qs, & + & temp_Dz, temp_P, temp_T, temp_xland, gridkm, & + & config_flags%insert_init_cloud, max_relh, & + & kts, kte-1, debug_flag) + + do k = kts, kte-1 + grid%cldfra(i,k,j) = temp_CF(k) + enddo + + if (debug_flag) then + do k = kts, kte-1 + write(*,*) ' DEBUG_column: ', temp_P(k), temp_T(k), temp_Qv(k), temp_Qc(k), temp_Qi(k), moist(i,k,j,P_QC), moist(i,k,j,P_QI) + enddo + endif + + if (config_flags%insert_init_cloud) then + do k = kts, kte-1 + moist(i,k,j,P_QV) = MAX(temp_Qv(k), moist(i,k,j,P_QV)) + moist(i,k,j,P_QC) = temp_Qc(k) + moist(i,k,j,P_QI) = temp_Qi(k) + enddo + endif + + enddo + enddo + + DEALLOCATE(temp_P) + DEALLOCATE(temp_Dz) + DEALLOCATE(temp_T) + DEALLOCATE(temp_R) + DEALLOCATE(temp_Qv) + DEALLOCATE(temp_Qc) + DEALLOCATE(temp_Nc) + DEALLOCATE(temp_Qi) + DEALLOCATE(temp_Ni) + DEALLOCATE(temp_Qs) + DEALLOCATE(temp_CF) + + endif + !+---+-----------------------------------------------------------------+ !..Let us ensure that double-moment microphysics variables have numbers !.. where there is mass. Currently doing this for Thompson-MP only, but diff --git a/phys/module_diag_nwp.F b/phys/module_diag_nwp.F index f502ef5515..811af567a7 100644 --- a/phys/module_diag_nwp.F +++ b/phys/module_diag_nwp.F @@ -606,8 +606,8 @@ SUBROUTINE diagnostic_output_nwp( & DO i=i_start(ij),i_end(ij) DO k=kme-1, kms, -1 if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE - zans1 = (3.5 + 2./7. * (ALOG10(temp_qg(i,k,j))+7.)) - zans1 = MAX(2., MIN(zans1, 7.)) + zans1 = 3.0 + 2./7.*(ALOG10(temp_qg(i,k,j))+8.) + zans1 = MAX(2., MIN(zans1, 6.)) N0exp = 10.**zans1 lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g) diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index fbc32a7420..b7bbcec4b3 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -107,8 +107,8 @@ MODULE module_mp_thompson !.. mixing ratio. Also, when mu_g is non-zero, these become equiv !.. y-intercept for an exponential distrib and proper values are !.. computed based on same mixing ratio and total number concentration. - REAL, PARAMETER, PRIVATE:: gonv_min = 1.E3 - REAL, PARAMETER, PRIVATE:: gonv_max = 9.E6 + REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2 + REAL, PARAMETER, PRIVATE:: gonv_max = 1.E7 !..Mass power law relations: mass = am*D**bm !.. Snow from Field et al. (2005), others assume spherical form. @@ -288,11 +288,11 @@ MODULE module_mp_thompson !..Lookup tables for graupel y-intercept parameter (/m**4). REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: & - N0g_exp = (/1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & - 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & - 1.e7/) + 1.e6/) !..Lookup tables for ice number concentration (/m**3). REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: & @@ -1577,7 +1577,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & enddo #endif -!..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments. +!..Bug fix (2016Jun15 and 2020Nov30), prevent use of uninitialized value(s). do k = kts, kte smo0(k) = 0. smo1(k) = 0. @@ -1587,6 +1587,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & smod(k) = 0. smoe(k) = 0. smof(k) = 0. + mvd_r(k) = 0. + mvd_c(k) = 0. enddo !+---+-----------------------------------------------------------------+ @@ -1848,8 +1850,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = (3.5 + 2./7. * (ygra1+7.)) - zans1 = MAX(2., MIN(zans1, 7.)) + zans1 = 3.0 + 2./7.*(ygra1+8.) + zans1 = MAX(2., MIN(zans1, 6.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg @@ -2909,8 +2911,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = (3.5 + 2./7. * (ygra1+7.)) - zans1 = MAX(2., MIN(zans1, 7.)) + zans1 = 3.0 + 2./7.*(ygra1+8.) + zans1 = MAX(2., MIN(zans1, 6.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg @@ -3578,10 +3580,10 @@ subroutine qr_acr_qg good = 0 IF ( wrf_dm_on_monitor() ) THEN - INQUIRE(FILE="qr_acr_qgV2.dat",EXIST=lexist) + INQUIRE(FILE="qr_acr_qgV3.dat",EXIST=lexist) IF ( lexist ) THEN - CALL wrf_message("ThompMP: read qr_acr_qgV2.dat instead of computing") - OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=1234) + CALL wrf_message("ThompMP: read qr_acr_qgV3.dat instead of computing") + OPEN(63,file="qr_acr_qgV3.dat",form="unformatted",err=1234) READ(63,err=1234) tcg_racg READ(63,err=1234) tmr_racg READ(63,err=1234) tcr_gacr @@ -3594,12 +3596,12 @@ subroutine qr_acr_qg INQUIRE(63,opened=lopen) IF (lopen) THEN IF( force_read_thompson ) THEN - CALL wrf_error_fatal("Error reading qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.") + CALL wrf_error_fatal("Error reading qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.") ENDIF CLOSE(63) ELSE IF( force_read_thompson ) THEN - CALL wrf_error_fatal("Error opening qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.") + CALL wrf_error_fatal("Error opening qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.") ENDIF ENDIF ELSE @@ -3610,7 +3612,7 @@ subroutine qr_acr_qg ENDIF ELSE IF( force_read_thompson ) THEN - CALL wrf_error_fatal("Non-existent qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.") + CALL wrf_error_fatal("Non-existent qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.") ENDIF ENDIF ENDIF @@ -3722,8 +3724,8 @@ subroutine qr_acr_qg #endif IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN - CALL wrf_message("Writing qr_acr_qgV2.dat in Thompson MP init") - OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=9234) + CALL wrf_message("Writing qr_acr_qgV3.dat in Thompson MP init") + OPEN(63,file="qr_acr_qgV3.dat",form="unformatted",err=9234) WRITE(63,err=9234) tcg_racg WRITE(63,err=9234) tmr_racg WRITE(63,err=9234) tcr_gacr @@ -3733,7 +3735,7 @@ subroutine qr_acr_qg CLOSE(63) RETURN ! ----- RETURN 9234 CONTINUE - CALL wrf_error_fatal("Error writing qr_acr_qgV2.dat") + CALL wrf_error_fatal("Error writing qr_acr_qgV3.dat") ENDIF ENDIF @@ -5294,8 +5296,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & if (ANY(L_qg .eqv. .true.)) then do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = (3.5 + 2./7. * (ygra1+7.)) - zans1 = MAX(2., MIN(zans1, 7.)) + zans1 = 3.0 + 2./7.*(ygra1+8.) + zans1 = MAX(2., MIN(zans1, 6.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index b8a6b06820..93243461a9 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -3765,6 +3765,17 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & qvs(k) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+35.) endif + if (modify_qvapor) then + if (qc(k).gt.1.E-8) then + qv(k) = MAX(qv(k), qvsw) + qvs(k) = qvsw + endif + if (qc(k).le.1.E-8 .and. qi(k).ge.1.E-9) then + qv(k) = MAX(qv(k), qvsi*1.005) !..To ensure a tiny bit ice supersaturation + qvs(k) = qvsi + endif + endif + rh(k) = MAX(0.01, qv(k)/qvs(k)) rhoa(k) = p(k)/(287.0*t(k)) ENDDO @@ -3779,14 +3790,13 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.65 + SQRT(1./(25.0+gridkm*gridkm*delz*0.01)) - RH_00O = 0.81 + SQRT(1./(50.0+gridkm*gridkm*delz*0.01)) + RH_00L = 0.53 + MIN(0.46,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00O = 0.86 + MIN(0.13,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).gt.1.E-7 .or. qi(k).ge.1.E-7 & - & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then + if (qc(k).gt.1.E-6 .or. qi(k).ge.1.E-6 & + & .or. (qs(k).gt.1.E-5 .and. t(k).lt.273.)) then CLDFRA(K) = 1.0 - qvs(k) = qv(k) else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -3798,36 +3808,36 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & tc = t(k) - 273.15 if (tc .lt. -12.0) RH_00 = RH_00L - if (tc .ge. 20.0) then + if (tc .ge. 29.0) then CLDFRA(K) = 0.0 elseif (tc .ge. -12.0) then RHUM = MIN(rh(k), 1.0) - CLDFRA(K) = MAX(0., 1.0-SQRT((1.005-RHUM)/(1.005-RH_00))) + CLDFRA(K) = MAX(0., 1.0-SQRT((1.001-RHUM)/(1.001-RH_00))) else if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then !..For HRRR model, the following look OK. RHUM = MIN(rh(k), 1.45) - RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+100.) + RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+112.) if (RH_00 .ge. 1.5) then WRITE (dbg_msg,*) ' FATAL: RH_00 too large (1.5): ', RH_00, RH_00L, tc CALL wrf_error_fatal (dbg_msg) endif - CLDFRA(K) = MAX(0., 1.0-SQRT((1.5-RHUM)/(1.5-RH_00))) + CLDFRA(K) = MAX(0., 1.0-SQRT((1.46-RHUM)/(1.46-RH_00))) else !..but for the GFS model, RH is way lower. RHUM = MIN(rh(k), 1.05) - RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+100.) + RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+112.) if (RH_00 .ge. 1.05) then WRITE (dbg_msg,*) ' FATAL: RH_00 too large (1.05): ', RH_00, RH_00L, tc CALL wrf_error_fatal (dbg_msg) endif - CLDFRA(K) = MAX(0., 1.0-SQRT((1.05-RHUM)/(1.05-RH_00))) + CLDFRA(K) = MAX(0., 1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) endif endif - if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.9)) + if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.95)) if (debug_flag) then - WRITE (dbg_msg,*) 'DEBUG-GT: cloud fraction: ', RH_00, RHUM, CLDFRA(K) + WRITE (dbg_msg,*) 'DEBUG-GT: cloud fraction (k,RH_00, RHUM, CF): ',k,RH_00,RHUM,CLDFRA(K) CALL wrf_debug (150, dbg_msg) endif @@ -3847,8 +3857,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & WRITE (dbg_msg,*) 'DEBUG-GT: Made-up fake profile of clouds' CALL wrf_debug (150, dbg_msg) do k = kte, kts, -1 - write(dbg_msg,'(f7.2, 2x, f7.2, 2x, f6.4, 2x, f7.3, 2x, f15.7, 2x, f15.7)') & - & T(k)-273.15, P(k)*0.01, rh(k), cldfra(k)*100., qc(k)*1000.,qi(k)*1000. + write(dbg_msg,'(f7.2,2x,f7.2,2x,f6.4,2x,f7.3,x,f15.7,x,f15.7,x,f15.7)') & + & T(k)-273.15, P(k)*0.01, rh(k), cldfra(k)*100., qc(k)*1000.,qi(k)*1000., qs(k)*1000. CALL wrf_debug (150, dbg_msg) enddo endif @@ -3857,7 +3867,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & if (modify_qvapor) then DO k = kts,kte if (cldfra(k).gt.0.20 .and. cldfra(k).lt.1.0) then - qv(k) = qvs(k) + qv(k) = MAX(qv(k),qvs(k)) endif ENDDO endif @@ -3884,16 +3894,18 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& !..Local vars. REAL, DIMENSION(kts:kte):: theta REAL:: theta1, theta2, delz - INTEGER:: k, k2, k_tropo, k_m12C, k_cldb, k_cldt, kbot + INTEGER:: k, k2, k_tropo, k_m12C, k_cldb, k_cldt, kbot, k_p150 LOGICAL:: in_cloud character*512 dbg_msg !+---+ k_m12C = 0 + k_p150 = 0 DO k = kte, kts, -1 theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.)) if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10100.0) k_m12C = MAX(k_m12C, k) + if (P1d(k).gt.14999.0 .and. k_p150.eq.0) k_p150 = k ENDDO if (k_m12C .le. kts) k_m12C = kts @@ -3917,18 +3929,25 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& !.. tropopause height, as would any other diagnostic, so ensure resulting !.. k_tropo level is above 700hPa. - DO k = kte-3, kts, -1 + if ( (kte-k_p150) .lt. 3) k_p150 = kte-3 + DO k = k_p150-2, kts, -1 theta1 = theta(k) theta2 = theta(k+2) - delz = dz1d(k) + dz1d(k+1) + dz1d(k+2) - if ( ((((theta2-theta1)/delz) .lt. 10./1500. ) .AND. & - & (P1d(k).gt.8500.)) .or. (P1d(k).gt.70000.) ) then - goto 86 - endif + delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2) + if ( (((theta2-theta1)/delz).lt.10./1500.) .OR. P1d(k).gt.70000.) EXIT ENDDO - 86 continue k_tropo = MAX(kts+2, MIN(k+2, kte-1)) + if (k_tropo .gt. k_p150) then + DO k = kte-3, k_p150-2, -1 + theta1 = theta(k) + theta2 = theta(k+2) + delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2) + if ( (((theta2-theta1)/delz).lt.10./1500.) .AND. P1d(k).gt.9500.) EXIT + ENDDO + k_tropo = MAX(k_p150-1, MIN(k+2, kte-1)) + endif + if (k_tropo.gt.kte-2) then WRITE (dbg_msg,*) 'DEBUG-GT: CAUTION, tropopause appears to be very high up: ', k_tropo CALL wrf_debug (150, dbg_msg) @@ -3948,19 +3967,20 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& endif ENDDO -!..We would like to prevent fractional clouds below LCL in idealized -!.. situation with deep well-mixed convective PBL, that otherwise is -!.. likely to get clouds in more realistic capping inversion layer. +!..Be a bit more conservative with lower cloud fraction in scenario with +!.. well-mixed convective boundary layer below LCL. - kbot = kts+2 + kbot = kts+1 DO k = kbot, k_m12C - if ( (theta(k)-theta(k-1)) .gt. 0.025E-3*Dz1d(k)) EXIT + if ( (theta(k)-theta(k-1)) .gt. 0.010E-3*Dz1d(k)) EXIT ENDDO kbot = MAX(kts+1, k-2) DO k = kts, kbot - if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) cfr1d(k) = 0. + if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) cfr1d(k) = MAX(0.01,0.5*cfr1d(k)) + ENDDO + DO k = kts,k_tropo + if (cfr1d(k).gt.0.0) kbot = MIN(k,kbot) ENDDO - !..Starting below tropo height, if cloud fraction greater than 1 percent, !.. compute an approximate total layer depth of cloud, determine a total @@ -4001,16 +4021,16 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& CALL wrf_debug (150, dbg_msg) endif if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) & - & qi1d(k_cldb)=0.05*qvs1d(k_cldb) + & qi1d(k_cldb)=qi1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb) k = k_cldb endif k = k - 1 ENDDO - k_cldb = k_m12C + 3 + k_cldb = k_m12C + 5 in_cloud = .false. - k = k_m12C + 2 + k = k_m12C + 4 DO WHILE (.not. in_cloud .AND. k.gt.kbot) k_cldt = 0 if (cfr1d(k).ge.0.01) then @@ -4037,7 +4057,7 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& k = k_cldb elseif ((k_cldt - k_cldb + 1) .eq. 1) then if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) & - & qc1d(k_cldb)=0.05*qvs1d(k_cldb) + & qc1d(k_cldb)=qc1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb) k = k_cldb endif k = k - 1 @@ -4067,9 +4087,9 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) ! print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz do k = k1, k2 - max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k))) + max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k))) enddo - max_iwc = MIN(1.E-3, max_iwc) + max_iwc = MIN(2.E-3, max_iwc) this_dz = 0.0 do k = k1, k2 @@ -4079,7 +4099,7 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz - iwc = MAX(1.E-6, this_iwc*(1.-entr)) + iwc = MAX(5.E-6, this_iwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then qi(k) = qi(k) + cfr(k)*cfr(k)*iwc endif @@ -4108,9 +4128,9 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz do k = k1, k2 - max_lwc = MAX(1.E-6, max_lwc - qc(k)) + max_lwc = MAX(1.E-5, max_lwc - qc(k)) enddo - max_lwc = MIN(1.E-3, max_lwc) + max_lwc = MIN(2.E-3, max_lwc) this_dz = 0.0 do k = k1, k2 @@ -4120,7 +4140,7 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_lwc = max_lwc*this_dz/tdz - lwc = MAX(1.E-6, this_lwc*(1.-entr)) + lwc = MAX(5.E-6, this_lwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc endif @@ -4146,14 +4166,14 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) lwp = 0. iwp = 0. do k = kts, kte - if (cfr(k).gt.0.0) then + if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then lwp = lwp + qc(k)*Rho(k)*dz(k) iwp = iwp + qi(k)*Rho(k)*dz(k) endif enddo - if (lwp .gt. 1.5) then - xfac = 1.5/lwp + if (lwp .gt. 1.0) then + xfac = 1.0/lwp do k = kts, kte if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then qc(k) = qc(k)*xfac @@ -4161,8 +4181,8 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) enddo endif - if (iwp .gt. 1.5) then - xfac = 1.5/iwp + if (iwp .gt. 1.0) then + xfac = 1.0/iwp do k = kts, kte if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then qi(k) = qi(k)*xfac