From 57f14bef0b9801c6efec6f0e3a94826425c1b4ea Mon Sep 17 00:00:00 2001 From: Jerold Herwehe Date: Thu, 3 Nov 2022 20:29:22 -0400 Subject: [PATCH 1/3] Add the Pleim surface layer scheme to MPAS-A This EPA_PXsfclay commit will add the Pleim surface layer scheme (Pleim, 2006) to MPAS-A as a new "config_sfclayer_scheme" option called "sf_pxsfclay". The Pleim surface layer scheme code (module_sf_pxsfclay.F) was migrated from WRFV4.4.1 with only minor changes and has been generalized to work with either MPAS or WRF. Note that this implementation of the Pleim surface layer in MPAS-A does not yet support fractional sea ice, in contrast to the other available surface layer schemes. The Pleim surface layer works with any LSM or PBL scheme, but is normally used in conjunction with the PX LSM and the ACM2 PBL. Refer to the description in the WRF-ARW technical document available at https://opensky.ucar.edu/islandora/object/opensky:2898 for an overview of the physics behind this option. New file: src/core_atmosphere/physics/physics_wrf/module_sf_pxsfclay.F Modified files: Makefile src/core_atmosphere/Registry.xml src/core_atmosphere/physics/mpas_atmphys_control.F src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F src/core_atmosphere/physics/physics_wrf/Makefile These EPA_PXsfclay code changes to MPAS-A are based on the 22 August 2022 "develop" branch of MPAS v7.3. Reference: Pleim, J. E., 2006: A simple, efficient solution of flux-profile relationships in the atmospheric surface layer, J. Appl. Meteor. Climatol., 45, 341347, https://doi.org/10.1175/JAM2339.1. --- Makefile | 12 +- src/core_atmosphere/Registry.xml | 2 +- .../physics/mpas_atmphys_control.F | 13 +- .../physics/mpas_atmphys_driver_sfclayer.F | 70 +- .../physics/physics_wrf/Makefile | 1 + .../physics/physics_wrf/module_sf_pxsfclay.F | 622 ++++++++++++++++++ 6 files changed, 707 insertions(+), 13 deletions(-) create mode 100644 src/core_atmosphere/physics/physics_wrf/module_sf_pxsfclay.F diff --git a/Makefile b/Makefile index 1d31bb2843..326b36789c 100644 --- a/Makefile +++ b/Makefile @@ -360,11 +360,11 @@ intel-mpi: # BUILDTARGET Intel compiler suite with Intel MPI library "CC_SERIAL = icc" \ "CXX_SERIAL = icpc" \ "FFLAGS_PROMOTION = -real-size 64" \ - "FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \ - "CFLAGS_OPT = -O3" \ - "CXXFLAGS_OPT = -O3" \ + "FFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback -no-wrap-margin -convert big_endian -free -align array64byte" \ + "CFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ + "CXXFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback" \ + "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback -no-wrap-margin" \ "CFLAGS_DEBUG = -g -traceback" \ "CXXFLAGS_DEBUG = -g -traceback" \ "LDFLAGS_DEBUG = -g -fpe0 -traceback" \ @@ -779,7 +779,7 @@ ifneq "$(LAPACK)" "" endif RM = rm -f -CPP = cpp -P -traditional +CPP = ${CXX_SERIAL} -E # Modified for use with the Intel C++ compiler RANLIB = ranlib ifdef CORE @@ -976,7 +976,7 @@ ifeq "$(findstring clean, $(MAKECMDGOALS))" "clean" # CHECK FOR CLEAN TARGET override AUTOCLEAN=false endif # END OF CLEAN TARGET CHECK -VER=$(shell git describe --dirty 2> /dev/null) +VER="v7.3.develop.EPA_PXsfclay" # Hard-coded specific version identifier #override CPPFLAGS += -DMPAS_GIT_VERSION=$(VER) ifeq "$(findstring v, $(VER))" "v" diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 6378596797..348bc96134 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -2212,7 +2212,7 @@ + possible_values="`suite',`sf_monin_obukhov',`sf_mynn',`sf_pxsfclay',`off'"/> 1 + REGIME(I) = 1.0 + ZOLL = BR(I) * GZ1OZ0(I) / (1.0 - RICRITI * RICUT(I)) + PSIM(I) = 1.0 - BETAM - ZOLL + PSIH(I) = 1.0 - BETAH - ZOLL + + ELSE IF (BR(I) .GE. 0.0) THEN +! -----CLASS 2; STABLE: for 1 > Z/L >0 + REGIME(I) = 2.0 + ZOLL = BR(I) * GZ1OZ0(I) / (1.0 - RICRITI * BR(I)) + PSIM(I) = -BETAM * ZOLL + PSIH(I) = -BETAH * ZOLL + + ELSE +! ----- CLASS 3 or 4; UNSTABLE: +! ----- CLASS 4 IS FOR ACM NON-LOCAL CONVECTION (H/L < -3) + REGIME(I) = 3.0 + AM = 0.031 + 0.276 * ALOG(GZ1OZ0(I)) + AH = 0.04 + 0.355 * ALOG(GZ1OZ0(I)) + SQLNZZ0 = SQRT(GZ1OZ0(I)) + PSIM(I) = AM * ALOG(1.0 - BM * SQLNZZ0 * BR(I)) + PSIH(I) = AH * ALOG(1.0 - BH * SQLNZZ0 * BR(I)) + + ENDIF + ENDDO + +! -------- COMPUTE THE FRICTIONAL VELOCITY AND SURFACE FLUXES: + DO I = its,ite + DTG = THETA1(I) - THETAG(I) + PSIX = GZ1OZ0(I) - PSIM(I) + UST(I) = 0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX + USTM(I) = UST(I) + +! ------- OVER WATER, ALTER ROUGHNESS LENGTH (Z0) ACCORDING TO WIND (UST). + IF ((XLAND(I)-1.5) .GE. 0.0) THEN + ZNT(I) = CZO * USTM(I) * USTM(I) / G + OZO + GZ1OZ0(I) = ALOG(ZA(I) / ZNT(I)) + PSIX = GZ1OZ0(I) - PSIM(I) + UST(I) = KARMAN * WSPD(I) / PSIX + USTM(I) = UST(I) + ENDIF + + RA(I) = PR0 * (GZ1OZ0(I) - PSIH(I)) / (KARMAN * UST(I)) + RBH = 5.0 / UST(I) + RBW = 4.503/UST(I) + CHS(I) = 1./(RA(I) + RBH) + CQS = 1./(RA(I) + RBW) + MOL(I) = DTG * CHS(I) / UST(I) + TSTV = (THETAV1(I) - TH0(I)) * CHS(I) / UST(I) + IF (ABS(TSTV) .LT. 1.E-5) TSTV = 1.E-5 + MOLENGTH(I) = THETAV1(I) * UST(I) * UST(I) / (KARMAN * G * TSTV) + +! ---Compute 2m surface exchange coefficients for heat and moisture + XMOL = MOLENGTH(I) + IF(MOLENGTH(I).GT.0.0) XMOL = AMAX1(MOLENGTH(I),2.0) + RMOL(I) = 1/XMOL + ZOL(I) = ZA(I)*RMOL(I) + ZOBOL = 1.5*RMOL(I) + Z10OL = 10.0*RMOL(I) + ZNTOL = ZNT(I)*RMOL(I) + IF(XMOL.LT.0.0) THEN + YNT = ( 1.0 - GAMAH * ZNTOL )**0.5 + YOB = ( 1.0 - GAMAH * ZOBOL )**0.5 + PSIH2 = 2. * ALOG((YOB+1.0)/(YNT+1.0)) + x1 = (1.0 - gamam * z10ol)**0.25 + x2 = (1.0 - gamam * zntol)**0.25 + psim10 = 2.0 * ALOG( (1.0+x1) / (1.0+x2) ) + & + ALOG( (1.0+x1*x1) / (1.0+x2*x2)) - & + 2.0 * ATAN(x1) + 2.0 * ATAN(x2) + ELSE + IF((ZOBOL-ZNTOL).LE.1.0) THEN + PSIH2 = -BETAH*(ZOBOL-ZNTOL) + ELSE + PSIH2 = 1.-BETAH-(ZOBOL-ZNTOL) + ENDIF + IF((Z10OL-ZNTOL).LE.1.0) THEN + PSIM10 = -BETAM*(Z10OL-ZNTOL) + ELSE + PSIM10 = 1.-BETAM-(Z10OL-ZNTOL) + ENDIF + ENDIF + G2OZ0 = ALOG(1.5 / ZNT(I)) + G10OZ0 = ALOG(10.0 / ZNT(I)) + RA2 = PR0 * (G2OZ0 - PSIH2) / (KARMAN * UST(I)) + CHS2(I) = 1.0/(RA2 + RBH) + CQS2(I) = 1.0/(RA2 + RBW) + U10(I) = US(I)*(G10OZ0-PSIM10)/PSIX + V10(I) = VS(I)*(G10OZ0-PSIM10)/PSIX + +! -----COMPUTE SURFACE HEAT AND MOIST FLUX: + FLHC(i) = CPM(I)*RHOX(I)*CHS(I) +#if defined(mpas) + FLQC(i) = RHOX(I)*CQS +#else + FLQC(i) = RHOX(I)*CQS*MAVAIL(I) +#endif + QFX(I) = FLQC(I)*(QSFC(I)-QV1D(I)) + QFX(I) = AMAX1(QFX(I),0.) + LH(I) = XLV*QFX(I) + IF(XLAND(I)-1.5.GT.0.)THEN + HFX(I)= -FLHC(I)*DTG + ELSEIF(XLAND(I)-1.5.LT.0.)THEN + HFX(I)= -FLHC(I)*DTG + HFX(I)= AMAX1(HFX(I),-250.) + ENDIF +#if defined(mpas) + TH2(I) = THETAG(I) - HFX(I) / (CPM(I)*RHOX(I)*CHS2(I)) + CPOT = (100./PSFC(I))**ROVCP + TA2(I) = TH2(I)/CPOT + QA2(I) = QSFC(I) - QFX(I) / (RHOX(I)*CQS2(I)) +#endif + ENDDO + + + END SUBROUTINE PXSFCLAY1D + +!==================================================================== + SUBROUTINE pxsfclayinit( allowed_to_read ) + + LOGICAL , INTENT(IN) :: allowed_to_read + INTEGER :: N + REAL :: ZOLN,X,Y + + + END SUBROUTINE pxsfclayinit + +!------------------------------------------------------------------- + +END MODULE module_sf_pxsfclay From d0165022e90b10bc600a2486989d19a733e72d36 Mon Sep 17 00:00:00 2001 From: Laura Fowler Date: Tue, 22 Oct 2024 19:22:06 -0600 Subject: [PATCH 2/3] * In the Pleim surface layer scheme sourcecode: -> restored Makefile to the default Makefile for v8.2.2. --- Makefile | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Makefile b/Makefile index 326b36789c..1d31bb2843 100644 --- a/Makefile +++ b/Makefile @@ -360,11 +360,11 @@ intel-mpi: # BUILDTARGET Intel compiler suite with Intel MPI library "CC_SERIAL = icc" \ "CXX_SERIAL = icpc" \ "FFLAGS_PROMOTION = -real-size 64" \ - "FFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback -no-wrap-margin -convert big_endian -free -align array64byte" \ - "CFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ - "CXXFLAGS_OPT = -O3 -xBROADWELL -fma -fp-model precise -traceback" \ + "FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \ + "CFLAGS_OPT = -O3" \ + "CXXFLAGS_OPT = -O3" \ "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback -no-wrap-margin" \ + "FFLAGS_DEBUG = -g -convert big_endian -free -CU -CB -check all -fpe0 -traceback" \ "CFLAGS_DEBUG = -g -traceback" \ "CXXFLAGS_DEBUG = -g -traceback" \ "LDFLAGS_DEBUG = -g -fpe0 -traceback" \ @@ -779,7 +779,7 @@ ifneq "$(LAPACK)" "" endif RM = rm -f -CPP = ${CXX_SERIAL} -E # Modified for use with the Intel C++ compiler +CPP = cpp -P -traditional RANLIB = ranlib ifdef CORE @@ -976,7 +976,7 @@ ifeq "$(findstring clean, $(MAKECMDGOALS))" "clean" # CHECK FOR CLEAN TARGET override AUTOCLEAN=false endif # END OF CLEAN TARGET CHECK -VER="v7.3.develop.EPA_PXsfclay" # Hard-coded specific version identifier +VER=$(shell git describe --dirty 2> /dev/null) #override CPPFLAGS += -DMPAS_GIT_VERSION=$(VER) ifeq "$(findstring v, $(VER))" "v" From 84a8e725c15c10202039b8142c121396fbee4194 Mon Sep 17 00:00:00 2001 From: Laura Fowler Date: Wed, 23 Oct 2024 10:23:57 -0600 Subject: [PATCH 3/3] * In the Pleim surface layer scheme sourcecode, added minor changes: -> in ./src/core_atmosphere/physics/mpas_atmphys_control.F, revised the "if" statement when checking the different options for config_sfclayer_scheme. -> in ./src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F, removed lines with empty "case("sf_pxsfclay")" options; reorganized calls to pxsfclay. --- .../physics/mpas_atmphys_control.F | 2 +- .../physics/mpas_atmphys_driver_sfclayer.F | 101 ++++++++---------- 2 files changed, 48 insertions(+), 55 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_control.F b/src/core_atmosphere/physics/mpas_atmphys_control.F index 72466d473e..0dcfa6030c 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_control.F +++ b/src/core_atmosphere/physics/mpas_atmphys_control.F @@ -278,7 +278,7 @@ subroutine physics_namelist_check(configs) write(mpas_err_message,'(A,A20)') 'illegal value for surface layer scheme: ', & trim(config_sfclayer_scheme) call physics_error_fatal(mpas_err_message) - else if(config_sfclayer_scheme .ne. 'sf_pxsfclay') then + else if(config_pbl_scheme == 'bl_mynn') then config_sfclayer_scheme = 'sf_mynn' elseif(config_pbl_scheme == 'bl_ysu') then diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F index febaf23ec6..9be30d1d93 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F @@ -224,8 +224,6 @@ subroutine allocate_sfclayer(configs) if(.not.allocated(ch_sea)) allocate(ch_sea(ims:ime,jms:jme)) endif - case("sf_pxsfclay") - case default end select sfclayer_select @@ -348,8 +346,6 @@ subroutine deallocate_sfclayer(configs) if(allocated(ch_sea)) deallocate(ch_sea) endif - case("sf_pxsfclay") - case default end select sfclayer_select @@ -593,8 +589,6 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) enddo enddo - case("sf_pxsfclay") - case default end select sfclayer_select @@ -799,8 +793,6 @@ subroutine sfclayer_to_MPAS(configs,sfc_input,diag_physics,its,ite) enddo endif - case("sf_pxsfclay") - case default end select sfclayer_select @@ -1092,57 +1084,58 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite call mpas_timer_stop('sf_mynn') case("sf_pxsfclay") - call mpas_timer_start('PX_sfclay') + call mpas_timer_start('sf_pxsfclay') call pxsfclay( & - p3d = pres_hyd_p , psfc = psfc_p , t3d = t_p , & - u3d = u_p , v3d = v_p , qv3d = qv_p , & - dz8w = dz_p , cp = cp , g = gravity , & - rovcp = rcp , R = R_d , xlv = xlv , & - chs = chs_p , chs2 = chs2_p , cqs2 = cqs2_p , & - cpm = cpm_p , znt = znt_p , ust = ust_p , & - pblh = hpbl_p , mavail = mavail_p , zol = zol_p , & - mol = mol_p , regime = regime_p , psim = psim_p , & - psih = psih_p , th3d = th_p , & - xland = xland_p , hfx = hfx_p , qfx = qfx_p , & - lh = lh_p , tsk = tsk_p , flhc = flhc_p , & - ta2 = t2m_p , th2 = th2m_p , qa2 = q2_p , & - flqc = flqc_p , qgh = qgh_p , qsfc = qsfc_p , & - rmol = rmol_p , u10 = u10_p , v10 = v10_p , & - gz1oz0 = gz1oz0_p , wspd = wspd_p , br = br_p , & - isfflx = isfflx , dx = dx_p , svp1 = svp1 , & - svp2 = svp2 , svp3 = svp3 , svpt0 = svpt0 , & - ep1 = ep_1 , ep2 = ep_2 , karman = karman , & + p3d = pres_hyd_p , psfc = psfc_p , t3d = t_p , & + u3d = u_p , v3d = v_p , qv3d = qv_p , & + dz8w = dz_p , cp = cp , g = gravity , & + rovcp = rcp , R = R_d , xlv = xlv , & + chs = chs_p , chs2 = chs2_p , cqs2 = cqs2_p , & + cpm = cpm_p , znt = znt_p , ust = ust_p , & + pblh = hpbl_p , mavail = mavail_p , zol = zol_p , & + mol = mol_p , regime = regime_p , psim = psim_p , & + psih = psih_p , th3d = th_p , xland = xland_p , & + hfx = hfx_p , qfx = qfx_p , lh = lh_p , & + tsk = tsk_p , flhc = flhc_p , ta2 = t2m_p , & + th2 = th2m_p , qa2 = q2_p , flqc = flqc_p , & + qgh = qgh_p , qsfc = qsfc_p , rmol = rmol_p , & + u10 = u10_p , v10 = v10_p , gz1oz0 = gz1oz0_p , & + wspd = wspd_p , br = br_p , isfflx = isfflx , & + dx = dx_p , svp1 = svp1 , svp2 = svp2 , & + svp3 = svp3 , svpt0 = svpt0 , ep1 = ep_1 , & + ep2 = ep_2 , karman = karman , & ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , & ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , & its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte & ) -! Note that fractional sea-ice is currently not supported by PX. -! if(config_frac_seaice) then -! call pxsfclay( & -! p3d = pres_hyd_p , psfc = psfc_p , t3d = t_p , & -! u3d = u_p , v3d = v_p , qv3d = qv_p , & -! dz8w = dz_p , cp = cp , g = gravity , & -! rovcp = rcp , R = R_d , xlv = xlv , & -! chs = chs_sea , chs2 = chs2_sea , cqs2 = cqs2_sea , & -! cpm = cpm_sea , znt = znt_sea , ust = ust_sea , & -! pblh = hpbl_p , mavail = mavail_sea , zol = zol_sea , & -! mol = mol_sea , regime = regime_sea , psim = psim_sea , & -! psih = psih_sea , th3d = th_p , & -! xland = xland_sea , hfx = hfx_sea , qfx = qfx_sea , & -! lh = lh_sea , tsk = tsk_sea , flhc = flhc_sea , & -! ta2 = t2m_sea , th2 = th2m_sea , qa2 = q2_sea , & -! flqc = flqc_sea , qgh = qgh_sea , qsfc = qsfc_sea , & -! rmol = rmol_sea , u10 = u10_sea , v10 = v10_sea , & -! gz1oz0 = gz1oz0_sea , wspd = wspd_sea , br = br_sea , & -! isfflx = isfflx , dx = dx_p , svp1 = svp1 , & -! svp2 = svp2 , svp3 = svp3 , svpt0 = svpt0 , & -! ep1 = ep_1 , ep2 = ep_2 , karman = karman , & -! ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , & -! ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , & -! its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte & -! ) -! endif - call mpas_timer_stop('PX_sfclay') + +! Note that fractional sea-ice is currently not supported by PX. +! if(config_frac_seaice) then +! call pxsfclay( & +! p3d = pres_hyd_p , psfc = psfc_p , t3d = t_p , & +! u3d = u_p , v3d = v_p , qv3d = qv_p , & +! dz8w = dz_p , cp = cp , g = gravity , & +! rovcp = rcp , R = R_d , xlv = xlv , & +! chs = chs_sea , chs2 = chs2_sea , cqs2 = cqs2_sea , & +! cpm = cpm_sea , znt = znt_sea , ust = ust_sea , & +! pblh = hpbl_p , mavail = mavail_sea , zol = zol_sea , & +! mol = mol_sea , regime = regime_sea , psim = psim_sea , & +! psih = psih_sea , th3d = th_p , xland = xland_sea , & +! hfx = hfx_sea , qfx = qfx_sea , lh = lh_sea , & +! tsk = tsk_sea , flhc = flhc_sea , ta2 = t2m_sea , & +! th2 = th2m_sea , qa2 = q2_sea , flqc = flqc_sea , & +! qgh = qgh_sea , qsfc = qsfc_sea , rmol = rmol_sea , & +! u10 = u10_sea , v10 = v10_sea , gz1oz0 = gz1oz0_sea , & +! wspd = wspd_sea , br = br_sea , isfflx = isfflx , & +! dx = dx_p , svp1 = svp1 , svp2 = svp2 , & +! svp3 = svp3 , svpt0 = svpt0 , ep1 = ep_1 , & +! ep2 = ep_2 , karman = karman , & +! ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , & +! ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , & +! its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte & +! ) +! endif + call mpas_timer_stop('sf_pxsfclay') case default