diff --git a/README b/README index d5b942d812..6e98338938 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -WRF Model Version 4.2.1 +WRF Model Version 4.2.2 http://www2.mmm.ucar.edu/wrf/users/ diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 869dff658c..9e7202afed 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -3063,7 +3063,7 @@ package icedepth_one seaice_thickness_opt==1 - state:icedept #Time series options for text output package notseries process_time_series==0 - - package tseries process_time_series==1 - state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_p_profile,ts_w_profile -package tseries_add_solar process_time_series==2 - state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_cldfrac2d,ts_wvp,ts_lwp,ts_iwp,ts_swp,ts_lwp_tot,ts_iwp_tot,ts_swp_tot,ts_re_qc,ts_re_qi,ts_re_qs,ts_re_qc_tot,ts_re_qi_tot,ts_re_qs_tot,ts_tau_qc,ts_tau_qi,ts_tau_qs,ts_tau_qc_tot,ts_tau_qi_tot,ts_tau_qs_tot,ts_cbaseht,ts_ctopht,ts_cbaseht_tot,ts_ctopht_tot,ts_clrnidx,ts_p_profile,ts_w_profile +package tseries_add_solar process_time_series==2 - state:ts_hour,ts_u,ts_v,ts_q,ts_t,ts_psfc,ts_glw,ts_gsw,ts_hfx,ts_lh,ts_tsk,ts_tslb,ts_clw,ts_rainc,ts_rainnc,ts_u_profile,ts_v_profile,ts_gph_profile,ts_th_profile,ts_cldfrac2d,ts_wvp,ts_lwp,ts_iwp,ts_swp,ts_lwp_tot,ts_iwp_tot,ts_swp_tot,ts_re_qc,ts_re_qi,ts_re_qs,ts_re_qc_tot,ts_re_qi_tot,ts_re_qs_tot,ts_tau_qc,ts_tau_qi,ts_tau_qs,ts_tau_qc_tot,ts_tau_qi_tot,ts_tau_qs_tot,ts_cbaseht,ts_ctopht,ts_cbaseht_tot,ts_ctopht_tot,ts_clrnidx,ts_p_profile,ts_w_profile,ts_swdown,ts_swddni,ts_swddif,ts_swdownc,ts_swddnic,ts_swdown2,ts_swddni2,ts_swddif2,ts_swdownc2,ts_swddnic2 # WRF-HAILCAST state real HAILCAST_DHAIL1 ij misc 1 - r "HAILCAST_DHAIL1" "WRF-HAILCAST Hail Diameter, 1st rank order" "mm" diff --git a/Registry/registry.solar_fields b/Registry/registry.solar_fields index 6cd55e6ba6..2855d4642b 100644 --- a/Registry/registry.solar_fields +++ b/Registry/registry.solar_fields @@ -61,6 +61,16 @@ state real ts_cbaseht_tot ?! misc - - - "TS_C state real ts_ctopht_tot ?! misc - - - "TS_CTOPHT_TOT" "CLOUD TOP HEIGHT RES + UNRES" state real ts_clrnidx ?! misc - - - "TS_CLRNIDX" "CLEARNESS INDEX" state real ts_sza ?! misc - - - "TS_SZA" "SOLAR ZENITH ANGLE" +state real ts_swdown ?! misc - - - "TS_SWDOWN" "DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE" +state real ts_swddni ?! misc - - - "TS_SWDDNI" "SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE" +state real ts_swddif ?! misc - - - "TS_SWDDIF" "SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE" +state real ts_swdownc ?! misc - - - "TS_SWDOWNC" "DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE" +state real ts_swddnic ?! misc - - - "TS_SWDDNIC" "CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE" +state real ts_swdown2 ?! misc - - - "TS_SWDOWN2" "DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE FROM FARMS" +state real ts_swddni2 ?! misc - - - "TS_SWDDNI2" "SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS" +state real ts_swddif2 ?! misc - - - "TS_SWDDIF2" "SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE FROM FARMS" +state real ts_swdownc2 ?! misc - - - "TS_SWDOWNC2" "DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE FROM FARMS" +state real ts_swddnic2 ?! misc - - - "TS_SWDDNIC2" "CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS" # Package declarations diff --git a/arch/configure.defaults b/arch/configure.defaults index 140ebf7b1c..692fbfe9c3 100644 --- a/arch/configure.defaults +++ b/arch/configure.defaults @@ -44,7 +44,7 @@ RLFLAGS = CC_TOOLS = cc ########################################################### -#ARCH Linux i486 i586 i686, gfortran compiler with gcc #serial smpar dmpar dm+sm +#ARCH Linux i486 i586 i686 armv7l aarch64, gfortran compiler with gcc #serial smpar dmpar dm+sm # DESCRIPTION = GNU ($SFC/$SCC) DMPARALLEL = # 1 @@ -1834,7 +1834,7 @@ RLFLAGS = CC_TOOLS = $(SCC) ########################################################### -#ARCH Linux x86_64 ppc64le i486 i586 i686 #serial smpar dmpar dm+sm +#ARCH Linux x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm # DESCRIPTION = INTEL ($SFC/$SCC): HSW/BDW DMPARALLEL = # 1 @@ -1878,7 +1878,7 @@ RLFLAGS = CC_TOOLS = $(SCC) ########################################################### -#ARCH Linux KNL x86_64 ppc64le i486 i586 i686 #serial smpar dmpar dm+sm +#ARCH Linux KNL x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm # DESCRIPTION = INTEL ($SFC/$SCC): KNL MIC DMPARALLEL = # 1 diff --git a/chem/chemics_init.F b/chem/chemics_init.F index ea27ae3c58..66bc2f5024 100755 --- a/chem/chemics_init.F +++ b/chem/chemics_init.F @@ -1944,7 +1944,11 @@ subroutine chem_init (id,chem,emis_ant,scalar,dt,bioemdt,photdt,chemdt,stepbioe, ! drydep_select: SELECT CASE(config_flags%gas_drydep_opt) CASE (WESELY) - CALL wrf_debug(15,'initializing dry dep (wesely)') + IF (numgas .eq. 0) THEN + CALL wrf_error_fatal("ERROR: numgas = 0, SELECTED CHEM OPT IS & + &NOT COMPATIBLE WITH WESELY DRY DEPOSITION") + ENDIF + CALL wrf_debug(15,'initializing dry dep (wesely)') call dep_init( id, config_flags, numgas, mminlu_loc, & its, ite, jts, jte, ide, jde ) diff --git a/chem/module_emissions_anthropogenics.F b/chem/module_emissions_anthropogenics.F index 6f168e67c9..52358038a2 100755 --- a/chem/module_emissions_anthropogenics.F +++ b/chem/module_emissions_anthropogenics.F @@ -241,7 +241,7 @@ subroutine add_anthropogenics(id,dtstep,dz8w,config_flags,rho_phy,alt, & if (config_flags%aircraft_emiss_opt == 1 ) then do j=jts,jte do k=kts,min(config_flags%kemit_aircraft,kte) - conv_rho(its:ite)=4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.) + conv_rho(its:ite)=4.828e-4/rho_phy(its:ite,k,j)*dtstep/(dz8w(its:ite,k,j)*60.) if( p_no >= param_first_scalar ) then chem(its:ite,k,j,p_no) = chem(its:ite,k,j,p_no) + emis_aircraft(its:ite,k,j,p_eac_no) *conv_rho(its:ite) endif diff --git a/chem/module_mosaic_addemiss.F b/chem/module_mosaic_addemiss.F index 874bdcde5a..723224c97d 100644 --- a/chem/module_mosaic_addemiss.F +++ b/chem/module_mosaic_addemiss.F @@ -601,7 +601,7 @@ subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, & do i = its, ite ! compute mass biomass burning emissions [(ug/m3)*m/s] for each species ! using the apportioning fractions - aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf) + IF ( p_ebu_sulf .gt. 1) aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf) aem_oc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_oc) aem_bc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_bc) ! Option to calculate OIN fraction of total PM for fire emissions diff --git a/dyn_em/module_big_step_utilities_em.F b/dyn_em/module_big_step_utilities_em.F index e31c0874db..3da86cdb99 100644 --- a/dyn_em/module_big_step_utilities_em.F +++ b/dyn_em/module_big_step_utilities_em.F @@ -1506,7 +1506,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & jtf=MIN(jte,jde-1) IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1 - IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2 + IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-2 DO j = j_start, jtf @@ -1539,7 +1539,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & jtf=MIN(jte,jde-1) IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1 - IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2 + IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = ite-2 DO j = j_start, jtf @@ -1574,7 +1574,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & jtf=MIN(jte,jde-1) IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2 - IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3 + IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-3 DO j = j_start, jtf @@ -1664,7 +1664,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & jtf=MIN(jte,jde-1) IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2 - IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3 + IF ( (config_flags%open_xe) .and. ite == ide ) itf = ite-3 DO j = j_start, jtf diff --git a/dyn_em/module_diffusion_em.F b/dyn_em/module_diffusion_em.F index 57fd621e7a..0c98282352 100644 --- a/dyn_em/module_diffusion_em.F +++ b/dyn_em/module_diffusion_em.F @@ -2160,9 +2160,18 @@ SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv, & ! tke_seed if the drag and flux are off. c_k = config_flags%c_k - tke_seed = tke_seed_value - if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. & - (config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0. + tke_seed = 0. + IF (config_flags%isfflx .eq. 0) THEN + IF ((config_flags%diff_opt .eq. 2) .and. (config_flags%bl_pbl_physics .eq. 0)) THEN + IF( (config_flags%tke_drag_coefficient .lt. epsilon) .and. & + (config_flags%tke_heat_flux .lt. epsilon) ) THEN + tke_seed = tke_seed_value + ENDIF + ELSE + !tke_drag_coefficient and tke_heat_flux are irrelevant here + tke_seed = tke_seed_value + ENDIF + ENDIF DO j = j_start, j_end DO k = kts+1, ktf-1 @@ -4282,11 +4291,13 @@ SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, & DO i = i_start, i_end cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV)) hfx(i,j)=heat_flux*cpm*rho(i,kts,j) ! provided for output only - rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & - -g*heat_flux*rho(i,kts,j)/dnw(kts) if(config_flags%use_theta_m == 1)THEN rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & + -g*heat_flux*(1.+rvovrd*moist(i,kts,j,P_QV))*rho(i,kts,j)/dnw(kts) & -g*1.61*theta(i,kts,j)*qfx(i,j)/dnw(kts) + ELSE + rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & + -g*heat_flux*rho(i,kts,j)/dnw(kts) ENDIF ENDDO ENDDO @@ -4297,11 +4308,13 @@ SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, & cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV)) heat_flux = hfx(i,j)/cpm - rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & - -g*heat_flux/dnw(kts) if(config_flags%use_theta_m == 1)THEN rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & + -g*heat_flux*(1.+rvovrd*moist(i,kts,j,P_QV))/dnw(kts) & -g*1.61*theta(i,kts,j)*qfx(i,j)/dnw(kts) + ELSE + rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & + -g*heat_flux/dnw(kts) ENDIF ENDDO @@ -5133,7 +5146,7 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, & DO i = i_start, i_end delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j)) pu1=pu(delxy,hpbl(i,j)) - IF(sfcflg(i,j).and.sflux(i,j).GT.0.0)THEN + IF((sfcflg(i,j).and.sflux(i,j).GT.0.0) .and. (hpbl(i,j) .GT. 0)) THEN !nonlocal momentum flux based on Brown and Grant (1997) brint = -sm*ust(i,j)*ust(i,j)*wstar3(i,j)/(hpbl(i,j)*wscale(i,j)**4) gamu(i,j) = pu1 * brint*u_phy(i,1,j)/wspd(i,j) @@ -5165,8 +5178,12 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, & deltaoh(i,j) = max(ezfac*deltaoh(i,j),hpbl(i,j)-za(i,kpbl(i,j)-1,j)-1.) deltaoh(i,j) = min(deltaoh(i,j), hpbl(i,j)) - rigs(i,j) = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.) - rigs(i,j) = max(min(rigs(i,j), rigsmax),rimin) + if ((du .ne. 0) .or. (dv .ne. 0)) then + rigs(i,j) = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.) + rigs(i,j) = max(min(rigs(i,j), rigsmax),rimin) + else + rigs(i,j) = rigsmax + endif enlfrac2(i,j) = max(min(wm3/wstar3(i,j)/(1.+cpent/rigs(i,j)),entfmax),entfmin) enlfrac2(i,j) = enlfrac2(i,j)*enlfrac ENDIF @@ -5182,44 +5199,41 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, & ENDDO ENDDO ENDDO - - DO j = j_start, j_end - DO i = i_start, i_end - deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j) - ENDDO - ENDDO DO j = j_start, j_end DO i = i_start, i_end - delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j)) - mlfrac = mltop-deltaoh(i,j) - ezfrac = mltop+deltaoh(i,j) - zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.) - sfcfracn = max(sfcfracn1,zfacmf(i,1,j)) -! - sflux0 = (a11+a12*sfcfracn)*sflux(i,j) - snlflux0 = nlfrac*sflux0 - amf1 = snlflux0/sfcfracn - amf2 = -snlflux0/(mlfrac-sfcfracn) - bmf2 = -mlfrac*amf2 - amf3 = snlflux0*enlfrac2(i,j)/deltaoh(i,j) - bmf3 = -amf3*mlfrac - hfxpbl(i,j) = amf3+bmf3 - pth1 = pthnl(delxy,hpbl(i,j)) - hfxpbl(i,j) = hfxpbl(i,j)*pth1 + IF (pblflg(i,j)) THEN + deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j) + delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j)) + mlfrac = mltop-deltaoh(i,j) + ezfrac = mltop+deltaoh(i,j) + zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.) + sfcfracn = max(sfcfracn1,zfacmf(i,1,j)) + ! + sflux0 = (a11+a12*sfcfracn)*sflux(i,j) + snlflux0 = nlfrac*sflux0 + amf1 = snlflux0/sfcfracn + amf2 = -snlflux0/(mlfrac-sfcfracn) + bmf2 = -mlfrac*amf2 + amf3 = snlflux0*enlfrac2(i,j)/deltaoh(i,j) + bmf3 = -amf3*mlfrac + hfxpbl(i,j) = amf3+bmf3 + pth1 = pthnl(delxy,hpbl(i,j)) + hfxpbl(i,j) = hfxpbl(i,j)*pth1 - DO k = kts, ktf - zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin) - IF(pblflg(i,j).and.k.LT.kpbl(i,j)) THEN - IF(zfacmf(i,k,j).LE.sfcfracn) THEN - nlflux(i,k,j) = amf1*zfacmf(i,k,j) - ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN - nlflux(i,k,j) = amf2*zfacmf(i,k,j)+bmf2 - ENDIF - nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j)) - nlflux(i,k,j) = nlflux(i,k,j)*pth1 - ENDIF - ENDDO + DO k = kts, ktf + zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin) + IF(k.LT.kpbl(i,j)) THEN + IF(zfacmf(i,k,j).LE.sfcfracn) THEN + nlflux(i,k,j) = amf1*zfacmf(i,k,j) + ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN + nlflux(i,k,j) = amf2*zfacmf(i,k,j)+bmf2 + ENDIF + nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j)) + nlflux(i,k,j) = nlflux(i,k,j)*pth1 + ENDIF + ENDDO + ENDIF ENDDO ENDDO END SUBROUTINE nonlocal_flux @@ -5237,10 +5251,15 @@ FUNCTION pthnl(d,h) REAL,PARAMETER :: b1 = 2.0, b2 = 0.875 real :: d,h,doh,num,den - doh = d/h - num = a1*(doh)**b1 + a2*(doh)**b2+a3 - den = a4*(doh)**b1 + a5*(doh)**b2+a6 - pthnl = a7*num/den + (1. - a7) + if (h .ne. 0) then + doh = d/h + num = a1*(doh)**b1 + a2*(doh)**b2+a3 + den = a4*(doh)**b1 + a5*(doh)**b2+a6 + pthnl = a7*num/den + (1. - a7) + else + pthnl = 1. + endif + pthnl = max(pthnl,pmin) pthnl = min(pthnl,pmax) @@ -5261,10 +5280,14 @@ FUNCTION pthl(d,h) REAL,PARAMETER :: b1 = 2.0, b2 = 0.5 REAL :: d,h,doh,num,den - doh = d/h - num = a1*(doh)**b1 + a2*(doh)**b2+a3 - den = a4*(doh)**b1 + a5*(doh)**b2+a6 - pthl = a7*num/den+(1. - a7) + if (h .ne. 0) then + doh = d/h + num = a1*(doh)**b1 + a2*(doh)**b2+a3 + den = a4*(doh)**b1 + a5*(doh)**b2+a6 + pthl = a7*num/den+(1. - a7) + else + pthl = 1. + endif pthl = max(pthl,pmin) pthl = min(pthl,pmax) @@ -5284,10 +5307,14 @@ FUNCTION pu(d,h) REAL,PARAMETER :: b1 = 2.0, b2 = 0.6666667 REAL :: d,h,doh,num,den - doh = d/h - num = a1*(doh)**b1 + a2*(doh)**b2 - den = a3*(doh)**b1 + a4*(doh)**b2+a5 - pu = num/den + if (h .ne. 0) then + doh = d/h + num = a1*(doh)**b1 + a2*(doh)**b2 + den = a3*(doh)**b1 + a4*(doh)**b2+a5 + pu = num/den + else + pu = 1. + endif pu = max(pu,pmin) pu = min(pu,pmax) @@ -7717,7 +7744,7 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,& INTEGER :: i_start, i_end, j_start, j_end REAL :: V0_u,V0_v,ustar,beta - REAL :: heat_flux, moist_flux + REAL :: heat_flux, moist_flux, qfac REAL :: cpm,rdz_w,rhoavg_u,rhoavg_v,rdz_z ! End declarations. !----------------------------------------------------------------------- @@ -8094,17 +8121,23 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,& - dt*rdz_w*nlflux_rho(i,k+1,j) & + dt*rdz_w*hfx(i,j)/cpm ELSE + qfac = 1.+rvovrd*moist(i,kts,j,P_QV) d(k) = var_mix(i,kts,j) & - - dt*rdz_w*nlflux_rho(i,k+1,j) & - + dt*rdz_w*(hfx(i,j)/cpm+1.61*theta(i,kts,j)*qfx(i,j)) + - dt*rdz_w*nlflux_rho(i,k+1,j)*qfac & + + dt*rdz_w*(qfac*hfx(i,j)/cpm+1.61*theta(i,kts,j)*qfx(i,j)) ENDIF + IF(config_flags%use_theta_m == 1)THEN + qfac = 1.+rvovrd*moist(i,kts,j,P_QV) + ELSE + qfac = 1. + ENDIF DO k = kts+1, ktf-1 rdz_w = -g/dnw(k)/(c1h(k)*mu(i,j) + c2h(k)) a(k) = -rdz_w*rdz(i,k,j)*xkxavg(i,k,j)*dt b(k) = 1.+rdz_w*(rdz(i,k+1,j)*xkxavg(i,k+1,j)+rdz(i,k,j)*xkxavg(i,k,j))*dt c(k) = -rdz_w*rdz(i,k+1,j)*xkxavg(i,k+1,j)*dt - d(k) = -rdz_w*(nlflux_rho(i,k+1,j)-nlflux_rho(i,k,j))*dt + var_mix(i,k,j) + d(k) = -rdz_w*(nlflux_rho(i,k+1,j)-nlflux_rho(i,k,j))*qfac*dt + var_mix(i,k,j) ENDDO a(ktf) = 0. @@ -8146,15 +8179,21 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,& DO i = i_start, i_end DO k = kts, ktf-1 rdz_w = - g/dnw(k)/(c1h(k)*mu(i,j)+c2h(k)) - beta = 1.5*sqrt(tke(i,k,j))/l_diss(i,k,j) + IF (l_diss(i,k,j) .ne. 0) THEN + beta = 1.5*sqrt(tke(i,k,j))/l_diss(i,k,j) + ELSE + beta = 0. + ENDIF a(k) = - 2.0*xkxavg(i,k,j)*dt*rdz_w*rdz(i,k,j) b(k) = 1.0 + 2.0*dt*rdz_w*(rdz(i,k,j)*xkxavg(i,k,j) & + rdz(i,k+1,j)*xkxavg(i,k+1,j)) & + dt*beta c(k) = - 2.0*xkxavg(i,k+1,j)*dt*rdz(i,k+1,j)*rdz_w - d(k) = tke(i,k,j) & - + 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j) + d(k) = tke(i,k,j) + IF (l_diss(i,k,j) .ne. 0) THEN + d(k) = d(k) + 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j) + ENDIF ENDDO a(ktf) = 0. !-1 diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index bcd346f49b..ae0a7f1388 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -1097,8 +1097,8 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_CHEM_E_5.inc" ELSE - WRITE(message,*)'solve_em: invalid h_mom_adv_order = ',& - & config_flags%h_mom_adv_order + WRITE(message,*)'solve_em: invalid h_sca_adv_order = ',& + & config_flags%h_sca_adv_order ENDIF ENDIF #endif @@ -1421,12 +1421,12 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & #ifdef DM_PARALLEL IF( config_flags%shcu_physics == CAMUWSHCUSCHEME ) THEN CALL wrf_debug ( 200 , ' call HALO CHEM AFTER SHALLOW CUMULUS' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_CHEM_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_CHEM_E_5.inc" ELSE - WRITE(message,*)'module_first_rk_step_part1: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(message,*)'module_first_rk_step_part1: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(message)) ENDIF ENDIF diff --git a/dyn_em/module_first_rk_step_part2.F b/dyn_em/module_first_rk_step_part2.F index 009b09fca6..518592858f 100644 --- a/dyn_em/module_first_rk_step_part2.F +++ b/dyn_em/module_first_rk_step_part2.F @@ -689,12 +689,12 @@ SUBROUTINE first_rk_step_part2 ( grid , config_flags & # include "HALO_EM_PHYS_DIFFUSION.inc" ENDIF - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_TKE_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_TKE_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF #endif diff --git a/dyn_em/module_initialize_ideal.F b/dyn_em/module_initialize_ideal.F index 22817bd550..058b299233 100644 --- a/dyn_em/module_initialize_ideal.F +++ b/dyn_em/module_initialize_ideal.F @@ -1120,10 +1120,10 @@ SUBROUTINE init_domain_rk ( grid & ! rebalance hydrostatically - DO k = 2,kte + DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) @@ -1165,10 +1165,10 @@ SUBROUTINE init_domain_rk ( grid & ! rebalance hydrostatically - DO k = 2,kte + DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) @@ -1212,8 +1212,8 @@ SUBROUTINE init_domain_rk ( grid & DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) @@ -1245,8 +1245,8 @@ SUBROUTINE init_domain_rk ( grid & DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) ENDDO @@ -1306,8 +1306,8 @@ SUBROUTINE init_domain_rk ( grid & DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) @@ -1355,8 +1355,8 @@ SUBROUTINE init_domain_rk ( grid & DO k = 2,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*( & - ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & - (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) + ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ & + (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index f2d6906f07..be68cc3756 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -4930,6 +4930,8 @@ END SUBROUTINE find_my_parent2 #ifdef VERT_UNIT +!gfortran -DVERT_UNIT -ffree-form -ffree-line-length-none module_initialize_real.F -o vert.exe + !This is a main program for a small unit test for the vertical interpolation. program vint @@ -4961,6 +4963,7 @@ program vint logical, parameter :: use_surface = .TRUE. ! .FALSE. ! .TRUE. real , parameter :: zap_close_levels = 500. ! 100. integer, parameter :: force_sfc_in_vinterp = 6 ! 0 ! 6 + integer, parameter :: id = 1 integer :: k @@ -4975,7 +4978,7 @@ program vint print *,'UNIT TEST FOR VERTICAL INTERPOLATION' print *,'------------------------------------' print *,' ' - do lagrange_order = 1 , 9 , 8 + do lagrange_order = 1 , 1 print *,' ' print *,'------------------------------------' print *,'Lagrange Order = ',lagrange_order @@ -5006,7 +5009,7 @@ program vint generic , 'T' , & interp_type , lagrange_order , extrap_type , & lowest_lev_from_sfc , use_levels_below_ground , use_surface , & - zap_close_levels , force_sfc_in_vinterp , & + zap_close_levels , force_sfc_in_vinterp , id , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) @@ -5052,8 +5055,6 @@ subroutine fillitup ( fo , po , fn , pn , & integer :: i , j , k - real , parameter :: piov2 = 3.14159265358 / 2. - k = 1 do j = jts , jte do i = its , ite @@ -5065,6 +5066,7 @@ subroutine fillitup ( fo , po , fn , pn , & do j = jts , jte do i = its , ite po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) ) +! po(i,k,j) = FILL IN YOUR INPUT PRESSURE LEVELS end do end do end do @@ -5074,7 +5076,7 @@ subroutine fillitup ( fo , po , fn , pn , & do j = jts , jte do i = its , ite fo(i,k,j) = po(i,k,j) -! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. ) +! fo(i,k,j) = FILL IN YOUR COLUMN OF PRESS_LEVEL FIELD end do end do end do @@ -5083,7 +5085,6 @@ subroutine fillitup ( fo , po , fn , pn , & do j = jts , jte do i = its , ite fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)*((102000.-po(i,k,j))/102000.))*102000. -! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. ) end do end do end do @@ -5095,6 +5096,7 @@ subroutine fillitup ( fo , po , fn , pn , & do j = jts , jte do i = its , ite pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) ) +! pn(i,k,j) = FILL IN A COLUMN OF KNOWN FULL-LEVEL PRESSURES ON ETA SURFACES end do end do end do @@ -5113,7 +5115,7 @@ subroutine fillitup ( fo , po , fn , pn , & do j = jts , jte do i = its , ite fn(i,k,j) = pn(i,k,j) -! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. ) +! fn(i,k,j) = FILL IN COLUMN OF HALF LEVEL FIELD end do end do end do @@ -5137,6 +5139,12 @@ function skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ skip_middle_points_t = .false. end function skip_middle_points_t +subroutine wrf_message(level,message) + character(len=*), intent(in) :: message + integer, intent(in) :: level + print *,trim(message) +end subroutine wrf_message + #endif !--------------------------------------------------------------------- diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 0fe7edebdd..6b2aca6064 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -320,7 +320,7 @@ SUBROUTINE solve_em ( grid , config_flags & ! see matching halo calls below for stencils !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_CHEM' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_CHEM_E_3.inc" IF( config_flags%progn > 0 ) THEN # include "HALO_EM_SCALAR_E_3.inc" @@ -328,7 +328,7 @@ SUBROUTINE solve_em ( grid , config_flags & IF( config_flags%cu_physics == CAMZMSCHEME ) THEN # include "HALO_EM_SCALAR_E_3.inc" ENDIF - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_CHEM_E_5.inc" IF( config_flags%cu_physics == CAMZMSCHEME ) THEN # include "HALO_EM_SCALAR_E_5.inc" @@ -337,7 +337,7 @@ SUBROUTINE solve_em ( grid , config_flags & # include "HALO_EM_SCALAR_E_5.inc" ENDIF ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF @@ -346,12 +346,12 @@ SUBROUTINE solve_em ( grid , config_flags & ! see matching halo calls below for stencils !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_tracer' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_TRACER_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_TRACER_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF @@ -2288,12 +2288,12 @@ SUBROUTINE solve_em ( grid , config_flags & BENCH_START(tke_adv_tim) TKE_advance: IF (config_flags%km_opt .eq. 2.or.config_flags%km_opt.eq.5) then ! XZ #ifdef DM_PARALLEL - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_TKE_ADVECT_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_TKE_ADVECT_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF #endif @@ -3116,12 +3116,13 @@ SUBROUTINE solve_em ( grid , config_flags & ! scalar x #ifdef DM_PARALLEL - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_mom_adv_order <= 4 .and. config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_D2_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_mom_adv_order <= 6 .and. config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_D2_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order or h_sca_adv_order = ', & + config_flags%h_mom_adv_order, config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF # include "PERIOD_BDY_EM_D.inc" @@ -3252,13 +3253,13 @@ SUBROUTINE solve_em ( grid , config_flags & ! moist, chem, scalar, tke x - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN IF ( (config_flags%tke_adv_opt /= ORIGINAL .and. config_flags%tke_adv_opt /= WENO_SCALAR) .and. (rk_step == rk_order-1) ) THEN # include "HALO_EM_TKE_5.inc" ELSE # include "HALO_EM_TKE_3.inc" ENDIF - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN IF ( (config_flags%tke_adv_opt /= ORIGINAL .and. config_flags%tke_adv_opt /= WENO_SCALAR) .and. (rk_step == rk_order-1) ) THEN # include "HALO_EM_TKE_7.inc" ELSE @@ -4264,12 +4265,13 @@ SUBROUTINE solve_em ( grid , config_flags & #ifdef DM_PARALLEL - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_mom_adv_order <= 4 .and. config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_D3_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_mom_adv_order <= 6 .and. config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_D3_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order or h_sca_adv_order = ', & + config_flags%h_mom_adv_order, config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF # include "PERIOD_BDY_EM_D3.inc" @@ -4609,12 +4611,13 @@ SUBROUTINE solve_em ( grid , config_flags & ! see above !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_E' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_mom_adv_order <= 4 .and. config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_mom_adv_order <= 6 .and. config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order or h_sca_adv_order = ', & + config_flags%h_mom_adv_order, config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF #endif @@ -4625,12 +4628,12 @@ SUBROUTINE solve_em ( grid , config_flags & ! see above !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_MOIST' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_MOIST_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_MOIST_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF @@ -4639,12 +4642,12 @@ SUBROUTINE solve_em ( grid , config_flags & ! see above !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_CHEM' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_CHEM_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_CHEM_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF @@ -4653,12 +4656,12 @@ SUBROUTINE solve_em ( grid , config_flags & ! see above !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_TRACER' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_TRACER_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_TRACER_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF @@ -4667,12 +4670,12 @@ SUBROUTINE solve_em ( grid , config_flags & ! see above !-------------------------------------------------------------- CALL wrf_debug ( 200 , ' call HALO_RK_SCALAR' ) - IF ( config_flags%h_mom_adv_order <= 4 ) THEN + IF ( config_flags%h_sca_adv_order <= 4 ) THEN # include "HALO_EM_SCALAR_E_3.inc" - ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN + ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN # include "HALO_EM_SCALAR_E_5.inc" ELSE - WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order + WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order CALL wrf_error_fatal(TRIM(wrf_err_message)) ENDIF ENDIF diff --git a/external/RSL_LITE/module_dm.F b/external/RSL_LITE/module_dm.F index de586ca32b..ddc09bc5ef 100644 --- a/external/RSL_LITE/module_dm.F +++ b/external/RSL_LITE/module_dm.F @@ -228,6 +228,9 @@ SUBROUTINE wrf_dm_initialize CALL nl_set_nproc_y ( 1, ntasks_y ) WRITE( wrf_err_message , * )'Ntasks in X ',ntasks_x,', ntasks in Y ',ntasks_y CALL wrf_message( wrf_err_message ) +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + CALL MPI_BARRIER( local_communicator, ierr ) +#endif RETURN END SUBROUTINE wrf_dm_initialize @@ -1295,7 +1298,7 @@ INTEGER function getrealmpitype() RETURN END FUNCTION getrealmpitype - REAL FUNCTION wrf_dm_max_int ( inval ) + INTEGER FUNCTION wrf_dm_max_int ( inval ) IMPLICIT NONE #ifndef STUBMPI INTEGER, intent(in) :: inval diff --git a/external/io_int/module_io_int_read.F90 b/external/io_int/module_io_int_read.F90 index 5e2a8f6e68..8bd9d7d1e1 100644 --- a/external/io_int/module_io_int_read.F90 +++ b/external/io_int/module_io_int_read.F90 @@ -87,7 +87,7 @@ end module module_io_int_read module module_io_int_read use module_io_int_idx, only: io_int_loc, r_info - use, intrinsic :: iso_c_binding, only: c_int32_t + use, intrinsic :: iso_c_binding, only: c_int32_t, c_int64_t implicit none @@ -127,9 +127,9 @@ subroutine read_i0(ifd, records, varname, dst, ierr) integer, intent(out) :: dst integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count - integer :: tmp + integer(c_int64_t) :: offset + integer(c_int32_t) :: count + integer :: tmp call io_int_loc(varname, records, offset, count, ierr) if (ierr .ne. 0) then @@ -161,8 +161,8 @@ subroutine read_i1(ifd, records, varname, dst, ierr) integer, intent(inout) :: dst(:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i integer :: its, ite @@ -219,8 +219,8 @@ subroutine read_i2(ifd, records, varname, dst, ierr) integer, intent(inout) :: dst(:,:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i, j integer :: its, ite, jts, jte @@ -279,8 +279,8 @@ subroutine read_i3(ifd, records, varname, dst, ierr) integer, intent(inout) :: dst(:,:,:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i, j, k integer :: its, ite, jts, jte, kts, kte @@ -343,8 +343,8 @@ subroutine read_r0(ifd, records, varname, dst, ierr) real, intent(out) :: dst integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: tmp call io_int_loc(varname, records, offset, count, ierr) @@ -377,8 +377,8 @@ subroutine read_r1(ifd, records, varname, dst, ierr) real, intent(inout) :: dst(:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i integer :: its, ite @@ -435,8 +435,8 @@ subroutine read_r2(ifd, records, varname, dst, ierr) real, intent(inout) :: dst(:,:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i, j integer :: its, ite, jts, jte @@ -495,8 +495,8 @@ subroutine read_r3(ifd, records, varname, dst, ierr) real, intent(inout) :: dst(:,:,:) integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i, j, k integer :: its, ite, jts, jte, kts, kte @@ -559,8 +559,8 @@ subroutine read_c1(ifd, records, varname, dst, ierr) character(len=*), intent(inout) :: dst integer, intent(out) :: ierr - integer(kind=mpi_offset_kind) :: offset - integer :: count + integer(c_int64_t) :: offset + integer(c_int32_t) :: count integer :: num integer :: i integer, allocatable, dimension(:) :: tmp diff --git a/inc/version_decl b/inc/version_decl index fc13f34730..0a849d69f4 100644 --- a/inc/version_decl +++ b/inc/version_decl @@ -1 +1 @@ - CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.2.1' + CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.2.2' diff --git a/phys/module_cu_mskf.F b/phys/module_cu_mskf.F index 1f107a9ebe..0292f99713 100644 --- a/phys/module_cu_mskf.F +++ b/phys/module_cu_mskf.F @@ -1762,8 +1762,8 @@ subroutine mskf_mphy(su, qu, mu, du, cmel, cmei, zf, pm, te, qe, ep mtime=deltat/900._r8 mtimec=deltat/900._r8 - mtime = AMAX1(1.0,mtime) !TWG remove time scale limitation from CAM5 - mtimec = AMAX1(1.0,mtimec) + mtime = max(1.0_r8,mtime) !TWG remove time scale limitation from CAM5 + mtimec = max(1.0_r8,mtimec) ! conservation of qc @@ -5160,19 +5160,19 @@ SUBROUTINE MSKF_eta_PARA (I, J, & JK = KX-KQ+1 ! print *,'kf qliq=', QLIQ(KQ) - QLIQ(KQ) = amax1(0.0,zmqliq(1,JK)) - QICE(KQ) = amax1(0.0,zmqice(1,JK)) + QLIQ(KQ) = max(0._r8,zmqliq(1,JK)) + QICE(KQ) = max(0._r8,zmqice(1,JK)) !TWG 06/14/16 - QRAIN(KQ) = amax1(0.0,zmqrain(1,JK)) - QSNOW(KQ) = amax1(0.0,zmqsnow(1,JK)) - NLIQ(KQ) = amax1(0.0,ncmp(1,JK)) - NICE(KQ) = amax1(0.0,nimp(1,JK)) - NRAIN(KQ) = amax1(0.0,nrmp(1,JK)) - NSNOW(KQ) = amax1(0.0,nsmp(1,JK)) - CCN(KQ) = amax1(0.0,zmccn(1,JK)) - EFFCH(KQ) = MAX(2.49, MIN(effc(1,JK), 50.)) - EFFIH(KQ) = MAX(4.99, MIN(effi(1,JK), 125.)) - EFFSH(KQ) = MAX(9.99, MIN(effs(1,JK), 999.)) + QRAIN(KQ) = max(0._r8,zmqrain(1,JK)) + QSNOW(KQ) = max(0._r8,zmqsnow(1,JK)) + NLIQ(KQ) = max(0._r8,ncmp(1,JK)) + NICE(KQ) = max(0._r8,nimp(1,JK)) + NRAIN(KQ) = max(0._r8,nrmp(1,JK)) + NSNOW(KQ) = max(0._r8,nsmp(1,JK)) + CCN(KQ) = max(0._r8,zmccn(1,JK)) + EFFCH(KQ) = max(2.49_r8, min(effc(1,JK), 50._r8)) + EFFIH(KQ) = max(4.99_r8, min(effi(1,JK), 125._r8)) + EFFSH(KQ) = max(9.99_r8, min(effs(1,JK), 999._r8)) ! END TWG DETLQ(KQ)= QLIQ(KQ)*UDR(KQ) DETIC(KQ)= QICE(KQ)*UDR(KQ) diff --git a/phys/module_cu_scalesas.F b/phys/module_cu_scalesas.F index 937ed93172..ec11287d73 100755 --- a/phys/module_cu_scalesas.F +++ b/phys/module_cu_scalesas.F @@ -4222,7 +4222,7 @@ subroutine mfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, & if(cnvflg(i)) then if (gdx(i) < dxcrt) then scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) - scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0._kind_phys) + scaldfunc(i) = max(min(scaldfunc(i), 1._kind_phys), 0._kind_phys) else scaldfunc(i) = 1.0 endif diff --git a/phys/module_diag_solar.F b/phys/module_diag_solar.F index 919e5bab14..4461aefa54 100644 --- a/phys/module_diag_solar.F +++ b/phys/module_diag_solar.F @@ -28,7 +28,8 @@ MODULE module_diag_solar ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, qc, qi, qs, & + SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, & + param_first_scalar, p_qc, p_qi, p_qs, qv, qc, qi, qs, & qc_tot, qi_tot, has_reqc, has_reqi, has_reqs, f_qv, f_qc, f_qi, f_qs, & re_cloud, re_ice, re_snow, clrnidx, sza, cldfrac2d, wvp2d, lwp2d, iwp2d, swp2d, & wp2d_sum, lwp2d_tot, iwp2d_tot, wp2d_tot_sum, re_cloud_path, re_ice_path, re_snow_path, & @@ -45,7 +46,7 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: coszen, swdnb, swdnt REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ph, phb, cldfrac3d, qv, qc, qi, qs, qc_tot, qi_tot, & re_cloud, re_ice, re_snow, rho, dz8w - INTEGER, INTENT(IN) :: has_reqc, has_reqi, has_reqs + INTEGER, INTENT(IN) :: param_first_scalar, p_qc, p_qi, p_qs, has_reqc, has_reqi, has_reqs LOGICAL, INTENT(IN) :: f_qv, f_qc, f_qi, f_qs REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: clrnidx, sza, cldfrac2d, wvp2d, lwp2d, iwp2d, swp2d, wp2d_sum, & lwp2d_tot, iwp2d_tot, wp2d_tot_sum, re_cloud_path, re_ice_path, re_snow_path, re_cloud_path_tot, re_ice_path_tot, & @@ -101,9 +102,9 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, end if !!! LIQUID WATER - if (f_qc) then + if (f_qc .and. p_qc > param_first_scalar) then !!! RESOLVED !!! - ! Calc liquid water pth + ! Calc liquid water path q_aux = integrate_1var (rhodz, qc(i, :, j), kms, kme, kts, kte) lwp = q_aux lwp2d(i, j) = SIGN( MAX( lwp, 0.0 ), 1.0 ) @@ -124,9 +125,6 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, else tau_qc(i, j) = 0.0 end if - else - re_cloud_path(i, j) = MISSING - tau_qc(i, j) = MISSING end if !!! TOTAL (RESOLVED + UNRESOLVED) !!! @@ -151,14 +149,18 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, else tau_qc_tot(i, j) = 0.0 end if - else - re_cloud_path_tot(i, j) = MISSING - tau_qc_tot(i, j) = MISSING end if + else + lwp2d(i, j) = MISSING + re_cloud_path(i, j) = MISSING + tau_qc(i, j) = MISSING + lwp2d_tot(i, j) = MISSING + re_cloud_path_tot(i, j) = MISSING + tau_qc_tot(i, j) = MISSING end if !!! ICE - if (f_qi) then + if (f_qi .and. p_qi > param_first_scalar) then !!! RESOLVED !!! ! Calc ice water path q_aux = integrate_1var (rhodz, qi(i, :, j), kms, kme, kts, kte) @@ -186,9 +188,6 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, else tau_qi(i, j) = 0.0 end if - else - re_ice_path(i, j) = MISSING - tau_qi(i, j) = MISSING end if !!! TOTAL (RESOLVED + UNRESOLVED) !!! @@ -218,14 +217,18 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, else tau_qi_tot(i, j) = 0.0 end if - else - re_ice_path_tot(i, j) = MISSING - tau_qi_tot(i, j) = MISSING end if + else + iwp2d(i, j) = MISSING + re_ice_path(i, j) = MISSING + tau_qi(i, j) = MISSING + iwp2d_tot(i, j) = MISSING + re_ice_path_tot(i, j) = MISSING + tau_qi_tot(i, j) = MISSING end if !!! SNOW - if (f_qs) then + if (f_qs .and. p_qs > param_first_scalar) then !!! RESOLVED !!! ! Calc effective radius snow q_aux = integrate_1var (rhodz, qs(i, :, j), kms, kme, kts, kte) @@ -252,78 +255,123 @@ SUBROUTINE solar_diag (rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, qv, else tau_qs(i, j) = 0.0 end if - else - re_snow_path(i, j) = MISSING - tau_qs(i, j) = MISSING end if + else + swp2d(i, j) = MISSING + re_snow_path(i, j) = MISSING + tau_qs(i, j) = MISSING end if - if (f_qc .or. f_qi .or. f_qs) then + if ( (f_qc .or. f_qi .or. f_qs) .and. & + (p_qc > param_first_scalar .or. & + p_qi > param_first_scalar .or. & + p_qs > param_first_scalar) ) then !!! RESOLVED !!! ! Sum water paths (cloud liquid + ice + snow) - wp2d_sum(i, j) = lwp2d(i, j) + iwp2d(i, j) + swp2d(i, j) + wp2d_sum(i, j) = MAX( lwp2d(i, j), 0.0 ) + MAX( iwp2d(i, j), 0.0 ) + MAX( swp2d(i, j), 0.0 ) !!! CLOUD BASE & TOP HEIGHTS ! Cloud base first - k = kts - wc = qc(i, k, j) + qi(i, k, j) + qs(i, k, j) - do while ( (wc < WC_MIN) .and. (k < kte) ) - k = k + 1 - wc = qc(i, k, j) + qi(i, k, j) + qs(i, k, j) - end do - - if (k == kte) then - cbase(i, j) = MISSING - else - cbase(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) - end if - - ! Cloud top second - k = kte - wc = qc(i, k, j) + qi(i, k, j) + qs(i, k, j) - do while ( (wc < WC_MIN) .and. (k > kts) ) - k = k - 1 - wc = qc(i, k, j) + qi(i, k, j) + qs(i, k, j) - end do + if (wp2d_sum(i, j) > Q_MIN) then + k = kts + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + do while ( (wc < WC_MIN) .and. (k < kte) ) + k = k + 1 + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + end do + + if (k == kte) then + cbase(i, j) = MISSING + else + cbase(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + end if - if (k == kts) then - ctop(i, j) = MISSING + ! Cloud top second + k = kte + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + do while ( (wc < WC_MIN) .and. (k > kts) ) + k = k - 1 + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + end do + + if (k == kts) then + ctop(i, j) = MISSING + else + ctop(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + end if else - ctop(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + cbase(i, j) = MISSING + ctop(i, j) = MISSING end if !!! TOTAL (RESOLVED + UNRESOLVED) !!! ! Sum water paths (cloud liquid + ice + snow) ! Note that snow is from resolved only - wp2d_tot_sum(i, j) = lwp2d_tot(i, j) + iwp2d_tot(i, j) + swp2d(i, j) + wp2d_tot_sum(i, j) = MAX( lwp2d_tot(i, j), 0.0 ) + MAX( iwp2d_tot(i, j), 0.0 ) + MAX( swp2d(i, j), 0.0 ) ! Cloud base first - k = kts - wc = qc_tot(i, k, j) + qi_tot(i, k, j) + qs(i, k, j) - do while ( (wc < WC_MIN) .and. (k < kte) ) - k = k + 1 - wc = qc_tot(i, k, j) + qi_tot(i, k, j) + qs(i, k, j) - end do - - if (k == kte) then - cbase_tot(i, j) = MISSING - else - cbase_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) - end if - - ! Cloud top second - k = kte - wc = qc_tot(i, k, j) + qi_tot(i, k, j) + qs(i, k, j) - do while ( (wc < WC_MIN) .and. (k > kts) ) - k = k - 1 - wc = qc_tot(i, k, j) + qi_tot(i, k, j) + qs(i, k, j) - end do + if (wp2d_tot_sum(i, j) > Q_MIN) then + k = kts + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + do while ( (wc < WC_MIN) .and. (k < kte) ) + k = k + 1 + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + end do + + if (k == kte) then + cbase_tot(i, j) = MISSING + else + cbase_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + end if - if (k == kts) then - ctop_tot(i, j) = MISSING + ! Cloud top second + k = kte + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + do while ( (wc < WC_MIN) .and. (k > kts) ) + k = k - 1 + wc = 0.0 + if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j) + if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j) + if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j) + end do + + if (k == kts) then + ctop_tot(i, j) = MISSING + else + ctop_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + end if else - ctop_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G) + cbase_tot(i, j) = MISSING + ctop_tot(i, j) = MISSING end if + else + wp2d_sum(i, j) = MISSING + cbase(i, j) = MISSING + ctop(i, j) = MISSING + wp2d_tot_sum(i, j) = MISSING + cbase_tot(i, j) = MISSING + ctop_tot(i, j) = MISSING end if ENDDO diff --git a/phys/module_diagnostics_driver.F b/phys/module_diagnostics_driver.F index 5b10eac981..9b45e20d70 100644 --- a/phys/module_diagnostics_driver.F +++ b/phys/module_diagnostics_driver.F @@ -39,6 +39,7 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & SKIP_PRESS_DIAGS, SKIP_Z_DIAGS, & DO_TRAD_FIELDS, & DO_SOLAR_OUTPUT, & + PARAM_FIRST_SCALAR, & P_QG, P_QH, P_QV, P_QC, P_QI, P_QS, & P_QNG, P_QH, P_QNH, P_QR, P_QNR, & F_QV, F_QC, F_QI, F_QS, & @@ -1084,7 +1085,8 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & CALL solar_diag ( & rho=grid%rho, dz8w=dz8w, ph=grid%ph_2, phb=grid%phb, & cldfrac3d=grid%cldfra, coszen=grid%coszen, swdnb=grid%swdnb, & - swdnt=grid%swdnt, qv=moist(ims,kms,jms,P_QV), & + swdnt=grid%swdnt, param_first_scalar=param_first_scalar, & + p_qc=p_qc, p_qi=p_qi, p_qs=p_qs, qv=moist(ims,kms,jms,P_QV), & qc=moist(ims,kms,jms,P_QC), qi=moist(ims,kms,jms,P_QI), & qs=moist(ims,kms,jms,P_QS), qc_tot=grid%qc_tot, & qi_tot=grid%qi_tot, & diff --git a/phys/module_dust_emis.F b/phys/module_dust_emis.F index f022eb5116..b171ff979f 100644 --- a/phys/module_dust_emis.F +++ b/phys/module_dust_emis.F @@ -142,7 +142,7 @@ subroutine bulk_dust_emis (ktau,dt,num_soil_layers,u_phy,v_phy, & ! increase erodability where the surface albedo is high to account better for real deserts if (erodin .gt. 1.E-8 .AND. albbck(i,j).gt.0.175 .and. vegfra(i,j).lt.12.5) then - erodin = MIN(0.5, erodin + 0.1*albbck(i,j)) + erodin = min(0.5d0, dble(erodin + 0.1*albbck(i,j))) endif ! volumetric soil moisture over porosity @@ -158,14 +158,14 @@ subroutine bulk_dust_emis (ktau,dt,num_soil_layers,u_phy,v_phy, & ! Case of surface dry enough to erode IF (gwet < 0.5) THEN ! Pete's modified value ! IF (gwet < 0.2) THEN - u_ts = MAX(0.0D+0,u_ts0*(1.2D+0+2.0D-1*LOG10(MAX(1.0D-3, gwet)))) + u_ts = max(0.0d0,dble(u_ts0*(1.2d0+2.0d-1*log10(max(1.0d-3, dble(gwet)))))) ELSE ! Case of wet surface, no erosion u_ts = 100.0 END IF srce = frac_s(n)*erodin*dxy ! (m2) ! srce = 1.1*erodin*dxy ! (m2) - dsrc = MAX(0.0, ch_dust(n,month)*srce*w10m*w10m *(w10m - u_ts)*dt) ! (kg) + dsrc = max(0.0d0, dble(ch_dust(n,month)*srce*w10m*w10m *(w10m - u_ts)*dt)) ! (kg) ! unit change from kg/timestep/cell to ug/m2/s ! totalemis=((dsrc)/dt)*converi/dxy diff --git a/phys/module_irrigation.F b/phys/module_irrigation.F index aae8d014af..660ec1b634 100644 --- a/phys/module_irrigation.F +++ b/phys/module_irrigation.F @@ -35,7 +35,8 @@ SUBROUTINE drip_irrigation( julian_in & INTEGER, INTENT(INOUT) :: irr_rand_field_val IRRIGATION_CHANNEL=0. IF(RAINBL.LE.0.01 .AND. IRRIGATION.GE.0.001) THEN - end_hour=irr_start_hour+irr_num_hours + end_hour=irr_start_hour+irr_num_hours + if(end_hour.gt.23) end_hour=end_hour-24 constants_irrigation=irr_freq*irr_daily_amount*0.000277778*0.01/irr_num_hours ! hours in second:1/3600=0.000277778 phase=0. timing=modulo((int(julian_in)-irr_start_julianday),irr_freq) @@ -44,7 +45,9 @@ SUBROUTINE drip_irrigation( julian_in & xt24=mod(xtime,1440.) tloc=floor(gmt+xt24/60.) if(tloc.lt.0) tloc=tloc+24 - IF ((tloc.GE.irr_start_hour .AND. tloc.LT.end_hour) .AND. (julian_in.GE.irr_start_julianday .AND. julian_in.LT.irr_end_julianday) .AND. timing.EQ.0. ) THEN + IF (((tloc.GE.irr_start_hour .AND. tloc.LT.24) .OR. ( tloc.GE.0 .AND. tloc.LT.end_hour) ) & + .AND. ((julian_in.GE.irr_start_julianday .AND. julian_in.LT.367) .OR. & + ( julian_in.GE.0 .AND. julian_in.LT.irr_end_julianday)) .AND. timing.EQ.0. ) THEN IF(irr_ph.EQ.0) THEN RAINBL =RAINBL +dt*IRRIGATION*constants_irrigation IRRIGATION_CHANNEL=0. @@ -81,6 +84,7 @@ SUBROUTINE channel_irrigation( julian_in & !ARI IRRIGATION_CHANNEL=0. IF(RAINBL.LE.0.01 .AND. IRRIGATION.GE.0.001) THEN end_hour=irr_start_hour+irr_num_hours + if(end_hour.gt.23) end_hour=end_hour-24 constants_irrigation=irr_freq*irr_daily_amount*0.000277778*0.01/irr_num_hours ! hours in second:1/3600=0.000277778 phase=0. timing=modulo((int(julian_in)-irr_start_julianday),irr_freq) @@ -89,7 +93,9 @@ SUBROUTINE channel_irrigation( julian_in & !ARI xt24=mod(xtime,1440.) tloc=floor(gmt+xt24/60.) if(tloc.lt.0) tloc=tloc+24 - IF ((tloc.GE.irr_start_hour .AND. tloc.LT.end_hour) .AND. (julian_in.GE.irr_start_julianday .AND. julian_in.LT.irr_end_julianday) .AND. timing.EQ.0. ) THEN + IF (((tloc.GE.irr_start_hour .AND. tloc.LT.24) .OR. ( tloc.GE.0 .AND. tloc.LT.end_hour) ) & + .AND. ((julian_in.GE.irr_start_julianday .AND. julian_in.LT.367) .OR. & + ( julian_in.GE.0 .AND. julian_in.LT.irr_end_julianday)) .AND. timing.EQ.0. ) THEN IF(irr_ph.EQ.0) THEN IRRIGATION_CHANNEL=dt*IRRIGATION*constants_irrigation ELSE @@ -137,7 +143,7 @@ SUBROUTINE sprinkler_irrigation( julian_in qr_curr REAL, INTENT(IN ) :: dt end_hour=irr_start_hour+irr_num_hours - + if(end_hour.gt.23) end_hour=end_hour-24 xt24=mod(xtime,1440.) tloc=floor(gmt+xt24/60.) if(tloc.lt.0) tloc=tloc+24 @@ -147,9 +153,9 @@ SUBROUTINE sprinkler_irrigation( julian_in DO b=jts,jte constants_irrigation=irr_freq*irr_daily_amount/(irr_num_hours*3600*rho(a,kms,b)*dz8w(a,kms,b)*100) - IF ( (tloc.GE.irr_start_hour .AND. tloc.LT.end_hour) .AND. & + IF ( ((tloc.GE.irr_start_hour .AND. tloc.LT.24) .OR. ( tloc.GE.0 .AND. tloc.LT.end_hour) ) .AND. & (irrigation(a,b).GE.0.1) .AND. & - (julian_in.GE.irr_start_julianday .AND. julian_in.LT.irr_end_julianday ) ) THEN + ((julian_in.GE.irr_start_julianday .AND. julian_in.LT.367) .OR. ( julian_in.GE.0 .AND. julian_in.LT.irr_end_julianday)) ) THEN CALL irr_calc_phase(irr_ph,phase,irr_rand_field_val(a,b),a,b,irrigation(a,b),irr_freq) IF(irr_ph.EQ.0) THEN qr_curr(a,kms,b)=qr_curr(a,kms,b)+irrigation(a,b)*constants_irrigation*dt diff --git a/phys/module_mixactivate.F b/phys/module_mixactivate.F index 8bdc8bbc87..a44c0cc279 100644 --- a/phys/module_mixactivate.F +++ b/phys/module_mixactivate.F @@ -413,7 +413,7 @@ subroutine mixactivate( msectional, & (/0.02,0.05,0.1,0.2,0.5,1.0/) real super(psat) ! supersaturation real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m) - real :: ccnfact(psat,maxd_asize, maxd_atype) + real :: ccnfact real :: amcube(maxd_asize, maxd_atype) ! cube of dry mode radius (m) real :: argfactor(maxd_asize, maxd_atype) real aten ! surface tension parameter @@ -432,6 +432,8 @@ subroutine mixactivate( msectional, & real :: colmass_worst( 0:maxd_acomp, maxd_asize, maxd_atype ) real :: colmass_maxworst_r real :: rhodz( kts:kte ), rhodzsum + real :: tmp_amcube, tmp_dpvolmean, tmp_npv, tmp_num_mr + real :: tmp_vol_mr( kts:kte ) !!$#if (defined AIX) !!$#define ERF erf @@ -1163,36 +1165,64 @@ subroutine mixactivate( msectional, & chem(i,kts:kte,j,lnumcw)= raercol(kts:kte,lnumcw,nnew)*scale chem(i,kts:kte,j,lnum)= raercol(kts:kte,lnum,nnew)*scale endif + tmp_vol_mr(kts:kte) = 0.0 do l=1,ncomp(n) lmass=massptr_aer(l,m,n,ai_phase) lmasscw=massptr_aer(l,m,n,cw_phase) -! scale = mwdry/mw_aer(l,n) scale = 1.e9 chem(i,kts:kte,j,lmasscw)=raercol(kts:kte,lmasscw,nnew)*scale ! ug/kg chem(i,kts:kte,j,lmass)=raercol(kts:kte,lmass,nnew)*scale ! ug/kg + tmp_vol_mr(kts:kte) = tmp_vol_mr(kts:kte) + & + (raercol(kts:kte,lmass,nnew) + raercol(kts:kte,lmasscw,nnew))/(1.0e-3*dens_aer(l,n)) + ! (kg_dmap/kg_air)/(kg_dmap/cm3_dvap) = (cm3_dvap/kg_air) + ! note: dmap (or dvap) means dry mass (or volume) of aerosol particles enddo lwater=waterptr_aer(m,n) if(lwater>0)chem(i,kts:kte,j,lwater)=raercol(kts:kte,lwater,nnew) ! don't convert units + + exp45logsig=exp(4.5*alogsig(m,n)*alogsig(m,n)) do k=kts,kte - sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n))) + if (lnum > 0) then + tmp_num_mr = raercol(k,lnum,nnew) + raercol(k,lnumcw,nnew) ! (num_ap/kg_air) + if (tmp_num_mr .lt. 1.0e-14) then ! this is about 1e-20 num_ap/cm3_air + sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n))) + else + ! rce 2020/09/24 - calculate sm using the actual dgnum (that varies in + ! space & time) rather than the default value in dgnum_aer(m,n) + tmp_dpvolmean = (1.90985*tmp_vol_mr(k)/tmp_num_mr)**0.3333333 ! (cm) + tmp_dpvolmean = max( dlo_sect(m,n), min( dhi_sect(m,n), tmp_dpvolmean ) ) + tmp_npv = 6./(pi*((0.01*tmp_dpvolmean)**3)) ! (num_ap/m3_dvap) + tmp_amcube = 3./(4.*pi*exp45logsig*tmp_npv) ! tmp_amcube = (0.5*dgnum)**3 in m3 + sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*tmp_amcube)) + ! sm = critical supersaturation for diameter = dgnum + endif + else + tmp_num_mr = (tmp_vol_mr(k)*1.0e-6)*npv(m,n) ! (num_ap/kg_air) + sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n))) + end if + +! calculate ccn concentrations (num_ap/cm3_air) as diagnostics +! assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles do l=1,psat arg=argfactor(m,n)*log(sm/super(l)) + ! since scrit is proportional to dry_diam**(-3/2) + ! arg = (log(dp_for_super_l) - log(dgnum))/(sqrt(2)*alogsig) + ! where dp_for_super_l is diameter at which scrit = super(l) if(arg<2)then if(arg<-2)then - ccnfact(l,m,n)=1.e-6 ! convert from #/m3 to #/cm3 + ccnfact =1.0 else - ccnfact(l,m,n)=1.e-6*0.5*ERFC_NUM_RECIPES(arg) + ccnfact = 0.5*ERFC_NUM_RECIPES(arg) ! fraction of particles in bin/mode with scrit < super(l) ! fraction of particles in bin/mode with scrit < super(l) endif else - ccnfact(l,m,n) = 0. + ccnfact = 0.0 endif -! ccn concentration as diagnostic -! assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles - ccn(k,l)=ccn(k,l)+(raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew))*cs(k)*ccnfact(l,m,n) - enddo - enddo - enddo - enddo + ccn(k,l) = ccn(k,l) + (tmp_num_mr*ccnfact)*cs(k)*1.0e-6 + enddo ! l + enddo ! k + + enddo ! m + enddo ! n do l=1,psat !wig, 22-Nov-2006: added vertical bounds to prevent out-of-bounds at top if(l.eq.1)ccn1(i,kts:kte,j)=ccn(:,l) diff --git a/phys/module_mp_fast_sbm.F b/phys/module_mp_fast_sbm.F index 233eea27d6..333c093f17 100644 --- a/phys/module_mp_fast_sbm.F +++ b/phys/module_mp_fast_sbm.F @@ -1978,7 +1978,7 @@ SUBROUTINE JERSUPSAT_KS (DEL1,DEL2,DEL1N,DEL2N, & IRI = 1 IPI = 1 - IF(AMAX1(RW,PW,RI,PI)<=RW_PW_RI_PI_MIN) THEN + IF(max(RW,PW,RI,PI)<=RW_PW_RI_PI_MIN) THEN RW = 0.0 IRW = 0 diff --git a/phys/module_mp_milbrandt2mom.F b/phys/module_mp_milbrandt2mom.F index a1304348ad..e7b1a33838 100644 --- a/phys/module_mp_milbrandt2mom.F +++ b/phys/module_mp_milbrandt2mom.F @@ -190,7 +190,7 @@ real FUNCTION gamma(xx) gammadp= exp(gammadp) !!GEM: - gamma = sngl(gammadp) + gamma = gammadp END FUNCTION gamma !======================================================================! @@ -260,7 +260,7 @@ real FUNCTION gammln(xx) enddo !! GEM: - gammln= sngl( tmp+log(stp*ser/x) ) + gammln=tmp+log(stp*ser/x) END FUNCTION gammln !======================================================================! @@ -3051,7 +3051,7 @@ SUBROUTINE mp_milbrandt2mom_main(WZ,T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,PS, LAMr = 1./iLAMr !note: The following coding of 'No_r=...' prevents overflow: !No_r = NR(i,k)*LAMr**(1.+alpha_r))*iGR31 - No_r = sngl(dble(NR(i,k))*dble(LAMr)**dble(1.+alpha_r))*iGR31 + No_r = dble(NR(i,k))*dble(LAMr)**dble(1.+alpha_r)*iGR31 !note: There is an error in MY05a_eq(8) for VENTx (corrected in code) VENTr = Avx*GR32*iLAMr**cexr5 + Bvx*ScTHRD*sqrt(gam*afr*iMUkin)*GR17*iLAMr**cexr6 ABw = CHLC**2/(Ka*RGASV*T(i,k)**2)+1./(DE(i,k)*(QSW(i,k))*Cdiff) diff --git a/phys/module_mp_morr_two_moment.F b/phys/module_mp_morr_two_moment.F index 0c1b077557..94ec70e9fc 100644 --- a/phys/module_mp_morr_two_moment.F +++ b/phys/module_mp_morr_two_moment.F @@ -76,6 +76,8 @@ ! 1) changes and cleanup of code comments ! 2) correction to universal gas constant (very small change) +! CHANGES FOR WRFV4.3 +! 1) fix to saturation vapor pressure polysvp to work at T < -80 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING @@ -1287,6 +1289,7 @@ SUBROUTINE MORR_TWO_MOMENT_MICRO(QC3DTEN,QI3DTEN,QNI3DTEN,QR3DTEN, & ! SATURATION VAPOR PRESSURE AND MIXING RATIO ! hm, add fix for low pressure, 5/12/10 + EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0)) ! PA EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1)) ! PA @@ -4080,11 +4083,18 @@ REAL FUNCTION POLYSVP (T,TYPE) ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* & ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ & ! LOG10(6.1071))*100. +! hm 11/16/20, use Goff-Gratch for T < 195.8 K and Flatau et al. equal or above 195.8 K + if (t.ge.195.8) then + dt=t-273.15 + polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt))))))) + polysvp = polysvp*100. + else + polysvp = 10.**(-9.09718*(273.16/t-1.)-3.56654* & + alog10(273.16/t)+0.876793*(1.-t/273.16)+ & + alog10(6.1071))*100. - dt = max(-80.,t-273.16) - polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt))))))) - polysvp = polysvp*100. + end if END IF @@ -4092,17 +4102,27 @@ REAL FUNCTION POLYSVP (T,TYPE) IF (TYPE.EQ.0) THEN - dt = max(-80.,t-273.16) - polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) - polysvp = polysvp*100. - ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ & ! 5.02808*LOG10(373.16/T)- & ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ & ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ & ! LOG10(1013.246))*100. +! hm 11/16/20, use Goff-Gratch for T < 202.0 K and Flatau et al. equal or above 202.0 K + if (t.ge.202.0) then + dt = t-273.15 + polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) + polysvp = polysvp*100. + else + +! note: uncertain below -70 C, but produces physical values (non-negative) unlike flatau + polysvp = 10.**(-7.90298*(373.16/t-1.)+ & + 5.02808*alog10(373.16/t)- & + 1.3816e-7*(10**(11.344*(1.-t/373.16))-1.)+ & + 8.1328e-3*(10**(-3.49149*(373.16/t-1.))-1.)+ & + alog10(1013.246))*100. + end if - END IF + END IF END FUNCTION POLYSVP diff --git a/phys/module_mp_morr_two_moment_aero.F b/phys/module_mp_morr_two_moment_aero.F index 0a27d2e2d8..71cc3f8009 100644 --- a/phys/module_mp_morr_two_moment_aero.F +++ b/phys/module_mp_morr_two_moment_aero.F @@ -76,6 +76,10 @@ ! 1) changes and cleanup of code comments ! 2) correction to universal gas constant (very small change) +! CHANGES FOR WRFV4.3 +! 1) fix to saturation vapor pressure polysvp to work at T < -80 C +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! TWG 2017 ! TWG = Timothy Glotfelty, EPA ! Adapted from WRFV3.8.1 Morrison Double Moment Scheme @@ -5031,11 +5035,18 @@ REAL FUNCTION POLYSVP (T,TYPE) ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* & ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ & ! LOG10(6.1071))*100. +! hm 11/16/20, use Goff-Gratch for T < 195.8 K and Flatau et al. equal or above 195.8 K + if (t.ge.195.8) then + dt=t-273.15 + polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt))))))) + polysvp = polysvp*100. + else + polysvp = 10.**(-9.09718*(273.16/t-1.)-3.56654* & + alog10(273.16/t)+0.876793*(1.-t/273.16)+ & + alog10(6.1071))*100. - dt = max(-80.,t-273.16) - polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt))))))) - polysvp = polysvp*100. + end if END IF @@ -5043,17 +5054,27 @@ REAL FUNCTION POLYSVP (T,TYPE) IF (TYPE.EQ.0) THEN - dt = max(-80.,t-273.16) - polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) - polysvp = polysvp*100. - ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ & ! 5.02808*LOG10(373.16/T)- & ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ & ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ & ! LOG10(1013.246))*100. +! hm 11/16/20, use Goff-Gratch for T < 202.0 K and Flatau et al. equal or above 202.0 K + if (t.ge.202.0) then + dt = t-273.15 + polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) + polysvp = polysvp*100. + else + +! note: uncertain below -70 C, but produces physical values (non-negative) unlike flatau + polysvp = 10.**(-7.90298*(373.16/t-1.)+ & + 5.02808*alog10(373.16/t)- & + 1.3816e-7*(10**(11.344*(1.-t/373.16))-1.)+ & + 8.1328e-3*(10**(-3.49149*(373.16/t-1.))-1.)+ & + alog10(1013.246))*100. + end if - END IF + END IF END FUNCTION POLYSVP diff --git a/phys/module_mp_p3.F b/phys/module_mp_p3.F index f0b7c6d6b4..b6968444c7 100644 --- a/phys/module_mp_p3.F +++ b/phys/module_mp_p3.F @@ -354,7 +354,7 @@ subroutine p3_init(lookup_file_dir,nCat,model,stat,abort_on_err) if (global_status == STATUS_ERROR) then if (err_abort) then print*,'Stopping in P3 init' - call flush(6) + flush(6) stop endif return @@ -430,7 +430,7 @@ subroutine p3_init(lookup_file_dir,nCat,model,stat,abort_on_err) if (global_status == STATUS_ERROR) then if (err_abort) then print*,'Stopping in P3 init' - call flush(6) + flush(6) stop endif return diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index fbc32a7420..4fb4e9ddeb 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -49,6 +49,7 @@ MODULE module_mp_thompson #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) USE module_dm, ONLY : wrf_dm_max_real #endif + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG IMPLICIT NONE @@ -1304,16 +1305,16 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte - re_qc1d(k) = 2.49E-6 - re_qi1d(k) = 4.99E-6 - re_qs1d(k) = 9.99E-6 + re_qc1d(k) = RE_QC_BG + re_qi1d(k) = RE_QI_BG + re_qs1d(k) = RE_QS_BG enddo call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & re_qc1d, re_qi1d, re_qs1d, kts, kte) do k = kts, kte - re_cloud(i,k,j) = MAX(2.49E-6, MIN(re_qc1d(k), 50.E-6)) - re_ice(i,k,j) = MAX(4.99E-6, MIN(re_qi1d(k), 125.E-6)) - re_snow(i,k,j) = MAX(9.99E-6, MIN(re_qs1d(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc1d(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi1d(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs1d(k), 999.E-6)) enddo ENDIF @@ -3008,7 +3009,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt) prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k)) pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt), & - -tnc_wev(idx_d, idx_c, idx_n)*orho*odt) + DBLE(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt)) endif else @@ -5065,7 +5066,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qc) then do k = kts, kte - re_qc1d(k) = 2.49E-6 + re_qc1d(k) = RE_QC_BG if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE if (nc(k).lt.100) then inu_c = 15 @@ -5081,16 +5082,16 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qi) then do k = kts, kte - re_qi1d(k) = 2.49E-6 + re_qi1d(k) = RE_QI_BG if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi - re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) + re_qi1d(k) = MAX(5.01E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) enddo endif if (has_qs) then do k = kts, kte - re_qs1d(k) = 4.99E-6 + re_qs1d(k) = RE_QS_BG if (rs(k).le.R1) CYCLE tc0 = MIN(-0.1, t1d(k)-273.15) smob = rs(k)*oams @@ -5125,7 +5126,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smo2**b_ - re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) + re_qs1d(k) = MAX(10.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) enddo endif diff --git a/phys/module_mp_wdm5.F b/phys/module_mp_wdm5.F index caafdfd77d..6e69b082e2 100644 --- a/phys/module_mp_wdm5.F +++ b/phys/module_mp_wdm5.F @@ -10,6 +10,7 @@ MODULE module_mp_wdm5 ! USE module_mp_radar + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -289,9 +290,9 @@ SUBROUTINE wdm5(th, q, qc, qr, qi, qs & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i = its, ite do k = kts, kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -304,9 +305,9 @@ SUBROUTINE wdm5(th, q, qc, qr, qi, qs & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) do k = kts, kte - re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6)) - re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6)) - re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6)) enddo ! k loop enddo ! i loop endif ! has_reqc, etc... diff --git a/phys/module_mp_wdm6.F b/phys/module_mp_wdm6.F index 039ebfa64d..4ff157910c 100644 --- a/phys/module_mp_wdm6.F +++ b/phys/module_mp_wdm6.F @@ -10,6 +10,7 @@ module module_mp_wdm6 !------------------------------------------------------------------------------- ! use module_mp_radar + use module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! ! !------------------------------------------------------------------------------- @@ -301,9 +302,9 @@ subroutine wdm6(th, q, qc, qr, qi, qs, qg, & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i = its,ite do k = kts,kte - re_qc(k) = 2.51e-6 - re_qi(k) = 10.01e-6 - re_qs(k) = 25.e-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG ! t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -318,9 +319,9 @@ subroutine wdm6(th, q, qc, qr, qi, qs, qg, & kts, kte, i, j) ! do k = kts,kte - re_cloud(i,k,j) = max(2.51e-6, min(re_qc(k), 50.e-6)) - re_ice(i,k,j) = max(10.01e-6, min(re_qi(k), 125.e-6)) - re_snow(i,k,j) = max(25.e-6, min(re_qs(k), 999.e-6)) + re_cloud(i,k,j) = max(RE_QC_BG, min(re_qc(k), 50.e-6)) + re_ice(i,k,j) = max(RE_QI_BG, min(re_qi(k), 125.e-6)) + re_snow(i,k,j) = max(RE_QS_BG, min(re_qs(k), 999.e-6)) enddo enddo endif diff --git a/phys/module_mp_wdm7.F b/phys/module_mp_wdm7.F index edb97b5e3d..e063038d78 100644 --- a/phys/module_mp_wdm7.F +++ b/phys/module_mp_wdm7.F @@ -9,6 +9,7 @@ MODULE module_mp_wdm7 ! ! USE module_mp_radar + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -334,9 +335,9 @@ SUBROUTINE wdm7(th, q, qc, qr, qi, qs, qg, qh, & IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN DO i=its,ite DO k=kts,kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -349,9 +350,9 @@ SUBROUTINE wdm7(th, q, qc, qr, qi, qs, qg, qh, & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) DO k=kts,kte - re_cloud(i,k,j) = max(2.51E-6, min(re_qc(k), 50.E-6)) - re_ice(i,k,j) = max(10.01E-6, min(re_qi(k), 125.E-6)) - re_snow(i,k,j) = max(25.E-6, min(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = max(RE_QC_BG, min(re_qc(k), 50.E-6)) + re_ice(i,k,j) = max(RE_QI_BG, min(re_qi(k), 125.E-6)) + re_snow(i,k,j) = max(RE_QS_BG, min(re_qs(k), 999.E-6)) ENDDO ENDDO ENDIF diff --git a/phys/module_mp_wsm3.F b/phys/module_mp_wsm3.F index 1d38089bb6..a71d3cbfa6 100644 --- a/phys/module_mp_wsm3.F +++ b/phys/module_mp_wsm3.F @@ -11,6 +11,7 @@ MODULE module_mp_wsm3 ! + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -162,9 +163,9 @@ SUBROUTINE wsm3(th, q, qci, qrs & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i=its,ite do k=kts,kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -182,9 +183,9 @@ SUBROUTINE wsm3(th, q, qci, qrs & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) do k=kts,kte - re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6)) - re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6)) - re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6)) enddo enddo endif ! has_reqc, etc... diff --git a/phys/module_mp_wsm5.F b/phys/module_mp_wsm5.F index faf3982b26..e081a7b6e4 100644 --- a/phys/module_mp_wsm5.F +++ b/phys/module_mp_wsm5.F @@ -13,6 +13,7 @@ MODULE module_mp_wsm5 ! USE module_mp_radar + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -218,9 +219,9 @@ SUBROUTINE wsm5(th, q, qc, qr, qi, qs & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i=its,ite do k=kts,kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -232,9 +233,9 @@ SUBROUTINE wsm5(th, q, qc, qr, qi, qs & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) do k=kts,kte - re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6)) - re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6)) - re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6)) enddo enddo endif ! has_reqc, etc... diff --git a/phys/module_mp_wsm6.F b/phys/module_mp_wsm6.F index 5c52d40f28..69a2d49592 100644 --- a/phys/module_mp_wsm6.F +++ b/phys/module_mp_wsm6.F @@ -9,6 +9,7 @@ MODULE module_mp_wsm6 ! USE module_mp_radar + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -237,9 +238,9 @@ SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i=its,ite do k=kts,kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -251,9 +252,9 @@ SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg & qmin, t0c, re_qc, re_qi, re_qs, & kts, kte, i, j) do k=kts,kte - re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6)) - re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6)) - re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6)) enddo enddo endif ! has_reqc, etc... diff --git a/phys/module_mp_wsm7.F b/phys/module_mp_wsm7.F index 31e5636559..8f757f08c5 100644 --- a/phys/module_mp_wsm7.F +++ b/phys/module_mp_wsm7.F @@ -9,6 +9,7 @@ MODULE module_mp_wsm7 ! USE module_mp_radar + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG ! REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain @@ -245,9 +246,9 @@ SUBROUTINE wsm7(th, q, qc, qr, qi, qs, qg, qh & if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then do i=its,ite do k=kts,kte - re_qc(k) = 2.51E-6 - re_qi(k) = 10.01E-6 - re_qs(k) = 25.E-6 + re_qc(k) = RE_QC_BG + re_qi(k) = RE_QI_BG + re_qs(k) = RE_QS_BG ! t1d(k) = th(i,k,j)*pii(i,k,j) den1d(k)= den(i,k,j) @@ -261,9 +262,9 @@ SUBROUTINE wsm7(th, q, qc, qr, qi, qs, qg, qh & kts, kte, i, j) ! do k=kts,kte - re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6)) - re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6)) - re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6)) + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6)) + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6)) enddo enddo endif ! if has_reqc ne 0 end diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index d9edbf6c53..f5bf1b4c68 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -942,15 +942,22 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & ! for cloud, ice, and snow. IF ( config_flags%swint_opt .eq. 2 ) THEN - IF ( ( has_reqc .EQ. 1 ) .AND. & - ( has_reqi .EQ. 1 ) .AND. & - ( has_reqs .EQ. 1 ) .AND. & + IF (( config_flags%mp_physics == THOMPSON .OR. & + config_flags%mp_physics == THOMPSONAERO .OR. & + config_flags%mp_physics == WSM3SCHEME .OR. & + config_flags%mp_physics == WSM5SCHEME .OR. & + config_flags%mp_physics == WSM6SCHEME .OR. & + config_flags%mp_physics == WSM7SCHEME .OR. & + config_flags%mp_physics == WDM5SCHEME .OR. & + config_flags%mp_physics == WDM6SCHEME .OR. & + config_flags%mp_physics == WDM7SCHEME ).OR. & + (( has_reqc .EQ. 0 .AND. has_reqi .EQ. 0 .and. has_reqs .EQ. 0 ) .AND. & ( f_qc ) .AND. & ( f_qi ) .AND. & - ( f_qs ) ) THEN + ( f_qs ))) THEN ! everything is A-OK for FARMS ELSE - CALL wrf_error_fatal ('--- ERROR: FARMS (swint_opt==2) requires a different MP scheme') + CALL wrf_error_fatal ('--- ERROR: FARMS (swint_opt==2) requires a different MP scheme (Please see the module_physics_init') END IF END IF @@ -1163,13 +1170,13 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & end do !..Fill initial starting values of radiative effective radii for -!.. cloud water (2.51 microns), cloud ice (5.01 microns), and -!.. snow (10.01 microns). +!.. cloud water (2.49 microns), cloud ice (4.99 microns), and +!.. snow (9.99 microns). if (has_reqc.ne.0) then do j=jts,jtf do k=kts,ktf do i=its,itf - re_cloud(i,k,j) = 2.51E-6 + re_cloud(i,k,j) = RE_QC_BG end do end do end do @@ -1178,7 +1185,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & do j=jts,jtf do k=kts,ktf do i=its,itf - re_ice(i,k,j) = 5.01E-6 + re_ice(i,k,j) = RE_QI_BG end do end do end do @@ -1187,7 +1194,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & do j=jts,jtf do k=kts,ktf do i=its,itf - re_snow(i,k,j) = 10.01E-6 + re_snow(i,k,j) = RE_QS_BG end do end do end do diff --git a/phys/module_ra_farms.F b/phys/module_ra_farms.F index d63926c19a..0db1e92319 100644 --- a/phys/module_ra_farms.F +++ b/phys/module_ra_farms.F @@ -27,26 +27,29 @@ module module_ra_farms real, parameter :: DE_CLOUD_MIN = 5.0, DE_CLOUD_MAX = 120.0 real, parameter :: TAU_MIN = 0.0001, TAU_MAX = 300.0 real, parameter :: AOD550_VAL = 0.12, ANGEXP_VAL = 1.3, AERSSA_VAL = 0.85, AERASY_VAL = 0.9 + real, parameter :: RE_CLOUD_CLIM = 8.E-6, RE_ICE_CLIM = 24.E-6, RE_SNOW_CLIM = 24.E-6 contains subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, kte, & p8w, rho, dz8w, albedo, aer_opt, aerssa2d, aerasy2d, aod5502d, angexp2d, & coszen_loc, qv, qi, qs, qc, re_cloud, re_ice, re_snow, & - julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2) + julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2, & + has_reqc, has_reqi, has_reqs, cldfra) implicit None integer, intent(in) :: ims, ime, jms, jme, its, ite, jts, jte, kms, kme, & kts, kte - integer, intent(in) :: aer_opt + integer, intent(in) :: aer_opt, has_reqc, has_reqi, has_reqs real, intent(in) :: julian real, dimension(ims:ime, jms:jme), intent(in) :: albedo, coszen_loc real, dimension(ims:ime, kms:kme, jms:jme ), intent(in) :: qv, qi, qs, qc, & - p8w, rho, dz8w, re_cloud, re_ice, re_snow + p8w, rho, dz8w, cldfra + real, dimension(ims:ime, kms:kme, jms:jme ), intent(in) :: re_cloud, re_ice, re_snow real, dimension(ims:ime, jms:jme), intent(inout) :: aerssa2d, aerasy2d, aod5502d, angexp2d real, dimension(ims:ime,jms:jme), intent(inout) :: swddir2, swdown2, & @@ -54,10 +57,10 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, swdownc2, swddnic2 ! Local - integer :: i, j + integer :: i, j, k real :: tau_qv, tau_qi, tau_qs, pmw, swp, iwp, lwp, beta - real :: re_cloud_path, re_ice_path, re_snow_path, q_aux - real, dimension(kms:kme) :: rhodz + real :: re_cloud_path, re_ice_path, re_snow_path, q_aux, cldfra_h + real, dimension(kms:kme) :: rhodz, re_cloud_k, re_ice_k, re_snow_k j_loop: do j = jts, jte @@ -71,6 +74,34 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, swddnic2(i, j) = 0.0 else rhodz(:) = rho(i, :, j) * dz8w(i, :, j) / (1. + qv(i, :, j)) + re_cloud_k(:) = re_cloud(i, :, j) + re_ice_k(:) = re_ice(i, :, j) + re_snow_k(:) = re_snow(i, :, j) + + if (has_reqc == 1) then + do k = kts, kte + if (CLDFRA (i, k, j) > 0.0 .and. re_cloud_k(k) < 2.5E-6) re_cloud_k(k) = RE_CLOUD_CLIM + end do + else + re_cloud_k(:) = RE_CLOUD_CLIM + end if + + if (has_reqi == 1) then + do k = kts, kte + if (cldfra(i, k, j) > 0.0 .and. re_ice_k(k) < 5.0E-6) re_ice_k(k) = RE_ICE_CLIM + end do + else + re_ice_k(:) = RE_ICE_CLIM + end if + + if (has_reqs == 1) then + do k = kts, kte + if (cldfra(i, k, j) > 0.0 .and. re_snow_k(k) < 10.0E-6) re_snow_k(k) = RE_SNOW_CLIM + end do + else + re_snow_k(:) = RE_SNOW_CLIM + end if + ! PMW pmw = integrate_1var (rhodz, qv(i, :, j), kms, kme, kts, kte) @@ -80,7 +111,7 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, if (q_aux > 0.0) then re_cloud_path = integrate_2var (rhodz, qc(i, :, j), & - re_cloud(i, :, j), kms, kme, kts, kte) + re_cloud_k, kms, kme, kts, kte) re_cloud_path = re_cloud_path / q_aux else re_cloud_path = 0.0 @@ -92,7 +123,7 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, if (q_aux > 0.0) then re_ice_path = integrate_2var (rhodz, qi(i, :, j), & - re_ice(i, :, j), kms, kme, kts, kte) + re_ice_k, kms, kme, kts, kte) re_ice_path = re_ice_path / q_aux else re_ice_path = 0.0 @@ -104,12 +135,22 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, if (q_aux > 0.0) then re_snow_path = integrate_2var (rhodz, qs(i, :, j), & - re_snow(i, :, j), kms, kme, kts, kte) + re_snow_k, kms, kme, kts, kte) re_snow_path = re_snow_path / q_aux else re_snow_path = 0.0 end if + ! Calculate horizontal cloud fraction + q_aux = integrate_1var (rhodz, qc(i, :, j) + qi(i, :, j) + qs(i, :, j), kms, kme, kts, kte) + if (q_aux > 0.0) then + cldfra_h = integrate_2var (rhodz, qc(i, :, j) + qi(i, :, j) + qs(i, :, j), & + cldfra(i, :, j), kms, kme, kts, kte) + cldfra_h = cldfra_h / q_aux + else + cldfra_h = 0.0 + end if + ! optical thickness water if (re_cloud_path > 0.0) then tau_qv = THREE_OVER_TWO * lwp / re_cloud_path / 1000.0 @@ -154,11 +195,11 @@ subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, end if beta = aod5502d(i, j) * (1000.0/ 550.0) ** (- angexp2d(i, j)) - - Call Farms (p8w(i, 1, j), albedo(i, j), aerssa2d(i, j), & - aerasy2d(i, j), coszen_loc(i, j), beta, & - angexp2d(i, j), pmw, tau_qv, tau_qi, tau_qs, & - re_cloud_path, re_ice_path, re_snow_path, int(julian), & + + Call Farms (p8w(i, 1, j), albedo(i, j), aerssa2d(i, j), & + aerasy2d(i, j), coszen_loc(i, j), beta, & + angexp2d(i, j), pmw, tau_qv, tau_qi, tau_qs, cldfra_h, & + re_cloud_path, re_ice_path, re_snow_path, int(julian), & swdown2(i, j), swddni2(i, j), swddif2(i, j), swddir2(i, j), & swdownc2(i, j), swddnic2(i, j)) @@ -210,7 +251,7 @@ end function Integrate_2var subroutine FARMS (p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, & - tau_qv, tau_qi, tau_qs, re_cloud_path_m, re_ice_path_m, re_snow_path_m, & + tau_qv, tau_qi, tau_qs, cldfra_h, re_cloud_path_m, re_ice_path_m, re_snow_path_m, & juday, ghi, dni, dif, dir, ghi_clear, dni_clear) !!!!!! This Fast All-sky Radiation Model for Solar applications (FARMS) was developed by @@ -248,7 +289,7 @@ subroutine FARMS (p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, & implicit none real, intent(in) :: p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, & - tau_qv, tau_qi, tau_qs, re_cloud_path_m, re_ice_path_m, re_snow_path_m + tau_qv, tau_qi, tau_qs, re_cloud_path_m, re_ice_path_m, re_snow_path_m, cldfra_h integer, intent(in) :: juday real, intent(out) :: ghi, dni, dir, dif, ghi_clear, dni_clear @@ -355,6 +396,9 @@ subroutine FARMS (p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, & ghi = Ftotal ghi_clear = solarangle * F0 * ((Tddclr + Tduclr) / (1.0 - albdo * Ruuclr)) + ghi = cldfra_h * ghi + (1.0 - cldfra_h) * ghi_clear + dni = cldfra_h * dni + (1.0 - cldfra_h) * dni_clear + dif = ghi - dir end subroutine farms diff --git a/phys/module_ra_goddard.F b/phys/module_ra_goddard.F index f4ebeb3204..54d0d12d3e 100644 --- a/phys/module_ra_goddard.F +++ b/phys/module_ra_goddard.F @@ -5254,7 +5254,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, & !dir$ vector aligned DO ic=1,irestrict if(lmask(ic) .eqv. .true.) then - df_sub(ic,k) = max(df(ic,k) - df(ic,k-1), 0.) !df for each layer (remove negative df_sub) + df_sub(ic,k) = max(df(ic,k) - df(ic,k-1), 0._fp_kind) !df for each layer (remove negative df_sub) endif ENDDO enddo @@ -5275,7 +5275,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, & !dir$ vector aligned DO ic=1,irestrict if(lmask(ic) .eqv. .true.) then - flc(ic,k)=max(flc(ic,k)-df_clr(ic,k),0.) !this filter is for small cosine zenith angle. + flc(ic,k)=max(flc(ic,k)-df_clr(ic,k),0._fp_kind) !this filter is for small cosine zenith angle. endif ENDDO enddo @@ -5316,7 +5316,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, & i_tau = min(max(int(cld_alb*10.)+1,1),10) !1~10 ratio = ratio_lut(i_tau,i_cos) else !use computed clear and cloudy flux ratio (not fast_overcast) - ratio = max(0.01, min(1.,(flx(ic,k)/flc(ic,k)))) + ratio = max(0.01_fp_kind, min(1._fp_kind,(flx(ic,k)/flc(ic,k)))) endif df_sub(ic,k) = df_sub(ic,k)*ratio !compute cloudy-sky df_sub enddo !k @@ -5340,8 +5340,8 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, & !dir$ vector aligned DO ic=1,irestrict if(lmask(ic) .eqv. .true.) then - flx(ic,k) = max(flx(ic,k)-df_cld(ic,k) , 0.) !this max is for small cosz - flxd(ic,k) = max(flxd(ic,k)-df_cld(ic,k), 0.) !this max is for small cosz + flx(ic,k) = max(flx(ic,k)-df_cld(ic,k) , 0._fp_kind) !this max is for small cosz + flxd(ic,k) = max(flxd(ic,k)-df_cld(ic,k), 0._fp_kind) !this max is for small cosz flxu(ic,k) = flx(ic,k)-flxd(ic,k) endif ENDDO @@ -5361,7 +5361,7 @@ subroutine swrad ( np,cosz, pl,ta,wa,oa, fcld,ict,icb, & if ( (fdirir(ic)-df_cld(ic,np+1)) >= 0. ) then ! normal fdirir(ic)=fdirir(ic)-df_cld(ic,np+1) ! updated else ! if negative, it also reduces diffuse component. - fdifir(ic) = MAX(0., fdifir(ic) + (fdirir(ic)-df_cld(ic,np+1)) ) + fdifir(ic) = max(0._fp_kind, fdifir(ic) + (fdirir(ic)-df_cld(ic,np+1)) ) fdirir(ic)=0. endif ! @@ -5717,7 +5717,7 @@ subroutine sw_uvpar (np,wh,oh,dp, & !-----for direct incident radiation ! the effective layer optical properties. eqs. (6.2)-(6.4) !modify max & min - tautob(ic)=tausto(ic) + max(tauclb(ic,k,ib),0.e0) + tautob(ic)=tausto(ic) + max(tauclb(ic,k,ib),0._fp_kind) ssatob(ic)=min(max((ssatau(ic)+ssacl(ic,k,ib)*tauclb(ic,k,ib))/tautob(ic) , ssa_min), ssa_max) !SSA of cloud is unity asytob(ic)=min(max((asysto(ic)+asycl(ic,k,ib)*ssacl(ic,k,ib)*tauclb(ic,k,ib)) & /( max(ssatob(ic)*tautob(ic),const_tiny) ),asy_min), asy_max) @@ -6176,7 +6176,7 @@ subroutine sw_ir (np,wh,dp, & !-----compute reflectance and transmittance of the cloudy portion of a layer !-----for direct incident radiation. eqs.(6.2)-(6.4) - tautob(ic,k)=tausto(ic,k)+max(tauclb(ic,k),0.e0) + tautob(ic,k)=tausto(ic,k)+max(tauclb(ic,k),0._fp_kind) ssatob(ic,k)=min(max((ssatau(ic,k)+ssacl(ic,k,iv)*tauclb(ic,k))/tautob(ic,k),ssa_min),ssa_max) asytob(ic,k)=min(max((asysto(ic,k)+asycl(ic,k,iv)*ssacl(ic,k,iv)*tauclb(ic,k)) & /( max(ssatob(ic,k)*tautob(ic,k),const_tiny) ),ssa_min),ssa_max) @@ -6187,7 +6187,7 @@ subroutine sw_ir (np,wh,dp, & ssatof(ic,k)=ssatob(ic,k) asytof(ic,k)=asytob(ic,k) else - tautof(ic,k)=tausto(ic,k)+max(tauclf(ic,k),0.e0) + tautof(ic,k)=tausto(ic,k)+max(tauclf(ic,k),0._fp_kind) ssatof(ic,k)=min(max((ssatau(ic,k)+ssacl(ic,k,iv)*tauclf(ic,k))/tautof(ic,k),ssa_min),ssa_max) asytof(ic,k)=min(max((asysto(ic,k)+asycl(ic,k,iv)*ssacl(ic,k,iv)*tauclf(ic,k)) & /( max(ssatof(ic,k)*tautof(ic,k),const_tiny) ),asy_min),asy_max) @@ -6529,7 +6529,7 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, & !-----normalize cloud cover following eq. (7.8) fa=fcld(ic,k)/cc(ic,kk) !-----table look-up - taux=min(taux,32.) + taux=min(taux,32._fp_kind) fm=cosz(ic)/dm ft=(log10(taux)-t1)/dt fa=fa/da @@ -6555,8 +6555,8 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, & xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ & caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa) xai= xai-2.*caib(im,it,ia) - xai=max(xai,0.0) - xai=min(xai,1.0) + xai=max(xai,0._fp_kind) + xai=min(xai,1._fp_kind) tauclb(ic,k) = taux*xai !-----scale cloud optical thickness for diffuse radiation following eq. (7.4) ! the scaling factor, xai, is a function of the cloud optical @@ -6566,8 +6566,8 @@ subroutine cloud_scale (np,cosz,fcld,taucld,ict,icb, & xai=xai+(-caif(it,ia-1)*(1.-fa)+ & caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa) xai= xai-caif(it,ia) - xai=max(xai,0.0) - xai=min(xai,1.0) + xai=max(xai,0._fp_kind) + xai=min(xai,1._fp_kind) tauclf(ic,k) = taux*xai endif endif @@ -10427,11 +10427,11 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,& !new dp (ic,k) = pl(ic,k+1)-pl(ic,k) dh2o (ic,k) = 1.02*wa(ic,k)*dp(ic,k) - dh2o (ic,k) = max(dh2o (ic,k),1.e-30) + dh2o (ic,k) = max(dh2o (ic,k),1.e-30_fp_kind) do3 (ic,k) = 476.*oa(ic,k)*dp(ic,k) - do3 (ic,k) = max(do3 (ic,k),1.e-30) + do3 (ic,k) = max(do3 (ic,k),1.e-30_fp_kind) dco2 (ic,k) = 789.*co2*dp(ic,k) - dco2 (ic,k) = max(dco2 (ic,k),1.e-30) + dco2 (ic,k) = max(dco2 (ic,k),1.e-30_fp_kind) dch4 (ic,k) = 789.*ch4*dp(ic,k) dn2o (ic,k) = 789.*n2o*dp(ic,k) df11 (ic,k) = 789.*cfc11*dp(ic,k) @@ -10576,7 +10576,7 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,& taux=taux*(1.-ww*ff) !-----compute cloud diffuse transmittance. it is approximated by using ! a diffusivity factor of 1.66. - tauxa=max(0.,1.66*taux) + tauxa=max(0._fp_kind,1.66*taux) tcldlyr(ic,k)=0. if(tauxa.lt.80.) tcldlyr(ic,k)=exp(-tauxa) endif @@ -10885,8 +10885,8 @@ subroutine lwrad ( np, emiss, tb, ts, ict, icb,& if (k2 .eq. k1+1) then !dir$ vector aligned DO ic=1,irestrict - yy=min(0.999,trant(ic)) - yy=max(0.001,yy) + yy=min(0.999_fp_kind,trant(ic)) + yy=max(0.001_fp_kind,yy) !-hmhj use log instead of alog for default intrinsic function xx=(blevel(ic,k1)-blevel(ic,k2))/ log(yy) bu=(blevel(ic,k1)-blevel(ic,k2)*yy)/(1.0-yy)+xx @@ -11652,9 +11652,9 @@ subroutine tablup(k1,k2,np,nx,nh,sabs,spre,stem,w1,p1, & !-----normalize we and pe pe=(log10(x2)-p1)/dpe !-----restrict the magnitudes of the normalized we and pe. - we=min(we,REAL(nh-1)) - pe=max(pe,0.0) - pe=min(pe,REAL(nx-1)) + we=min(we,real(nh-1,kind=fp_kind)) + pe=max(pe,0._fp_kind) + pe=min(pe,real(nx-1,kind=fp_kind)) !-----assign iw and ip and compute the distance of we and pe ! from iw and ip. iw=int(we+1.0) @@ -11684,8 +11684,8 @@ subroutine tablup(k1,k2,np,nx,nh,sabs,spre,stem,w1,p1, & t2 = ca*(1.-fw) + cb*fw !-----update the total transmittance between levels k1 and k2 tran(ic)= (ax + (t1+t2*x3) * x3)*tran(ic) - tran(ic)=min(tran(ic),0.9999999) - tran(ic)=max(tran(ic),0.0000001) + tran(ic)=min(tran(ic),0.9999999_fp_kind) + tran(ic)=max(tran(ic),0.0000001_fp_kind) else tran(ic)=0.9999999 endif diff --git a/phys/module_ra_rrtmg_swk.F b/phys/module_ra_rrtmg_swk.F index 46e5fbc694..a18146b9dc 100644 --- a/phys/module_ra_rrtmg_swk.F +++ b/phys/module_ra_rrtmg_swk.F @@ -1732,7 +1732,7 @@ subroutine cldprmc_sw(nlayers, inflag, iceflag, liqflag, cldfmc, & 'ERROR: SNOW GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS' factor = (radsno - 2._rb)/3._rb index = int(factor) - if(index.eq.167) index = 166 + if(index.eq.46) index = 45 fint = factor-real(index) ib = ngb(ig) extcosno(ig) = extice3(index,ib)+fint* & @@ -4308,7 +4308,7 @@ subroutine swdatinit(cpdair) 7700._rb, 8050._rb,12850._rb,16000._rb,22650._rb,29000._rb, & 38000._rb, 820._rb/) wavenum2(:) = (/3250._rb, 4000._rb, 4650._rb, 5150._rb, 6150._rb, 7700._rb, & - 8050._rb, 2850._rb,16000._rb,22650._rb,29000._rb,38000._rb, & + 8050._rb,12850._rb,16000._rb,22650._rb,29000._rb,38000._rb, & 50000._rb, 2600._rb/) delwave(:) = (/ 650._rb, 750._rb, 650._rb, 500._rb, 1000._rb, 1550._rb, & 350._rb, 4800._rb, 3150._rb, 6650._rb, 6350._rb, 9000._rb, & diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index b8a6b06820..5df9182e9f 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -1984,10 +1984,10 @@ SUBROUTINE radiation_driver ( & lwupflx, lwupflxc, lwdnflx, lwdnflxc, & swupflx, swupflxc, swdnflx, swdnflxc, & lwupt, lwuptc, lwdnt, lwdntc, & - lwupb, lwupbc, lwdnb, lwdnb, & + lwupb, lwupbc, lwdnb, lwdnbc, & glw, olr, lwcf, & swupt, swuptc, swdnt, swdntc, & - swupb, swupbc, swdnb, swdnb, & + swupb, swupbc, swdnb, swdnbc, & gsw, swcf, & coszen, solcon, albedo, emiss, & t,t8w, tsk, rho, p, p8w, cldfra, & @@ -2797,8 +2797,9 @@ SUBROUTINE radiation_driver ( & jte = j_end(ij) call farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, kte, & p8w, rho, dz8w, albedo, aer_opt, aerssa2d, aerasy2d, aod5502d, angexp2d, & - coszen_loc, qv, qi, qs, qc, re_cloud, re_ice, re_snow, & - julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2) + coszen_loc, qv, qi, qs, qc, re_cloud, re_ice, re_snow, & + julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2, & + has_reqc, has_reqi, has_reqs, CLDFRA) enddo !$OMP END PARALLEL DO end if @@ -3247,6 +3248,7 @@ SUBROUTINE calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, & xxlat=xlat(i,j)*degrad coszen(i,j)=sin(xxlat)*sin(declin) & +cos(xxlat)*cos(declin) *cos(hrang(i,j)) + coszen(i, j) = min (max (coszen(i, j), -1.0), 1.0) enddo enddo END SUBROUTINE calc_coszen diff --git a/phys/module_sf_noahmpdrv.F b/phys/module_sf_noahmpdrv.F index f77b959845..f2d0f67467 100644 --- a/phys/module_sf_noahmpdrv.F +++ b/phys/module_sf_noahmpdrv.F @@ -1161,6 +1161,17 @@ SUBROUTINE TRANSFER_MP_PARAMETERS(VEGTYPE,SOILTYPE,SLOPETYPE,SOILCOLOR,CROPTYPE, parameters%GDDS3 = GDDS3_TABLE(CROPTYPE) ! GDD from seeding to post vegetative parameters%GDDS4 = GDDS4_TABLE(CROPTYPE) ! GDD from seeding to intial reproductive parameters%GDDS5 = GDDS5_TABLE(CROPTYPE) ! GDD from seeding to pysical maturity + parameters%C3PSN = C3PSNI_TABLE(CROPTYPE) ! parameters from stomata ! Zhe Zhang 2020-07-13 + parameters%KC25 = KC25I_TABLE(CROPTYPE) + parameters%AKC = AKCI_TABLE(CROPTYPE) + parameters%KO25 = KO25I_TABLE(CROPTYPE) + parameters%AKO = AKOI_TABLE(CROPTYPE) + parameters%AVCMX = AVCMXI_TABLE(CROPTYPE) + parameters%VCMX25 = VCMX25I_TABLE(CROPTYPE) + parameters%BP = BPI_TABLE(CROPTYPE) + parameters%MP = MPI_TABLE(CROPTYPE) + parameters%FOLNMX = FOLNMXI_TABLE(CROPTYPE) + parameters%QE25 = QE25I_TABLE(CROPTYPE) ! ends here parameters%C3C4 = C3C4_TABLE(CROPTYPE) ! photosynthetic pathway: 1. = c3 2. = c4 parameters%AREF = AREF_TABLE(CROPTYPE) ! reference maximum CO2 assimulation rate parameters%PSNRF = PSNRF_TABLE(CROPTYPE) ! CO2 assimulation reduction factor(0-1) (caused by non-modeling part,e.g.pest,weeds) @@ -1187,6 +1198,9 @@ SUBROUTINE TRANSFER_MP_PARAMETERS(VEGTYPE,SOILTYPE,SLOPETYPE,SOILCOLOR,CROPTYPE, parameters%STPT = STPT_TABLE(CROPTYPE,:) ! fraction of carbohydrate flux to stem parameters%RTPT = RTPT_TABLE(CROPTYPE,:) ! fraction of carbohydrate flux to root parameters%GRAINPT = GRAINPT_TABLE(CROPTYPE,:) ! fraction of carbohydrate flux to grain + parameters%LFCT = LFCT_TABLE(CROPTYPE,:) ! fraction of translocation to grain ! Zhe Zhang 2020-07-13 + parameters%STCT = STCT_TABLE(CROPTYPE,:) ! fraction of translocation to grain + parameters%RTCT = RTCT_TABLE(CROPTYPE,:) ! fraction of translocation to grain parameters%BIO2LAI = BIO2LAI_TABLE(CROPTYPE) ! leaf are per living leaf biomass [m^2/kg] END IF @@ -1728,7 +1742,7 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT croptype(i,1,j) > croptype(i,4,j) ) then ! choose corn cropcat (i,j) = 1 - lfmassxy(i,j) = lai(i,j)/0.035 + lfmassxy(i,j) = lai(i,j)/0.015 ! Initialize lfmass Zhe Zhang 2020-07-13 stmassxy(i,j) = xsaixy(i,j)/0.003 elseif(croptype(i,2,j) > croptype(i,1,j) .and. & @@ -1736,7 +1750,7 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT croptype(i,2,j) > croptype(i,4,j) ) then ! choose soybean cropcat (i,j) = 2 - lfmassxy(i,j) = lai(i,j)/0.015 + lfmassxy(i,j) = lai(i,j)/0.030 ! Initialize lfmass Zhe Zhang 2020-07-13 stmassxy(i,j) = xsaixy(i,j)/0.003 else @@ -3193,7 +3207,7 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, lh(i,j) = qfx(i,j)*xlv hfx(i,j) = hfx_urb(i,j) + (1-frc_urb2d(i,j))*hfx_rural(i,j) ![W/m/m] sh_urb2d(i,j) = hfx_urb(i,j)/frc_urb2d(i,j) - lh_urb2d(i,j) = qfx_urb(i,j)*xlv + lh_urb2d(i,j) = qfx_urb(i,j)*xlv/frc_urb2d(i,j) g_urb2d(i,j) = grdflx_urb(i,j) rn_urb2d(i,j) = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j) ust(i,j) = (umom**2.+vmom**2.)**.25 diff --git a/phys/module_sf_noahmplsm.F b/phys/module_sf_noahmplsm.F index 9686a403c6..d0d557b55e 100644 --- a/phys/module_sf_noahmplsm.F +++ b/phys/module_sf_noahmplsm.F @@ -346,6 +346,9 @@ MODULE MODULE_SF_NOAHMPLSM REAL :: LFPT(NSTAGE) ! fraction of carbohydrate flux to leaf REAL :: STPT(NSTAGE) ! fraction of carbohydrate flux to stem REAL :: RTPT(NSTAGE) ! fraction of carbohydrate flux to root + REAL :: LFCT(NSTAGE) ! fraction of carbohydrate flux transallocate from leaf to grain ! Zhe Zhang 2020-07-13 + REAL :: STCT(NSTAGE) ! fraction of carbohydrate flux transallocate from stem to grain + REAL :: RTCT(NSTAGE) ! fraction of carbohydrate flux transallocate from root to grain REAL :: GRAINPT(NSTAGE) ! fraction of carbohydrate flux to grain REAL :: BIO2LAI ! leaf are per living leaf biomass [m^2/kg] @@ -8589,6 +8592,7 @@ SUBROUTINE CO2FLUX_CROP (parameters, REAL :: STMSMN REAL :: SAPM !stem area per unit mass (m2/g) REAL :: DIEST + REAL :: LFCONVERT !leaf to grain conversion ! Zhe Zhang 2020-07-13 REAL :: STCONVERT !stem to grain conversion [g/m2/s] REAL :: RTCONVERT !root to grain conversion [g/m2/s] ! -------------------------- constants ------------------------------- @@ -8691,15 +8695,23 @@ SUBROUTINE CO2FLUX_CROP (parameters, GPP = CBHYDRAFX* 0.4 !!g/m2/s C 0.4=12/30, CH20 to C + LFCONVERT = 0.0 ! Zhe Zhang 2020-07-13 STCONVERT = 0.0 RTCONVERT = 0.0 - IF(PGS==6) THEN - STCONVERT = STMASS*(0.00005*DT/3600.0) - STMASS = STMASS - STCONVERT - RTCONVERT = RTMASS*(0.0005*DT/3600.0) - RTMASS = RTMASS - RTCONVERT - GRAIN = GRAIN + STCONVERT + RTCONVERT - END IF + LFCONVERT = LFMASS*(parameters%LFCT(PGS)*DT/3600.0) + STCONVERT = STMASS*(parameters%STCT(PGS)*DT/3600.0) + RTCONVERT = RTMASS*(parameters%RTCT(PGS)*DT/3600.0) + LFMASS = LFMASS - LFCONVERT + STMASS = STMASS - STCONVERT + RTMASS = RTMASS - RTCONVERT + GRAIN = GRAIN + STCONVERT + RTCONVERT + LFCONVERT + !IF(PGS==6) THEN + ! STCONVERT = STMASS*(0.00005*DT/3600.0) + ! STMASS = STMASS - STCONVERT + ! RTCONVERT = RTMASS*(0.0005*DT/3600.0) + ! RTMASS = RTMASS - RTCONVERT + ! GRAIN = GRAIN + STCONVERT + RTCONVERT + !END IF IF(RTMASS.LT.0.0) THEN RTTOVR = NPPR @@ -9300,6 +9312,18 @@ MODULE NOAHMP_TABLES REAL :: GDDS4_TABLE(NCROP) ! GDD from seeding to intial reproductive REAL :: GDDS5_TABLE(NCROP) ! GDD from seeding to pysical maturity + REAL :: C3PSNI_TABLE(NCROP) !photosynthetic pathway: 0. = c4, 1. = c3 ! Zhe Zhang 2020-07-03 + REAL :: KC25I_TABLE(NCROP) !co2 michaelis-menten constant at 25c (pa) + REAL :: AKCI_TABLE(NCROP) !q10 for kc25 + REAL :: KO25I_TABLE(NCROP) !o2 michaelis-menten constant at 25c (pa) + REAL :: AKOI_TABLE(NCROP) !q10 for ko25 + REAL :: VCMX25I_TABLE(NCROP) !maximum rate of carboxylation at 25c (umol co2/m**2/s) + REAL :: AVCMXI_TABLE(NCROP) !q10 for vcmx25 + REAL :: BPI_TABLE(NCROP) !minimum leaf conductance (umol/m**2/s) + REAL :: MPI_TABLE(NCROP) !slope of conductance-to-photosynthesis relationship + REAL :: QE25I_TABLE(NCROP) !quantum efficiency at 25c (umol co2 / umol photon) + REAL :: FOLNMXI_TABLE(NCROP) !foliage nitrogen concentration when + INTEGER :: C3C4_TABLE(NCROP) ! photosynthetic pathway: 1. = c3 2. = c4 REAL :: AREF_TABLE(NCROP) ! reference maximum CO2 assimulation rate REAL :: PSNRF_TABLE(NCROP) ! CO2 assimulation reduction factor(0-1) (caused by non-modeling part,e.g.pest,weeds) @@ -9330,6 +9354,9 @@ MODULE NOAHMP_TABLES REAL :: STPT_TABLE(NCROP,NSTAGE) ! fraction of carbohydrate flux to stem REAL :: RTPT_TABLE(NCROP,NSTAGE) ! fraction of carbohydrate flux to root REAL :: GRAINPT_TABLE(NCROP,NSTAGE) ! fraction of carbohydrate flux to grain + REAL :: LFCT_TABLE(NCROP,NSTAGE) ! fraction of carbohydrate translocation from leaf to grain ! Zhe Zhang 2020-07-13 + REAL :: STCT_TABLE(NCROP,NSTAGE) ! stem to grain + REAL :: RTCT_TABLE(NCROP,NSTAGE) ! root to grain REAL :: BIO2LAI_TABLE(NCROP) ! leaf are per living leaf biomass [m^2/kg] ! MPTABLE.TBL optional parameters @@ -9868,6 +9895,17 @@ subroutine read_mp_crop_parameters() REAL, DIMENSION(NCROP) :: GDDS3 REAL, DIMENSION(NCROP) :: GDDS4 REAL, DIMENSION(NCROP) :: GDDS5 + REAL, DIMENSION(NCROP) :: C3PSN ! this session copied from stomata parameters Zhe Zhang 2020-07-13 + REAL, DIMENSION(NCROP) :: KC25 + REAL, DIMENSION(NCROP) :: AKC + REAL, DIMENSION(NCROP) :: KO25 + REAL, DIMENSION(NCROP) :: AKO + REAL, DIMENSION(NCROP) :: AVCMX + REAL, DIMENSION(NCROP) :: VCMX25 + REAL, DIMENSION(NCROP) :: BP + REAL, DIMENSION(NCROP) :: MP + REAL, DIMENSION(NCROP) :: FOLNMX + REAL, DIMENSION(NCROP) :: QE25 ! until here INTEGER, DIMENSION(NCROP) :: C3C4 REAL, DIMENSION(NCROP) :: AREF REAL, DIMENSION(NCROP) :: PSNRF @@ -9894,12 +9932,20 @@ subroutine read_mp_crop_parameters() REAL, DIMENSION(NCROP) :: STPT_S1,STPT_S2,STPT_S3,STPT_S4,STPT_S5,STPT_S6,STPT_S7,STPT_S8 REAL, DIMENSION(NCROP) :: RTPT_S1,RTPT_S2,RTPT_S3,RTPT_S4,RTPT_S5,RTPT_S6,RTPT_S7,RTPT_S8 REAL, DIMENSION(NCROP) :: GRAINPT_S1,GRAINPT_S2,GRAINPT_S3,GRAINPT_S4,GRAINPT_S5,GRAINPT_S6,GRAINPT_S7,GRAINPT_S8 + REAL, DIMENSION(NCROP) :: LFCT_S1,LFCT_S2,LFCT_S3,LFCT_S4,LFCT_S5,LFCT_S6,LFCT_S7,LFCT_S8 + REAL, DIMENSION(NCROP) :: STCT_S1,STCT_S2,STCT_S3,STCT_S4,STCT_S5,STCT_S6,STCT_S7,STCT_S8 + REAL, DIMENSION(NCROP) :: RTCT_S1,RTCT_S2,RTCT_S3,RTCT_S4,RTCT_S5,RTCT_S6,RTCT_S7,RTCT_S8 REAL, DIMENSION(NCROP) :: BIO2LAI - NAMELIST / noahmp_crop_parameters /DEFAULT_CROP, PLTDAY, HSDAY, PLANTPOP, IRRI, GDDTBASE, GDDTCUT, GDDS1, GDDS2, & - GDDS3, GDDS4, GDDS5, C3C4, AREF, PSNRF, I2PAR, TASSIM0, & - TASSIM1, TASSIM2, K, EPSI, Q10MR, FOLN_MX, LEFREEZ, & +! NAMELIST / noahmp_crop_parameters /DEFAULT_CROP, PLTDAY, HSDAY, PLANTPOP, IRRI, GDDTBASE, GDDTCUT, GDDS1, GDDS2, & +! GDDS3, GDDS4, GDDS5, C3C4, AREF, PSNRF, I2PAR, TASSIM0, & +! TASSIM1, TASSIM2, K, EPSI, Q10MR, FOLN_MX, LEFREEZ, & +! Zhe Zhang 2020-07-13 + NAMELIST / noahmp_crop_parameters /DEFAULT_CROP, PLTDAY, HSDAY, PLANTPOP, IRRI, GDDTBASE, GDDTCUT, GDDS1, GDDS2, GDDS3, GDDS4, GDDS5, & ! + C3PSN, KC25, AKC, KO25, AKO, AVCMX, VCMX25, BP, MP, FOLNMX, QE25, & ! parameters added from stomata + C3C4, AREF, PSNRF, I2PAR, TASSIM0, & + TASSIM1, TASSIM2, K, EPSI, Q10MR, FOLN_MX, LEFREEZ, & DILE_FC_S1,DILE_FC_S2,DILE_FC_S3,DILE_FC_S4,DILE_FC_S5,DILE_FC_S6,DILE_FC_S7,DILE_FC_S8, & DILE_FW_S1,DILE_FW_S2,DILE_FW_S3,DILE_FW_S4,DILE_FW_S5,DILE_FW_S6,DILE_FW_S7,DILE_FW_S8, & FRA_GR, & @@ -9911,6 +9957,9 @@ subroutine read_mp_crop_parameters() STPT_S1, STPT_S2, STPT_S3, STPT_S4, STPT_S5, STPT_S6, STPT_S7, STPT_S8, & RTPT_S1, RTPT_S2, RTPT_S3, RTPT_S4, RTPT_S5, RTPT_S6, RTPT_S7, RTPT_S8, & GRAINPT_S1,GRAINPT_S2,GRAINPT_S3,GRAINPT_S4,GRAINPT_S5,GRAINPT_S6,GRAINPT_S7,GRAINPT_S8, & + LFCT_S1,LFCT_S2,LFCT_S3,LFCT_S4,LFCT_S5,LFCT_S6,LFCT_S7,LFCT_S8, & + STCT_S1,STCT_S2,STCT_S3,STCT_S4,STCT_S5,STCT_S6,STCT_S7,STCT_S8, & + RTCT_S1,RTCT_S2,RTCT_S3,RTCT_S4,RTCT_S5,RTCT_S6,RTCT_S7,RTCT_S8, & BIO2LAI @@ -9927,6 +9976,17 @@ subroutine read_mp_crop_parameters() GDDS3_TABLE = -1.E36 GDDS4_TABLE = -1.E36 GDDS5_TABLE = -1.E36 + C3PSNI_TABLE = -1.E36 ! parameter from PSN copied from stomata ! Zhe Zhang 2020-07-13 + KC25I_TABLE = -1.E36 + AKCI_TABLE = -1.E36 + KO25I_TABLE = -1.E36 + AKOI_TABLE = -1.E36 + AVCMXI_TABLE = -1.E36 + VCMX25I_TABLE = -1.E36 + BPI_TABLE = -1.E36 + MPI_TABLE = -1.E36 + FOLNMXI_TABLE = -1.E36 + QE25I_TABLE = -1.E36 ! ends here C3C4_TABLE = -99999 AREF_TABLE = -1.E36 PSNRF_TABLE = -1.E36 @@ -9953,6 +10013,9 @@ subroutine read_mp_crop_parameters() STPT_TABLE = -1.E36 RTPT_TABLE = -1.E36 GRAINPT_TABLE = -1.E36 + LFCT_TABLE = -1.E36 ! convert start + STCT_TABLE = -1.E36 + RTCT_TABLE = -1.E36 ! convert end BIO2LAI_TABLE = -1.E36 @@ -9983,6 +10046,17 @@ subroutine read_mp_crop_parameters() GDDS3_TABLE = GDDS3 GDDS4_TABLE = GDDS4 GDDS5_TABLE = GDDS5 + C3PSNI_TABLE(1:5) = C3PSN(1:5) ! parameters from stomata ! Zhe Zhang 2020-07-13 + KC25I_TABLE(1:5) = KC25(1:5) + AKCI_TABLE(1:5) = AKC(1:5) + KO25I_TABLE(1:5) = KO25(1:5) + AKOI_TABLE(1:5) = AKO(1:5) + AVCMXI_TABLE(1:5) = AVCMX(1:5) + VCMX25I_TABLE(1:5) = VCMX25(1:5) + BPI_TABLE(1:5) = BP(1:5) + MPI_TABLE(1:5) = MP(1:5) + FOLNMXI_TABLE(1:5) = FOLNMX(1:5) + QE25I_TABLE(1:5) = QE25(1:5) ! ends here C3C4_TABLE = C3C4 AREF_TABLE = AREF PSNRF_TABLE = PSNRF @@ -10072,6 +10146,30 @@ subroutine read_mp_crop_parameters() GRAINPT_TABLE(:,6) = GRAINPT_S6 GRAINPT_TABLE(:,7) = GRAINPT_S7 GRAINPT_TABLE(:,8) = GRAINPT_S8 + LFCT_TABLE(:,1) = LFCT_S1 + LFCT_TABLE(:,2) = LFCT_S2 + LFCT_TABLE(:,3) = LFCT_S3 + LFCT_TABLE(:,4) = LFCT_S4 + LFCT_TABLE(:,5) = LFCT_S5 + LFCT_TABLE(:,6) = LFCT_S6 + LFCT_TABLE(:,7) = LFCT_S7 + LFCT_TABLE(:,8) = LFCT_S8 + STCT_TABLE(:,1) = STCT_S1 + STCT_TABLE(:,2) = STCT_S2 + STCT_TABLE(:,3) = STCT_S3 + STCT_TABLE(:,4) = STCT_S4 + STCT_TABLE(:,5) = STCT_S5 + STCT_TABLE(:,6) = STCT_S6 + STCT_TABLE(:,7) = STCT_S7 + STCT_TABLE(:,8) = STCT_S8 + RTCT_TABLE(:,1) = RTCT_S1 + RTCT_TABLE(:,2) = RTCT_S2 + RTCT_TABLE(:,3) = RTCT_S3 + RTCT_TABLE(:,4) = RTCT_S4 + RTCT_TABLE(:,5) = RTCT_S5 + RTCT_TABLE(:,6) = RTCT_S6 + RTCT_TABLE(:,7) = RTCT_S7 + RTCT_TABLE(:,8) = RTCT_S8 BIO2LAI_TABLE = BIO2LAI end subroutine read_mp_crop_parameters diff --git a/phys/module_sf_sfclayrev.F b/phys/module_sf_sfclayrev.F index 76574cc115..6f9c30f9c5 100644 --- a/phys/module_sf_sfclayrev.F +++ b/phys/module_sf_sfclayrev.F @@ -1202,7 +1202,7 @@ function psim_stable(zolf) real :: rzol nzol = int(zolf*100.) rzol = zolf*100. - nzol - if(nzol+1 .le. 1000)then + if(nzol+1 .lt. 1000)then psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol)) else psim_stable = psim_stable_full(zolf) @@ -1215,7 +1215,7 @@ function psih_stable(zolf) real :: rzol nzol = int(zolf*100.) rzol = zolf*100. - nzol - if(nzol+1 .le. 1000)then + if(nzol+1 .lt. 1000)then psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol)) else psih_stable = psih_stable_full(zolf) @@ -1228,7 +1228,7 @@ function psim_unstable(zolf) real :: rzol nzol = int(-zolf*100.) rzol = -zolf*100. - nzol - if(nzol+1 .le. 1000)then + if(nzol+1 .lt. 1000)then psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol)) else psim_unstable = psim_unstable_full(zolf) @@ -1241,7 +1241,7 @@ function psih_unstable(zolf) real :: rzol nzol = int(-zolf*100.) rzol = -zolf*100. - nzol - if(nzol+1 .le. 1000)then + if(nzol+1 .lt. 1000)then psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol)) else psih_unstable = psih_unstable_full(zolf) diff --git a/phys/module_shcu_deng.F b/phys/module_shcu_deng.F index 408135ebd2..edba0e5ce6 100644 --- a/phys/module_shcu_deng.F +++ b/phys/module_shcu_deng.F @@ -255,8 +255,7 @@ SUBROUTINE deng_shcu_driver( & CA_RAD(I,K,J)=(1.0-CLDAREAA(I,K,J))*CS+CLDAREAA(I,K,J) CLDFRA_SH(I,K,J)= CA_RAD(I,K,J) -! CW_RAD(I,K,J)=CLDLIQA(I,K,J)+QC(I,K,J) - CW_RAD(I,K,J)=CLDLIQA(I,K,J)*CLDAREAA(I,K,J) + QC(I,K,J) + CW_RAD(I,K,J)=CLDLIQA(I,K,J)*CLDAREAA(I,K,J) ENDDO ENDDO ENDDO @@ -414,12 +413,9 @@ SUBROUTINE deng_shcu_driver( & ENDIF IF(CLDDPTHB(I,J) .LE. 0.0) THEN ! No active updraft - IF( CLDAREAB(I,K,J) .LE. 1.0e-3 .OR. CLDLIQB(I,K,J) .LE. 1.0e-17 ) THEN - ! QC(I,K,J)=QC(I,K,J)+CLDAREAB(I,K,J)*CLDLIQB(I,K,J) RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)+CLDAREAB(I,K,J)*CLDLIQB(I,K,J)/DT CLDAREAB(I,K,J) = 0.0 CLDLIQB(I,K,J) = 0.0 - ENDIF ENDIF ENDDO ENDDO @@ -679,8 +675,8 @@ SUBROUTINE deng_shcu(I,J, & ! in DCLDTOP(I,J)=0.0 DCLDBASE(I,J)=0.0 CLDBMFLX(I,J)=0.0 - RADIUSC(I,J)=0.0 #endif + RADIUSC(I,J)=0.0 ! !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED FROM THE !...BOTTOM-UP IN THE SHALLOW CONVECTION SCHEME... @@ -2099,6 +2095,7 @@ SUBROUTINE deng_shcu(I,J, & ! in ! IF(CAPEI(I,J) .LT. 100.0) THEN ! averaging KF number of clouds NT = NINT(15.0*60.0/(0.5*DT))-1 ! for 15 min. + NT = MIN(NT, 100) ELSE NT = 0 ENDIF @@ -2799,7 +2796,7 @@ SUBROUTINE deng_shcu(I,J, & ! in ! ! CALCULATE THE CONVECTIVE RAINFALL ! - RAINSHV(I,J)=.1*.5*DT2*PPTFLX/DXSQ + RAINSHV(I,J)=.5*DT2*PPTFLX/DXSQ IF ( wrf_dm_on_monitor()) THEN CALL get_wrf_debug_level( dbg_level ) @@ -4478,7 +4475,7 @@ SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0,ktau,i,j,nk,tra IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN write(98,'(a,6i5,2f9.2,i3)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, & p/100.0, thes, tracker - call flush(98) + flush(98) ENDIF ! t00=ttab(ithtb ,iptb ) @@ -4796,7 +4793,7 @@ SUBROUTINE TPMIX2DD(p,thes,ts,qs,ktau,i,j,nk) ithtb=int(tth)+1 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN write(97,'(a,6i5,2f9.2)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, p/100., thes - call flush(97) + flush(97) ENDIF ! t00=ttab(ithtb ,iptb ) diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index d082f547ad..5bb4ba964c 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -1392,8 +1392,8 @@ SUBROUTINE surface_driver( & REAL, PARAMETER :: PI_GRECO=3.14159 INTEGER :: end_hour, irr_start,xt24,irr_day REAL :: constants_irrigation - REAL,DIMENSION( ims:ime, jms:jme ) :: IRRIGATION_CHANNEL - REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN),OPTIONAL:: IRRIGATION + REAL, DIMENSION( ims:ime, jms:jme ) :: IRRIGATION_CHANNEL + REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) , OPTIONAL:: IRRIGATION REAL, INTENT(IN),OPTIONAL:: irr_daily_amount INTEGER :: phase INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT),OPTIONAL :: irr_rand_field @@ -1499,6 +1499,7 @@ SUBROUTINE surface_driver( & IF ( PRESENT( rainshv ))RAINBL(i,j) = RAINBL(i,j) + RAINSHV(i,j) RAINBL(i,j) = MAX (RAINBL(i,j), 0.0) #if (EM_CORE==1) + IRRIGATION_CHANNEL(i,j) = 0. sf_surf_irr: SELECT CASE(sf_surf_irr_scheme) CASE(DRIP) CALL drip_irrigation( & @@ -2740,7 +2741,7 @@ SUBROUTINE surface_driver( & dl_u_bep,sf_bep,vl_bep & !O multi-layer urban ,sfcheadrt,INFXSRT, soldrain & !hydro ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & ! fasdas - ,RS,XLAIDYN) + ,RS,XLAIDYN,IRRIGATION_CHANNEL) ELSE CALL wrf_error_fatal('Lack arguments to call lsm_mosaic') ENDIF @@ -2860,7 +2861,7 @@ SUBROUTINE surface_driver( & ! ! END FASDAS ! - ,RS,XLAIDYN) + ,RS,XLAIDYN,IRRIGATION_CHANNEL) ENDIF call seaice_noah( SEAICE_ALBEDO_OPT, SEAICE_ALBEDO_DEFAULT, SEAICE_THICKNESS_OPT, & @@ -4369,7 +4370,7 @@ subroutine myjsfc_seaice_wrapper(ITIMESTEP,HT,DZ, & ! QSFC_SEA calculation as done in myjsfc for open water points PSFC = PINT(I,LOWLYR(I,J),J) QSFC_SEA(i,j) = PQ0SEA/PSFC*EXP(A2S*(TSK(i,j)-A3S)/(TSK(i,j)-A4S)) - QSFC(i,j) = QSFC(i,j) - (1.0-XICE(i,j)) * QSFC_SEA(i,j) / XICE(i,j) + QSFC(i,j) = (QSFC(i,j) - (1.0-XICE(i,j)) * QSFC_SEA(i,j)) / XICE(i,j) ! HFX_SEA(i,j) = HFX(i,j) QFX_SEA(i,j) = QFX(i,j) @@ -4736,7 +4737,7 @@ subroutine qnsesfc_seaice_wrapper(ITIMESTEP,HT,DZ, & ! QSFC_SEA calculation as done in qnsesfc for open water points PSFC = PINT(I,LOWLYR(I,J),J) QSFC_SEA(i,j) = PQ0SEA/PSFC*EXP(A2S*(TSK(i,j)-A3S)/(TSK(i,j)-A4S)) - QSFC(i,j) = QSFC(i,j) - (1.0-XICE(i,j)) * QSFC_SEA(i,j) / XICE(i,j) + QSFC(i,j) = (QSFC(i,j) - (1.0-XICE(i,j)) * QSFC_SEA(i,j)) / XICE(i,j) ! HFX_SEA(i,j) = HFX(i,j) QFX_SEA(i,j) = QFX(i,j) diff --git a/run/MPTABLE.TBL b/run/MPTABLE.TBL index be26e9b248..03afe68d7a 100644 --- a/run/MPTABLE.TBL +++ b/run/MPTABLE.TBL @@ -360,21 +360,31 @@ DEFAULT_CROP = 0 ! The default crop type(1- ! 1 2 3 4 5 !---------------------------------------------------------- -PLTDAY = 130, 111, 111, 111, 111, ! Planting date -HSDAY = 280, 300, 300, 300, 300, ! Harvest date +PLTDAY = 111, 131, 111, 111, 111, ! Planting date +HSDAY = 300, 280, 300, 300, 300, ! Harvest date PLANTPOP = 78.0, 78.0, 78.0, 78.0, 78.0, ! Plant density [per ha] - used? IRRI = 0.0, 0.0, 0.0, 0.0, 0.0, ! Irrigation strategy 0= non-irrigation 1=irrigation (no water-stress) GDDTBASE = 10.0, 10.0, 10.0, 10.0, 10.0, ! Base temperature for GDD accumulation [C] GDDTCUT = 30.0, 30.0, 30.0, 30.0, 30.0, ! Upper temperature for GDD accumulation [C] -GDDS1 = 60.0, 50.0, 50.0, 50.0, 50.0, ! GDD from seeding to emergence -GDDS2 = 675.0, 718.0, 718.0, 718.0, 718.0, ! GDD from seeding to initial vegetative -GDDS3 = 1183.0, 933.0, 933.0, 933.0, 933.0, ! GDD from seeding to post vegetative -GDDS4 = 1253.0, 1103.0, 1103.0, 1103.0, 1103.0, ! GDD from seeding to intial reproductive -GDDS5 = 1605.0, 1555.0, 1555.0, 1555.0, 1555.0, ! GDD from seeding to pysical maturity - +GDDS1 = 50.0, 60.0, 50.0, 50.0, 50.0, ! GDD from seeding to emergence +GDDS2 = 625.0, 675.0, 718.0, 718.0, 718.0, ! GDD from seeding to initial vegetative +GDDS3 = 933.0, 1183.0, 933.0, 933.0, 933.0, ! GDD from seeding to post vegetative +GDDS4 = 1103.0, 1253.0, 1103.0, 1103.0, 1103.0, ! GDD from seeding to intial reproductive +GDDS5 = 1555.0, 1605.0, 1555.0, 1555.0, 1555.0, ! GDD from seeding to pysical maturity +C3PSN = 0.0, 1.0, 1.0, 1.0, 1.0, ! transfer crop-specific photosynthetic parameters +KC25 = 30.0, 30.0, 30.0, 30.0, 30.0, ! Zhe Zhang +AKC = 2.1, 2.1, 2.1, 2.1, 2.1, ! 2020-02-05 +KO25 = 3.E4, 3.E4, 3.E4, 3.E4, 3.E4, ! +AKO = 1.2, 1.2, 1.2, 1.2, 1.2, ! +AVCMX = 2.4, 2.4, 2.4, 2.4, 2.4, ! +VCMX25 = 60.0, 80.0, 60.0, 60.0, 55.0, ! +BP = 4.E4, 1.E4, 2.E3, 2.E3, 2.E3, ! +MP = 4., 9., 6., 9., 9., ! +FOLNMX = 1.5, 1.5, 1.5, 1.5, 1.5, ! +QE25 = 0.05, 0.06, 0.06, 0.06, 0.06, ! C3C4 = 2, 1, 2, 2, 2, ! photosynthetic pathway: 1. = c3 2. = c4 -Aref = 7.0, 7.0, 7.0, 7.0, 7.0, ! reference maximum CO2 assimulation rate +Aref = 7.0, 7.0, 7.0, 7.0, 7.0, ! reference maximum CO2 assimulation rate PSNRF = 0.85, 0.85, 0.85, 0.85, 0.85, ! CO2 assimulation reduction factor(0-1) (caused by non-modeling part,e.g.pest,weeds) I2PAR = 0.5, 0.5, 0.5, 0.5, 0.5, ! Fraction of incoming solar radiation to photosynthetically active radiation TASSIM0 = 8.0, 8.0, 8.0, 8.0, 8.0, ! Minimum temperature for CO2 assimulation [C] @@ -390,7 +400,7 @@ LEFREEZ = 268, 268, 268, 268, 268, ! characteristic T for lea DILE_FC_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! coeficient for temperature leaf stress death [1/s] DILE_FC_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages DILE_FC_S3 = 0.0, 0.0, 0.0, 0.0, 0.0, -DILE_FC_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, +DILE_FC_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, DILE_FC_S5 = 0.5, 0.5, 0.5, 0.5, 0.5, DILE_FC_S6 = 0.5, 0.5, 0.5, 0.5, 0.5, DILE_FC_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, @@ -411,8 +421,8 @@ LF_OVRC_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of leaf turnove LF_OVRC_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages LF_OVRC_S3 = 0.0, 0.0, 0.0, 0.0, 0.0, LF_OVRC_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, -LF_OVRC_S5 = 0.2, 0.48, 0.48, 0.48, 0.48, -LF_OVRC_S6 = 0.3, 0.48, 0.48, 0.48, 0.48, +LF_OVRC_S5 = 0.2, 0.2, 0.48, 0.48, 0.48, +LF_OVRC_S6 = 0.3, 0.3, 0.48, 0.48, 0.48, LF_OVRC_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, LF_OVRC_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, @@ -420,8 +430,8 @@ ST_OVRC_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of stem turnove ST_OVRC_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages ST_OVRC_S3 = 0.0, 0.0, 0.0, 0.0, 0.0, ST_OVRC_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, -ST_OVRC_S5 = 0.12, 0.12, 0.12, 0.12, 0.12, -ST_OVRC_S6 = 0.06, 0.06, 0.06, 0.06, 0.06, +ST_OVRC_S5 = 0.2, 0.12, 0.12, 0.12, 0.12, +ST_OVRC_S6 = 0.3, 0.06, 0.06, 0.06, 0.06, ST_OVRC_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, ST_OVRC_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, @@ -429,21 +439,20 @@ RT_OVRC_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of root tunrove RT_OVRC_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages RT_OVRC_S3 = 0.0, 0.0, 0.0, 0.0, 0.0, RT_OVRC_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, -RT_OVRC_S5 = 0.12, 0.12, 0.12, 0.12, 0.12, -RT_OVRC_S6 = 0.06, 0.06, 0.06, 0.06, 0.06, +RT_OVRC_S5 = 0.12, 0.12, 0.12, 0.12, 0.12, +RT_OVRC_S6 = 0.06, 0.06, 0.06, 0.06, 0.06, RT_OVRC_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, RT_OVRC_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, - -LFMR25 = 1.0, 1.0, 1.0, 1.0, 1.0, ! leaf maintenance respiration at 25C [umol CO2/m**2 /s] -STMR25 = 0.05, 0.1, 0.1, 0.1, 0.1, ! stem maintenance respiration at 25C [umol CO2/kg bio/s] -RTMR25 = 0.05, 0.0, 0.0, 0.0, 0.0, ! root maintenance respiration at 25C [umol CO2/kg bio/s] -GRAINMR25 = 0.0, 0.1, 0.1, 0.1, 0.1, ! grain maintenance respiration at 25C [umol CO2/kg bio/s] +LFMR25 = 0.8, 1.0, 1.0, 1.0, 1.0, ! leaf maintenance respiration at 25C [umol CO2/m**2 /s] +STMR25 = 0.05, 0.05, 0.1, 0.1, 0.1, ! stem maintenance respiration at 25C [umol CO2/kg bio/s] +RTMR25 = 0.05, 0.05, 0.0, 0.0, 0.0, ! root maintenance respiration at 25C [umol CO2/kg bio/s] +GRAINMR25 = 0.0, 0.0, 0.1, 0.1, 0.1, ! grain maintenance respiration at 25C [umol CO2/kg bio/s] LFPT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of carbohydrate flux to leaf LFPT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages -LFPT_S3 = 0.4, 0.4, 0.4, 0.4, 0.4, -LFPT_S4 = 0.2, 0.2, 0.2, 0.2, 0.2, +LFPT_S3 = 0.36, 0.4, 0.4, 0.4, 0.4, +LFPT_S4 = 0.1, 0.2, 0.2, 0.2, 0.2, LFPT_S5 = 0.0, 0.0, 0.0, 0.0, 0.0, LFPT_S6 = 0.0, 0.0, 0.0, 0.0, 0.0, LFPT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, @@ -451,33 +460,60 @@ LFPT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, STPT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of carbohydrate flux to stem STPT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages -STPT_S3 = 0.2, 0.2, 0.2, 0.2, 0.2, -STPT_S4 = 0.5, 0.5, 0.5, 0.5, 0.5, -STPT_S5 = 0.0, 0.15, 0.15, 0.15, 0.15, -STPT_S6 = 0.0, 0.05, 0.05, 0.05, 0.05, +STPT_S3 = 0.24, 0.2, 0.2, 0.2, 0.2, +STPT_S4 = 0.6, 0.5, 0.5, 0.5, 0.5, +STPT_S5 = 0.0, 0.0, 0.15, 0.15, 0.15, +STPT_S6 = 0.0, 0.0, 0.05, 0.05, 0.05, STPT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, -STPT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, +STPT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, RTPT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of carbohydrate flux to root RTPT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages -RTPT_S3 = 0.34, 0.4, 0.4, 0.4, 0.4, +RTPT_S3 = 0.4, 0.4, 0.4, 0.4, 0.4, RTPT_S4 = 0.3, 0.3, 0.3, 0.3, 0.3, -RTPT_S5 = 0.05, 0.05, 0.05, 0.05, 0.05, -RTPT_S6 = 0.0, 0.05, 0.05, 0.05, 0.05, +RTPT_S5 = 0.05, 0.05, 0.05, 0.05, 0.05, +RTPT_S6 = 0.0, 0.0, 0.05, 0.05, 0.05, RTPT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, RTPT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, - + GRAINPT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! fraction of carbohydrate flux to grain GRAINPT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, ! One row for each of 8 stages GRAINPT_S3 = 0.0, 0.0, 0.0, 0.0, 0.0, GRAINPT_S4 = 0.0, 0.0, 0.0, 0.0, 0.0, -GRAINPT_S5 = 0.95, 0.8, 0.8, 0.8, 0.8, -GRAINPT_S6 = 1.0, 0.9, 0.9, 0.9, 0.9, +GRAINPT_S5 = 0.95, 0.95, 0.8, 0.8, 0.8, +GRAINPT_S6 = 1.0, 1.0, 0.9, 0.9, 0.9, GRAINPT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, GRAINPT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, +LFCT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! carbohydrate translocation +LFCT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, +LFCT_S3 = 0.0, 0., 0.4, 0.4, 0.4, +LFCT_S4 = 0.0, 0., 0.3, 0.3, 0.3, +LFCT_S5 = 0.0, 0.0, 0.05, 0.05, 0.05, +LFCT_S6 = 0.0, 0.0, 0.05, 0.05, 0.05, +LFCT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, +LFCT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, + +STCT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! carbohydrate translocation +STCT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, +STCT_S3 = 0.0, 0., 0.4, 0.4, 0.4, +STCT_S4 = 0.0, 0., 0.3, 0.3, 0.3, +STCT_S5 = 0.0, 0., 0.05, 0.05, 0.05, +STCT_S6 = 0.0, 0.0, 0.05, 0.05, 0.05, +STCT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, +STCT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, + +RTCT_S1 = 0.0, 0.0, 0.0, 0.0, 0.0, ! carbohydrate translocation +RTCT_S2 = 0.0, 0.0, 0.0, 0.0, 0.0, +RTCT_S3 = 0.0, 0., 0.4, 0.4, 0.4, +RTCT_S4 = 0.0, 0., 0.3, 0.3, 0.3, +RTCT_S5 = 0.0, 0., 0.05, 0.05, 0.05, +RTCT_S6 = 0.0, 0.0, 0.05, 0.05, 0.05, +RTCT_S7 = 0.0, 0.0, 0.0, 0.0, 0.0, +RTCT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, + +BIO2LAI = 0.015, 0.030, 0.015, 0.015, 0.015, ! leaf are per living leaf biomass [m^2/kg] -BIO2LAI = 0.035, 0.015, 0.015, 0.015, 0.015, ! leaf are per living leaf biomass [m^2/kg] / diff --git a/run/README.tslist b/run/README.tslist index e65e46ad7c..1e912e6f5c 100644 --- a/run/README.tslist +++ b/run/README.tslist @@ -125,6 +125,16 @@ cbaseht_tot: CLOUD BASE HEIGHT RES + UNRES (m) ctopht_tot: CLOUD TOP HEIGHT RES + UNRES (m) clrnidx: CLEARNESS INDEX () sza: SOLAR ZENITH ANGLE (deg) +swdown: DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE (W m-2) +swddni: SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE (W m-2) +swddif: SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE (W m-2) +swdownc: DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE (W m-2) +swddnic: CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE (W m-2) +swdown2: DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE FROM FARMS (W m-2) +swddni2: SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS (W m-2) +swddif2: SHORTWAVE SURFACE DOWNWARD DIFFUSE IRRADIANCE FROM FARMS (W m-2) +swdownc2: DOWNWARD CLEAR-SKY SHORTWAVE FLUX AT GROUND SURFACE FROM FARMS (W m-2) +swddnic2: CLEAR-SKY SHORTWAVE SURFACE DOWNWARD DIRECT NORMAL IRRADIANCE FROM FARMS (W m-2) Format of the files of vertical profile: diff --git a/share/module_bc.F b/share/module_bc.F index ea2c5b2382..c0b071893a 100644 --- a/share/module_bc.F +++ b/share/module_bc.F @@ -680,7 +680,8 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) :: dat TYPE( grid_config_rec_type ) config_flags - INTEGER :: i, j, k, istag, jstag, itime, k_end + INTEGER :: i, j, k, istag, jstag, itime, k_end, & + i_start, i_end LOGICAL :: debug, open_bc_copy @@ -779,7 +780,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & symmetry_xs: IF( ( config_flags%symmetric_xs ) .and. & ( its == ids ) ) THEN - IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN + IF ( istag == -1 ) THEN DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag) DO k = kts, k_end @@ -823,7 +824,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & symmetry_xe: IF( ( config_flags%symmetric_xe ) .and. & ( ite == ide ) ) THEN - IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN + IF ( istag == -1 ) THEN DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag) DO k = kts, k_end @@ -927,13 +928,27 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & ! same procedure in y +! Set the starting and ending loop indexes in the 'i' direction, so that +! halo cells on the edge of the domain are also updated. Begin with a default +! start and end index for inner tiles, and then modify if the tile is on the +! edge of the domain. + + i_start = MAX(ids, its-1) + i_end = MIN(ite+1, ide+istag) + IF ( its .eq. ids) THEN + i_start = ims + END IF + IF ( ite .eq. ide) THEN + i_end = ime + END IF + periodicity_y: IF( ( config_flags%periodic_y ) ) THEN IF ( ( jds == jps ) .and. ( jde == jpe ) ) THEN ! test if both north and south on processor IF( jts == jds ) then DO j = 0, -(bdyzone-1), -1 DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jds+j-1) = dat(i,k,jde+j-1) ENDDO ENDDO @@ -945,7 +960,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & DO j = -jstag, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde+j+jstag) = dat(i,k,jds+j+jstag) ENDDO ENDDO @@ -960,11 +975,11 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & symmetry_ys: IF( ( config_flags%symmetric_ys ) .and. & ( jts == jds) ) THEN - IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN + IF ( jstag == -1 ) THEN DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jds-j) = dat(i,k,jds+j-1) ENDDO ENDDO @@ -976,7 +991,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jds-j) = - dat(i,k,jds+j) ENDDO ENDDO @@ -986,7 +1001,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jds-j) = dat(i,k,jds+j) ENDDO ENDDO @@ -1003,11 +1018,11 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & symmetry_ye: IF( ( config_flags%symmetric_ye ) .and. & ( jte == jde ) ) THEN - IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN + IF ( jstag == -1 ) THEN DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde+j-1) = dat(i,k,jde-j) ENDDO ENDDO @@ -1019,7 +1034,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde+j) = - dat(i,k,jde-j) ENDDO ENDDO @@ -1029,7 +1044,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & DO j = 1, bdyzone DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde+j) = dat(i,k,jde-j) ENDDO ENDDO @@ -1050,7 +1065,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & ( jts == jds) .and. open_bc_copy ) THEN DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jds-1) = dat(i,k,jds) dat(i,k,jds-2) = dat(i,k,jds) dat(i,k,jds-3) = dat(i,k,jds) @@ -1070,7 +1085,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & IF (variable /= 'v' .and. variable /= 'y' ) THEN DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde ) = dat(i,k,jde-1) dat(i,k,jde+1) = dat(i,k,jde-1) dat(i,k,jde+2) = dat(i,k,jde-1) @@ -1080,7 +1095,7 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & ELSE DO k = kts, k_end - DO i = MAX(ids,its-1), MIN(ite+1,ide+istag) + DO i = i_start, i_end dat(i,k,jde+1) = dat(i,k,jde) dat(i,k,jde+2) = dat(i,k,jde) dat(i,k,jde+3) = dat(i,k,jde) @@ -1095,54 +1110,6 @@ SUBROUTINE set_physical_bc3d( dat, variable_in, & END IF periodicity_y -! fix corners for doubly periodic domains - - IF ( config_flags%periodic_x .and. config_flags%periodic_y & - .and. (ids == ips) .and. (ide == ipe) & - .and. (jds == jps) .and. (jde == jpe) ) THEN - - IF ( (its == ids) .and. (jts == jds) ) THEN ! lower left corner fill - DO j = 0, -(bdyzone-1), -1 - DO k = kts, k_end - DO i = 0, -(bdyzone-1), -1 - dat(ids+i-1,k,jds+j-1) = dat(ide+i-1,k,jde+j-1) - ENDDO - ENDDO - ENDDO - END IF - - IF ( (ite == ide) .and. (jts == jds) ) THEN ! lower right corner fill - DO j = 0, -(bdyzone-1), -1 - DO k = kts, k_end - DO i = 1, bdyzone - dat(ide+i+istag,k,jds+j-1) = dat(ids+i+istag,k,jde+j-1) - ENDDO - ENDDO - ENDDO - END IF - - IF ( (ite == ide) .and. (jte == jde) ) THEN ! upper right corner fill - DO j = 1, bdyzone - DO k = kts, k_end - DO i = 1, bdyzone - dat(ide+i+istag,k,jde+j+jstag) = dat(ids+i+istag,k,jds+j+jstag) - ENDDO - ENDDO - ENDDO - END IF - - IF ( (its == ids) .and. (jte == jde) ) THEN ! upper left corner fill - DO j = 1, bdyzone - DO k = kts, k_end - DO i = 0, -(bdyzone-1), -1 - dat(ids+i-1,k,jde+j+jstag) = dat(ide+i-1,k,jds+j+jstag) - ENDDO - ENDDO - ENDDO - END IF - - END IF - END SUBROUTINE set_physical_bc3d SUBROUTINE init_module_bc diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index db28b315a3..69cca36211 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -576,24 +576,6 @@ END FUNCTION bep_bem_nbui_max END IF ENDDO -!----------------------------------------------------------------------- -! The FARMS radiation option (swint_opt==2) requires both effective radii -! and mass for cloud, ice, and snow. A run-time option is available to -! disable the use of effective radii in the MP schemes. These two options -! may not be used together. -!----------------------------------------------------------------------- - oops = 0 - IF ( ( model_config_rec%swint_opt .EQ. 2 ) .AND. & - ( model_config_rec%use_mp_re .NE. 1 ) ) THEN - oops = oops + 1 - END IF - - IF ( oops .GT. 0 ) THEN - wrf_err_message = '--- ERROR: FARMS (swint_opt=2) requires effective radii (use_mp_re=1)' - CALL wrf_message ( wrf_err_message ) - count_fatal_error = count_fatal_error + 1 - END IF - #endif !----------------------------------------------------------------------- diff --git a/share/module_model_constants.F b/share/module_model_constants.F index d9a31823c7..e4a7b9eb77 100644 --- a/share/module_model_constants.F +++ b/share/module_model_constants.F @@ -63,6 +63,10 @@ MODULE module_model_constants REAL , PARAMETER :: rhowater = 1000. ! density of liquid water at 0^oC (kg m^-3) REAL , PARAMETER :: rhosnow = 100. ! density of snow (kg m^-3) REAL , PARAMETER :: rhoair0 = 1.28 ! density of dry air at 0^oC and 1000mb pressure (kg m^-3) + + REAL , PARAMETER :: RE_QC_BG = 2.49E-6 ! effective radius of cloud for background (m) + REAL , PARAMETER :: RE_QI_BG = 4.99E-6 ! effective radius of ice for background (m) + REAL , PARAMETER :: RE_QS_BG = 9.99E-6 ! effective radius of snow for background (m) ! ! Now namelist-specified parameter: ccn_conc - RAS ! REAL , PARAMETER :: n_ccn0 = 1.0E8 diff --git a/share/wrf_timeseries.F b/share/wrf_timeseries.F index c40e38a7f6..226a070284 100644 --- a/share/wrf_timeseries.F +++ b/share/wrf_timeseries.F @@ -351,6 +351,7 @@ SUBROUTINE calc_ts( grid ) INTEGER :: i, k, mm, n, ix, iy, rc REAL :: earth_u, earth_v, & output_t, output_q, clw, xtime_minutes + REAL, PARAMETER :: MISSING = -999.0 REAL, ALLOCATABLE, DIMENSION(:) :: p8w REAL, ALLOCATABLE, DIMENSION(:) :: earth_u_profile, earth_v_profile TYPE (grid_config_rec_type) :: config_flags @@ -516,6 +517,37 @@ SUBROUTINE calc_ts( grid ) grid%ts_ctopht_tot(n,i) = grid%ctopht_tot(ix,iy) grid%ts_clrnidx(n,i) = grid%clrnidx(ix,iy) grid%ts_sza(n,i) = grid%sza(ix,iy) + grid%ts_swdown(n,i) = grid%swdown(ix,iy) + grid%ts_swddni(n,i) = grid%swddni(ix,iy) + grid%ts_swddif(n,i) = grid%swddif(ix,iy) + ! Calc extra solar variables if FARMS or RRTMG/RRTMG FAST + if ( config_flags%swint_opt == 2 .or. & + config_flags%ra_sw_physics == RRTMG_SWSCHEME .or. & + config_flags%ra_sw_physics == RRTMG_SWSCHEME_FAST ) then + grid%ts_swdownc(n,i) = grid%swdownc(ix,iy) + grid%ts_swddnic(n,i) = grid%swddnic(ix,iy) + if ( config_flags%swint_opt == 2 ) then ! FARMS + grid%ts_swdown2(n,i) = grid%swdown2(ix,iy) + grid%ts_swddni2(n,i) = grid%swddni2(ix,iy) + grid%ts_swddif2(n,i) = grid%swddif2(ix,iy) + grid%ts_swdownc2(n,i) = grid%swdownc2(ix,iy) + grid%ts_swddnic2(n,i) = grid%swddnic2(ix,iy) + else + grid%ts_swdown2(n,i) = MISSING + grid%ts_swddni2(n,i) = MISSING + grid%ts_swddif2(n,i) = MISSING + grid%ts_swdownc2(n,i) = MISSING + grid%ts_swddnic2(n,i) = MISSING + end if + else + grid%ts_swdownc(n,i) = MISSING + grid%ts_swddnic(n,i) = MISSING + grid%ts_swdown2(n,i) = MISSING + grid%ts_swddni2(n,i) = MISSING + grid%ts_swddif2(n,i) = MISSING + grid%ts_swdownc2(n,i) = MISSING + grid%ts_swddnic2(n,i) = MISSING + end if END IF #else grid%ts_tsk(n,i) = grid%nmm_tsk(ix,iy) @@ -580,6 +612,16 @@ SUBROUTINE calc_ts( grid ) grid%ts_ctopht_tot(n,i) = 1.E30 grid%ts_clrnidx(n,i) = 1.E30 grid%ts_sza(n,i) = 1.E30 + grid%ts_swdown(n,i) = 1.E30 + grid%ts_swddni(n,i) = 1.E30 + grid%ts_swddif(n,i) = 1.E30 + grid%ts_swdownc(n,i) = 1.E30 + grid%ts_swddnic(n,i) = 1.E30 + grid%ts_swdown2(n,i) = 1.E30 + grid%ts_swddni2(n,i) = 1.E30 + grid%ts_swddif2(n,i) = 1.E30 + grid%ts_swdownc2(n,i) = 1.E30 + grid%ts_swddnic2(n,i) = 1.E30 END IF #endif grid%ts_tsk(n,i) = 1.E30 @@ -787,6 +829,36 @@ SUBROUTINE write_ts( grid ) ts_buf(:,:) = grid%ts_sza(:,:) CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_sza(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swdown(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdown(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddni(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddni(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddif(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddif(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swdownc(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdownc(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddnic(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddnic(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swdown2(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdown2(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddni2(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddni2(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddif2(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddif2(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swdownc2(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdownc2(:,:),grid%ts_buf_size*grid%max_ts_locs) + + ts_buf(:,:) = grid%ts_swddnic2(:,:) + CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddnic2(:,:),grid%ts_buf_size*grid%max_ts_locs) END IF #endif @@ -836,7 +908,7 @@ SUBROUTINE write_ts( grid ) grid%ts_clw(n,i) ELSE !!! WRF-Solar diagnostics - WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,39(f13.5,1x))') & + WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,49(f13.5,1x))') & grid%id, grid%ts_hour(n,i), & grid%id_tsloc(i), ix, iy, & grid%ts_t(n,i), & @@ -877,7 +949,17 @@ SUBROUTINE write_ts( grid ) grid%ts_cbaseht_tot(n,i), & grid%ts_ctopht_tot(n,i), & grid%ts_clrnidx(n,i), & - grid%ts_sza(n,i) + grid%ts_sza(n,i), & + grid%ts_swdown(n,i), & + grid%ts_swddni(n,i), & + grid%ts_swddif(n,i), & + grid%ts_swdownc(n,i), & + grid%ts_swddnic(n,i), & + grid%ts_swdown2(n,i), & + grid%ts_swddni2(n,i), & + grid%ts_swddif2(n,i), & + grid%ts_swdownc2(n,i), & + grid%ts_swddnic2(n,i) END IF #else WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,7(f13.5,1x))') & diff --git a/test/em_real/examples.namelist b/test/em_real/examples.namelist index db55bd6e8c..c32f6043ed 100755 --- a/test/em_real/examples.namelist +++ b/test/em_real/examples.namelist @@ -492,13 +492,48 @@ Price, J. F., T. B. Sanford, and G. Z. Forristal, 1994: Forced stage response to -** Using U. Miami Forward Lagrangian trajectory calculation - (add it in namelist record &physics): -&domain +** Using the ACOM Forward Lagrangian trajectory calculation: + +&domains num_traj = 25, &physics traj_opt = 1, + dm_has_traj = .true., ..true., .true. + + +For domain #1, the file must "wrfinput_traj_d01" exist in the working directory. Similarly for domain 2, 3, etc. Each domain +has a separate file for a namelist. + +&traj_default + traj_def%start_time = '2000-01-24_12:00:00', + traj_def%stop_time = '2000-01-25_12:00:00', + traj_def%dyn_name(1:6) = 'p', 'T', 'z', 'u', 'v', 'w', + traj_def%hyd_name(1) = 'QVAPOR', +/ + +&traj_spec + traj_type%start_time = '2000-01-24_12:00:00', '2000-01-24_12:00:00', + '2000-01-24_12:00:00', '2000-01-24_12:00:00', + '2000-01-24_12:00:00', '2000-01-24_12:00:00', + '2000-01-24_12:00:00', '2000-01-24_12:00:00', + '2000-01-24_12:00:00', '2000-01-24_12:00:00', + '2000-01-24_12:00:00', + traj_type%stop_time = '2000-01-25_12:00:00', '2000-01-25_12:00:00', + '2000-01-25_12:00:00', '2000-01-25_12:00:00', + '2000-01-25_12:00:00', '2000-01-25_12:00:00', + '2000-01-25_12:00:00', '2000-01-25_12:00:00', + '2000-01-25_12:00:00', '2000-01-25_12:00:00', + '2000-01-25_12:00:00', + traj_type%lev = 60., 60., 60., 60., 60., 60., 60., 60., 60., 60., 60., + traj_type%lon = -79.88470, -79.74551, -79.60422, -79.46072, + -79.31503, -79.16708, -79.01682, -78.86417, + -78.70911, -78.55151, -78.39142, + traj_type%lat = 29.18063, 29.70515, 30.23069, 30.75718, + 31.28461, 31.81292, 32.34208, 32.87204, + 33.40276, 33.93421, 34.46631, +/ + ** Vertical nesting diff --git a/var/da/da_minimisation/da_get_innov_vector.inc b/var/da/da_minimisation/da_get_innov_vector.inc index 792c3856dd..d3e60d2ca5 100644 --- a/var/da/da_minimisation/da_get_innov_vector.inc +++ b/var/da/da_minimisation/da_get_innov_vector.inc @@ -285,7 +285,7 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags) ! [3] Optionally rescale observation errors: !------------------------------------------------------------------------ - if (use_obs_errfac) call da_use_obs_errfac( iv) + if (use_obs_errfac .and. it == 1) call da_use_obs_errfac(iv) !------------------------------------------------------------------------ ! [4] Optionally add Gaussian noise to O, O-B: diff --git a/var/da/da_setup_structures/da_setup_flow_predictors_ep_format2.inc b/var/da/da_setup_structures/da_setup_flow_predictors_ep_format2.inc index f26c8b537e..7905516743 100644 --- a/var/da/da_setup_structures/da_setup_flow_predictors_ep_format2.inc +++ b/var/da/da_setup_structures/da_setup_flow_predictors_ep_format2.inc @@ -359,6 +359,7 @@ subroutine da_setup_flow_predictors_ep_format2( ix, jy, kz, ne, ep, its, ite, jt read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz) temp3d = temp3d_r4 end if + call wrf_dm_bcast_real(temp3d, ijk) te = ie + (it-1)*nens ep % ci(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * & temp3d(its:ite,jts:jte,kts:kte) @@ -393,6 +394,7 @@ subroutine da_setup_flow_predictors_ep_format2( ix, jy, kz, ne, ep, its, ite, jt read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz) temp3d = temp3d_r4 end if + call wrf_dm_bcast_real(temp3d, ijk) te = ie + (it-1)*nens ep % sn(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * & temp3d(its:ite,jts:jte,kts:kte) @@ -427,6 +429,7 @@ subroutine da_setup_flow_predictors_ep_format2( ix, jy, kz, ne, ep, its, ite, jt read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz) temp3d = temp3d_r4 end if + call wrf_dm_bcast_real(temp3d, ijk) te = ie + (it-1)*nens ep % gr(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * & temp3d(its:ite,jts:jte,kts:kte) diff --git a/wrftladj/module_diffusion_em_ad.F b/wrftladj/module_diffusion_em_ad.F index 10f3dbc3c1..3e9d9e404d 100644 --- a/wrftladj/module_diffusion_em_ad.F +++ b/wrftladj/module_diffusion_em_ad.F @@ -6711,11 +6711,20 @@ SUBROUTINE a_tke_km(config_flags,xkmh,a_xkmh,xkmv,a_xkmv,xkhh,a_xkhh,xkhv, & !LPB[12] c_k = config_flags%c_k - tke_seed = tke_seed_value + tke_seed = 0. !LPB[13] - if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. & - (config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0. + IF (config_flags%isfflx .eq. 0) THEN + IF ((config_flags%diff_opt .eq. 2) .and. (config_flags%bl_pbl_physics .eq. 0)) THEN + IF( (config_flags%tke_drag_coefficient .lt. epsilon) .and. & + (config_flags%tke_heat_flux .lt. epsilon) ) THEN + tke_seed = tke_seed_value + ENDIF + ELSE + !tke_drag_coefficient and tke_heat_flux are irrelevant here + tke_seed = tke_seed_value + ENDIF + ENDIF !LPB[14] DO j = j_start, j_end diff --git a/wrftladj/module_diffusion_em_tl.F b/wrftladj/module_diffusion_em_tl.F index b07a832579..ff4014ffa2 100644 --- a/wrftladj/module_diffusion_em_tl.F +++ b/wrftladj/module_diffusion_em_tl.F @@ -2645,10 +2645,19 @@ SUBROUTINE g_tke_km(config_flags,xkmh,g_xkmh,xkmv,g_xkmv,xkhh,g_xkhh, & g_c_k =0.0 c_k =config_flags%c_k - tke_seed =tke_seed_value - - if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. & - (config_flags%tke_heat_flux .gt. epsilon) ) tke_seed =0. + tke_seed = 0. + + IF (config_flags%isfflx .eq. 0) THEN + IF ((config_flags%diff_opt .eq. 2) .and. (config_flags%bl_pbl_physics .eq. 0)) THEN + IF( (config_flags%tke_drag_coefficient .lt. epsilon) .and. & + (config_flags%tke_heat_flux .lt. epsilon) ) THEN + tke_seed = tke_seed_value + ENDIF + ELSE + !tke_drag_coefficient and tke_heat_flux are irrelevant here + tke_seed = tke_seed_value + ENDIF + ENDIF DO j =j_start,j_end DO k =kts+1,ktf-1