From 90f4b65627e06a4bcc5d33385a01eb3dc0afb4f3 Mon Sep 17 00:00:00 2001 From: joeolson42 Date: Sat, 13 May 2023 01:47:52 +0000 Subject: [PATCH 1/6] temporary patch for snow mixing which causes numerical instabilities --- physics/mynnedmf_wrapper.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/physics/mynnedmf_wrapper.F90 b/physics/mynnedmf_wrapper.F90 index 83a73e6b3..3c7de235f 100644 --- a/physics/mynnedmf_wrapper.F90 +++ b/physics/mynnedmf_wrapper.F90 @@ -418,7 +418,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .true. + FLAG_QS = .false. FLAG_QNC= .true. FLAG_QNWFA= .true. FLAG_QNIFA= .true. @@ -428,7 +428,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = qgrs_snow(i,k) + sqs(i,k) = 0. !qgrs_snow(i,k) qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -441,7 +441,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .true. + FLAG_QS = .false. FLAG_QNC= .true. FLAG_QNWFA= .false. FLAG_QNIFA= .false. @@ -451,7 +451,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = qgrs_snow(i,k) + sqs(i,k) = 0. !qgrs_snow(i,k) qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -464,7 +464,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .true. + FLAG_QS = .false. FLAG_QNC= .false. FLAG_QNWFA= .false. FLAG_QNIFA= .false. @@ -474,7 +474,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = qgrs_snow(i,k) + sqs(i,k) = 0. !qgrs_snow(i,k) qnc(i,k) = 0. qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -834,7 +834,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) !dqdt_ozone(i,k) = 0.0 dqdt_water_aer_num_conc(i,k) = RQNWFABLTEN(i,k) dqdt_ice_aer_num_conc(i,k) = RQNIFABLTEN(i,k) @@ -869,7 +869,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) enddo enddo if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then @@ -887,7 +887,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_liquid_cloud(i,k) = RQCBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) !dqdt_ozone(i,k) = 0.0 enddo enddo @@ -896,7 +896,7 @@ SUBROUTINE mynnedmf_wrapper_run( & call dtend_helper(100+ntcw,RQCBLTEN) call dtend_helper(100+ntiw,RQIBLTEN) call dtend_helper(100+ntinc,RQNIBLTEN) - call dtend_helper(100+ntsw,RQSBLTEN) + !call dtend_helper(100+ntsw,RQSBLTEN) endif !do k=1,levs ! do i=1,im From 17fc149be79d3046ae061198602ae9e2184ea54c Mon Sep 17 00:00:00 2001 From: joeolson42 Date: Sat, 13 May 2023 01:50:02 +0000 Subject: [PATCH 2/6] small change to table values of Leaf Area Index to match them better with the LAI monthly climatology that is used in HRRR --- physics/set_soilveg_ruc.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/set_soilveg_ruc.F90 b/physics/set_soilveg_ruc.F90 index db56fb8a4..8ce6023ff 100644 --- a/physics/set_soilveg_ruc.F90 +++ b/physics/set_soilveg_ruc.F90 @@ -190,9 +190,9 @@ subroutine set_soilveg_ruc(me,isot,ivet,nlunit,errmsg,errflg) ifortbl =(/1, 2, 4, 3, 2, 4, 4, 5, 5, 5, 4, 7, 9, 7, & & 9, 9, 9, 5, 5, 5, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0/) - laitbl =(/6.40, 6.48, 5.16, 3.31, 5.50, 3.66, 2.60, & - & 3.66, 3.66, 2.90, 5.72, 5.68, 1.00, 4.29, & - & 0.01, 0.75, 0.01, 3.35, 3.35, 3.35, 0.01, & + laitbl =(/2.80, 5.18, 4.16, 4.81, 4.20, 1.16, 0.90, & + & 3.00, 3.00, 1.10, 1.72, 2.58, 1.00, 2.29, & + & 0.01, 0.75, 0.01, 1.00, 1.00, 1.00, 0.01, & & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, & & 0.00, 0.00/) From b766e8102a8f441e0485e5059d810a7fd4a78f4c Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 19 May 2023 11:27:32 -0400 Subject: [PATCH 3/6] fix Cray compiler error having to do with binary/unary operator precedence standards --- physics/module_bl_mynn.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index 51a906faf..dcfdc1011 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -7521,7 +7521,7 @@ FUNCTION phim(zet) dummy_0=(1.-am_unst*zet) ! parentesis arg dummy_1=dummy_0**0.333333 ! y - dummy_11=-0.33333*am_unst*dummy_0**-0.6666667 ! dy/dzet + dummy_11=-0.33333*am_unst*dummy_0**(-0.6666667) ! dy/dzet dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.) ! f dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.) ! df/dzet dummy_3 = 0.57735*(2.*dummy_1+1.) ! g @@ -7573,7 +7573,7 @@ FUNCTION phih(zet) dummy_0=(1.-ah_unst*zet) ! parentesis arg dummy_1=dummy_0**0.333333 ! y - dummy_11=-0.33333*ah_unst*dummy_0**-0.6666667 ! dy/dzet + dummy_11=-0.33333*ah_unst*dummy_0**(-0.6666667) ! dy/dzet dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.) ! f dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.) ! df/dzet dummy_3 = 0.57735*(2.*dummy_1+1.) ! g From 1591a273fa085c443a50214e7aaf7cd673633c58 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 19 May 2023 11:53:15 -0400 Subject: [PATCH 4/6] add Anders Jensen to Thompson MP CODEOWNERS list --- CODEOWNERS | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CODEOWNERS b/CODEOWNERS index 4b7e45310..189fabd95 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -116,7 +116,7 @@ physics/module_gfdl_cloud_microphys.* @RuiyuSun physics/module_MP_FER_HIRES.* @ericaligo-NOAA @grantfirl @Qingfu-Liu @dustinswales physics/module_mp_nssl_2mom.F90 @grantfirl @Qingfu-Liu @dustinswales physics/module_mp_radar.* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales -physics/module_mp_thompson* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales +physics/module_mp_thompson* @gthompsnWRF @RuiyuSun @AndersJensen-NOAA @grantfirl @Qingfu-Liu @dustinswales physics/module_nst* @XuLi-NOAA @grantfirl @Qingfu-Liu @dustinswales physics/module_sf_exchcoef.f90 @grantfirl @Qingfu-Liu @dustinswales physics/module_SF_JSFC.F90 @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales @@ -126,7 +126,7 @@ physics/module_soil_pre.* @tanyasmirnova physics/moninshoc.* @SMoorthi-emc @grantfirl @Qingfu-Liu @dustinswales physics/mp_fer_hires.* @ericaligo-NOAA @grantfirl @Qingfu-Liu @dustinswales physics/mp_nssl.* @grantfirl @Qingfu-Liu @dustinswales -physics/mp_thompson* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales +physics/mp_thompson* @gthompsnWRF @RuiyuSun @AndersJensen-NOAA @grantfirl @Qingfu-Liu @dustinswales physics/multi_gases.F90 @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales physics/myjpbl_wrapper.* @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales physics/myjsfc_wrapper.* @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales From 5ffa961ae46c8901658ff76790c24fc3de10f238 Mon Sep 17 00:00:00 2001 From: joeolson42 Date: Fri, 26 May 2023 16:34:20 +0000 Subject: [PATCH 5/6] remove benign bug from smoke/dust mixing when nchem /= ndvel --- physics/module_bl_mynn.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index dcfdc1011..ec6b5700d 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -1003,9 +1003,8 @@ SUBROUTINE mynn_bl_driver( & if ( mix_chem ) then do ic = 1,ndvel vd1(ic) = vdep(i,ic) ! dry deposition velocity - chem1(kts,ic) = chem3d(i,kts,ic) enddo - do k = kts+1,kte + do k = kts,kte do ic = 1,nchem chem1(k,ic) = chem3d(i,k,ic) enddo @@ -1013,9 +1012,8 @@ SUBROUTINE mynn_bl_driver( & else do ic = 1,ndvel vd1(ic) = 0. ! dry deposition velocity - chem1(kts,ic) = 0. enddo - do k = kts+1,kte + do k = kts,kte do ic = 1,nchem chem1(k,ic) = 0. enddo From 3a1c48dfaf227ffd061cf5ccc81a7cccef1bb9d6 Mon Sep 17 00:00:00 2001 From: joeolson42 Date: Fri, 26 May 2023 18:01:46 +0000 Subject: [PATCH 6/6] changes to the irrigation algorithm in RUC to use real-time vegetation fraction instead of LAI to define the growing phase of the crops --- physics/module_sf_ruclsm.F90 | 96 +++++++++++++++--------------------- 1 file changed, 41 insertions(+), 55 deletions(-) diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index 850e3ee5e..6294bc068 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -415,7 +415,7 @@ SUBROUTINE LSMRUC(xlat,xlon, & curat, & INFILTRP real (kind_phys) :: cq,r61,r273,arp,brp,x,evs,eis - real (kind_phys) :: cropsm + real (kind_phys) :: cropfr, cropsm, newsm, factor real (kind_phys) :: meltfactor, ac,as, wb,rovcp INTEGER :: NROOT @@ -445,8 +445,8 @@ SUBROUTINE LSMRUC(xlat,xlon, & NDDZS=2*(nzs-2) !-- - testptlat = 48.7074_kind_phys !39.958 !42.05 !39.0 !74.12 !29.5 - testptlon = 289.03_kind_phys !271.622 !286.75 !280.6 !164.0 !283.0 + testptlat = 35.55 !48.7074_kind_phys !39.958 !42.05 !39.0 !74.12 !29.5 + testptlon = 278.66 !289.03_kind_phys !271.622 !286.75 !280.6 !164.0 !283.0 !-- @@ -983,63 +983,49 @@ SUBROUTINE LSMRUC(xlat,xlon, & ! Fraction of cropland category in the grid box should not have soil moisture below ! wilting point during the growing season. -! Let's keep soil moisture 20% above wilting point for the fraction of grid box under -! croplands. +! Let's keep soil moisture 5% above wilting point for the crop fraction of grid box. ! This change violates LSM moisture budget, but ! can be considered as a compensation for irrigation not included into LSM. -!tgs - "irrigation" uses fractional landuse, therefore mosaic_lu=1. +! "Irigation" could be applied when landuse fractional information +! is available and mosaic_lu=1. if(mosaic_lu == 1) then - IF (lufrac(crop) > zero .and. lai(i,j) > 1.1_kind_phys) THEN - ! cropland - do k=1,nroot - cropsm=1.1_kind_phys*wilt - qmin - if(soilm1d(k) < cropsm*lufrac(crop)) then - IF (debug_print ) THEN - if (abs(xlat-testptlat).lt.0.2 .and. & - abs(xlon-testptlon).lt.0.2)then -print * ,'Soil moisture is below wilting in cropland category at time step',ktau - print*,' lat,lon=',xlat,xlon & - ,'lufrac(crop),k,soilm1d(k),wilt,cropsm', & - lufrac(crop),k,soilm1d(k),wilt,cropsm - endif - ENDIF - soilm1d(k) = cropsm*lufrac(crop) - IF (debug_print ) THEN - if (abs(xlat-testptlat).lt.0.2 .and. & - abs(xlon-testptlon).lt.0.2)then - print*,' lat,lon=',xlat,xlon - print * ,'Added soil water to cropland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k) - endif - ENDIF - endif - enddo + ! greenness factor: between 0 for min greenness and 1 for max greenness. + factor = max(zero,min(one,(vegfra(i,j)-shdmin(i,j))/max(one,(shdmax(i,j)-shdmin(i,j))))) + if (abs(xlat-testptlat).lt.0.1 .and. & + abs(xlon-testptlon).lt.0.1)then + print *,' lat,lon=',xlat,xlon,' factor=',factor + endif - ELSEIF (ivgtyp(i,j) == natural .and. lai(i,j) > 0.7) THEN - ! grassland: assume that 40% of grassland is irrigated cropland - do k=1,nroot - cropsm=1.2_kind_phys*wilt - qmin - if(soilm1d(k) < cropsm*lufrac(natural)*0.4) then - IF (debug_print ) THEN - if (abs(xlat-testptlat).lt.0.2 .and. & - abs(xlon-testptlon).lt.0.2)then -print * ,'Soil moisture is below wilting in mixed grassland/cropland category at time step',ktau - print*,' lat,lon=',xlat,xlon, & - 'lufrac(natural),k,soilm1d(k),wilt', & - lufrac(natural),k,soilm1d(k),wilt - endif - ENDIF - soilm1d(k) = cropsm * lufrac(natural)*0.4_kind_phys + if((ivgtyp(i,j) == natural .or. ivgtyp(i,j) == crop) .and. factor > 0.75) then + ! cropland or grassland, apply irrigation during the growing seaspon when fraction + ! of greenness is > 0.75. - IF (debug_print ) THEN - if (abs(xlat-testptlat).lt.0.2 .and. & - abs(xlon-testptlon).lt.0.2)then - print*,' lat,lon=',xlat,xlon - print * ,'Added soil water to grassland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k) - endif - ENDIF - endif + do k=1,nroot + cropsm=1.05_kind_phys*wilt - qmin + cropfr = min(one,lufrac(crop) + 0.4*lufrac(natural)) ! assume that 40% of natural is cropland + newsm = cropsm*cropfr + (1.-cropfr)*soilm1d(k) + if(soilm1d(k) < newsm) then + IF (debug_print ) THEN + if (abs(xlat-testptlat).lt.0.1 .and. & + abs(xlon-testptlon).lt.0.1)then + print * ,'Soil moisture is below wilting in cropland areas at time step',ktau + print * ,' lat,lon=',xlat,xlon + print * ,' lufrac=',lufrac,'factor=',factor & + ,'lai,ivgtyp,lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm,', & + lai(i,j),ivgtyp(i,j),lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm + endif + ENDIF + soilm1d(k) = newsm + IF (debug_print ) THEN + if (abs(xlat-testptlat).lt.0.1 .and. & + abs(xlon-testptlon).lt.0.1)then + print*,' lat,lon=',xlat,xlon + print * ,'Added soil water to cropland areas, k,soilm1d(k)',k,soilm1d(k) + endif + ENDIF + endif ! < cropsm enddo - ENDIF + endif ! crop endif ! mosaic_lu !*** DIAGNOSTICS @@ -6599,7 +6585,7 @@ SUBROUTINE TRANSF( debug_print, & TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID ENDIF !-- uncomment next line for non-linear root distribution -! TRANF(1)=part(1) + TRANF(1)=part(1) DO K=2,NROOT totliq=soiliqw(k)+qmin