diff --git a/.gitmodules b/.gitmodules index f6c56f6c3f..9adff1253d 100644 --- a/.gitmodules +++ b/.gitmodules @@ -16,3 +16,6 @@ [submodule "src/core_atmosphere/physics/physics_noaa/SMOKE"] path = src/core_atmosphere/physics/physics_noaa/SMOKE url = https://github.com/NOAA-GSL/SMOKE.git +[submodule "src/core_atmosphere/physics/physics_noaa/MYNN-SFC"] + path = src/core_atmosphere/physics/physics_noaa/MYNN-SFC + url = https://github.com/NCAR/MYNN-SFC.git diff --git a/src/core_atmosphere/CMakeLists.txt b/src/core_atmosphere/CMakeLists.txt index ea4f57448e..4441d9732d 100644 --- a/src/core_atmosphere/CMakeLists.txt +++ b/src/core_atmosphere/CMakeLists.txt @@ -99,9 +99,6 @@ set(ATMOSPHERE_CORE_PHYSICS_WRF_SOURCES module_sf_mynn.F sf_mynn_pre.F sf_sfclayrev_pre.F - module_mynnsfclay_driver.F90 - module_mynnsfclay.F90 - module_mynnsfclay_pre.F90 ) list(TRANSFORM ATMOSPHERE_CORE_PHYSICS_WRF_SOURCES PREPEND physics/physics_wrf/) @@ -248,6 +245,31 @@ set(ATMOSPHERE_CORE_PHYSICS_MYNN-EDMF_SOURCES ) list(TRANSFORM ATMOSPHERE_CORE_PHYSICS_MYNN-EDMF_SOURCES PREPEND physics/physics_noaa/MYNN-EDMF/) +set(ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/physics/physics_noaa/MYNN-SFC) + +if(NOT EXISTS ${ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_DIR}) + set(PHYSICS_MYNN-SFC_REPO_URL "https://github.com/NCAR/MYNN-SFC.git") + execute_process(COMMAND git clone ${PHYSICS_MYNN-SFC_REPO_URL} ${ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_DIR} + RESULT_VARIABLE GIT_CLONE_RESULT + OUTPUT_VARIABLE GIT_CLONE_OUTPUT + ERROR_VARIABLE GIT_CLONE_ERROR) + if(NOT GIT_CLONE_RESULT EQUAL 0) + message(FATAL_ERROR "Git clone failed with error: ${GIT_CLONE_ERROR}") + endif() + +else() + message(STATUS "Directory ${DIR_TO_CHECK} already exists, skipping clone") +endif() + +set(ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_SOURCES + module_sf_mynnsfc_land.F90 + module_sf_mynnsfc_water.F90 + module_sf_mynnsfc_ice.F90 + MPAS/module_sf_mynnsfc_common.F90 + MPAS/module_sf_mynnsfc_driver.F90 +) +list(TRANSFORM ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_SOURCES PREPEND physics/physics_noaa/MYNN-SFC/) + set(ATMOSPHERE_CORE_PHYSICS_RUCLSM_DIR ${CMAKE_CURRENT_SOURCE_DIR}/physics/physics_noaa/RUCLSM) if(NOT EXISTS ${ATMOSPHERE_CORE_PHYSICS_RUCLSM_DIR}) @@ -489,6 +511,7 @@ add_library(core_atmosphere ${ATMOSPHERE_CORE_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_UGWP_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_TEMPO_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_MYNN-EDMF_SOURCES} + ${ATMOSPHERE_CORE_PHYSICS_MYNN-SFC_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_RUCLSM_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_GFL_SOURCES} ${ATMOSPHERE_CORE_PHYSICS_SMOKE_SOURCES} diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index d7d7f4bb5f..5397e8960d 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -2820,6 +2820,16 @@ description="configuration for surface layer-scheme" possible_values="`suite',`sf_monin_obukhov','sf_monin_obukhov_rev',`sf_mynn',`sf_mynnsfclay',`sf_myj', `off'"/> + + + + diff --git a/src/core_atmosphere/physics/Makefile b/src/core_atmosphere/physics/Makefile index 2fe1019b08..5f8e20c440 100644 --- a/src/core_atmosphere/physics/Makefile +++ b/src/core_atmosphere/physics/Makefile @@ -6,7 +6,7 @@ endif all: ./../tools/manage_externals/checkout_externals --externals ./../Externals.cfg - $(MAKE) lookup_tables core_physics_init core_physics_mmm core_microphysics core_UGWP core_RUCLSM core_mynnedmf core_GFL core_SMOKE core_physics_wrf core_physics_noahmp core_physics + $(MAKE) lookup_tables core_physics_init core_physics_mmm core_microphysics core_UGWP core_RUCLSM core_mynnedmf core_mynnsfc core_GFL core_SMOKE core_physics_wrf core_physics_noahmp core_physics dummy: echo "****** compiling physics ******" @@ -73,6 +73,9 @@ core_SMOKE: core_physics_init core_mynnedmf: core_physics_init core_physics_mmm (cd physics_noaa/MYNN-EDMF; cp ./MPAS/Makefile .; cp ./MPAS/module_bl_mynnedmf_driver.F90 .; cp ./MPAS/module_bl_mynnedmf_common.F90 .; $(MAKE) all) +core_mynnsfc: + (cd physics_noaa/MYNN-SFC; cp ./MPAS/Makefile .; cp ./MPAS/module_sf_mynnsfc_driver.F90 .; cp ./MPAS/module_sf_mynnsfc_common.F90 .; $(MAKE) all) + core_RUCLSM: core_physics_init core_physics_mmm (cd physics_noaa/RUCLSM; cp ./MPAS/Makefile .; cp ./MPAS/module_ruc_land.F90 .; cp ./MPAS/module_ruc_ice.F90 .; cp ./MPAS/module_sf_sfcdiags_ruclsm.F .; $(MAKE) all COREDEF="$(COREDEF)") @@ -92,7 +95,7 @@ core_physics_noahmp: core_physics_init: $(OBJS_init) -core_physics: core_microphysics core_mynnedmf core_RUCLSM core_GFL core_SMOKE core_physics_wrf core_physics_noahmp +core_physics: core_microphysics core_mynnedmf core_mynnsfc core_RUCLSM core_GFL core_SMOKE core_physics_wrf core_physics_noahmp ($(MAKE) phys_interface COREDEF="$(COREDEF)") ar -ru libphys.a $(OBJS_init) $(OBJS) ($(MAKE) -C ./physics_mmm -f Makefile.mpas physics_mmm_lib) @@ -290,6 +293,7 @@ clean: ( if [ -d physics_noaa/UGWP ]; then cd physics_noaa/UGWP; $(MAKE) clean; fi ) ( if [ -d physics_noaa/RUCLSM ]; then cd physics_noaa/RUCLSM; $(MAKE) -f MPAS/Makefile clean; fi ) ( if [ -d physics_noaa/MYNN-EDMF ]; then cd physics_noaa/MYNN-EDMF; cp ./MPAS/Makefile .; $(MAKE) clean; fi ) + ( if [ -d physics_noaa/MYNN-SFC ]; then cd physics_noaa/MYNN-SFC; cp ./MPAS/Makefile .; $(MAKE) clean; fi ) ( if [ -d physics_noaa/GFL ]; then cd physics_noaa/GFL; cp ./MPAS/Makefile .; $(MAKE) clean; fi ) ( if [ -d physics_noaa/SMOKE ]; then cd physics_noaa/SMOKE; cp ./MPAS/Makefile .; $(MAKE) clean; fi ) @# Certain systems with intel compilers generate *.i files @@ -300,16 +304,16 @@ clean: $(RM) $@ $*.mod ifeq "$(GEN_F90)" "true" $(CPP) $(CPPFLAGS) $(COREDEF) $(HYDROSTATIC) $(CPPINCLUDES) $< > $*.f90 - $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noahmp -I./physics_noahmp/utility -I./physics_noahmp/drivers/mpas -I./physics_noahmp/src -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 + $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noahmp -I./physics_noahmp/utility -I./physics_noahmp/drivers/mpas -I./physics_noahmp/src -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/MYNN-SFC -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 else - $(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noahmp -I./physics_noahmp/utility -I./physics_noahmp/drivers/mpas -I./physics_noahmp/src -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 + $(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noahmp -I./physics_noahmp/utility -I./physics_noahmp/drivers/mpas -I./physics_noahmp/src -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/MYNN-SFC -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 endif .F90.o: $(RM) $@ $*.mod ifeq "$(GEN_F90)" "true" $(CPP) $(CPPFLAGS) $(COREDEF) $(HYDROSTATIC) $(CPPINCLUDES) $< > $*.f90 - $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 + $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/MYNN-SFC -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 else - $(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 + $(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./physics_noaa/TEMPO -I./physics_noaa/UGWP -I./physics_noaa/MYNN-EDMF -I./physics_noaa/MYNN-SFC -I./physics_noaa/RUCLSM -I./physics_noaa/GFL -I./physics_noaa/SMOKE -I.. -I../../framework -I../../external/esmf_time_f90 endif diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F b/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F index 8c9b951474..98dbb20cda 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F @@ -590,7 +590,11 @@ subroutine seaice_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) hfx_p(i,j) = xice_p(i,j)*hfx_p(i,j) + (1._RKIND-xice_p(i,j))*hfx_sea(i,j) lh_p(i,j) = xice_p(i,j)*lh_p(i,j) + (1._RKIND-xice_p(i,j))*lh_sea(i,j) qfx_p(i,j) = xice_p(i,j)*qfx_p(i,j) + (1._RKIND-xice_p(i,j))*qfx_sea(i,j) - qgh_p(i,j) = xice_p(i,j)*qgh_p(i,j) + (1._RKIND-xice_p(i,j))*qgh_sea(i,j) + if(trim(sfclayer_scheme) .ne. "sf_mynnsfclay") then + qgh_p(i,j) = xice_p(i,j)*qgh_p(i,j) + (1._RKIND-xice_p(i,j))*qgh_sea(i,j) + else + qgh_p(i,j) = qgh_p(i,j) + endif qsfc_p(i,j) = xice_p(i,j)*qsfc_p(i,j) + (1._RKIND-xice_p(i,j))*qsfc_sea(i,j) tsk_p(i,j) = xice_p(i,j)*tsk_p(i,j) + (1._RKIND-xice_p(i,j))*tsk_sea(i,j) sfc_albedo_p(i,j) = xice_p(i,j)*sfc_albedo_p(i,j) + (1._RKIND-xice_p(i,j))*0.08_RKIND diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F index 7e9aa600d8..1d828e5024 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F @@ -15,11 +15,11 @@ module mpas_atmphys_driver_sfclayer use mpas_atmphys_vars use module_sf_mynn,only: sfclay_mynn - use module_mynnsfclay_driver,only:mynnsfclay_driver + use module_sf_mynnsfc_driver,only: mynnsfc_driver use module_sf_sfclay use module_sf_sfclayrev,only: sfclayrev use sf_mynn,only: sf_mynn_init - use module_mynnsfclay,only: mynnsfclay_init + use module_sf_mynnsfc_driver,only: mynnsfc_init use sf_sfclayrev,only: sf_sfclayrev_init use module_sf_myjsfc @@ -678,7 +678,6 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) case("sf_mynnsfclay") !input variables: - call mpas_pool_get_array(diag_physics,'qcg' ,qcg ) call mpas_pool_get_array(sfc_input ,'snowh',snowh) !inout variables: call mpas_pool_get_array(diag_physics,'ch' ,ch ) @@ -688,7 +687,6 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) do i = its,ite !input variables: snowh_p(i,j) = snowh(i) - qcg_p(i,j) = qcg(i) !inout variables: ch_p(i,j) = ch(i) cqs_p(i,j) = cqs(i) @@ -1097,7 +1095,7 @@ subroutine init_sfclayer(configs,diag_physics,sfc_input) call sf_mynn_init(errmsg,errflg) case("sf_mynnsfclay") - call mynnsfclay_init(errmsg,errflg) + call mynnsfc_init(allowed_to_read,errmsg,errflg) case("sf_myj") call mpas_pool_get_config(configs,'config_do_restart' ,config_do_restart ) @@ -1164,6 +1162,8 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite logical,pointer:: config_do_restart,config_frac_seaice logical,pointer:: config_do_DAcycling character(len=StrKIND),pointer:: sfclayer_scheme + character(len=StrKIND),pointer:: lsm_scheme + integer,pointer:: sf_mynn_sfcflux_land,sf_mynn_sfcflux_water real(kind=RKIND),dimension(:),pointer:: areaCell !local variables: @@ -1390,32 +1390,33 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite case("sf_mynnsfclay") call mpas_timer_start('sf_mynnsfclay') - call mpas_log_write('--- enter subroutine mynnsfclay_driver:') - call mynnsfclay_driver( & - p3d = pres_hyd_p , pi3d = pi_p , psfcpa = psfc_p , & - th3d = th_p , t3d = t_p , u3d = u_p , & - v3d = v_p , qv3d = qv_p , qc3d = qc_p , & - rho3d = rho_p , dz8w = dz_p , cp = cp , & - g = gravity , rovcp = rcp , R = R_d , & - xlv = xlv , chs = chs_p , chs2 = chs2_p , & + call mpas_log_write('--- enter subroutine mynnsfc_driver:') + call mpas_pool_get_config(configs,'config_mynn_sfcflux_land' ,sf_mynn_sfcflux_land ) + call mpas_pool_get_config(configs,'config_mynn_sfcflux_water',sf_mynn_sfcflux_water) + call mpas_pool_get_config(configs,'config_lsm_scheme' ,lsm_scheme ) + call mynnsfc_driver( & + p3d = pres_hyd_p , th3d = th_p , t3d = t_p , & + u3d = u_p , v3d = v_p , qv3d = qv_p , & + rho3d = rho_p , dz8w = dz_p , & + chs = chs_p , chs2 = chs2_p , psfcpa = psfc_p , & cqs = cqs_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 , & + zol = zol_p , mol = mol_p , & psim = psim_p , psih = psih_p , xland = xland_p , & hfx = hfx_p , qfx = qfx_p , lh = lh_p , & tsk = tsk_p , flhc = flhc_p , flqc = flqc_p , & - qgh = qgh_p , qsfc = qsfc_p , rmol = rmol_p , & + qsfc = qsfc_p , rmol = rmol_p , & u10 = u10_p , v10 = v10_p , th2 = th2m_p , & t2 = t2m_p , q2 = q2_p , snowh = snowh_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 , & + isfflx = isfflx , dx = dx_p , & ustm = ustm_p , ck = ck_p , cka = cka_p , & cd = cd_p , cda = cda_p , ch = ch_p , & - qcg = qcg_p , spp_pbl = spp_pbl , isftcflx = isftcflx , & - iz0tlnd = iz0tlnd , & + spp_pbl = spp_pbl , & + sf_mynn_sfcflux_land = sf_mynn_sfcflux_land , & + sf_mynn_sfcflux_water = sf_mynn_sfcflux_water , & + flagc_lsm = lsm_scheme , & itimestep= itimestep , initflag = initflag , & restart = config_do_restart , & cycling = config_do_DAcycling , & @@ -1426,31 +1427,29 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite ) if(config_frac_seaice) then - call mynnsfclay_driver( & - p3d = pres_hyd_p , pi3d = pi_p , psfcpa = psfc_p , & - th3d = th_p , t3d = t_p , u3d = u_p , & - v3d = v_p , qv3d = qv_p , qc3d = qc_p , & - rho3d = rho_p , dz8w = dz_p , cp = cp , & - g = gravity , rovcp = rcp , R = R_d , & - xlv = xlv , chs = chs_sea , chs2 = chs2_sea , & + call mynnsfc_driver( & + p3d = pres_hyd_p , th3d = th_p , t3d = t_p , & + u3d = u_p , v3d = v_p , qv3d = qv_p , & + rho3d = rho_p , dz8w = dz_p , & + chs = chs_sea , chs2 = chs2_sea , psfcpa = psfc_p , & cqs = cqs_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 , & + zol = zol_sea , mol = mol_sea , & psim = psim_sea , psih = psih_sea , xland = xland_sea , & hfx = hfx_sea , qfx = qfx_sea , lh = lh_sea , & tsk = tsk_sea , flhc = flhc_sea , flqc = flqc_sea , & - qgh = qgh_sea , qsfc = qsfc_sea , rmol = rmol_sea , & + qsfc = qsfc_sea , rmol = rmol_sea , & u10 = u10_sea , v10 = v10_sea , th2 = th2m_sea , & t2 = t2m_sea , q2 = q2_sea , snowh = snowh_p , & 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 , & + isfflx = isfflx , dx = dx_p , & ustm = ustm_sea , ck = ck_sea , cka = cka_sea , & cd = cd_sea , cda = cda_sea , ch = ch_sea , & - qcg = qcg_p , spp_pbl = spp_pbl , isftcflx = isftcflx , & - iz0tlnd = iz0tlnd , & + spp_pbl = spp_pbl , & + sf_mynn_sfcflux_land = sf_mynn_sfcflux_land , & + sf_mynn_sfcflux_water = sf_mynn_sfcflux_water , & + flagc_lsm = lsm_scheme , & itimestep= itimestep , initflag = initflag , & restart = config_do_restart , & cycling = config_do_DAcycling , & diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index b379086243..eff23474d5 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -758,7 +758,7 @@ module mpas_atmphys_vars lakedepth_p, &!depth of lakes needed to run the lake model physics. lakemask_p !mask needed to detect the location of lakes to run the lake model physics. -!... arrays only in mynn surface layer scheme (module_sf_mynn.F): +!... arrays only in mynn surface layer scheme (module_sf_mynn.F): (note qcg is not needed for mynn sfclayer submodule) real(kind=RKIND),dimension(:,:),allocatable:: & ch_p, &!surface exchange coeff for heat [m/s] qcg_p !cloud water mixing ratio at the ground surface [kg/kg] diff --git a/src/core_atmosphere/physics/physics_noaa/MYNN-SFC b/src/core_atmosphere/physics/physics_noaa/MYNN-SFC new file mode 160000 index 0000000000..6f267c89b5 --- /dev/null +++ b/src/core_atmosphere/physics/physics_noaa/MYNN-SFC @@ -0,0 +1 @@ +Subproject commit 6f267c89b59041aeacf36d1d04a10b79e74cc781 diff --git a/src/core_atmosphere/physics/physics_wrf/Makefile b/src/core_atmosphere/physics/physics_wrf/Makefile index d8e17263d2..0f39483dbe 100644 --- a/src/core_atmosphere/physics/physics_wrf/Makefile +++ b/src/core_atmosphere/physics/physics_wrf/Makefile @@ -52,9 +52,6 @@ OBJS = \ module_sf_urban.o \ bl_mynn_post.o \ bl_mynn_pre.o \ - module_mynnsfclay.o \ - module_mynnsfclay_pre.o \ - module_mynnsfclay_driver.o \ cu_ntiedtke_post.o \ cu_ntiedtke_pre.o \ sf_mynn_pre.o \ @@ -71,10 +68,6 @@ module_bl_mynn.o: \ bl_mynn_post.o \ bl_mynn_pre.o -module_mynnsfclay_driver.o: \ - module_mynnsfclay.o \ - module_mynnsfclay_pre.o - module_cam_support.o: \ module_cam_shr_kind_mod.o diff --git a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay.F90 b/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay.F90 deleted file mode 100644 index 60d3c67686..0000000000 --- a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay.F90 +++ /dev/null @@ -1,2279 +0,0 @@ -!================================================================================================================= - module module_mynnsfclay - -!------------------------------------------------------------------- -!Modifications implemented by Joseph Olson NOAA/GSL -!The following overviews the current state of this scheme:: -! -! BOTH LAND AND WATER: -!1) Calculation of stability parameter (z/L) taken from Li et al. (2010 BLM) -! for first iteration of first time step; afterwards, exact calculation -! using a brute force iterative method described in Olson et al. (2021 NOAA -! Tech memorandum). This method replaces the iterative technique used in -! module_sf_sfclayrev.F (Jimenez et al. 2013) with mods. Either technique -! gives about the same result. The former technique is retained in this -! module (commented out) for potential subsequent benchmarking. -!2) Fixed isflux=0 option to turn off scalar fluxes, but keep momentum -! fluxes for idealized studies (credit: Anna Fitch). -!3) Kinematic viscosity varies with temperature according to Andreas (1989). -!4) Uses the blended Monin-Obukhov flux-profile relationships COARE (Fairall -! et al 2003) for the unstable regime (a blended mix of Dyer-Hicks 1974 and -! Grachev et al (2000). Uses Cheng and Brutsaert (2005) for stable conditions. -!5) The following overviews the namelist variables that control the -! aerodynamic roughness lengths (over water) and the thermal and moisture -! roughness lengths (defaults are recommended): -! -! LAND only: -! "iz0tlnd" namelist option is used to select the following options: -! (default) =0: Zilitinkevich (1995); Czil now set to 0.095 -! =1: Czil_new (modified according to Chen & Zhang 2008) -! =2: Modified Yang et al (2002, 2008) - generalized for all landuse -! =3: constant zt = z0/7.4 (original form; Garratt 1992) -! -! WATER only: -! "isftcflx" namelist option is used to select the following options: -! (default) =0: z0, zt, and zq from the COARE algorithm. Set COARE_OPT (below) to -! 3.0 (Fairall et al. 2003, default) -! 3.5 (Edson et al 2013) - now with bug fix (Edson et al. 2014, JPO) -! =1: z0 from Davis et al (2008), zt & zq from COARE 3.0/3.5 -! =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) -! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5 -! -! SNOW/ICE only: -! Andreas (2002) snow/ice parameterization for thermal and -! moisture roughness is used over all gridpoints with snow deeper than -! 0.1 m. This algorithm calculates a z0 for snow (Andreas et al. 2005, BLM), -! which is only used as part of the thermal and moisture roughness -! length calculation, not to directly impact the surface winds. -! -! Misc: -!1) Added a more elaborate diagnostic for u10 & V10 for high vertical resolution -! model configurations but for most model configurations with depth of -! the lowest half-model level near 10 m, a neutral-log diagnostic is used. -! -!2) Option to activate stochastic parameter perturbations (SPP), which -! perturb z0, zt, and zq, along with many other parameters in the MYNN- -! EDMF scheme. -! -!NOTE: This code was primarily tested in combination with the RUC LSM. -! Performance with the Noah (or other) LSM is relatively unknown. -!------------------------------------------------------------------- - use ccpp_kind_types,only: kind_phys -! use mynn_shared,only: esat_blend,qsat_blend,xl_blend - - implicit none - private - public:: mynnsfclay, & - mynnsfclay_init, & - mynnsfclay_finalize - - - logical,parameter:: debug_code = .false. - integer,parameter:: psi_opt = 0 ! 0 = stable: Cheng and Brustaert - ! unstable: blended COARE - ! 1 = GFS - real,parameter:: wmin = 0.1 - real,parameter:: vconvc = 1.25 - real,parameter:: snowz0 = 0.011 - real,parameter:: coare_opt = 3.5 ! 3.0 or 3.5 - !For debugging purposes: - - real,dimension(0:1000),save:: psim_stab,psim_unstab, & - psih_stab,psih_unstab - - - contains - - -!================================================================================================================= -!>\section arg_table_sf_mynn_init -!!\html\include sf_mynn_init.html -!! - subroutine mynnsfclay_init(errmsg,errflg) -!================================================================================================================= - -!--- output arguments: - character(len=*),intent(out):: errmsg - integer,intent(out):: errflg - -!----------------------------------------------------------------------------------------------------------------- - - call psi_init(psi_opt) - - errmsg = ' ' - errflg = 0 - - end subroutine mynnsfclay_init - -!================================================================================================================= -!>\section arg_table_sf_mynn_finalize -!!\html\include sf_mynn_finalize.html -!! - subroutine mynnsfclay_finalize(errmsg,errflg) -!================================================================================================================= - -!--- output arguments: - character(len=*),intent(out):: errmsg - integer,intent(out):: errflg - -!----------------------------------------------------------------------------------------------------------------- - - errmsg = ' ' - errflg = 0 - - end subroutine mynnsfclay_finalize - -!================================================================================================================= -!>\section arg_table_sf_mynn_run -!!\html\include sf_mynn_run.html -!! - subroutine mynnsfclay( & - u1d,v1d,t1d,qv1d,p1d,dz8w1d,rho1d, & - u1d2,v1d2,dz2w1d,cp,g,rovcp,r,xlv, & - psfcpa,chs,chs2,cqs,cqs2, & - cpm,pblh,rmol, & - znt,ust,mavail,zol,mol,regime,psim, & - psih,xland,hfx,qfx,tsk,u10,v10,th2, & - t2,q2,flhc,flqc,snowh,qgh,qsfc,lh, & - gz1oz0,wspd,br,isfflx,dx,svp1,svp2, & - svp3,svpt0,ep1,ep2,karman,ch,qcg, & - itimestep,wstar,qstar,ustm,ck,cka, & - cd,cda,spp_pbl,rstoch1d,isftcflx, & - iz0tlnd,its,ite,restart_or_cycle, & - errmsg,errflg & - ) - implicit none -!================================================================================================================= - -!----------------------------- -! scalars: -!----------------------------- - integer,intent(in):: its,ite - integer,intent(in):: itimestep - - real(kind=kind_phys),intent(in):: svp1,svp2,svp3,svpt0,ep1,ep2 - real(kind=kind_phys),intent(in):: karman,cp,g,rovcp,r,xlv - - real(kind=kind_phys),parameter:: prt=1. !prandlt number - real(kind=kind_phys),parameter:: xka=2.4e-5 !molecular diffusivity - - character(len=*),intent(out):: errmsg - integer,intent(out):: errflg - -!----------------------------- -! namelist options -!----------------------------- - logical,intent(in):: spp_pbl,restart_or_cycle - - integer,intent(in):: isfflx - integer,intent(in),optional:: isftcflx,iz0tlnd - -!----------------------------- -! 1d arrays -!----------------------------- - real(kind=kind_phys),intent(in),dimension(its:ite):: mavail, & - pblh, & - xland, & - tsk, & - psfcpa, & - qcg, & - snowh, & - dx - - real(kind=kind_phys),intent(in),dimension(its:ite):: u1d, & - v1d, & - u1d2, & - v1d2, & - qv1d, & - p1d, & - t1d, & - dz8w1d, & - dz2w1d, & - rho1d - real(kind=kind_phys),intent(in),dimension(its:ite):: & - rstoch1d - - - real(kind=kind_phys),intent(inout),dimension(its:ite):: & - regime, & - hfx, & - qfx, & - lh, & - mol, & - rmol, & - qgh, & - qsfc, & - znt, & - zol, & - ust, & - cpm, & - chs2, & - cqs2, & - cqs, & - chs, & - ch, & - flhc, & - flqc, & - gz1oz0, & - wspd, & - br, & - psim, & - psih - -!----------------------------- -! diagnostic outputs: -!----------------------------- - real(kind=kind_phys),intent(out),dimension(its:ite):: & - u10, & - v10, & - th2, & - t2, & - q2 - - real(kind=kind_phys),intent(out),dimension(its:ite):: & - wstar, & - qstar - - real(kind=kind_phys),intent(out),dimension(its:ite),optional:: & - ck, & - cka, & - cd, & - cda, & - ustm - -!----------------------------- -! local variables -!----------------------------- - integer:: n,i,k,l,yesno - - real(kind=kind_phys):: ep3 - real(kind=kind_phys):: pl,thcon,tvcon,e1 - real(kind=kind_phys):: dthvdz,dthvm,vconv,zol2,zol10,zolza,zolz0,zolzt - real(kind=kind_phys):: dtg,psix,dtthx,dthdz,psix10,psit,psit2, & - psiq,psiq2,psiq10,dzdt,psix_min - real(kind=kind_phys):: fluxc,vsgd - real(kind=kind_phys):: restar,visc,dqg,oldust,oldtst - - real(kind=kind_phys),dimension(its:ite) :: & - za, & !height of lowest 1/2 sigma level(m) - za2, & !height of 2nd lowest 1/2 sigma level(m) - thv1d, & !theta-v at lowest 1/2 sigma (K) - th1d, & !theta at lowest 1/2 sigma (K) - tc1d, & !t at lowest 1/2 sigma (Celsius) - tv1d, & !tv at lowest 1/2 sigma (K) - qvsh, & !qv at lowest 1/2 sigma (spec humidity) - psih2, & !m-o stability functions at z=2 m - psim10, & !m-o stability functions at z=10 m - psih10, & !m-o stability functions at z=10 m - wspdi, & - z_q, & !moisture roughness length - z_t, & !thermalroughness length - ZNTstoch, & - govrth, & !g/theta - thgb, & !theta at ground - thvgb, & !theta-v at ground - psfc, & !press at surface (Pa/1000) - qsfcmr, & !qv at surface (mixing ratio, kg/kg) - gz2oz0, & !log((2.0+znt(i))/znt(i)) - gz10oz0, & !log((10.+znt(i))/znt(i)) - gz2ozt, & !log((2.0+z_t(i))/z_t(i)) - gz10ozt, & !log((10.+z_t(i))/z_t(i)) - gz1ozt, & !log((za(i)+z_t(i))/z_t(i)) - zratio !z0/zt - -!----------------------------------------------------------------------------------------------------------------- - - ep3 = 1.-ep2 - if (itimestep == 1) then - if (debug_code) then - print*,"======= printing of constants:" - print*,"cp=", cp," g=", g - print*,"Rd=", R," ep1=", ep1 - print*,"xlv=", XLV," karman=", karman - print*,"ep2=", ep2," ep3=", ep3 - endif - endif - - do i=its,ite - !convert ground & lowest layer temperature to potential temperature: - !psfc cmb - psfc(i)=psfcpa(i)/1000. - thgb(i)=tsk(i)*(100./psfc(i))**rovcp !(K) - !PL cmb - pl=p1d(i)/1000. - thcon=(100./pl)**rovcp - th1d(i)=t1d(i)*thcon !(Theta, K) - tc1d(i)=t1d(i)-273.15 !(T, Celsius) - - !convert to virtual temperature - qvsh(i)=qv1d(i)/(1.+qv1d(i)) !convert to spec hum (kg/kg) - tvcon=(1.+ep1*qvsh(i)) - thv1d(i)=th1d(i)*tvcon !(K) - tv1d(i)=t1d(i)*tvcon !(K) - - !rho1d(i)=psfcpa(i)/(r*tv1d(i)) !now using value calculated in sfc driver - za(i)=0.5*dz8w1d(i) !height of first half-sigma level - za2(i)=dz8w1d(i) + 0.5*dz2w1d(i) !height of 2nd half-sigma level - govrth(i)=g/th1d(i) - enddo - - do i=its,ite - if (tsk(i) .lt. 273.15) then - !saturation vapor pressure wrt ice (svp1=.6112; 10*mb) - e1=svp1*exp(4648*(1./273.15 - 1./tsk(i)) - & - & 11.64*log(273.15/tsk(i)) + 0.02265*(273.15 - tsk(i))) - else - !saturation vapor pressure wrt water (Bolton 1980) - e1=svp1*exp(svp2*(tsk(i)-svpt0)/(tsk(i)-svp3)) - endif - !for land points, qsfc can come from lsm, only recompute over water - if (xland(i).gt.1.5 .or. qsfc(i).le.0.0) then !water - qsfc(i)=ep2*e1/(psfc(i)-ep3*e1) !specific humidity - qsfcmr(i)=ep2*e1/(psfc(i)-e1) !mixing ratio - else !land - qsfcmr(i)=qsfc(i)/(1.-qsfc(i)) - endif - - !qgh changed to use lowest-level air temp consistent with myjsfc change - !q2sat = qgh in LSM - if (tsk(i) .lt. 273.15) then - !saturation vapor pressure wrt ice - e1=svp1*exp(4648*(1./273.15 - 1./t1d(i)) - & - & 11.64*log(273.15/t1d(i)) + 0.02265*(273.15 - t1d(i))) - else - !saturation vapor pressure wrt water (Bolton 1980) - e1=svp1*exp(svp2*(t1d(i)-svpt0)/(t1d(i)-svp3)) - endif - pl=p1d(i)/1000. - !qgh(i)=ep2*e1/(pl-ep_3*e1) !specific humidity - qgh(i)=ep2*e1/(pl-e1) !mixing ratio - cpm(i)=cp*(1.+0.84*qv1d(i)) - enddo - - do i=its,ite - wspd(i)=sqrt(u1d(i)*u1d(i)+v1d(i)*v1d(i)) - - !tgs:thvgb(i)=thgb(i)*(1.+ep1*qsfc(i)*mavail(i)) - thvgb(i)=thgb(i)*(1.+ep1*qsfc(i)) - - dthdz=(th1d(i)-thgb(i)) - dthvdz=(thv1d(i)-thvgb(i)) - - !-------------------------------------------------------- - ! Calculate the convective velocity scale (WSTAR) and - ! subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) - ! and Mahrt and Sun (1995, MWR) low-res correction, respectively - ! (vsgd; for coeff 0.32, 13 km ~ 0.37 m/s; for 3 km == 0 m/s) - !------------------------------------------------------- - !Use Beljaars over land and water - fluxc = max(hfx(i)/rho1d(i)/cp & - & + ep1*thvgb(i)*qfx(i)/rho1d(i),0.) - !wstar(i) = vconvc*(g/tsk(i)*pblh(i)*fluxc)**.33 - if (xland(i).gt.1.5 .or. qsfc(i).le.0.0) then !water - wstar(i) = vconvc*(g/tsk(i)*pblh(i)*fluxc)**.33 - vsgd = 0.20 * (max(dx(i)/5000.-1.,0.))**.33 - else !land - !increase height scale, assuming that the non-local transoport - !from the mass-flux (plume) mixing exceedsd the pblh. - wstar(i) = vconvc*(g/tsk(i)*min(1.5*pblh(i),4000.)*fluxc)**.33 - vsgd = 0.20 * (max(dx(i)/5000.-1.,0.))**.33 - endif - wspd(i)=sqrt(wspd(i)*wspd(i)+wstar(i)*wstar(i)+vsgd*vsgd) - wspd(i)=max(wspd(i),wmin) - - !-------------------------------------------------------- - ! calculate the bulk richardson number of surface layer, - ! according to Akb(1976), Eq(12). - !-------------------------------------------------------- - br(i)=govrth(i)*za(i)*dthvdz/(wspd(i)*wspd(i)) - if (.not. restart_or_cycle .and. itimestep == 1) then - !set limits according to Li et al. (2010) boundary-layer meteorol (p.158) - br(i)=max(br(i),-2.0) - br(i)=min(br(i),2.0) - else - br(i)=max(br(i),-4.0) - br(i)=min(br(i), 4.0) - endif - - ! if previously unstable, do not let into regimes 1 and 2 (stable) - ! if (itimestep .gt. 1) then - ! if(mol(i).lt.0.)br(i)=min(br(i),0.0) - !endif - - enddo - - 1006 format(a,f7.3,a,f9.4,a,f9.5,a,f9.4) - 1007 format(a,f2.0,a,f6.2,a,f7.3,a,f7.2) - - if (debug_code) then - do i=its,ite - print*,"===important input to mynnsfclay: i=",i," ITIME=",ITIMESTEP - print*,"xland=",xland(i)," pblh=",pblh(i)," tsk=",tsk(i), & - " th1=", th1d(i)," qsfc=", qsfc(i)," znt=", znt(i), & - " ust=", ust(i)," snowh=", snowh(i)," psfcpa=",PSFCPA(i), & - " fluxc=",fluxc," qfx=",qfx(i)," hfx=",hfx(i), & - " u1=",u1d(i)," v1=",v1d(i)," qv1=",qvsh(i), & - " rho=",rho1d(i)," w*=",wstar(i)," dz1=",za(i), & - " br=",br(i)," THVGB=",THVGB(i)," DTHVDZ=",DTHVDZ, & - " psim_stab=",psim_stab(1)," psim_unstab=",psim_stab(1), & - " psih_stab=",psih_stab(1)," psih_unstab=",psih_unstab(1) - enddo - endif - -!-------------------------------------------------------------------- -!-------------------------------------------------------------------- -!--- begin i-loop -!-------------------------------------------------------------------- -!-------------------------------------------------------------------- - - do i=its,ite - - !compute kinematic viscosity (m2/s) Andreas (1989) CRREL Rep. 89-11 - !valid between -173 and 277 degrees C. - visc=1.326e-5*(1. + 6.542e-3*tc1d(i) + 8.301e-6*tc1d(i)*tc1d(i) & - - 4.84e-9*tc1d(i)*tc1d(i)*tc1d(i)) - - if ((xland(i)-1.5).ge.0) then - !-------------------------------------- - ! water - !-------------------------------------- - ! calculate z0 (znt) - !-------------------------------------- - if ( present(isftcflx) ) then - if ( isftcflx .eq. 0 ) then - if (coare_opt .eq. 3.0) then - !COARE 3.0 (misleading subroutine name) - call charnock_1955(znt(i),ust(i),wspd(i),visc,za(i)) - else - !COARE 3.5 - call edson_etal_2013(znt(i),ust(i),wspd(i),visc,za(i)) - endif - elseif ( isftcflx .eq. 1 .or. isftcflx .eq. 2 ) then - call davis_etal_2008(znt(i),ust(i)) - elseif ( isftcflx .eq. 3 ) then - call taylor_yelland_2001(znt(i),ust(i),wspd(i)) - elseif ( isftcflx .eq. 4 ) then - if (coare_opt .eq. 3.0) then - !COARE 3.0 (MISLEADING SUBROUTINE NAME) - call charnock_1955(znt(i),ust(i),wspd(i),visc,za(i)) - else - !COARE 3.5 - call edson_etal_2013(znt(i),ust(i),wspd(i),visc,za(i)) - endif - endif - else - !default to COARE 3.0/3.5 - if (coare_opt .eq. 3.0) then - !COARE 3.0 - call charnock_1955(znt(i),ust(i),wspd(i),visc,za(i)) - else - !COARE 3.5 - call edson_etal_2013(znt(i),ust(i),wspd(i),visc,za(i)) - endif - endif - - !add stochastic perturbaction of ZNT - if (spp_pbl) then - zntstoch(i) = max(znt(i) + znt(i)*1.0*rstoch1d(i), 1e-6) - else - zntstoch(i) = znt(i) - endif - - !compute roughness reynolds number (restar) using new znt - ! AHW: Garrattt formula: Calculate roughness Reynolds number - ! Kinematic viscosity of air (linear approx to - ! temp dependence at sea level) - restar=max(ust(i)*zntstoch(i)/visc, 0.1) - - !-------------------------------------- - !calculate z_t and z_q - !-------------------------------------- - if ( present(isftcflx) ) then - if ( isftcflx .eq. 0 ) then - if (coare_opt .eq. 3.0) then - call fairall_etal_2003(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - else - !presumably, this will be published soon, but hasn't yet - call fairall_etal_2014(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - endif - elseif ( isftcflx .eq. 1 ) then - if (coare_opt .eq. 3.0) then - call fairall_etal_2003(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - else - call fairall_etal_2014(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - endif - elseif ( isftcflx .eq. 2 ) then - call garratt_1992(z_t(i),z_q(i),zntstoch(i),restar,xland(i)) - elseif ( isftcflx .eq. 3 ) then - if (coare_opt .eq. 3.0) then - call fairall_etal_2003(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - else - call fairall_etal_2014(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - endif - endif - else - !default to COARE 3.0/3.5 - if (coare_opt .eq. 3.0) then - call fairall_etal_2003(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - else - call fairall_etal_2014(z_t(i),z_q(i),restar,ust(i),visc,rstoch1d(i),spp_pbl) - endif - endif - - else - - !add stochastic perturbaction of znt - if (spp_pbl) then - zntstoch(i) = max(znt(i) + znt(i)*1.0*rstoch1d(i), 1e-6) - else - zntstoch(i) = znt(i) - endif - - !-------------------------------------- - ! land - !-------------------------------------- - !compute roughness reynolds number (restar) using default znt - restar=max(ust(i)*zntstoch(i)/visc, 0.1) - - !-------------------------------------- - ! get z_t and z_q - !-------------------------------------- - !check for snow/ice points over land - if ( snowh(i) .ge. 0.1) then - call andreas_2002(zntstoch(i),visc,ust(i),z_t(i),z_q(i)) - else - if ( present(iz0tlnd) ) then - if ( iz0tlnd .le. 1 ) then - call zilitinkevich_1995(zntstoch(i),z_t(i),z_q(i),restar,& - ust(i),karman,xland(i),iz0tlnd,spp_pbl,rstoch1d(i)) - elseif ( iz0tlnd .eq. 2 ) then - call yang_2008(zntstoch(i),z_t(i),z_q(i),ust(i),mol(i),& - qstar(i),restar,visc,xland(i)) - elseif ( iz0tlnd .eq. 3 ) then - !original mynn in wrf-arw used this form: - call garratt_1992(z_t(i),z_q(i),zntstoch(i),restar,xland(i)) - endif - else - !default to zilitinkevich - call zilitinkevich_1995(zntstoch(i),z_t(i),z_q(i),restar,& - ust(i),karman,xland(i),0,spp_pbl,rstoch1d(i)) - endif - endif - - endif - zratio(i)=zntstoch(i)/z_t(i) !needed for Li et al. - - gz1oz0(i)= log((za(i)+zntstoch(i))/zntstoch(i)) - gz1ozt(i)= log((za(i)+zntstoch(i))/z_t(i)) - gz2oz0(i)= log((2.0+zntstoch(i))/zntstoch(i)) - gz2ozt(i)= log((2.0+zntstoch(i))/z_t(i)) - gz10oz0(i)=log((10.+zntstoch(i))/zntstoch(i)) - gz10ozt(i)=log((10.+zntstoch(i))/z_t(i)) - - !-------------------------------------------------------------------- - !--- DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATE STABILITY CLASS: - ! - ! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.). - ! - ! CRITERIA FOR THE CLASSES ARE AS FOLLOWS: - ! - ! 1. BR .GE. 0.2; - ! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1), - ! - ! 2. BR .LT. 0.2 .AND. BR .GT. 0.0; - ! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS - ! (REGIME=2), - ! - ! 3. BR .EQ. 0.0 - ! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3), - ! - ! 4. BR .LT. 0.0 - ! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4). - ! - !-------------------------------------------------------------------- - if (br(i) .gt. 0.0) then - if (br(i) .gt. 0.2) then - !---class 1; stable (nighttime) conditions: - regime(i)=1. - else - !---class 2; damped mechanical turbulence: - regime(i)=2. - endif - - if (.not. restart_or_cycle .and. itimestep .le. 1) then - !compute z/l first guess: - call li_etal_2010(zol(i),br(i),za(i)/zntstoch(i),zratio(i)) - else - zol(i)=za(i)*karman*g*mol(i)/(th1d(i)*max(ust(i)*ust(i),0.0001)) - zol(i)=max(zol(i),0.0) - zol(i)=min(zol(i),20.) - endif - - !Use Pedros iterative function to find z/L - !zol(i)=zolri(br(i),za(i),zntstoch(i),z_t(i),zol(i),psi_opt) - !Use brute-force method - zol(i)=zolrib(br(i),za(i),zntstoch(i),z_t(i),gz1oz0(i),gz1ozt(i),zol(i),psi_opt) - zol(i)=max(zol(i),0.0) - zol(i)=min(zol(i),20.) - - zolzt = zol(i)*z_t(i)/za(i) ! zt/l - zolz0 = zol(i)*zntstoch(i)/za(i) ! z0/l - zolza = zol(i)*(za(i)+zntstoch(i))/za(i) ! (z+z0/l - zol10 = zol(i)*(10.+zntstoch(i))/za(i) ! (10+z0)/l - zol2 = zol(i)*(2.+zntstoch(i))/za(i) ! (2+z0)/l - - !compute psim and psih - if ((xland(i)-1.5).ge.0) then - ! water - !call psi_suselj_sood_2010(psim(i),psih(i),zol(i)) - !call psi_beljaars_holtslag_1991(psim(i),psih(i),zol(i)) - !call psi_businger_1971(psim(i),psih(i),zol(i)) - !call psi_dyerhicks(psim(i),psih(i),zol(i),z_t(i),zntstoch(i),za(i)) - !call psi_cb2005(psim(i),psih(i),zolza,zolz0) - ! or use tables - psim(i)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt) - psih(i)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt) - psim10(i)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt) - psih10(i)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt) - psih2(i)=psih_stable(zol2,psi_opt)-psih_stable(zolz0,psi_opt) - else - ! land - !call psi_beljaars_holtslag_1991(psim(i),psih(i),zol(i)) - !call psi_businger_1971(psim(i),psih(i),zol(i)) - !call psi_zilitinkevich_esau_2007(psim(i),psih(i),zol(i)) - !call psi_dyerhicks(psim(i),psih(i),zol(i),z_t(i),zntstoch(i),za(i)) - !call psi_cb2005(psim(i),psih(i),zolza,zolz0) - ! or use tables - psim(i)=psim_stable(zolza,psi_opt)-psim_stable(zolz0,psi_opt) - psih(i)=psih_stable(zolza,psi_opt)-psih_stable(zolzt,psi_opt) - psim10(i)=psim_stable(zol10,psi_opt)-psim_stable(zolz0,psi_opt) - psih10(i)=psih_stable(zol10,psi_opt)-psih_stable(zolz0,psi_opt) - psih2(i)=psih_stable(zol2,psi_opt)-psih_stable(zolz0,psi_opt) - endif - - !psim10(i)=10./za(i)*psim(i) - !psih10(i)=10./za(i)*psih(i) - !psim2(i)=2./za(i)*psim(i) - !psih2(i)=2./za(i)*psih(i) - - ! 1.0 over monin-obukhov length - rmol(i)= zol(i)/za(i) - - elseif(br(i) .eq. 0.) then - !========================================================= - !-----class 3; forced convection/neutral: - !========================================================= - regime(i)=3. - - psim(i)=0.0 - psih(i)=psim(i) - psim10(i)=0. - psih10(i)=0. - psih2(i)=0. - - zol(i)=0. - !if (ust(i) .lt. 0.01) then - ! zol(i)=br(i)*gz1oz0(i) - !else - ! zol(i)=karman*govrth(i)*za(i)*mol(i)/(max(ust(i)*ust(i),0.001)) - !endif - rmol(i) = zol(i)/za(i) - - elseif(br(i) .lt. 0.)then - !========================================================== - !-----class 4; free convection: - !========================================================== - regime(i)=4. - - if (.not. restart_or_cycle .and. itimestep .le. 1) then - !compute z/l first guess: - call li_etal_2010(zol(i),br(i),za(i)/zntstoch(i),zratio(i)) - else - zol(i)=za(i)*karman*g*mol(i)/(th1d(i)*max(ust(i)*ust(i),0.001)) - zol(i)=max(zol(i),-20.0) - zol(i)=min(zol(i),0.0) - endif - - !Use Pedros iterative function to find z/L - !zol(I)=zolri(br(I),ZA(I),ZNTstoch(I),z_t(I),ZOL(I),psi_opt) - !Use alternative method - zol(i)=zolrib(br(i),za(i),zntstoch(i),z_t(i),gz1oz0(i),gz1ozt(i),zol(i),psi_opt) - zol(i)=max(zol(i),-20.0) - zol(i)=min(zol(i),0.0) - - zolzt = zol(i)*z_t(i)/za(i) ! zt/l - zolz0 = zol(i)*zntstoch(i)/za(i) ! z0/l - zolza = zol(i)*(za(i)+zntstoch(i))/za(i) ! (z+z0/l - zol10 = zol(i)*(10.+zntstoch(i))/za(i) ! (10+z0)/l - zol2 = zol(i)*(2.+zntstoch(i))/za(i) ! (2+z0)/l - - !compute psim and psih - if ((xland(i)-1.5).ge.0) then - ! water - !call psi_suselj_sood_2010(psim(i),psih(i),zol(i)) - !call psi_hogstrom_1996(psim(i),psih(i),zol(i), z_t(i), zntstoch(i), za(i)) - !call psi_businger_1971(psim(i),psih(i),zol(i)) - !call psi_dyerhicks(psim(i),psih(i),zol(i),z_t(i),zntstoch(i),za(i)) - ! use tables - psim(i)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt) - psih(i)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt) - psim10(i)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt) - psih10(i)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt) - psih2(i)=psih_unstable(zol2,psi_opt)-psih_unstable(zolz0,psi_opt) - else - ! land - !call psi_hogstrom_1996(psim(i),psih(i),zol(i), z_t(i), zntstoch(i), za(i)) - !call psi_businger_1971(psim(i),psih(i),zol(i)) - !call psi_dyerhicks(psim(i),psih(i),zol(i),z_t(i),zntstoch(i),za(i)) - ! use tables - psim(i)=psim_unstable(zolza,psi_opt)-psim_unstable(zolz0,psi_opt) - psih(i)=psih_unstable(zolza,psi_opt)-psih_unstable(zolzt,psi_opt) - psim10(i)=psim_unstable(zol10,psi_opt)-psim_unstable(zolz0,psi_opt) - psih10(i)=psih_unstable(zol10,psi_opt)-psih_unstable(zolz0,psi_opt) - psih2(i)=psih_unstable(zol2,psi_opt)-psih_unstable(zolz0,psi_opt) - endif - - !psim10(i)=10./za(i)*psim(i) - !psih2(i)=2./za(i)*psih(i) - - !---limit psih and psim in the case of thin layers and - !---high roughness. this prevents denominator in fluxes - !---from getting too small - psih(i)=min(psih(i),0.9*gz1ozt(i)) - psim(i)=min(psim(i),0.9*gz1oz0(i)) - psih2(i)=min(psih2(i),0.9*gz2ozt(i)) - psim10(i)=min(psim10(i),0.9*gz10oz0(i)) - psih10(i)=min(psih10(i),0.9*gz10ozt(i)) - - rmol(i) = zol(i)/za(i) - - endif - - !------------------------------------------------------------ - !-----compute the frictional velocity: - !------------------------------------------------------------ - ! over water, apply maximum cd (=3.1e-3), which corresponds to - ! psix10 = 7.1: - if ((xland(i)-1.5).lt.0.) then !land - psix=gz1oz0(i)-psim(i) - psix10=gz10oz0(i)-psim10(i) - else !water - psix10=max(gz10oz0(i)-psim10(i),7.1) - psix_min=7.1*log(za(i)/zntstoch(i))/log(10./zntstoch(i)) - psix=max(gz1oz0(i)-psim(i), psix_min) - endif - ! to prevent oscillations average with old value - oldust = ust(i) - ust(i)=0.5*ust(i)+0.5*karman*wspd(i)/psix - !non-averaged: ust(i)=karman*wspd(i)/psix - - ! compute u* without vconv for use in hfx calc when isftcflx > 0 - wspdi(i)=max(sqrt(u1d(i)*u1d(i)+v1d(i)*v1d(i)), wmin) - if ( present(ustm) ) then - ustm(i)=0.5*ustm(i)+0.5*karman*wspdi(i)/psix - endif - - if ((xland(i)-1.5).lt.0.) then !land - ust(i)=max(ust(i),0.005) !further relaxing this limit - no need to go lower - !keep ustm = ust over land. - if ( present(ustm) ) ustm(i)=ust(i) - endif - - !------------------------------------------------------------ - !-----compute the thermal and moisture resistance (psiq and psit): - !------------------------------------------------------------ - ! lower limit added to prevent large flhc in soil model - ! activates in unstable conditions with thin layers or high z0 - gz1ozt(i)= log((za(i)+zntstoch(i))/z_t(i)) - gz2ozt(i)= log((2.0+zntstoch(i))/z_t(i)) - - psit =max(gz1ozt(i)-psih(i) ,1.) - psit2=max(gz2ozt(i)-psih2(i),1.) - - psiq=max(log((za(i)+zntstoch(i))/z_q(i))-psih(i) ,1.0) - psiq2=max(log((2.0+zntstoch(i))/z_q(i))-psih2(i) ,1.0) - psiq10=max(log((10.0+zntstoch(i))/z_q(i))-psih10(i) ,1.0) - !---------------------------------------------------- - !compute the temperature scale (or friction temperature, T*) - !---------------------------------------------------- - dtg=thv1d(i)-thvgb(i) - oldtst=mol(i) - mol(i)=karman*dtg/psit/prt - !t_star(i) = -hfx(i)/(ust(i)*cpm(i)*rho1d(i)) - !t_star(i) = mol(i) - !---------------------------------------------------- - !compute the moisture scale (or q*) - dqg=(qvsh(i)-qsfc(i))*1000. !(kg/kg -> g/kg) - qstar(i)=karman*dqg/psiq/prt - - !if () then - ! write(*,1001)"regime:",regime(i)," z/l:",zol(i)," u*:",ust(i)," tstar:",mol(i) - ! write(*,1002)"psim:",psim(i)," psih:",psih(i)," w*:",wstar(i)," dthv:",thv1d(i)-thvgb(i) - ! write(*,1003)"cpm:",cpm(i)," rho1d:",rho1d(i)," l:",zol(i)/za(i)," dth:",th1d(i)-thgb(i) - ! write(*,1004)"z0/zt:",zratio(i)," z0:",zntstoch(i)," zt:",z_t(i)," za:",za(i) - ! write(*,1005)"re:",restar," mavail:",mavail(i)," qsfc(i):",qsfc(i)," qvsh(i):",qvsh(i) - ! print*,"visc=",visc," z0:",zntstoch(i)," t1d(i):",t1d(i) - ! write(*,*)"=============================================" - !endif - - enddo ! end i-loop - - 1000 format(a,f6.1, a,f6.1, a,f5.1, a,f7.1) - 1001 format(a,f2.0, a,f10.4,a,f5.3, a,f11.5) - 1002 format(a,f7.2, a,f7.2, a,f7.2, a,f10.3) - 1003 format(a,f7.2, a,f7.2, a,f10.3,a,f10.3) - 1004 format(a,f11.3,a,f9.7, a,f9.7, a,f6.2, a,f10.3) - 1005 format(a,f9.2,a,f6.4,a,f7.4,a,f7.4) - - !---------------------------------------------------------- - ! compute surface heat and moisture fluxes - !---------------------------------------------------------- - do i=its,ite - - !For computing the diagnostics and fluxes (below), whether the fluxes - !are turned off or on, we need the following: - ! over water, apply maximum cd (=3.1e-3), which corresponds to - ! psix = 7.1 (at 10 m): - if ((xland(i)-1.5).lt.0.) then !land - psix=gz1oz0(i)-psim(i) - psix10=gz10oz0(i)-psim10(i) - else !water - psix10=max(gz10oz0(i)-psim10(i),7.1) - psix_min=7.1*log(za(i)/zntstoch(i))/log(10./zntstoch(i)) - psix=max(gz1oz0(i)-psim(i), psix_min) - endif - - psit =max(gz1ozt(i)-psih(i), 1.0) - psit2=max(gz2ozt(i)-psih2(i),1.0) - - psiq=max(log((za(i)+z_q(i))/z_q(i))-psih(i) ,1.0) - psiq2=max(log((2.0+z_q(i))/z_q(i))-psih2(i) ,1.0) - psiq10=max(log((10.0+z_q(i))/z_q(i))-psih10(i) ,1.0) - - if (isfflx .lt. 1) then - - qfx(i) = 0. - hfx(i) = 0. - flhc(i) = 0. - flqc(i) = 0. - lh(i) = 0. - chs(i) = 0. - ch(i) = 0. - cqs(i) = 0. - chs2(i) = 0. - cqs2(i) = 0. - if(present(ck) .and. present(cd) .and. & - &present(cka) .and. present(cda)) then - ck(i) = 0. - cd(i) = 0. - cka(i)= 0. - cda(i)= 0. - endif - else - - !------------------------------------------ - ! calculate the exchange coefficients for heat (flhc) - ! and moisture (flqc) - !------------------------------------------ - flqc(i)=rho1d(i)*mavail(i)*ust(i)*karman/psiq - flhc(i)=rho1d(i)*cpm(i)*ust(i)*karman/psit - - !---------------------------------- - ! compute surface moisture flux: - !---------------------------------- - qfx(i)=flqc(i)*(qsfcmr(i)-qv1d(i)) - !joe: qfx(i)=max(qfx(i),0.) !originally did not allow neg qfx - qfx(i)=max(qfx(i),-0.02) !allows small neg qfx, like myj - lh(i)=xlv*qfx(i) - - !---------------------------------- - ! compute surface heat flux: - !---------------------------------- - if(xland(i)-1.5.gt.0.)then !water - hfx(i)=flhc(i)*(thgb(i)-th1d(i)) - if ( present(isftcflx) ) then - if ( isftcflx.ne.0 ) then - ! ahw: add dissipative heating term - hfx(i)=hfx(i)+rho1d(i)*ustm(i)*ustm(i)*wspdi(i) - endif - endif - elseif(xland(i)-1.5.lt.0.)then !land - hfx(i)=flhc(i)*(thgb(i)-th1d(i)) - hfx(i)=max(hfx(i),-250.) - endif - - !chs(i)=ust(i)*karman/(alog(karman*ust(i)*za(i) & - ! /xka+za(i)/zl)-psih(i)) - - chs(i)=ust(i)*karman/psit - - ! the exchange coefficient for cloud water is assumed to be the - ! same as that for heat. ch is multiplied by wspd. - - !ch(i)=chs(i) - ch(i)=flhc(i)/( cpm(i)*rho1d(i) ) - - !these are used for 2-m diagnostics only - cqs(i) =ust(i)*karman/psiq - cqs2(i)=ust(i)*karman/psiq2 - chs2(i)=ust(i)*karman/psit2 - - if(present(ck) .and. present(cd) .and. & - &present(cka) .and. present(cda)) then - ck(i)=(karman/psix10)*(karman/psiq10) - cd(i)=(karman/psix10)*(karman/psix10) - cka(i)=(karman/psix)*(karman/psiq) - cda(i)=(karman/psix)*(karman/psix) - endif - - endif !end isfflx option - - !----------------------------------------------------- - !compute diagnostics - !----------------------------------------------------- - !compute 10 m wnds - !----------------------------------------------------- - ! If the lowest model level is close to 10-m, use it - ! instead of the flux-based diagnostic formula. - if (za(i) .le. 7.0) then - ! high vertical resolution - if(za2(i) .gt. 7.0 .and. za2(i) .lt. 13.0) then - !use 2nd model level - u10(i)=u1d2(i) - v10(i)=v1d2(i) - else - u10(i)=u1d(i)*log(10./zntstoch(i))/log(za(i)/zntstoch(i)) - v10(i)=v1d(i)*log(10./zntstoch(i))/log(za(i)/zntstoch(i)) - endif - elseif(za(i) .gt. 7.0 .and. za(i) .lt. 13.0) then - !moderate vertical resolution - !u10(i)=u1d(i)*psix10/psix - !v10(i)=v1d(i)*psix10/psix - !use neutral-log: - u10(i)=u1d(i)*log(10./zntstoch(i))/log(za(i)/zntstoch(i)) - v10(i)=v1d(i)*log(10./zntstoch(i))/log(za(i)/zntstoch(i)) - else - ! very coarse vertical resolution - u10(i)=u1d(i)*psix10/psix - v10(i)=v1d(i)*psix10/psix - endif - - !----------------------------------------------------- - !compute 2m t, th, and q - !these will be overwritten for land points in the lsm - !----------------------------------------------------- - dtg=th1d(i)-thgb(i) - th2(i)=thgb(i)+dtg*psit2/psit - !*** be certain that the 2-m theta is bracketed by - !*** the values at the surface and lowest model level. - if ((th1d(i)>thgb(i) .and. (th2(i)th1d(i))) .or. & - (th1d(i)thgb(i) .or. th2(i) 1200. .or. hfx(i) < -700.)then - print*,"suspicious values in mynn sfclayer",& - i, "hfx: ",hfx(i) - yesno = 1 - endif - if (lh(i) > 1200. .or. lh(i) < -700.)then - print*,"suspicious values in mynn sfclayer",& - i, "lh: ",lh(i) - yesno = 1 - endif - if (ust(i) < 0.0 .or. ust(i) > 4.0 )then - print*,"suspicious values in mynn sfclayer",& - i, "ust: ",ust(i) - yesno = 1 - endif - if (wstar(i)<0.0 .or. wstar(i) > 6.0)then - print*,"suspicious values in mynn sfclayer",& - i, "wstar: ",wstar(i) - yesno = 1 - endif - if (rho1d(i)<0.0 .or. rho1d(i) > 1.6 )then - print*,"suspicious values in mynn sfclayer",& - i, "rho: ",rho1d(i) - yesno = 1 - endif - if (qsfc(i)*1000. <0.0 .or. qsfc(i)*1000. >40.)then - print*,"suspicious values in mynn sfclayer",& - i, "qsfc: ",qsfc(i) - yesno = 1 - endif - if (pblh(i)<0. .or. pblh(i)>6000.)then - print*,"suspicious values in mynn sfclayer",& - i, "pblh: ",pblh(i) - yesno = 1 - endif - - if (yesno == 1) then - print*," other info:" - write(*,1001)"regime:",regime(i)," z/l:",zol(i)," u*:",ust(i),& - " tstar:",mol(i) - write(*,1002)"psim:",psim(i)," psih:",psih(i)," w*:",wstar(i),& - " dthv:",thv1d(i)-thvgb(i) - write(*,1003)"cpm:",cpm(i)," rho1d:",rho1d(i)," l:",& - zol(i)/za(i)," dth:",th1d(i)-thgb(i) - write(*,*)" z0:",zntstoch(i)," zt:",z_t(i)," za:",za(i) - write(*,1005)"re:",restar," mavail:",mavail(i)," qsfc(i):",& - qsfc(i)," qvsh(i):",qvsh(i) - print*,"psix=",psix," z0:",zntstoch(i)," t1d(i):",t1d(i) - write(*,*)"=============================================" - endif - endif - - enddo !end i-loop - - errmsg = ' ' - errflg = 0 - - end subroutine mynnsfclay - -!================================================================================================================= - subroutine zilitinkevich_1995(z_0,zt,zq,restar,ustar,karman,landsea,iz0tlnd2,spp_pbl,rstoch) -!this subroutine returns the thermal and moisture roughness lengths -!from Zilitinkevich (1995) and Zilitinkevich et al. (2001) over -!land and water, respectively. -! -!MODS: -!20120705 : added IZ0TLND option. Note: This option was designed -! to work with the Noah LSM and may be specific for that -! LSM only. Tests with RUC LSM showed no improvements. - implicit none -!================================================================================================================= - -!--- input arguments: - logical,intent(in):: spp_pbl - integer,optional,intent(in):: iz0tlnd2 - - real(kind=kind_phys),intent(in):: rstoch - real(kind=kind_phys),intent(in):: z_0,restar,ustar,karman,landsea - -!--- output arguments: - real(kind=kind_phys),intent(out):: zt,zq - -!--- local variables: - real(kind=kind_phys):: czil !=0.100 in Chen et al. (1997) - !=0.075 in Zilitinkevich (1995) - !=0.500 in Lemone et al. (2008) - -!----------------------------------------------------------------------------------------------------------------- - - if (landsea-1.5 .gt. 0) then !water - -!this is based on Zilitinkevich, Grachev, and Fairall (2001): -!their equations 15 and 16). - if (restar .lt. 0.1) then - zt = z_0*exp(karman*2.0) - zt = min( zt, 6.0e-5) - zt = max( zt, 2.0e-9) - zq = z_0*exp(karman*3.0) - zq = min( zq, 6.0e-5) - zq = max( zq, 2.0e-9) - else - zt = z_0*exp(-karman*(4.0*sqrt(restar)-3.2)) - zt = min( zt, 6.0e-5) - zt = max( zt, 2.0e-9) - zq = z_0*exp(-karman*(4.0*sqrt(restar)-4.2)) - zq = min( zt, 6.0e-5) - zq = max( zt, 2.0e-9) - endif - - else !land - -!option to modify czil according to Chen & Zhang (2009): - if ( iz0tlnd2 .eq. 1 ) then - czil = 10.0 ** ( -0.40 * ( z_0 / 0.07 ) ) - else - czil = 0.095 !0.075 !0.10 - end if - - zt = z_0*exp(-karman*czil*sqrt(restar)) - zt = min( zt, 0.75*z_0) - - zq = z_0*exp(-karman*czil*sqrt(restar)) - zq = min( zq, 0.75*z_0) - -!stochastically perturb thermal and moisture roughness length. -!currently set to half the amplitude: - if (spp_pbl) then - zt = zt + zt * 0.5 * rstoch - zt = max(zt, 0.0001) - zq = zt - endif - - endif - - end subroutine zilitinkevich_1995 - -!================================================================================================================= - subroutine davis_etal_2008(Z_0,ustar) -!a.k.a. : Donelan et al. (2004) -!this formulation for roughness length was designed to match -!the labratory experiments of Donelan et al. (2004). -!this is an update version from Davis et al. 2008, which -!corrects a small-bias in Z_0 (AHW real-time 2012). - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: ustar - -!--- output arguments: - real(kind=kind_phys),intent(out):: z_0 - -!--- local variables: - real(kind=kind_phys):: zw, zn1, zn2 - real(kind=kind_phys),parameter:: g=9.81,ozo=1.59e-5 - -!----------------------------------------------------------------------------------------------------------------- - -!old form: z_0 = 10.*exp(-10./(ustar**(1./3.))) -!new form: - - zw = min((ustar/1.06)**(0.3),1.0) - zn1 = 0.011*ustar*ustar/g + ozo - zn2 = 10.*exp(-9.5*ustar**(-.3333)) + & - 0.11*1.5e-5/amax1(ustar,0.01) - z_0 = (1.0-zw) * zn1 + zw * zn2 - - z_0 = max( z_0, 1.27e-7) !these max/mins were suggested by - z_0 = min( z_0, 2.85e-3) !Davis et al. (2008) - - end subroutine davis_etal_2008 - -!================================================================================================================= - subroutine taylor_yelland_2001(z_0,ustar,wsp10) -!this formulation for roughness length was designed account for -!wave steepness. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: ustar,wsp10 - -!--- output arguments: - real(kind=kind_phys),intent(out):: z_0 - -!--- local variables: - real(kind=kind_phys),parameter:: g=9.81, pi=3.14159265 - real(kind=kind_phys):: hs, tp, lp - -!----------------------------------------------------------------------------------------------------------------- - -!hs is the significant wave height - hs = 0.0248*(wsp10**2.) -!Tp dominant wave period - tp = 0.729*max(wsp10,0.1) -!lp is the wavelength of the dominant wave - lp = g*tp**2/(2*pi) - - z_0 = 1200.*hs*(hs/lp)**4.5 - z_0 = max( z_0, 1.27e-7) !these max/mins were suggested by - z_0 = min( z_0, 2.85e-3) !Davis et al. (2008) - - end subroutine taylor_yelland_2001 - -!================================================================================================================= - subroutine charnock_1955(Z_0,ustar,wsp10,visc,zu) -!This version of Charnock's relation employs a varying -!Charnock parameter, similar to COARE3.0 [Fairall et al. (2003)]. -!The Charnock parameter CZC is varied from .011 to .018 -!between 10-m wsp = 10 and 18. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: ustar,visc,wsp10,zu - -!--- output arguments: - real(kind=kind_phys),intent(out):: z_0 - -!--- local variables: - real(kind=kind_phys),parameter:: G=9.81, CZO2=0.011 - real(kind=kind_phys):: czc !variable charnock "constant" - real(kind=kind_phys):: wsp10m ! logarithmically calculated 10 m - -!----------------------------------------------------------------------------------------------------------------- - - wsp10m = wsp10*log(10./1e-4)/log(zu/1e-4) - czc = czo2 + 0.007*min(max((wsp10m-10.)/8., 0.), 1.0) - - z_0 = czc*ustar*ustar/g + (0.11*visc/max(ustar,0.05)) - z_0 = max( z_0, 1.27e-7) !these max/mins were suggested by - z_0 = min( z_0, 2.85e-3) !Davis et al. (2008) - - end subroutine charnock_1955 - -!================================================================================================================= - subroutine edson_etal_2013(z_0,ustar,wsp10,visc,zu) -!This version of Charnock's relation employs a varying -!Charnock parameter, taken from COARE 3.5 [Edson et al. (2001, JPO)]. -!The Charnock parameter CZC is varied from about .005 to .028 -!between 10-m wind speeds of 6 and 19 m/s. -!11 Nov 2021: Note that this was finally fixed according to the -! Edson et al (2014) corrigendum, where "m" was corrected. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: ustar,visc,wsp10,zu - -!--- output arguments: - real(kind=kind_phys),intent(out):: z_0 - -!--- local variables: - real(kind=kind_phys),parameter:: g=9.81 - real(kind=kind_phys),parameter:: m=0.0017, b=-0.005 - real(kind=kind_phys):: czc ! variable charnock "constant" - real(kind=kind_phys):: wsp10m ! logarithmically calculated 10 m - -!----------------------------------------------------------------------------------------------------------------- - - wsp10m = wsp10*log(10/1e-4)/log(zu/1e-4) - wsp10m = min(19.,wsp10m) - czc = m*wsp10m + b - czc = max(czc, 0.0) - - z_0 = czc*ustar*ustar/g + (0.11*visc/max(ustar,0.07)) - z_0 = max( z_0, 1.27e-7) !These max/mins were suggested by - z_0 = min( z_0, 2.85e-3) !Davis et al. (2008) - - end subroutine edson_etal_2013 - -!================================================================================================================= - subroutine garratt_1992(zt,zq,z_0,ren,landsea) -!This formulation for the thermal and moisture roughness lengths -!(Zt and Zq) relates them to Z0 via the roughness Reynolds number (Ren). -!This formula comes from Fairall et al. (2003). It is modified from -!the original Garratt-Brutsaert model to better fit the COARE/HEXMAX -!data. The formula for land uses a constant ratio (Z_0/7.4) taken -!from Garratt (1992). - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: ren, z_0,landsea - -!--- output arguments: - real(kind=kind_phys),intent(out):: zt,zq - -!--- local variables: - real(kind=kind_phys):: rq - real(kind=kind_phys),parameter:: e=2.71828183 - -!----------------------------------------------------------------------------------------------------------------- - - if (landsea-1.5 .gt. 0) then !water - - zt = z_0*exp(2.0 - (2.48*(ren**0.25))) - zq = z_0*exp(2.0 - (2.28*(ren**0.25))) - - zq = min( zq, 5.5e-5) - zq = max( zq, 2.0e-9) - zt = min( zt, 5.5e-5) - zt = max( zt, 2.0e-9) !same lower limit as ecmwf - - else !land - - zq = z_0/(e**2.) !taken from Garratt (1980,1992) - zt = zq - - endif - - end subroutine garratt_1992 - -!================================================================================================================= - subroutine fairall_etal_2003(zt,zq,ren,ustar,visc,rstoch,spp_pbl) -!This formulation for thermal and moisture roughness length (Zt and Zq) -!as a function of the roughness Reynolds number (Ren) comes from the -!COARE3.0 formulation, empirically derived from COARE and HEXMAX data -![Fairall et al. (2003)]. Edson et al. (2004; JGR) suspected that this -!relationship overestimated the scalar roughness lengths for low Reynolds -!number flows, so an optional smooth flow relationship, taken from Garratt -!(1992, p. 102), is available for flows with Ren < 2. -! -!This is for use over water only. - implicit none -!================================================================================================================= - -!--- input arguments: - logical,intent(in):: spp_pbl - real(kind=kind_phys),intent(in):: ren,ustar,visc,rstoch - -!--- output arguments: - real(kind=kind_phys),intent(out):: zt,zq - -!----------------------------------------------------------------------------------------------------------------- - - if (ren .le. 2.) then - - zt = (5.5e-5)*(ren**(-0.60)) - zq = zt - !for smooth seas, can use Garratt - !zq = 0.2*visc/max(ustar,0.1) - !zq = 0.3*visc/max(ustar,0.1) - - else - - !for rough seas, use coare - zt = (5.5e-5)*(ren**(-0.60)) - zq = zt - - endif - - if (spp_pbl) then - zt = zt + zt * 0.5 * rstoch - zq = zt - endif - - zt = min(zt,1.0e-4) - zt = max(zt,2.0e-9) - - zq = min(zt,1.0e-4) - zq = max(zt,2.0e-9) - - end subroutine fairall_etal_2003 - -!================================================================================================================= - subroutine fairall_etal_2014(zt,zq,ren,ustar,visc,rstoch,spp_pbl) -!This formulation for thermal and moisture roughness length (Zt and Zq) -!as a function of the roughness Reynolds number (Ren) comes from the -!COARE 3.5/4.0 formulation, empirically derived from COARE and HEXMAX data -![Fairall et al. (2014? coming soon, not yet published as of July 2014)]. -!This is for use over water only. - implicit none -!================================================================================================================= - -!--- input arguments: - logical,intent(in):: spp_pbl - real(kind=kind_phys),intent(in):: ren,ustar,visc,rstoch - -!--- output arguments: - real(kind=kind_phys),intent(out):: Zt,Zq - -!----------------------------------------------------------------------------------------------------------------- - -!zt = (5.5e-5)*(ren**(-0.60)) - zt = min(1.6e-4, 5.8e-5/(ren**0.72)) - zq = zt - - if (spp_pbl) then - zt = max(zt + zt*0.5*rstoch,2.0e-9) - zq = max(zt + zt*0.5*rstoch,2.0e-9) - else - zt = max(zt,2.0e-9) - zq = max(zt,2.0e-9) - endif - - - end subroutine fairall_etal_2014 - -!================================================================================================================= - subroutine yang_2008(z_0,zt,zq,ustar,tstar,qst,ren,visc,landsea) -!This is a modified version of Yang et al (2002 QJRMS, 2008 JAMC) -!and Chen et al (2010, J of Hydromet). Although it was originally -!designed for arid regions with bare soil, it is modified -!here to perform over a broader spectrum of vegetation. -! -!The original formulation relates the thermal roughness length (Zt) -!to u* and T*: -! -! Zt = ht * EXP(-beta*(ustar**0.5)*(ABS(tstar)**0.25)) -! -!where ht = Renc*visc/ustar and the critical Reynolds number -!(Renc) = 70. Beta was originally = 10 (2002 paper) but was revised -!to 7.2 (in 2008 paper). Their form typically varies the -!ratio Z0/Zt by a few orders of magnitude (1-1E4). -! -!This modified form uses beta = 1.5 and a variable Renc (function of Z_0), -!so zt generally varies similarly to the Zilitinkevich form (with Czil ~ 0.1) -!for very small or negative surface heat fluxes but can become close to the -!Zilitinkevich with Czil = 0.2 for very large HFX (large negative T*). -!Also, the exponent (0.25) on tstar was changed to 1.0, since we found -!Zt was reduced too much for low-moderate positive heat fluxes. -! -!This should only be used over land! - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: z_0,ren,ustar,tstar,qst,visc,landsea - -!--- output arguments: - real(kind=kind_phys),intent(out):: zt,zq - -!--- local variables: - real(kind=kind_phys):: ht, &! roughness height at critical Reynolds number - tstar2, &! bounded T*, forced to be non-positive - qstar2, &! bounded q*, forced to be non-positive - z_02, &! bounded Z_0 for variable Renc2 calc - renc2 ! variable Renc, function of Z_0 - - real(kind=kind_phys),parameter:: renc=300., & !old constant Renc - beta=1.5, & !important for diurnal variation - m=170., & !slope for Renc2 function - b=691. !y-intercept for Renc2 function - -!----------------------------------------------------------------------------------------------------------------- - - z_02 = min(z_0,0.5) - z_02 = max(z_02,0.04) - renc2= b + m*log(z_02) - ht = renc2*visc/max(ustar,0.01) - tstar2 = min(tstar, 0.0) - qstar2 = min(qst,0.0) - - zt = ht * exp(-beta*(ustar**0.5)*(abs(tstar2)**1.0)) - zq = ht * exp(-beta*(ustar**0.5)*(abs(qstar2)**1.0)) -!zq = zt - - zt = min(zt, z_0/2.0) - zq = min(zq, z_0/2.0) - - end subroutine yang_2008 - -!================================================================================================================= - subroutine andreas_2002(z_0,bvisc,ustar,zt,zq) -! This is taken from Andreas (2002; J. of Hydromet) and -! Andreas et al. (2005; BLM). -! -! This should only be used over snow/ice! - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: z_0,bvisc,ustar - -!--- output arguments: - real(kind=kind_phys),intent(out):: zt, zq - -!--- local variables: - real(kind=kind_phys):: ren2,zntsno - - real(kind=kind_phys),parameter:: bt0_s=1.25, bt0_t=0.149, bt0_r=0.317, & - bt1_s=0.0, bt1_t=-0.55, bt1_r=-0.565, & - bt2_s=0.0, bt2_t=0.0, bt2_r=-0.183 - - real(kind=kind_phys),parameter:: bq0_s=1.61, bq0_t=0.351, bq0_r=0.396, & - bq1_s=0.0, bq1_t=-0.628, bq1_r=-0.512, & - bq2_s=0.0, bq2_t=0.0, bq2_r=-0.180 - -!----------------------------------------------------------------------------------------------------------------- - -!calculate zo for snow (Andreas et al. 2005, BLM): - zntsno = 0.135*bvisc/ustar + & - (0.035*(ustar*ustar)/9.8) * & - (5.*exp(-1.*(((ustar - 0.18)/0.1)*((ustar - 0.18)/0.1))) + 1.) - ren2 = ustar*zntsno/bvisc - -!Make sure that Re is not outside of the range of validity -!for using their equations - if (ren2 .gt. 1000.) ren2 = 1000. - - if (ren2 .le. 0.135) then - - zt = zntsno*exp(bt0_s + bt1_s*log(ren2) + bt2_s*log(ren2)**2) - zq = zntsno*exp(bq0_s + bq1_s*log(ren2) + bq2_s*log(ren2)**2) - - else if (ren2 .gt. 0.135 .and. ren2 .lt. 2.5) then - - zt = zntsno*exp(bt0_t + bt1_t*log(ren2) + bt2_t*log(ren2)**2) - zq = zntsno*exp(bq0_t + bq1_t*log(ren2) + bq2_t*log(ren2)**2) - - else - - zt = zntsno*exp(bt0_r + bt1_r*log(ren2) + bt2_r*log(ren2)**2) - zq = zntsno*exp(bq0_r + bq1_r*log(ren2) + bq2_r*log(ren2)**2) - - endif - - end subroutine andreas_2002 - -!================================================================================================================= - subroutine psi_hogstrom_1996(psi_m,psi_h,zl,zt,z_0,za) -!this subroutine returns the stability functions based off -!of hogstrom (1996). - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl,zt,z_0,za - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m,psi_h - -!--- local variables: - real(kind=kind_phys):: x,x0,y,y0,zml,zhl - -!----------------------------------------------------------------------------------------------------------------- - - zml = z_0*zl/za - zhl = zt*zl/za - - if (zl .gt. 0.) then !stable (not well tested - seem large) - - psi_m = -5.3*(zl - zml) - psi_h = -8.0*(zl - zhl) - - else !unstable - - x = (1.-19.0*zl)**0.25 - x0= (1.-19.0*zml)**0.25 - y = (1.-11.6*zl)**0.5 - y0= (1.-11.6*zhl)**0.5 - - psi_m = 2.*log((1.+x)/(1.+x0)) + & - &log((1.+x**2.)/(1.+x0**2.)) - & - &2.0*atan(x) + 2.0*atan(x0) - psi_h = 2.*log((1.+y)/(1.+y0)) - - endif - - end subroutine psi_hogstrom_1996 - -!================================================================================================================= - subroutine psi_dyerhicks(psi_m,psi_h,zl,zt,z_0,za) -!This subroutine returns the stability functions based off -!of Hogstrom (1996), but with different constants compatible -!with Dyer and Hicks (1970/74?). This formulation is used for -!testing/development by Nakanishi (personal communication). - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl,zt,z_0,za - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m,psi_h - -!--- local variables: - real(kind=kind_phys):: x,x0,y,y0,zml,zhl - -!----------------------------------------------------------------------------------------------------------------- - - zml = z_0*zl/za !zo/l - zhl = zt*zl/za !zt/l - - if (zl .gt. 0.) then !stable - - psi_m = -5.0*(zl - zml) - psi_h = -5.0*(zl - zhl) - - else !unstable - - x = (1.-16.*zl)**0.25 - x0= (1.-16.*zml)**0.25 - - y = (1.-16.*zl)**0.5 - y0= (1.-16.*zhl)**0.5 - - psi_m = 2.*log((1.+x)/(1.+x0)) + & - &log((1.+x**2.)/(1.+x0**2.)) - & - &2.0*atan(x) + 2.0*atan(x0) - psi_h = 2.*log((1.+y)/(1.+y0)) - - endif - - end subroutine psi_dyerhicks - -!================================================================================================================= - subroutine psi_beljaars_holtslag_1991(psi_m,psi_h,zl) -!this subroutine returns the stability functions based off -!of Beljaar and Holtslag 1991, which is an extension of Holtslag -!and Debruin 1989. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m, psi_h - -!--- local variables: - real(kind=kind_phys):: a=1.,b=0.666,c=5.,d=0.35 - -!----------------------------------------------------------------------------------------------------------------- - - if (zl .lt. 0.) then !unstable - - write(*,*)"WARNING: Universal stability functions from" - write(*,*)" Beljaars and Holtslag (1991) should only" - write(*,*)" be used in the stable regime!" - psi_m = 0. - psi_h = 0. - - else !stable - - psi_m = -(a*zl + b*(zl -(c/d))*exp(-d*zl) + (b*c/d)) - psi_h = -((1.+.666*a*zl)**1.5 + & - b*(zl - (c/d))*exp(-d*zl) + (b*c/d) -1.) - - endif - - end subroutine psi_beljaars_holtslag_1991 - -!================================================================================================================= - subroutine psi_zilitinkevich_esau_2007(psi_m,psi_h,zl) -!this subroutine returns the stability functions come from -!Zilitinkevich and Esau (2007, BM), which are formulatioed from the -!"generalized similarity theory" and tuned to the LES DATABASE64 -!to determine their dependence on z/L. - IMPLICIT NONE -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m, psi_h - -!--- local variables: - real(kind=kind_phys),parameter:: cm=3.0,ct=2.5 - -!----------------------------------------------------------------------------------------------------------------- - - if (zl .lt. 0.) then !unstable - -! write(*,*)"WARNING: Universal stability function from" -! write(*,*)" Zilitinkevich and Esau (2007) should only" -! write(*,*)" be used in the stable regime!" - psi_m = 0. - psi_h = 0. - - else !stable - - psi_m = -cm*(zl**(5./6.)) - psi_h = -ct*(zl**(4./5.)) - - endif - - end subroutine psi_zilitinkevich_esau_2007 - -!================================================================================================================= - subroutine psi_businger_1971(psi_m,psi_h,zl) -!this subroutine returns the flux-profile relationships -!of Businger el al. 1971. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m, psi_h - -!--- local variables: - real(kind=kind_phys):: x, y - real(kind=kind_phys),parameter:: pi180 = 3.14159265/180. - -!----------------------------------------------------------------------------------------------------------------- - - if (zl .lt. 0.) then !unstable - - x = (1. - 15.0*zl)**0.25 - y = (1. - 9.0*zl)**0.5 - - psi_m = log(((1.+x)/2.)**2.) + & - & log((1.+x**2.)/2.) - & - & 2.0*atan(x) + pi180*90. - psi_h = 2.*log((1.+y)/2.) - - else !stable - - psi_m = -4.7*zl - psi_h = -(4.7/0.74)*zl - - endif - - end subroutine psi_businger_1971 - -!================================================================================================================= - subroutine psi_suselj_sood_2010(psi_m,psi_h,zl) -!this subroutine returns flux-profile relatioships based off -!of Lobocki (1993), which is derived from the MY-level 2 model. -!Suselj and Sood (2010) applied the surface layer length scales -!from Nakanishi (2001) to get this new relationship. These functions -!are more agressive (larger magnitude) than most formulations. They -!showed improvement over water, but untested over land. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl - -!--- output arguments: - real(kind=kind_phys),intent(out):: psi_m, psi_h - -!--- local variables: - real(kind=kind_phys),parameter:: rfc=0.19, ric=0.183, phit=0.8 - -!----------------------------------------------------------------------------------------------------------------- - - if (zl .gt. 0.) then !stable - - psi_m = -(zl/rfc + 1.1223*exp(1.-1.6666/zl)) - !psi_h = -zl*ric/((rfc**2.)*phit) + 8.209*(zl**1.1091) - !their eq for psi_h crashes the model and does not match - !their fig 1. this eq (below) matches their fig 1 better: - psi_h = -(zl*ric/((rfc**2.)*5.) + 7.09*(zl**1.1091)) - - else !unstable - - psi_m = 0.9904*log(1. - 14.264*zl) - psi_h = 1.0103*log(1. - 16.3066*zl) - - endif - - end subroutine psi_suselj_sood_2010 - -!================================================================================================================= - subroutine psi_cb2005(psim1,psih1,zl,z0l) -!this subroutine returns the stability functions based off -!of Cheng and Brutseart (2005, BLM), for use in stable conditions only. -!the returned values are the combination of psi((za+zo)/L) - psi(z0/L) - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: zl,z0l - -!--- output arguments: - real(kind=kind_phys),intent(out):: psim1,psih1 - -!----------------------------------------------------------------------------------------------------------------- - - psim1 = -6.1*log(zl + (1.+ zl **2.5)**0.4) & - -6.1*log(z0l + (1.+ z0l**2.5)**0.4) - psih1 = -5.5*log(zl + (1.+ zl **1.1)**0.90909090909) & - -5.5*log(z0l + (1.+ z0l**1.1)**0.90909090909) - - end subroutine psi_cb2005 - -!================================================================================================================= - subroutine li_etal_2010(zl,rib,zaz0,z0zt) -!this subroutine returns a more robust z/l that best matches -!the z/l from hogstrom (1996) for unstable conditions and beljaars -!and holtslag (1991) for stable conditions. - implicit none -!================================================================================================================= - -!--- input arguments: - real(kind=kind_phys),intent(in):: rib,zaz0,z0zt - -!--- output arguments: - real(kind=kind_phys),intent(out):: zl - -!--- local variables: - real(kind=kind_phys):: alfa,beta,zaz02,z0zt2 - - real(kind=kind_phys),parameter:: au11=0.045, bu11=0.003, bu12=0.0059, & - &bu21=-0.0828, bu22=0.8845, bu31=0.1739, & - &bu32=-0.9213, bu33=-0.1057 - real(kind=kind_phys),parameter:: aw11=0.5738, aw12=-0.4399, aw21=-4.901, & - &aw22=52.50, bw11=-0.0539, bw12=1.540, & - &bw21=-0.669, bw22=-3.282 - real(kind=kind_phys),parameter:: as11=0.7529, as21=14.94, bs11=0.1569, & - &bs21=-0.3091, bs22=-1.303 - -!----------------------------------------------------------------------------------------------------------------- - -!set limits according to Li et al (2010), p 157. - zaz02=zaz0 - if (zaz0 .lt. 100.0) zaz02=100. - if (zaz0 .gt. 100000.0) zaz02=100000. - -!set more limits according to Li et al (2010) - z0zt2=z0zt - if (z0zt .lt. 0.5) z0zt2=0.5 - if (z0zt .gt. 100.0) z0zt2=100. - - alfa = log(zaz02) - beta = log(z0zt2) - - if (rib .le. 0.0) then - zl = au11*alfa*rib**2 + ( & - & (bu11*beta + bu12)*alfa**2 + & - & (bu21*beta + bu22)*alfa + & - & (bu31*beta**2 + bu32*beta + bu33))*rib - - !if(zL .LT. -15 .OR. zl .GT. 0.)print*,"VIOLATION Rib<0:",zL - zl = max(zl,-15.) !limits set according to Li et al (2010) - zl = min(zl,0.) !Figure 1. - elseif (rib .gt. 0.0 .and. rib .le. 0.2) then - zl = ((aw11*beta + aw12)*alfa + & - & (aw21*beta + aw22))*rib**2 + & - & ((bw11*beta + bw12)*alfa + & - & (bw21*beta + bw22))*rib - - !if(zl .lt. 0 .or. zl .gt. 4)print*,"violation 00.2:",zl - zl = min(zl,20.) !limits according to Li et al (2010), their Figure 1c. - zl = max(zl,1.) - endif - - end subroutine li_etal_2010 - -!================================================================================================================= - real(kind=kind_phys) function zolri(ri,za,z0,zt,zol1,psi_opt) - implicit none -! This iterative algorithm is a two-point secant method taken from the revised -! surface layer scheme in WRF-ARW, written by Pedro Jimenez and Jimy Dudhia and -! summarized in Jimenez et al. (2012, MWR). This function was adapted -! to input the thermal roughness length, zt, (as well as z0) and use initial -! estimate of z/L. -!================================================================================================================= - -!--- input arguments: - integer, intent(in):: psi_opt - real(kind=kind_phys),intent(in):: ri,za,z0,zt,zol1 - -!--- local variables and arrays: - integer:: n - integer,parameter:: nmax = 20 - real(kind=kind_phys):: x1,x2,fx1,fx2 - -!----------------------------------------------------------------------------------------------------------------- - - if (ri.lt.0.)then - x1=zol1 - 0.02 !-5. - x2=0. - else - x1=0. - x2=zol1 + 0.02 !5. - endif - - n=0 - fx1=zolri2(x1,ri,za,z0,zt,psi_opt) - fx2=zolri2(x2,ri,za,z0,zt,psi_opt) - - do while (abs(x1 - x2) > 0.01 .and. n < nmax) - if(abs(fx2) .lt. abs(fx1))then - x1=x1-fx1/(fx2-fx1)*(x2-x1) - fx1=zolri2(x1,ri,za,z0,zt,psi_opt) - zolri=x1 - else - x2=x2-fx2/(fx2-fx1)*(x2-x1) - fx2=zolri2(x2,ri,za,z0,zt,psi_opt) - zolri=x2 - endif - n=n+1 - enddo - - if (n==nmax .and. abs(x1 - x2) >= 0.01) then - !if convergence fails, use approximate values: - call li_etal_2010(zolri, ri, za/z0, z0/zt) - !print*,"failed, n=",n," ri=",ri," zt=",zt - else - !print*,"success,n=",n," ri=",ri," z/l=",zolri - endif - - end function zolri - -!================================================================================================================= - real(kind=kind_phys) function zolri2(zol2,ri2,za,z0,zt,psi_opt) - implicit none -! input: ================================= -! zol2 - estimated z/l -! ri2 - calculated bulk richardson number -! za - 1/2 depth of first model layer -! z0 - aerodynamic roughness length -! zt - thermal roughness length -! output: ================================ -! zolri2 - delta ri -!================================================================================================================= - -!--- input arguments: - integer,intent(in):: psi_opt - real(kind=kind_phys),intent(in):: ri2,za,z0,zt - -!--- inout arguments: - real(kind=kind_phys),intent(inout):: zol2 - -!--- local variables and arrays: - real(kind=kind_phys):: zol20,zol3,psim1,psih1,psix2,psit2,zolt - -!----------------------------------------------------------------------------------------------------------------- - - if(zol2*ri2 .lt. 0.) then - !print*,"wrong quadrants: z/l=",zol2," ri=",ri2 - zol2=0. - endif - - zol20=zol2*z0/za ! z0/l - zol3=zol2+zol20 ! (z+z0)/l - zolt=zol2*zt/za ! zt/l - - if (ri2.lt.0) then - psit2=max(log((za+z0)/zt)-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0) - psix2=max(log((za+z0)/z0)-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)),1.0) - else - psit2=max(log((za+z0)/zt)-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0) - psix2=max(log((za+z0)/z0)-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)),1.0) - endif - - zolri2=zol2*psit2/psix2**2 - ri2 -!print*," target ri=",ri2," est ri=",zol2*psit2/psix2**2 - - end function zolri2 - -!================================================================================================================= - real(kind=kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt) - implicit none -!this iterative algorithm to compute z/L from bulk-Ri -!================================================================================================================= - -!--- input arguments: - integer,intent(in):: psi_opt - real(kind=kind_phys),intent(in):: ri,za,z0,zt,logz0,logzt - -!--- inout arguments: - real(kind=kind_phys),intent(inout):: zol1 - -!--- local variables and arrays: - integer:: n - integer,parameter :: nmax = 20 - real(kind=kind_phys):: zol20,zol3,zolt,zolold - real(kind=kind_phys):: psit2,psix2 -!real(kind=kind_phys),dimension(nmax):: zlhux - -!----------------------------------------------------------------------------------------------------------------- - - if(zol1*ri .lt. 0.) then -! print*,"WRONG QUADRANTS: z/L=",zol1," ri=",ri - zol1=0. - endif - - if (ri .lt. 0.) then - zolold=-99999. - zolrib=-66666. - else - zolold=99999. - zolrib=66666. - endif - - n=1 - do while (abs(zolold - zolrib) > 0.01 .and. n < nmax) - - if(n==1)then - zolold=zol1 - else - zolold=zolrib - endif - zol20=zolold*z0/za ! z0/L - zol3=zolold+zol20 ! (z+z0)/L - zolt=zolold*zt/za ! zt/L - - if (ri.lt.0) then - psit2=MAX(logzt-(psih_unstable(zol3,psi_opt)-psih_unstable(zolt,psi_opt)), 1.0) - psix2=MAX(logz0-(psim_unstable(zol3,psi_opt)-psim_unstable(zol20,psi_opt)), 1.0) - else - psit2=MAX(logzt-(psih_stable(zol3,psi_opt)-psih_stable(zolt,psi_opt)), 1.0) - psix2=MAX(logz0-(psim_stable(zol3,psi_opt)-psim_stable(zol20,psi_opt)), 1.0) - endif - - zolrib=ri*psix2**2/psit2 - !zLhux(n)=zolrib - n=n+1 - enddo - - if (n==nmax .and. abs(zolold - zolrib) > 0.01 ) then - !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri - !if convergence fails, use approximate values: - call li_etal_2010(zolrib,ri,za/z0,z0/zt) - !zLhux(n)=zolri - !print*,"FAILED, n=",n," Ri=",ri," zt=",zt - !print*,"z/L=",zLhux(1:nmax) - else - !print*,"SUCCESS,n=",n," Ri=",ri," z/L=",zolrib - endif - - end function zolrib - -!================================================================================================================= - subroutine psi_init(psi_opt) - implicit none -!define tables from -10 <= z/L <= 10 -!================================================================================================================= - - integer,intent(in):: psi_opt - integer:: n - real(kind=kind_phys):: zolf - -!----------------------------------------------------------------------------------------------------------------- - - if (psi_opt == 0) then - do n = 0,1000 - !stable function tables - zolf = float(n)*0.01 - psim_stab(n)=psim_stable_full(zolf) - psih_stab(n)=psih_stable_full(zolf) - - !unstable function tables - zolf = -float(n)*0.01 - psim_unstab(n)=psim_unstable_full(zolf) - psih_unstab(n)=psih_unstable_full(zolf) - enddo - else - do n = 0,1000 - !stable function tables - zolf = float(n)*0.01 - psim_stab(n)=psim_stable_full_gfs(zolf) - psih_stab(n)=psih_stable_full_gfs(zolf) - - !unstable function tables - zolf = -float(n)*0.01 - psim_unstab(n)=psim_unstable_full_gfs(zolf) - psih_unstab(n)=psih_unstable_full_gfs(zolf) - enddo - endif - - end subroutine psi_init - -!================================================================================================================= -! ... Full equations for the integrated similarity functions ... -!================================================================================================================= - real(kind=kind_phys) function psim_stable_full(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - - psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5)) - - end function psim_stable_full - -!================================================================================================================= - real(kind=kind_phys) function psih_stable_full(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - - psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1)) - - end function psih_stable_full - -!================================================================================================================= - real(kind=kind_phys) function psim_unstable_full(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: x,ym,psimc,psimk - - x=(1.-16.*zolf)**.25 - psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.) - - ym=(1.-10.*zolf)**0.33 - psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) - - psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.) - - end function psim_unstable_full - -!================================================================================================================= - real(kind=kind_phys) function psih_unstable_full(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: y,yh,psihc,psihk - - y=(1.-16.*zolf)**.5 - psihk=2.*log((1+y)/2.) - - yh=(1.-34.*zolf)**0.33 - psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.) - - psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2.) - - end function psih_unstable_full - -!================================================================================================================= -! ... integrated similarity functions from GFS... -! -!================================================================================================================= - real(kind=kind_phys) function psim_stable_full_gfs(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: aa - real(kind=kind_phys),parameter:: alpha4 = 20. - - aa = sqrt(1. + alpha4 * zolf) - psim_stable_full_gfs = -1.*aa + log(aa + 1.) - - end function psim_stable_full_gfs - -!================================================================================================================= - real(kind=kind_phys) function psih_stable_full_gfs(zolf) - implicit none - - real(kind=kind_phys):: zolf - real(kind=kind_phys):: bb - real(kind=kind_phys),parameter:: alpha4 = 20. - - bb = sqrt(1. + alpha4 * zolf) - psih_stable_full_gfs = -1.*bb + log(bb + 1.) - - end function psih_stable_full_gfs - -!================================================================================================================= - real(kind=kind_phys) function psim_unstable_full_gfs(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: hl1,tem1 - real(kind=kind_phys),parameter:: a0=-3.975, a1=12.32, & - b1=-7.755, b2=6.041 - - if (zolf .ge. -0.5) then - hl1 = zolf - psim_unstable_full_gfs = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) - else - hl1 = -zolf - tem1 = 1.0 / sqrt(hl1) - psim_unstable_full_gfs = log(hl1) + 2. * sqrt(tem1) - .8776 - end if - - end function psim_unstable_full_gfs - -!================================================================================================================= - real(kind=kind_phys) function psih_unstable_full_gfs(zolf) - implicit none - - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: hl1,tem1 - real(kind=kind_phys),parameter:: a0p=-7.941, a1p=24.75, & - b1p=-8.705, b2p=7.899 - - if (zolf .ge. -0.5) then - hl1 = zolf - psih_unstable_full_gfs = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) - else - hl1 = -zolf - tem1 = 1.0 / sqrt(hl1) - psih_unstable_full_gfs = log(hl1) + .5 * tem1 + 1.386 - end if - - end function psih_unstable_full_gfs - -!================================================================================================================= -! These functions use the look-up table functions when |z/L| <= 10 -! but default to the full equations when |z/L| > 10. -!================================================================================================================= - real(kind=kind_phys) function psim_stable(zolf,psi_opt) - implicit none - - integer,intent(in):: psi_opt - integer:: nzol - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: rzol - -!----------------------------------------------------------------------------------------------------------------- - - nzol = int(zolf*100.) - rzol = zolf*100. - nzol - if(nzol+1 .le. 1000)then - psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol)) - else - if (psi_opt == 0) then - psim_stable = psim_stable_full(zolf) - else - psim_stable = psim_stable_full_gfs(zolf) - endif - endif - - end function psim_stable - -!================================================================================================================= - real(kind=kind_phys) function psih_stable(zolf,psi_opt) - implicit none - - integer,intent(in):: psi_opt - integer:: nzol - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: rzol - -!----------------------------------------------------------------------------------------------------------------- - - nzol = int(zolf*100.) - rzol = zolf*100. - nzol - if(nzol+1 .le. 1000)then - psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol)) - else - if (psi_opt == 0) then - psih_stable = psih_stable_full(zolf) - else - psih_stable = psih_stable_full_gfs(zolf) - endif - endif - - end function psih_stable - -!================================================================================================================= - real(kind=kind_phys) function psim_unstable(zolf,psi_opt) - implicit none - - integer,intent(in):: psi_opt - integer:: nzol - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: rzol - -!----------------------------------------------------------------------------------------------------------------- - - nzol = int(-zolf*100.) - rzol = -zolf*100. - nzol - if(nzol+1 .le. 1000)then - psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol)) - else - if (psi_opt == 0) then - psim_unstable = psim_unstable_full(zolf) - else - psim_unstable = psim_unstable_full_gfs(zolf) - endif - endif - - end function psim_unstable - -!================================================================================================================= - real(kind=kind_phys) function psih_unstable(zolf,psi_opt) - implicit none - - integer,intent(in):: psi_opt - integer:: nzol - real(kind=kind_phys),intent(in):: zolf - real(kind=kind_phys):: rzol - -!----------------------------------------------------------------------------------------------------------------- - - nzol = int(-zolf*100.) - rzol = -zolf*100. - nzol - if(nzol+1 .le. 1000)then - psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol)) - else - if (psi_opt == 0) then - psih_unstable = psih_unstable_full(zolf) - else - psih_unstable = psih_unstable_full_gfs(zolf) - endif - endif - - end function psih_unstable - -!================================================================================================================= - end module module_mynnsfclay -!================================================================================================================= - diff --git a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_driver.F90 b/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_driver.F90 deleted file mode 100644 index 1fdee508c7..0000000000 --- a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_driver.F90 +++ /dev/null @@ -1,471 +0,0 @@ -!================================================================================================================= - module module_mynnsfclay_driver -!================================================================================================================= - use mpas_kind_types,only: RKIND,StrKIND - - use module_mynnsfclay,only: mynnsfclay - use module_mynnsfclay_pre,only: mynnsfclay_pre - - - implicit none - private - public:: mynnsfclay_driver - - - contains - - -!================================================================================================================= - subroutine mynnsfclay_driver( & - u3d,v3d,t3d,qv3d,p3d,dz8w, & - cp,g,rovcp,r,xlv,psfcpa,chs,chs2,cqs,cqs2,cpm,& - znt,ust,pblh,mavail,zol,mol,regime,psim,psih, & - xland,hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc,rmol, & - u10,v10,th2,t2,q2,snowh, & - gz1oz0,wspd,br,isfflx,dx, & - svp1,svp2,svp3,svpt0,ep1,ep2, & - karman,itimestep,ch,th3d,pi3d,qc3d,rho3d,qcg, & - spp_pbl,pattern_spp_pbl, & - ustm,ck,cka,cd,cda,isftcflx,iz0tlnd, & - ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte, & - restart,cycling,initflag,errmsg,errflg & - ) -!------------------------------------------------------------------- -!-- u3d 3d u-velocity interpolated to theta points (m/s) -!-- v3d 3d v-velocity interpolated to theta points (m/s) -!-- t3d 3d temperature (k) -!-- qv3d 3d water vapor mixing ratio (kg/kg) -!-- p3d 3d pressure (pa) -!-- rho3d 3d density (kg/m3) -!-- dz8w 3d dz between full levels (m) -!-- cp heat capacity at constant pressure for dry air (j/kg/k) -!-- g acceleration due to gravity (m/s^2) -!-- rovcp r/cp -!-- r gas constant for dry air (j/kg/k) -!-- xlv latent heat of vaporization for water (j/kg) -!-- psfcpa surface pressure (pa) -!-- znt roughness length (m) -!-- ust u* in similarity theory (m/s) -!-- ustm u* in similarity theory (m/s) w* added to wspd. this is -! used to couple with tke scheme but not in mynn. -! (as of now, ustm = ust in this version) -!-- pblh pbl height from previous time (m) -!-- mavail surface moisture availability (between 0 and 1) -!-- zol z/l height over monin-obukhov length -!-- mol t* (similarity theory) (k) -!-- rmol reciprocal of m-o length (/m) -!-- regime flag indicating pbl regime (stable, unstable, etc.) -!-- psim similarity stability function for momentum -!-- psih similarity stability function for heat -!-- xland land mask (1 for land, 2 for water) -!-- hfx upward heat flux at the surface (w/m^2) -!-- qfx upward moisture flux at the surface (kg/m^2/s) -!-- lh net upward latent heat flux at surface (w/m^2) -!-- tsk surface temperature (k) -!-- flhc exchange coefficient for heat (w/m^2/k) -!-- flqc exchange coefficient for moisture (kg/m^2/s) -!-- chs heat/moisture exchange coefficient for lsm (m/s) -!-- qgh lowest-level saturated mixing ratio -!-- qsfc qv (specific humidity) at the surface -!-- qsfcmr qv (mixing ratio) at the surface -!-- u10 diagnostic 10m u wind -!-- v10 diagnostic 10m v wind -!-- th2 diagnostic 2m theta (k) -!-- t2 diagnostic 2m temperature (k) -!-- q2 diagnostic 2m mixing ratio (kg/kg) -!-- snowh snow height (m) -!-- gz1oz0 log((z1+znt)/znt) where znt is roughness length -!-- wspd wind speed at lowest model level (m/s) -!-- br bulk richardson number in surface layer -!-- isfflx isfflx=1 for surface heat and moisture fluxes -!-- dx horizontal grid size (m) -!-- svp1 constant for saturation vapor pressure (=0.6112 kpa) -!-- svp2 constant for saturation vapor pressure (=17.67 dimensionless) -!-- svp3 constant for saturation vapor pressure (=29.65 k) -!-- svpt0 constant for saturation vapor pressure (=273.15 k) -!-- ep1 constant for virtual temperature (rv/rd - 1) (dimensionless) -!-- ep2 constant for spec. hum. calc (rd/rv = 0.622) (dimensionless) -!-- ep3 constant for spec. hum. calc (1 - rd/rv = 0.378 ) (dimensionless) -!-- karman von karman constant -!-- ck enthalpy exchange coeff at 10 meters -!-- cd momentum exchange coeff at 10 meters -!-- cka enthalpy exchange coeff at the lowest model level -!-- cda momentum exchange coeff at the lowest model level -!-- isftcflx =0: z0, zt, and zq from coare3.0/3.5 (fairall et al 2003/edson et al 2013) -! (water =1: z0 from davis et al (2008), zt & zq from coare3.0/3.5 -! only) =2: z0 from davis et al (2008), zt & zq from garratt (1992) -! =3: z0 from taylor and yelland (2004), zt and zq from coare 3.0/3.5 -!-- iz0tlnd =0: zilitinkevich (1995) with czil=0.085, -! (land =1: czil_new (modified according to chen & zhang 2008) -! only) =2: modified yang et al (2002, 2008) - generalized for all landuse -! =3: constant zt = z0/7.4 (garratt 1992) -!-- u3d 3d u-velocity interpolated to theta points (m/s) -!-- v3d 3d v-velocity interpolated to theta points (m/s) -!-- t3d 3d temperature (k) -!-- qv3d 3d water vapor mixing ratio (kg/kg) -!-- p3d 3d pressure (pa) -!-- rho3d 3d density (kg/m3) -!-- dz8w 3d dz between full levels (m) -!-- cp heat capacity at constant pressure for dry air (j/kg/k) -!-- g acceleration due to gravity (m/s^2) -!-- rovcp r/cp -!-- r gas constant for dry air (j/kg/k) -!-- xlv latent heat of vaporization for water (j/kg) -!-- psfcpa surface pressure (pa) -!-- znt roughness length (m) -!-- ust u* in similarity theory (m/s) -!-- ustm u* in similarity theory (m/s) w* added to wspd. this is -! used to couple with tke scheme but not in mynn. -! (as of now, ustm = ust in this version) -!-- pblh pbl height from previous time (m) -!-- mavail surface moisture availability (between 0 and 1) -!-- zol z/l height over monin-obukhov length -!-- mol t* (similarity theory) (k) -!-- rmol reciprocal of m-o length (/m) -!-- regime flag indicating pbl regime (stable, unstable, etc.) -!-- psim similarity stability function for momentum -!-- psih similarity stability function for heat -!-- xland land mask (1 for land, 2 for water) -!-- hfx upward heat flux at the surface (w/m^2) -!-- qfx upward moisture flux at the surface (kg/m^2/s) -!-- lh net upward latent heat flux at surface (w/m^2) -!-- tsk surface temperature (k) -!-- flhc exchange coefficient for heat (w/m^2/k) -!-- flqc exchange coefficient for moisture (kg/m^2/s) -!-- chs heat/moisture exchange coefficient for lsm (m/s) -!-- qgh lowest-level saturated mixing ratio -!-- qsfc qv (specific humidity) at the surface -!-- qsfcmr qv (mixing ratio) at the surface -!-- u10 diagnostic 10m u wind -!-- v10 diagnostic 10m v wind -!-- th2 diagnostic 2m theta (k) -!-- t2 diagnostic 2m temperature (k) -!-- q2 diagnostic 2m mixing ratio (kg/kg) -!-- snowh snow height (m) -!-- gz1oz0 log((z1+znt)/znt) where znt is roughness length -!-- wspd wind speed at lowest model level (m/s) -!-- br bulk richardson number in surface layer -!-- isfflx isfflx=1 for surface heat and moisture fluxes -!-- dx horizontal grid size (m) -!-- svp1 constant for saturation vapor pressure (=0.6112 kpa) -!-- svp2 constant for saturation vapor pressure (=17.67 dimensionless) -!-- svp3 constant for saturation vapor pressure (=29.65 k) -!-- svpt0 constant for saturation vapor pressure (=273.15 k) -!-- ep1 constant for virtual temperature (rv/rd - 1) (dimensionless) -!-- ep2 constant for spec. hum. calc (rd/rv = 0.622) (dimensionless) -!-- ep3 constant for spec. hum. calc (1 - rd/rv = 0.378 ) (dimensionless) -!-- karman von karman constant -!-- ck enthalpy exchange coeff at 10 meters -!-- cd momentum exchange coeff at 10 meters -!-- cka enthalpy exchange coeff at the lowest model level -!-- cda momentum exchange coeff at the lowest model level -!-- isftcflx =0: z0, zt, and zq from coare3.0/3.5 (fairall et al 2003/edson et al 2013) -! (water =1: z0 from davis et al (2008), zt & zq from coare3.0/3.5 -! only) =2: z0 from davis et al (2008), zt & zq from garratt (1992) -! =3: z0 from taylor and yelland (2004), zt and zq from coare 3.0/3.5 -!-- iz0tlnd =0: zilitinkevich (1995) with czil=0.085, -! (land =1: czil_new (modified according to chen & zhang 2008) -! only) =2: modified yang et al (2002, 2008) - generalized for all landuse -! =3: constant zt = z0/7.4 (garratt 1992) -! -!-- ids start index for i in domain -!-- ide end index for i in domain -!-- jds start index for j in domain -!-- jde end index for j in domain -!-- kds start index for k in domain -!-- kde end index for k in domain -!-- ims start index for i in memory -!-- ime end index for i in memory -!-- jms start index for j in memory -!-- jme end index for j in memory -!-- kms start index for k in memory -!-- kme end index for k in memory -!-- its start index for i in tile -!-- ite end index for i in tile -!-- jts start index for j in tile -!-- jte end index for j in tile -!-- kts start index for k in tile -!-- kte end index for k in tile -!================================================================================================================= - -!--- input arguments: - integer,intent(in):: ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte - logical,intent(in):: restart,cycling - integer,intent(in):: itimestep - integer,intent(in):: initflag - integer,intent(in):: isfflx - integer,intent(in),optional:: isftcflx, iz0tlnd - integer,intent(in),optional:: spp_pbl - - real(kind=RKIND),intent(in):: svp1,svp2,svp3,svpt0 - real(kind=RKIND),intent(in):: ep1,ep2,karman - real(kind=RKIND),intent(in):: cp,g,rovcp,r,xlv - - real(kind=RKIND),intent(in),dimension(ims:ime,kms:kme,jms:jme):: & - dz8w, & - qv3d, & - p3d, & - t3d, & - qc3d, & - u3d, & - v3d, & - rho3d, & - th3d, & - pi3d - - real(kind=RKIND),intent(in),dimension(ims:ime,kms:kme,jms:jme),optional:: & - pattern_spp_pbl - - real(kind=RKIND),intent(in),dimension(ims:ime,jms:jme):: & - mavail, & - pblh, & - xland, & - tsk, & - qcg, & - psfcpa, & - snowh, & - dx - -!--- output arguments: - character(len=StrKIND),intent(out):: errmsg - integer,intent(out):: errflg - - real(kind=RKIND),intent(out),dimension(ims:ime,jms:jme):: & - u10, & - v10, & - th2, & - t2, & - q2 - - real(kind=RKIND),intent(out),dimension(ims:ime,jms:jme),optional:: & - ck, & - cka, & - cd, & - cda, & - ustm - -!--- inout arguments: - real(kind=RKIND),intent(inout),dimension(ims:ime,jms:jme):: & - regime, & - hfx, & - qfx, & - lh, & - mol, & - rmol, & - qsfc, & - qgh, & - znt, & - zol, & - ust, & - cpm, & - chs2, & - cqs2, & - cqs, & - chs, & - ch, & - flhc, & - flqc, & - gz1oz0, & - wspd, & - br, & - psim, & - psih - -!--- local variables and arrays: - integer:: i,j,k - logical:: restart_or_cycle - - real(kind=RKIND),dimension(its:ite):: & - u1d,v1d,u1d2,v1d2,qv1d,p1d,t1d,qc1d,rho1d,dz8w1d,dz2w1d - - real(kind=RKIND),dimension(its:ite):: rstoch1d - - real(kind=RKIND),dimension(ims:ime,jms:jme):: qstar,wstar - -!intermediate variables and arrays to accomodate the CCPP-compliant sourcecode: - logical:: f_spp - - real(kind=RKIND),dimension(its:ite):: mavail_hv,pblh_hv,xland_hv,tsk_hv,psfcpa_hv, & - qcg_hv,snowh_hv,dx_hv - real(kind=RKIND),dimension(its:ite):: regime_hv,hfx_hv,qfx_hv,lh_hv,mol_hv,rmol_hv, & - qgh_hv,qsfc_hv,znt_hv,zol_hv,ust_hv,cpm_hv,chs2_hv, & - cqs_hv,cqs2_hv,chs_hv,ch_hv,flhc_hv,flqc_hv,gz1oz0_hv, & - wspd_hv,br_hv,psim_hv,psih_hv - real(kind=RKIND),dimension(its:ite):: u10_hv,v10_hv,th2_hv,t2_hv,q2_hv,wstar_hv,qstar_hv - real(kind=RKIND),dimension(its:ite):: cd_hv,cda_hv,ck_hv,cka_hv,ustm_hv - -!----------------------------------------------------------------------------------------------------------------- - - f_spp = .false. - if(spp_pbl==1 .and. present(pattern_spp_pbl)) f_spp = .true. - restart_or_cycle = .false. - if(restart .or. cycling)restart_or_cycle=.true. - - errmsg = ' ' - errflg = 0 - - do j = jts,jte - - !initialization of arrays ust,mol,qsfc,and qstar that are initialized if itimestep equals 1: - do i = its,ite - ust_hv(i) = ust(i,j) - mol_hv(i) = mol(i,j) - qsfc_hv(i) = qsfc(i,j) - qstar_hv(i) = qstar(i,j) - enddo - - call mynnsfclay_pre(its,ite,kte,initflag,dz8w,u3d,v3d,p3d,t3d,rho3d,qv3d,qc3d,f_spp, & - pattern_spp_pbl,ust_hv,mol_hv,qsfc_hv,qstar_hv,dz8w1d,u1d,v1d,p1d,t1d,rho1d,& - qv1d,qc1d,rstoch1d,dz2w1d,u1d2,v1d2,errmsg,errflg) - - !input arguments: - do i = its,ite - mavail_hv(i) = mavail(i,j) - pblh_hv(i) = pblh(i,j) - xland_hv(i) = xland(i,j) - tsk_hv(i) = tsk(i,j) - psfcpa_hv(i) = psfcpa(i,j) - qcg_hv(i) = qcg(i,j) - snowh_hv(i) = snowh(i,j) - dx_hv(i) = dx(i,j) - enddo - - !inout arguments: - do i = its,ite - regime_hv(i) = regime(i,j) - hfx_hv(i) = hfx(i,j) - qfx_hv(i) = qfx(i,j) - lh_hv(i) = lh(i,j) - rmol_hv(i) = rmol(i,j) - qgh_hv(i) = qgh(i,j) - znt_hv(i) = znt(i,j) - zol_hv(i) = zol(i,j) - cpm_hv(i) = cpm(i,j) - chs2_hv(i) = chs2(i,j) - cqs2_hv(i) = cqs2(i,j) - cqs_hv(i) = cqs(i,j) - chs_hv(i) = chs(i,j) - ch_hv(i) = ch(i,j) - flhc_hv(i) = flhc(i,j) - flqc_hv(i) = flqc(i,j) - gz1oz0_hv(i) = gz1oz0(i,j) - wspd_hv(i) = wspd(i,j) - br_hv(i) = br(i,j) - psim_hv(i) = psim(i,j) - psih_hv(i) = psih(i,j) - enddo - - !output arguments: - do i = its,ite - u10_hv(i) = 0. - v10_hv(i) = 0. - th2_hv(i) = 0. - t2_hv(i) = 0. - q2_hv(i) = 0. - wstar_hv(i) = 0. - enddo - - !optional output arguments: - if(present(ck) .and. present(cka) .and. present(cd) .and. present(cda)) then - do i = its,ite - ck_hv(i) = 0. - cka_hv(i) = 0. - cd_hv(i) = 0. - cda_hv(i) = 0. - enddo - endif - if(present(ustm)) then - do i = its,ite - ustm_hv(i) = ustm(i,j) - enddo - endif - - call mynnsfclay( & - u1d = u1d , v1d = v1d , t1d = t1d , qv1d = qv1d , & - p1d = p1d , dz8w1d = dz8w1d , rho1d = rho1d , u1d2 = u1d2 , & - v1d2 = v1d2 , dz2w1d = dz2w1d , cp = cp , g = g , & - rovcp = rovcp , r = r , xlv = xlv , psfcpa = psfcpa_hv , & - cpm = cpm_hv , & - chs = chs_hv , chs2 = chs2_hv , cqs2 = cqs2_hv , cqs = cqs_hv , & - pblh = pblh_hv , rmol = rmol_hv , znt = znt_hv , ust = ust_hv , & - mavail = mavail_hv , zol = zol_hv , mol = mol_hv , regime = regime_hv , & - psim = psim_hv , psih = psih_hv , xland = xland_hv , hfx = hfx_hv , & - qfx = qfx_hv , tsk = tsk_hv , u10 = u10_hv , v10 = v10_hv , & - th2 = th2_hv , t2 = t2_hv , q2 = q2_hv , flhc = flhc_hv , & - flqc = flqc_hv , snowh = snowh_hv , qgh = qgh_hv , qsfc = qsfc_hv , & - lh = lh_hv , gz1oz0 = gz1oz0_hv , wspd = wspd_hv , br = br_hv , & - isfflx = isfflx , dx = dx_hv , svp1 = svp1 , svp2 = svp2 , & - svp3 = svp3 , svpt0 = svpt0 , ep1 = ep1 , ep2 = ep2 , & - karman = karman , ch = ch_hv , qcg = qcg_hv , itimestep = itimestep , & - wstar = wstar_hv , qstar = qstar_hv , ustm = ustm_hv , ck = ck_hv , & - cka = cka_hv , cd = cd_hv , cda = cda_hv , spp_pbl = f_spp , & - rstoch1d = rstoch1d , isftcflx = isftcflx , iz0tlnd = iz0tlnd , & - its = its , ite = ite , restart_or_cycle = restart_or_cycle , & - errmsg = errmsg , errflg = errflg & - ) - - !inout arguments: - do i = its,ite - regime(i,j) = regime_hv(i) - hfx(i,j) = hfx_hv(i) - qfx(i,j) = qfx_hv(i) - lh(i,j) = lh_hv(i) - mol(i,j) = mol_hv(i) - rmol(i,j) = rmol_hv(i) - qgh(i,j) = qgh_hv(i) - qsfc(i,j) = qsfc_hv(i) - znt(i,j) = znt_hv(i) - zol(i,j) = zol_hv(i) - ust(i,j) = ust_hv(i) - cpm(i,j) = cpm_hv(i) - chs2(i,j) = chs2_hv(i) - cqs2(i,j) = cqs2_hv(i) - cqs(i,j) = cqs_hv(i) - chs(i,j) = chs_hv(i) - ch(i,j) = ch_hv(i) - flhc(i,j) = flhc_hv(i) - flqc(i,j) = flqc_hv(i) - gz1oz0(i,j) = gz1oz0_hv(i) - wspd(i,j) = wspd_hv(i) - br(i,j) = br_hv(i) - psim(i,j) = psim_hv(i) - psih(i,j) = psih_hv(i) - enddo - - !output arguments: - do i = its,ite - u10(i,j) = u10_hv(i) - v10(i,j) = v10_hv(i) - th2(i,j) = th2_hv(i) - t2(i,j) = t2_hv(i) - q2(i,j) = q2_hv(i) - wstar(i,j) = wstar_hv(i) - qstar(i,j) = qstar_hv(i) - enddo - - !optional output arguments: - if(present(ck) .and. present(cka) .and. present(cd) .and. present(cda)) then - do i = its,ite - ck(i,j) = ck_hv(i) - cka(i,j) = cka_hv(i) - cd(i,j) = cd_hv(i) - cda(i,j) = cda_hv(i) - enddo - endif - if(present(ustm)) then - do i = its,ite - ustm(i,j) = ustm_hv(i) - enddo - endif - - enddo - - end subroutine mynnsfclay_driver - -!================================================================================================================= - end module module_mynnsfclay_driver -!================================================================================================================= diff --git a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_pre.F90 b/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_pre.F90 deleted file mode 100644 index be8eae81fb..0000000000 --- a/src/core_atmosphere/physics/physics_wrf/module_mynnsfclay_pre.F90 +++ /dev/null @@ -1,174 +0,0 @@ -!================================================================================================================= - module module_mynnsfclay_pre - use ccpp_kind_types,only: kind_phys - - implicit none - private - public:: mynnsfclay_pre_init, & - mynnsfclay_pre_finalize, & - mynnsfclay_pre - - - contains - -!================================================================================================================= -!>\section arg_table_sf_mynn_pre_init -!!\html\include sf_mynn_pre_init.html -!! - subroutine mynnsfclay_pre_init(errmsg,errflg) -!================================================================================================================= - -!--- output arguments: - character(len=*),intent(out):: & - errmsg ! output error message (-). - - integer,intent(out):: & - errflg ! output error flag (-). - -!----------------------------------------------------------------------------------------------------------------- - -!--- output error flag and message: - errflg = 0 - errmsg = " " - - end subroutine mynnsfclay_pre_init - -!================================================================================================================= -!>\section arg_table_sf_mynn_pre_finalize -!!\html\include sf_mynn_pre_finalize.html -!! - subroutine mynnsfclay_pre_finalize(errmsg,errflg) -!================================================================================================================= - -!--- output arguments: - character(len=*),intent(out):: & - errmsg ! output error message (-). - - integer,intent(out):: & - errflg ! output error flag (-). - -!----------------------------------------------------------------------------------------------------------------- - -!--- output error flag and message: - errflg = 0 - errmsg = " " - - end subroutine mynnsfclay_pre_finalize - -!================================================================================================================= -!>\section arg_table_sf_mynn_pre_run -!!\html\include sf_mynn_pre_run.html -!! - subroutine mynnsfclay_pre(its,ite,kte,initflag,dz3d,u3d,v3d,p3d,t3d,rho3d,qv3d,qc3d,f_spp,pattern_spp, & - ust,mol,qsfc,qstar,dz8w1d,u1d,v1d,p1d,t1d,rho1d,qv1d,qc1d,rstoch1d,dz2w1d,u1d2, & - v1d2,errmsg,errflg) -!================================================================================================================= - -!--- input arguments: - logical,intent(in):: f_spp - - integer,intent(in):: its,ite - integer,intent(in):: kte - integer,intent(in):: initflag - - real(kind=kind_phys),intent(in),dimension(its:ite,1:kte):: & - dz3d, &! - u3d, &! - v3d, &! - qv3d, &! - qc3d, &! - p3d, &! - t3d, &! - rho3d ! - - real(kind=kind_phys),intent(in),dimension(its:ite,1:kte):: & - pattern_spp ! - - -!--- inout arguments: - real(kind=kind_phys),intent(inout),dimension(its:ite):: & - ust, &! - mol, &! - qsfc, &! - qstar ! - - -!--- output arguments: - character(len=*),intent(out):: & - errmsg ! output error message (-). - - integer,intent(out):: & - errflg ! output error flag (-). - - real(kind=kind_phys),intent(out),dimension(its:ite):: & - dz8w1d, &! - u1d, &! - v1d, &! - qv1d, &! - qc1d, &! - p1d, &! - t1d, &! - rho1d, &! - rstoch1d ! - - real(kind=kind_phys),intent(out),dimension(its:ite):: & - dz2w1d, &! - u1d2, &! - v1d2 ! - - -!--- local variables: - integer:: i,kts - -!----------------------------------------------------------------------------------------------------------------- - - kts = 1 - - do i = its,ite - dz8w1d(i) = dz3d(i,kts) - u1d(i) = u3d(i,kts) - v1d(i) = v3d(i,kts) - qv1d(i) = qv3d(i,kts) - qc1d(i) = qc3d(i,kts) - p1d(i) = p3d(i,kts) - t1d(i) = t3d(i,kts) - rho1d(i) = rho3d(i,kts) - !--- 2nd model level winds - for diags with high-resolution grids: - dz2w1d(i) = dz3d(i,kts+1) - u1d2(i) = u3d(i,kts+1) - v1d2(i) = v3d(i,kts+1) - enddo - - if(f_spp) then - do i = its,ite - rstoch1d(i) = pattern_spp(i,kts) - enddo - else - do i = its,ite - rstoch1d(i)=0._kind_phys - enddo - endif - - if(initflag == 1) then - do i = its,ite - if (ust(i) .lt. 1e-4 .or. ust(i) .gt. 3.0) then - ust(i) = max(0.04*sqrt(u1d(i)*u1d(i) + v1d(i)*v1d(i)),0.001) - endif - !qfx(i) = 0._kind_phy - !hfx(i) = 0._kind_phy - mol(i) = 0._kind_phys - !qsfc(i) = qv1d(i)/(1.+qv1d(i)) - qsfc(i) = qsfc(i) !should be available - qstar(i) = 0._kind_phys - enddo - endif - -!--- output message and error flags: - errmsg = 'sf_mynn_mpas_run OK' - errflg = 0 - - end subroutine mynnsfclay_pre - -!================================================================================================================= - end module module_mynnsfclay_pre -!=================================================================================================================