From 7df983d896c36575cd6f1934104d63888febb112 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 24 Sep 2020 14:07:48 -0600 Subject: [PATCH] Coupling CTSM with WRF (#1275) TYPE: New Feature KEYWORDS: WRF-CTSM coupling, Community Terrestrial Systems Model (CTSM), Community Land Model (CLM), coupled land-atmosphere models SOURCE: Negin Sobhani 1 , Dave Lawrence 1 , Bill Sacks 1, Sam Levis 2, Mariana Vertenstein 1, Dave Gill 1, Mike Barlage 1,3, and Fei Chen 1 1 National Center for Atmospheric Research (NCAR) 2 SLevis Consulting 3 National Centers for Environmental Prediction (NCEP) DESCRIPTION OF CHANGES: This PR includes the capability to couple WRF with Community Terrestrial Systems Model (CTSM) via LILAC (Lightweight Infrastructure for Land-Atmosphere Coupling). CTSM can be used with WRF as one of WRF's land surface model schemes by setting `sf_surface_physics` in the WRF namelist.input file: ``` &physics sf_surface_physics = 6, 6, 6, / ``` For instructions on how to run WRF with CTSM please check [instructions on using CTSM with WRF.](https://escomp.github.io/ctsm-docs/versions/master/html/lilac/specific-atm-models/index.html) LIST OF MODIFIED FILES: M Registry/Registry.EM_COMMON M arch/Config.pl M arch/configure.defaults M arch/postamble M arch/preamble A doc/README.CTSM M dyn_em/module_first_rk_step_part1.F M dyn_em/module_initialize_real.F M dyn_em/solve_em.F M main/depend.common M phys/Makefile M phys/module_physics_init.F A phys/module_sf_ctsm.F M phys/module_surface_driver.F M share/input_wrf.F M share/module_check_a_mundo.F M wrftladj/solve_em_ad.F TESTS CONDUCTED: 1. The jenkins tests are all passing. 2. Initial assessment of the performance of WRF-CTSM5(NWP) compared to WRF-Noah and WRF-NoahMP has been completed over the CONUS domain. Overall, we are finding that WRF-CTSM5(NWP) is performing reasonably well, in some cases even showing improved performance improvement. For example, the following plot shows the bias for maximum daily temperature from NOAH, NOAH-MP, and CTSM for June 2013 over METAR observation sites: ![CTSM_PR_2013-04_Tmin_Tmax](https://user-images.githubusercontent.com/17344536/90686072-05acc780-e228-11ea-91fa-4d239e548d1b.png) RELEASE NOTE: The capability to couple CTSM with WRF is added and activated by setting namelist option `sf_surface_physics` to 6. This is the initial beta release of WRF-CTSM coupling capability. For instructions on how to run WRF with CTSM please check [instructions on using CTSM with WRF.](https://escomp.github.io/ctsm-docs/versions/master/html/lilac/specific-atm-models/index.html) --- Registry/Registry.EM_COMMON | 1 + arch/Config.pl | 34 ++ arch/configure.defaults | 80 +-- arch/postamble | 4 +- arch/preamble | 3 + doc/README.CTSM | 28 + dyn_em/module_first_rk_step_part1.F | 7 +- dyn_em/module_initialize_real.F | 5 +- dyn_em/solve_em.F | 10 + main/depend.common | 6 + phys/Makefile | 1 + phys/module_physics_init.F | 19 +- phys/module_sf_ctsm.F | 779 ++++++++++++++++++++++++++++ phys/module_surface_driver.F | 186 ++++++- share/input_wrf.F | 4 + share/module_check_a_mundo.F | 29 ++ wrftladj/solve_em_ad.F | 9 + 17 files changed, 1156 insertions(+), 49 deletions(-) create mode 100644 doc/README.CTSM create mode 100644 phys/module_sf_ctsm.F diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 869dff658c..a6c2d8f771 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -2909,6 +2909,7 @@ package lsmscheme sf_surface_physics==2 - state:flx4,fv package ruclsmscheme sf_surface_physics==3 - state:smfr3d,keepfr3dflag,soilt1,rhosnf,snowfallac,precipfr,acrunoff package noahmpscheme sf_surface_physics==4 - state:isnowxy,tvxy,tgxy,canliqxy,canicexy,eahxy,tahxy,cmxy,chxy,fwetxy,sneqvoxy,alboldxy,qsnowxy,wslakexy,zwtxy,waxy,wtxy,tsnoxy,zsnsoxy,snicexy,snliqxy,lfmassxy,rtmassxy,stmassxy,woodxy,stblcpxy,fastcpxy,xsaixy,taussxy,t2mvxy,t2mbxy,q2mvxy,q2mbxy,tradxy,neexy,gppxy,nppxy,fvegxy,qinxy,runsfxy,runsbxy,ecanxy,edirxy,etranxy,fsaxy,firaxy,aparxy,psnxy,savxy,sagxy,rssunxy,rsshaxy,bgapxy,wgapxy,tgvxy,tgbxy,chvxy,chbxy,shgxy,shcxy,shbxy,evgxy,evbxy,ghvxy,ghbxy,irgxy,ircxy,irbxy,trxy,evcxy,chleafxy,chucxy,chv2xy,chb2xy,chstarxy,smoiseq,smcwtdxy,rechxy,deeprechxy,fdepthxy,areaxy,rivercondxy,riverbedxy,eqzwt,pexpxy,qrfxy,qrfsxy,qspringxy,qspringsxy,qslatxy,stepwtd,rechclim,gddxy,grainxy,croptype,planting,harvest,season_gdd,cropcat,pgsxy,soilcomp,soilcl1,soilcl2,soilcl3,soilcl4 package clmscheme sf_surface_physics==5 - state:numc,nump,sabv,sabg,lwup,lhsoi,lhveg,lhtran,snl,snowdp,wtc,wtp,h2osno,t_grnd,t_veg,h2ocan,h2ocan_col,t2m_max,t2m_min,t2clm,t_ref2m,h2osoi_liq_s1,h2osoi_liq_s2,h2osoi_liq_s3,h2osoi_liq_s4,h2osoi_liq_s5,h2osoi_liq1,h2osoi_liq2,h2osoi_liq3,h2osoi_liq4,h2osoi_liq5,h2osoi_liq6,h2osoi_liq7,h2osoi_liq8,h2osoi_liq9,h2osoi_liq10,h2osoi_ice_s1,h2osoi_ice_s2,h2osoi_ice_s3,h2osoi_ice_s4,h2osoi_ice_s5,h2osoi_ice1,h2osoi_ice2,h2osoi_ice3,h2osoi_ice4,h2osoi_ice5,h2osoi_ice6,h2osoi_ice7,h2osoi_ice8,h2osoi_ice9,h2osoi_ice10,t_soisno_s1,t_soisno_s2,t_soisno_s3,t_soisno_s4,t_soisno_s5,t_soisno1,t_soisno2,t_soisno3,t_soisno4,t_soisno5,t_soisno6,t_soisno7,t_soisno8,t_soisno9,t_soisno10,dzsnow1,dzsnow2,dzsnow3,dzsnow4,dzsnow5,snowrds1,snowrds2,snowrds3,snowrds4,snowrds5,t_lake1,t_lake2,t_lake3,t_lake4,t_lake5,t_lake6,t_lake7,t_lake8,t_lake9,t_lake10,h2osoi_vol1,h2osoi_vol2,h2osoi_vol3,h2osoi_vol4,h2osoi_vol5,h2osoi_vol6,h2osoi_vol7,h2osoi_vol8,h2osoi_vol9,h2osoi_vol10,albedosubgrid,lhsubgrid,hfxsubgrid,lwupsubgrid,q2subgrid,sabvsubgrid,sabgsubgrid,nrasubgrid,swupsubgrid +package ctsmscheme sf_surface_physics==6 - - package pxlsmscheme sf_surface_physics==7 - state:t2_ndg_new,q2_ndg_new,t2_ndg_old,q2_ndg_old,t2obs,q2obs,vegf_px,imperv,canfra,lai_px,wwlt_px,wfc_px,wsat_px,clay_px,csand_px,fmsand_px package ssibscheme sf_surface_physics==8 - state:ssib_fm,ssib_fh,ssib_cm,ssibxdd,ssib_br,ssib_lhf,ssib_shf,ssib_ghf,ssib_egs,ssib_eci,ssib_ect,ssib_egi,ssib_egt,ssib_sdn,ssib_sup,ssib_ldn,ssib_lup,ssib_wat,ssib_shc,ssib_shg,ssib_lai,ssib_vcf,ssib_z00,ssib_veg,isnow,swe,snowden,snowdepth,tkair,dzo1,wo1,tssn1,tssno1,bwo1,bto1,cto1,fio1,flo1,bio1,blo1,ho1,dzo2,wo2,tssn2,tssno2,bwo2,bto2,cto2,fio2,flo2,bio2,blo2,ho2,dzo3,wo3,tssn3,tssno3,bwo3,bto3,cto3,fio3,flo3,bio3,blo3,ho3,dzo4,wo4,tssn4,tssno4,bwo4,bto4,cto4,fio4,flo4,bio4,blo4,ho4 diff --git a/arch/Config.pl b/arch/Config.pl index b6f45dbe97..ecccb670cd 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -15,6 +15,7 @@ $sw_jasperinc_path=""; $sw_esmflib_path=""; $sw_esmfinc_path=""; +$sw_ctsm_mkfile_path=""; $sw_ldflags=""; $sw_compileflags=""; $sw_opt_level=""; @@ -321,6 +322,20 @@ $sw_esmf_ldflag = "yes" ; } +# A separately-installed CTSM library and its dependencies are required +# to build WRF with the CTSM land surface model (the next generation of +# the CLM model, coupled via LILAC). The user must set the environment +# variable WRF_CTSM_MKFILE, which should be a path to a make-formatted +# file containing settings of CTSM_INCLUDES and CTSM_LIBS. When this +# environment variable is set, this also triggers adding -DWRF_USE_CTSM; +# when this env var is not set, then we instead use -DWRF_USE_CLM, which +# builds an old version of CLM that is included in the WRF source code. +# (We currently cannot build with both at once because of namespace +# collisions at link time.) +if ( $ENV{WRF_CTSM_MKFILE} ) { + $sw_ctsm_mkfile_path = $ENV{WRF_CTSM_MKFILE}; +} + # parse the configure.wrf file $validresponse = 0 ; @@ -584,6 +599,19 @@ $_ =~ s:ESMFIOEXTLIB:-L\$\(WRF_SRC_ROOT_DIR\)/external/esmf_time_f90 -lesmf_time:g ; } } + + # CTSM substitutions in configure.defaults and postamble + if ( $sw_ctsm_mkfile_path ) { + $_ =~ s:CONFIGURE_D_CTSM:-DWRF_USE_CTSM:g ; + $_ =~ s:CONFIGURE_CTSM_INC:\$\(CTSM_INCLUDES\):g ; + $_ =~ s:CONFIGURE_CTSM_LIB:\$\(CTSM_LIBS\):g ; + } + else { + $_ =~ s:CONFIGURE_D_CTSM:-DWRF_USE_CLM:g ; + $_ =~ s:CONFIGURE_CTSM_INC::g ; + $_ =~ s:CONFIGURE_CTSM_LIB::g ; + } + if ( $ENV{HWRF} ) { $_ =~ s:CONFIGURE_ATMOCN_LIB:-L\$\(WRF_SRC_ROOT_DIR\)/external/atm_ocn -latm_ocn:g ; @@ -779,6 +807,12 @@ $_ =~ s:ESMFIODEFS::g ; $_ =~ s:ESMFTARGET:esmf_time:g ; } + + # CTSM substitutions in preamble + if ( $sw_ctsm_mkfile_path ) { + $_ =~ s:\# CTSMINCLUDEGOESHERE:include $sw_ctsm_mkfile_path: ; + } + if ( $ENV{HWRF} ) { $_ =~ s:CONFIGURE_ATMOCN_LIB:-L\$\(WRF_SRC_ROOT_DIR\)/external/atm_ocn -latm_ocn:g ; diff --git a/arch/configure.defaults b/arch/configure.defaults index 90983753e8..ab58009d14 100644 --- a/arch/configure.defaults +++ b/arch/configure.defaults @@ -16,7 +16,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNEC -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNEC -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -c # -DRSL0_ONLY #-DNCARIBM_NOC99 -Xa -Kc99 LDFLAGS_LOCAL = -Wl,-h nodefs @@ -61,7 +61,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -104,7 +104,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -DF2CSTYLE # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -147,7 +147,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -190,7 +190,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = -L$(MPI_ROOT)/lib -lmpi CPLUSPLUSLIB = @@ -233,7 +233,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -D_ACCEL -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -D_ACCEL CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -308,7 +308,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars # -DRSL0_ONLY LDFLAGS_LOCAL = -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common CPLUSPLUSLIB = @@ -353,7 +353,7 @@ CC = $(DM_CC) LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DCHUNK=16 -DXEON_OPTIMIZED_WSM5 -DXEON_SIMD -DOPTIMIZE_CFL_TEST -DFSEEKO64_OK -DINTEL_YSU_KLUDGE -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DCHUNK=16 -DXEON_OPTIMIZED_WSM5 -DXEON_SIMD -DOPTIMIZE_CFL_TEST -DFSEEKO64_OK -DINTEL_YSU_KLUDGE CONFIGURE_D_CTSM OPTNOSIMD = OPTKNC = -fimf-precision=low -fimf-domain-exclusion=15 -opt-assume-safe-padding -opt-streaming-stores always -opt-streaming-cache-evict=0 -mP2OPT_hlo_pref_use_outer_strategy=F CFLAGS_LOCAL = -w -O3 $(OPTKNC) # -DRSL0_ONLY @@ -400,7 +400,7 @@ CC = $(DM_CC) LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DCHUNK=64 -DXEON_OPTIMIZED_WSM5 -DOPTIMIZE_CFL_TEST -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DCHUNK=64 -DXEON_OPTIMIZED_WSM5 -DOPTIMIZE_CFL_TEST CONFIGURE_D_CTSM OPTNOSIMD = OPTAVX = -xAVX CFLAGS_LOCAL = -w -O3 $(OPTAVX) # -DRSL0_ONLY @@ -473,7 +473,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars # -DRSL0_ONLY LDFLAGS_LOCAL = -ip -lmpi #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common CPLUSPLUSLIB = @@ -522,7 +522,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars # -DRSL0_ONLY LDFLAGS_LOCAL = -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common CPLUSPLUSLIB = @@ -601,7 +601,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars # -DRSL0_ONLY LDFLAGS_LOCAL = -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common CPLUSPLUSLIB = @@ -684,7 +684,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars # -DRSL0_ONLY LDFLAGS_LOCAL = -ip #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common CPLUSPLUSLIB = @@ -730,7 +730,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM -D__PATHSCALE__ +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM -D__PATHSCALE__ CFLAGS_LOCAL = # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -773,7 +773,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -816,7 +816,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -DMACOS # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -859,7 +859,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip -DMACOS #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars -DMACOS # -DRSL0_ONLY # increase stack size; also note that for OpenMP, set environment OMP_STACKSIZE 4G or greater LDFLAGS_LOCAL = -ip -Wl,-stack_addr,0xF10000000 -Wl,-stack_size,0x64000000 #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common @@ -905,7 +905,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -DMACOS #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars -DMACOS # -DRSL0_ONLY # increase stack size; also note that for OpenMP, set environment OMP_STACKSIZE 4G or greater LDFLAGS_LOCAL = -Wl,-stack_addr,0xF10000000 -Wl,-stack_size,0x64000000 #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common @@ -950,7 +950,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DG95 -DMACOS -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR -DRCONFIG_CHARLEN=64 -DWRF_USE_CLM +ARCH_LOCAL = -DG95 -DMACOS -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR -DRCONFIG_CHARLEN=64 CONFIGURE_D_CTSM CFLAGS_LOCAL = -DMACOS -DF2CSTYLE # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -994,7 +994,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c -DMACOS # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1037,7 +1037,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c -DMACOS # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1080,7 +1080,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 -ARCH_LOCAL = -DMAC_KLUDGE -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DMAC_KLUDGE -DF2CSTYLE -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -DMACOS -DF2CSTYLE # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1126,7 +1126,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DNATIVE_MASSV -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DNATIVE_MASSV CONFIGURE_D_CTSM CFLAGS_LOCAL = -DNOUNDERSCORE # -DRSL0_ONLY LDFLAGS_LOCAL = -lmass -lmassv -bnoquiet # print diagnostic messages CPLUSPLUSLIB = -lC @@ -1174,7 +1174,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DNATIVE_MASSV -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DNATIVE_MASSV CONFIGURE_D_CTSM CFLAGS_LOCAL = -DNOUNDERSCORE # -DRSL0_ONLY LDFLAGS_LOCAL = -lmass_64 -lmassvp7_64 -q64 -bnoquiet # linking diagnostics CPLUSPLUSLIB = -lC @@ -1290,7 +1290,7 @@ CC = $(DM_CC) LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -s integer32 -s real`expr 8 \* $(RWORDSIZE)` -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = # uncomment this for wrfda build @@ -1337,7 +1337,7 @@ CC = $(DM_CC) LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM OPTNOSIMD = # set this to override Cray 'craype' module setting #OPTAVX = -xAVX @@ -1390,7 +1390,7 @@ LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 # If system has even more processors, set VERY_LARGE_MAXPROC to that number -ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 -DBLUEGENE -DWRF_USE_CLM +ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 -DBLUEGENE CONFIGURE_D_CTSM CFLAGS_LOCAL = -DNOUNDERSCORE -DNCARIBM_NOC99 $(MPI_INC) # -DRSL0_ONLY LIB_LOCAL = $(MPI_LIB) LDFLAGS_LOCAL = -Wl,--allow-multiple-definition -qstatic @@ -1439,7 +1439,7 @@ LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 # If system has even more processors, set VERY_LARGE_MAXPROC to that number -ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 -DBLUEGENE -DWRF_USE_CLM +ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 -DBLUEGENE CONFIGURE_D_CTSM CFLAGS_LOCAL = -DNOUNDERSCORE # -DRSL0_ONLY LIB_LOCAL = LDFLAGS_LOCAL = -Wl,--allow-multiple-definition,--relax -qstatic @@ -1485,7 +1485,7 @@ LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -qrealsize=$(RWORDSIZE) -qintsize=4 # If system has even more processors, set VERY_LARGE_MAXPROC to that number -ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 -DWRF_USE_CLM +ARCH_LOCAL = -DMOVE_NL_OUTSIDE_MODULE_CONFIGURE -DNONSTANDARD_SYSTEM_SUBR -DVERY_LARGE_MAXPROC=36768 CONFIGURE_D_CTSM CFLAGS_LOCAL = -DNOUNDERSCORE # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = -lC @@ -1528,7 +1528,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1571,7 +1571,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -D_WIN32 -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -D_WIN32 CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -DMEMCPY_FOR_BCOPY # -DRSL0_ONLY LDFLAGS_LOCAL = Ws2_32.lib # -lnetcdff CPLUSPLUSLIB = @@ -1628,7 +1628,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1671,7 +1671,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -DMACOS # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1714,7 +1714,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DMACOS -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip -DMACOS #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars -DMACOS # -DRSL0_ONLY # increase stack size; also note that for OpenMP, set environment OMP_STACKSIZE 4G or greater LDFLAGS_LOCAL = -ip -Wl,-stack_addr,0xF10000000 -Wl,-stack_size,0x64000000 #-xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common @@ -1760,7 +1760,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DMACOS CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c -DMACOS # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1803,7 +1803,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -r$(RWORDSIZE) -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1846,7 +1846,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip -xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars -xCORE-AVX2 # -DRSL0_ONLY LDFLAGS_LOCAL = -ip -xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common -xCORE-AVX2 CPLUSPLUSLIB = @@ -1890,7 +1890,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -real-size `expr 8 \* $(RWORDSIZE)` -i4 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_FUNC CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -ip -xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -no-multibyte-chars -xMIC-AVX512 # -DRSL0_ONLY LDFLAGS_LOCAL = -ip -xHost -fp-model fast=2 -no-prec-div -no-prec-sqrt -ftz -align all -fno-alias -fno-common -xMIC-AVX512 CPLUSPLUSLIB = @@ -1934,7 +1934,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = #-fdefault-real-8 -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -w -O3 -c LDFLAGS_LOCAL = CPLUSPLUSLIB = @@ -1994,7 +1994,7 @@ CC = CONFIGURE_CC LD = $(FC) RWORDSIZE = CONFIGURE_RWORDSIZE PROMOTION = -CcdRR$(RWORDSIZE) -ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR -DWRF_USE_CLM +ARCH_LOCAL = -DNONSTANDARD_SYSTEM_SUBR CONFIGURE_D_CTSM CFLAGS_LOCAL = -Kfast -Xg -DSUN # -DRSL0_ONLY LDFLAGS_LOCAL = CPLUSPLUSLIB = diff --git a/arch/postamble b/arch/postamble index c8552595e2..152010fe3f 100644 --- a/arch/postamble +++ b/arch/postamble @@ -53,11 +53,11 @@ INCLUDE_MODULES = $(MODULE_SRCH_FLAG) \ -I$(WRF_SRC_ROOT_DIR)/wrftladj \ -I$(WRF_SRC_ROOT_DIR)/chem -I$(WRF_SRC_ROOT_DIR)/inc \ -I$(NETCDFPATH)/include \ - CONFIGURE_RTTOV_INC + CONFIGURE_RTTOV_INC CONFIGURE_CTSM_INC REGISTRY = Registry CC_TOOLS_CFLAGS = CONFIGURE_NMM_CORE -LIB = $(LIB_BUNDLED) $(LIB_EXTERNAL) $(LIB_LOCAL) $(LIB_WRF_HYDRO) +LIB = $(LIB_BUNDLED) $(LIB_EXTERNAL) $(LIB_LOCAL) $(LIB_WRF_HYDRO) CONFIGURE_CTSM_LIB LDFLAGS = $(OMP) $(FCFLAGS) $(LDFLAGS_LOCAL) CONFIGURE_LDFLAGS ENVCOMPDEFS = CONFIGURE_COMPILEFLAGS CPPFLAGS = $(ARCHFLAGS) $(ENVCOMPDEFS) -I$(LIBINCLUDE) $(TRADFLAG) CONFIGURE_COMMS_INCLUDE diff --git a/arch/preamble b/arch/preamble index 38a0773c74..ecefa1455b 100644 --- a/arch/preamble +++ b/arch/preamble @@ -100,6 +100,9 @@ NETCDF4_DEP_LIB = $(DEP_LIB_PATH) $(HDF5) $(ZLIB) $(GPFS) $(CURL) # NETCDF4INCLUDEGOESHERE +#### CTSM pieces + +# CTSMINCLUDEGOESHERE ############################################################################## diff --git a/doc/README.CTSM b/doc/README.CTSM new file mode 100644 index 0000000000..5608519ca9 --- /dev/null +++ b/doc/README.CTSM @@ -0,0 +1,28 @@ +===================== + Using CTSM with WRF +===================== + +CTSM has now been coupled to WRF. +CTSM can be used as the land surface model scheme in WRF by choosing sf_surface_physics = 6 +in the WRF namelist.input file. + +For comprehensive instructions on how to run WRF with CTSM please check: +https://escomp.github.io/ctsm-docs/versions/master/html/lilac/specific-atm-models/wrf.html + +In summary, the procedure for running WRF-CTSM is similar to the workflow for +running WRF real cases, except that it requires additional steps to: + +1. Clone the CTSM repository, +2. Build CTSM libraries and LILAC, +3. Set environments to build WRF with CTSM and LILAC, +4. Define namelist options required for CTSM, +5. Choosing sf_surface_physics = 6 in the WRF namelist.input. + +Currently nesting is not available for WRF-CTSM simulations. + +Please note that if you build WRF with CTSM, the previous WRF CLM option (sf_surface_physics=5) would not be available. + +** +For further assistance please contact Dr. Negin Sobhani (negins@ucar.edu) or Dr. David Lawrence (dlawren@ucar.edu). +** + diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index bcd346f49b..51d77ec230 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -33,6 +33,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & , k_start , k_end & , f_flux & , aerocu & + , restart_flag & ) USE module_state_description USE module_model_constants @@ -116,6 +117,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux + LOGICAL, INTENT(IN) :: restart_flag ! Local real :: HYDRO_dt REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA @@ -578,12 +580,13 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,tkdry3d=grid%tkdry3d, tksatu3d=grid%tksatu3d & & ,LakeModel=grid%sf_lake_physics, lake_min_elev=grid%lake_min_elev & !lake #if ( EM_CORE == 1) - & ,LakeMask=grid%LakeMask & ! lake + & ,LakeMask=grid%LakeMask & !lake + & ,restart_flag=restart_flag & !flag showing if is a restart timestep #endif ! CLM Varaibles & ,NUMC=grid%numc,NUMP=grid%nump,SABV=grid%sabv,SABG=grid%sabg, & & LWUP=grid%lwup,SNL=grid%snl, & - & HISTORY_INTERVAL=config_flags%history_interval , &!ylu add hist inverval for accumulation T max/min + & HISTORY_INTERVAL=config_flags%history_interval , & !ylu add hist inverval for accumulation T max/min & SNOWDP=grid%snowdp, WTC=grid%wtc,WTP=grid%wtp, H2OSNO=grid%h2osno, & & T_GRND=grid%t_grnd,T_VEG=grid%t_veg, & & H2OCAN=grid%h2ocan, H2OCAN_COL=grid%h2ocan_col,T2M_MAX=grid%t2m_max, & diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index f2d6906f07..00a559912d 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -3081,7 +3081,7 @@ SUBROUTINE init_domain_rk ( grid & fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) - CASE ( SLABSCHEME , LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME,CLMSCHEME, SSIBSCHEME ) + CASE ( SLABSCHEME , LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME,CLMSCHEME, CTSMSCHEME, SSIBSCHEME ) DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE @@ -3359,6 +3359,7 @@ SUBROUTINE init_domain_rk ( grid & ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ).AND. & ( model_config_rec%sf_surface_physics(grid%id) .NE. SSIBSCHEME ).AND. & !fds ( model_config_rec%sf_surface_physics(grid%id) .NE. CLMSCHEME ).AND. & + ( model_config_rec%sf_surface_physics(grid%id) .NE. CTSMSCHEME ).AND. & ( model_config_rec%sf_surface_physics(grid%id) .NE. PXLSMSCHEME ) ) THEN print *,'error in the grid%tslb' print *,'i,j=',i,j @@ -3380,7 +3381,7 @@ SUBROUTINE init_domain_rk ( grid & grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + & grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0) END DO - CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME,CLMSCHEME,SSIBSCHEME ) + CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME,CLMSCHEME, CTSMSCHEME, SSIBSCHEME ) ! CALL wrf_error_fatal ( 'Assigned constant soil moisture to 0.3, stopping') DO ns = 1 , model_config_rec%num_soil_layers grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + & diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 0fe7edebdd..b9937538a6 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -128,6 +128,7 @@ SUBROUTINE solve_em ( grid , config_flags & ! Flag for producing diagnostic fields (e.g., radar reflectivity) LOGICAL :: diag_flag INTEGER :: ke_diag ! tells reflectivity calculation whether to do full depth or only k=1 + LOGICAL :: restart_flag ! tells if it is a restart timestep to write restart files. #if (WRF_CHEM == 1) ! Index cross-referencing array for tendency accumulation @@ -277,6 +278,7 @@ SUBROUTINE solve_em ( grid , config_flags & curr_secs2 = real_time(tmpTimeInterval2) old_dt = grid%dt ! store old time step for flux averaging code at end of RK loop + !----------------------------------------------------------------------------- ! Adaptive time step: Added by T. Hutchinson, WSI 3/5/07 ! In this call, we do the time-step adaptation and set time-dependent lateral @@ -293,6 +295,13 @@ SUBROUTINE solve_em ( grid , config_flags & ! End of adaptive time step modifications !----------------------------------------------------------------------------- ! +! Set restart flag value history output time +!----------------------------------------------------------------------------- + restart_flag = .false. + if ( Is_alarm_tstep(grid%domain_clock, grid%alarms(restart_alarm)) ) then + restart_flag = .true. + endif +! ! Set diagnostic flag value history output time !----------------------------------------------------------------------------- @@ -756,6 +765,7 @@ SUBROUTINE solve_em ( grid , config_flags & , k_start , k_end & , f_flux & , aerocu & + , restart_flag & ) #ifdef DM_PARALLEL diff --git a/main/depend.common b/main/depend.common index d58517de04..e5a58206c4 100644 --- a/main/depend.common +++ b/main/depend.common @@ -449,6 +449,10 @@ module_sf_clm.o: module_cam_shr_kind_mod.o \ ../frame/module_wrf_error.o \ ../frame/module_configure.o +module_sf_ctsm.o: ../frame/module_dm.o \ + ../frame/module_configure.o \ + ../frame/module_wrf_error.o + module_sf_ssib.o: ../share/module_model_constants.o module_sf_noah_seaice_drv.o: module_sf_noah_seaice.o @@ -530,6 +534,7 @@ module_physics_init.o : \ module_sf_noahlsm.o \ module_sf_noahdrv.o \ module_sf_clm.o \ + module_sf_ctsm.o \ module_sf_ssib.o \ module_sf_noahmplsm.o \ module_sf_noahmpdrv.o \ @@ -734,6 +739,7 @@ module_surface_driver.o: \ module_sf_noahmp_groundwater.o \ module_sf_noahdrv.o \ module_sf_clm.o \ + module_sf_ctsm.o \ module_sf_ssib.o \ module_sf_noahmpdrv.o \ module_sf_ruclsm.o \ diff --git a/phys/Makefile b/phys/Makefile index 0833835900..ee46f37441 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -132,6 +132,7 @@ MODULES = \ module_sf_noahdrv.o \ module_sf_noahlsm.o \ module_sf_clm.o \ + module_sf_ctsm.o \ module_sf_ssib.o \ module_sf_noah_seaice.o \ module_sf_noah_seaice_drv.o \ diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index d9edbf6c53..cc37d864b8 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -1340,7 +1340,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & QKE, & SFCEVP,GRDFLX, & TRIM (MMINLU_LOC), & - XLAT, & + XLAT, XLONG, & ISNOWXY, ZSNSOXY, TSNOXY, & SNICEXY, SNLIQXY, TVXY, TGXY, CANICEXY, & CANLIQXY, EAHXY, TAHXY, CMXY, & @@ -2347,7 +2347,7 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & #endif QKE, SFCEVP,GRDFLX, & MMINLU, & - XLAT, & + XLAT, XLONG, & ISNOWXY, ZSNSOXY, TSNOXY, & SNICEXY, SNLIQXY, TVXY, TGXY, CANICEXY, & CANLIQXY, EAHXY, TAHXY, CMXY, & @@ -2497,6 +2497,9 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & USE noahmp_tables, ONLY: LOW_DENSITY_RESIDENTIAL_TABLE, HIGH_DENSITY_RESIDENTIAL_TABLE, HIGH_INTENSITY_INDUSTRIAL_TABLE #ifdef WRF_USE_CLM USE module_sf_clm, only : clminit +#endif +#ifdef WRF_USE_CTSM + USE module_sf_ctsm, only : ctsm_init #endif USE module_sf_urban USE module_sf_bep !BEP @@ -2625,6 +2628,7 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & !Noah-MP REAL, DIMENSION(ims:ime, jms:jme) :: XLAT + REAL, DIMENSION(ims:ime, jms:jme) :: XLONG INTEGER, OPTIONAL, DIMENSION(ims:ime,jms:jme) :: ISNOWXY REAL, OPTIONAL, DIMENSION(ims:ime,-2:num_soil_layers, jms:jme) :: ZSNSOXY REAL, OPTIONAL, DIMENSION(ims:ime,-2:0, jms:jme) :: TSNOXY @@ -3361,6 +3365,17 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & CALL wrf_error_fatal ( "SF CLM CURRENTLY ifdef'ed out, set -DWRF_USE_CLM in configure.wrf in ARCH_LOCAL" ) #endif +#ifdef WRF_USE_CTSM + CASE (CTSMSCHEME) + + CALL ctsm_init( & + ids=ids, ide=ide, jds=jds, jde=jde, & + ims=ims, ime=ime, jms=jms, jme=jme, & + its=its, ite=ite, jts=jts, jte=jte, & + dt=DT, xlat=XLAT, xlong=XLONG, & + atm_restart = restart) +#endif + CASE DEFAULT END SELECT sfc_select diff --git a/phys/module_sf_ctsm.F b/phys/module_sf_ctsm.F new file mode 100644 index 0000000000..590122ef27 --- /dev/null +++ b/phys/module_sf_ctsm.F @@ -0,0 +1,779 @@ +module module_sf_ctsm + + implicit none + private + +#ifdef WRF_USE_CTSM + + public :: ctsm_init + public :: ctsm_run + + ! FIXME(wjs, 2020-01-01) Introduce a lilac_kinds module, and get this from there + integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real + +contains + + !----------------------------------------------------------------------- + subroutine get_num_points(ide, jde, its, ite, jts, jte, & + num_points, ite_limited, jte_limited) + ! Return the number of points owned by this task + ! + ! If ite_limited and/or jte_limited are provided, then return those values, too. + ! These are the task end indices on the mass point grid. + integer, intent(in) :: ide ! domain end index, i + integer, intent(in) :: jde ! domain end index, j + integer, intent(in) :: its ! task start index, i + integer, intent(in) :: ite ! task end index, i + integer, intent(in) :: jts ! task start index, j + integer, intent(in) :: jte ! task end index, j + integer, intent(out) :: num_points ! number of points owned by this task + integer, optional, intent(out) :: ite_limited ! task end index on the mass point grid, i + integer, optional, intent(out) :: jte_limited ! task end index on the mass point grid, j + + integer :: my_ite_limited ! task end index on the mass point grid, i + integer :: my_jte_limited ! task end index on the mass point grid, j + + ! The very last index in both row & column space is just used on the momentum grid. + ! Here we are just working with the mass point grid, so we need to ignore that last + ! index. + my_ite_limited = min(ite, ide-1) + my_jte_limited = min(jte, jde-1) + + num_points = ((my_ite_limited - its + 1) * (my_jte_limited - jts + 1)) + + if (present(ite_limited)) then + ite_limited = my_ite_limited + end if + if (present(jte_limited)) then + jte_limited = my_jte_limited + end if + + end subroutine get_num_points + + subroutine create_gindex(ide, jde, its, ite, jts, jte, gindex) + ! Create a gindex array on each task. This gives the list of global indices owned by + ! each processor, on the mass point grid. + integer, intent(in) :: ide ! domain end index, i + integer, intent(in) :: jde ! domain end index, j + integer, intent(in) :: its ! task start index, i + integer, intent(in) :: ite ! task end index, i + integer, intent(in) :: jts ! task start index, j + integer, intent(in) :: jte ! task end index, j + integer, allocatable, intent(out) :: gindex(:) + + integer :: ite_limited ! task end index on the mass point grid, i + integer :: jte_limited ! task end index on the mass point grid, j + integer :: num_points + integer :: i, j, n + + call get_num_points( & + ide=ide, jde=jde, & + its=its, ite=ite, & + jts=jts, jte=jte, & + num_points = num_points, & + ite_limited = ite_limited, & + jte_limited = jte_limited) + + allocate(gindex(num_points)) + + n = 0 + do j = jts, jte_limited + do i = its, ite_limited + n = n + 1 + ! In the following, note that we use ide-1 rather than ide for the same reason + ! that we need ite_limited: ide gives the domain end index on the momentum grid, + ! but here we're just dealing with the mass point grid, which has one less point + ! in each direction. + gindex(n) = (j-1)*(ide-1) + i + end do + end do + end subroutine create_gindex + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + logical function is_land(xland, xice, xice_threshold) + ! Returns true if the given point is land, false if ocean or sea ice + real , intent(in) :: xland ! land mask (1 for land, 2 for water) + real , intent(in) :: xice ! fraction of grid that is seaice + real , intent(in) :: xice_threshold ! fraction of grid determining seaice + + if (xland >= 1.5) then + ! ocean + is_land = .false. + else if (xice >= xice_threshold) then + ! sea ice + is_land = .false. + else + is_land = .true. + end if + end function is_land + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + subroutine convert_2d_to_1d (ide, jde, ims, ime, jms, jme, its, ite, jts, jte, var_2d, var_1d) + ! Convert a 2-d native WRF array to a 1-d array appropriate for LILAC + ! + ! Allocates var_1d here + + use module_wrf_error + + ! input/output variables + integer , intent(in) :: ide ! domain end index, i + integer , intent(in) :: jde ! domain end index, j + integer , intent(in) :: ims ! memory start index (includes halo cells), i + integer , intent(in) :: ime ! memory end index (includes halo cells), i + integer , intent(in) :: jms ! memory start index (includes halo cells), j + integer , intent(in) :: jme ! memory end index (includes halo cells), j + integer , intent(in) :: its ! task start index, i + integer , intent(in) :: ite ! task end index, i + integer , intent(in) :: jts ! task start index, j + integer , intent(in) :: jte ! task end index, j + real , intent(in) :: var_2d(ims: , jms: ) + real(r8), allocatable, intent(out) :: var_1d(:) + + ! local variables + integer :: ite_limited ! task end index on the mass point grid, i + integer :: jte_limited ! task end index on the mass point grid, j + integer :: num_points + integer :: i, j, n + + if (.not. all(ubound(var_2d) == [ime, jme])) then + call wrf_error_fatal('convert_2d_to_1d: incorrect bounds for var_2d') + end if + + call get_num_points( & + ide=ide, jde=jde, & + its=its, ite=ite, & + jts=jts, jte=jte, & + num_points = num_points, & + ite_limited = ite_limited, & + jte_limited = jte_limited) + + allocate (var_1d (num_points)) + + n = 0 + + do j = jts, jte_limited + do i = its, ite_limited + n = n + 1 + var_1d(n) = var_2d(i,j) + end do + end do + + end subroutine convert_2d_to_1d + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + subroutine convert_1d_to_2d (ide, jde, ims, ime, jms, jme, its, ite, jts, jte, & + xland, xice, xice_threshold, & + var_1d, var_2d) + ! Convert a 1-d array from lilac to a 2-d native WRF array + ! + ! The output array is only set for land points (as determined by xland, xice, xice_threshold) + + use module_wrf_error + + ! input/output variables + integer , intent(in) :: ide ! domain end index, i + integer , intent(in) :: jde ! domain end index, j + integer , intent(in) :: ims ! memory start index (includes halo cells), i + integer , intent(in) :: ime ! memory end index (includes halo cells), i + integer , intent(in) :: jms ! memory start index (includes halo cells), j + integer , intent(in) :: jme ! memory end index (includes halo cells), j + integer , intent(in) :: its ! task start index, i + integer , intent(in) :: ite ! task end index, i + integer , intent(in) :: jts ! task start index, j + integer , intent(in) :: jte ! task end index, j + real , intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water) + real , intent(in) :: xice(ims: , jms: ) ! fraction of grid that is seaice + real , intent(in) :: xice_threshold ! fraction of grid determining seaice + real(r8) , intent(in) :: var_1d(:) + real , intent(inout) :: var_2d(ims: , jms: ) + + ! local variables + integer :: ite_limited ! task end index on the mass point grid, i + integer :: jte_limited ! task end index on the mass point grid, j + integer :: num_points + integer :: i, j, n + + if (.not. all(ubound(var_2d) == [ime, jme])) then + call wrf_error_fatal('convert_1d_to_2d: incorrect bounds for var_2d') + end if + + call get_num_points( & + ide=ide, jde=jde, & + its=its, ite=ite, & + jts=jts, jte=jte, & + num_points = num_points, & + ite_limited = ite_limited, & + jte_limited = jte_limited) + + if (.not. (ubound(var_1d, 1) == num_points)) then + call wrf_error_fatal('convert_1d_to_2d: incorrect size of var_1d') + end if + + n = 0 + + do j = jts, jte_limited + do i = its, ite_limited + n = n + 1 + if (is_land(xland(i,j), xice(i,j), xice_threshold)) then + var_2d(i,j) = var_1d(n) + end if + end do + end do + + end subroutine convert_1d_to_2d + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + subroutine ctsm_init( & + ids, ide, jds, jde, & + ims, ime, jms, jme, & + its, ite, jts, jte, & + dt, xlat, xlong, & + atm_restart) + ! Initialize CTSM via LILAC + + use module_wrf_error + use lilac_mod, only : lilac_init1, lilac_init2 + use ctsm_LilacCouplingFieldIndices + + integer, intent(in) :: ids + integer, intent(in) :: ide + integer, intent(in) :: jds + integer, intent(in) :: jde + integer, intent(in) :: ims + integer, intent(in) :: ime + integer, intent(in) :: jms + integer, intent(in) :: jme + integer, intent(in) :: its + integer, intent(in) :: ite + integer, intent(in) :: jts + integer, intent(in) :: jte + + real, intent(in) :: dt ! time step (seconds) + real, intent(in) :: xlat(ims: , jms: ) ! latitudes (deg) + real, intent(in) :: xlong(ims: , jms: ) ! longitudes (deg) + + logical, intent(in) :: atm_restart + integer :: comp_comm + integer , allocatable :: atm_global_index(:) + real :: xlong_0360(ims:ime, jms:jme) + real(r8) , allocatable :: atm_lons(:) + real(r8) , allocatable :: atm_lats(:) + + character(len=128) :: atm_calendar + character(len=128) :: atm_starttype + integer :: atm_timestep + integer :: atm_start_year ! (yyyy) + integer :: atm_start_mon ! (mm) + integer :: atm_start_day + integer :: atm_start_hour + integer :: atm_start_minute + integer :: atm_start_second ! seconds after the minute + integer :: atm_start_secs_since_midnight ! total seconds since midnight + + integer :: atm_global_nx + integer :: atm_global_ny + + character(len=512) :: message + + ! TODO(wjs, 2019-12-31) Is this the correct way to get the communicator? + ! (See https://github.com/ESCOMP/CTSM/issues/1078) + call wrf_get_dm_communicator(comp_comm) + + call create_gindex(ide, jde, its, ite, jts, jte, atm_global_index) + + ! Convert longitude from -180..180 to 0..360 + xlong_0360 = xlong + where (xlong_0360 < 0) + xlong_0360 = 360. + xlong_0360 + end where + + ! reshape lats and lons to 1d for lilac + call convert_2d_to_1d ( & + ide=ide, jde=jde, & + ims=ims, ime=ime, & + jms=jms, jme=jme, & + its=its, ite=ite, & + jts=jts, jte=jte, & + var_2d=xlong_0360, var_1d=atm_lons) + call convert_2d_to_1d ( & + ide=ide, jde=jde, & + ims=ims, ime=ime, & + jms=jms, jme=jme, & + its=its, ite=ite, & + jts=jts, jte=jte, & + var_2d=xlat, var_1d=atm_lats) + + atm_global_nx = ide-ids + atm_global_ny = jde-jds + + ! calendar stuff + + ! TODO(wjs, 2019-12-31) Is this the appropriate way to get the start and end times? + ! + ! Some specific questions: + ! - Should we use some passed-in argument rather than getting the namelist values + ! directly here? + ! - nl_get_start_year vs. nl_get_simulation_start_year (and similar for other units): + ! should we be using nl_get_simulation_* here? + ! - Is it correct to form the total seconds from hour, minute & second like this? + ! - I'm using 1 for the first argument; I think this gives the domain id; is this an + ! okay thing to do for now, or can we somehow determine the actual domain id? + ! + ! (See https://github.com/ESCOMP/CTSM/issues/1078) + call nl_get_start_year(1, atm_start_year) + call nl_get_start_month(1, atm_start_mon) + call nl_get_start_day(1, atm_start_day) + call nl_get_start_hour(1, atm_start_hour) + call nl_get_start_minute(1, atm_start_minute) + call nl_get_start_second(1, atm_start_second) + atm_start_secs_since_midnight = 60*(60*atm_start_hour + atm_start_minute) + atm_start_second + write(message, '("CTSM start time: ", I4, 1X, I2, 1X, I2, 1X, I5)') & + atm_start_year, atm_start_mon, atm_start_day, atm_start_secs_since_midnight + call wrf_message(message) + + atm_calendar = 'GREGORIAN' + + atm_starttype = 'startup' + if (atm_restart) then + atm_starttype = 'continue' + endif + + + ! TODO(wjs, 2019-12-31) Is there a way to directly get dt as an integer, rather than + ! relying on converting this real-valued dt to an integer? + ! (See https://github.com/ESCOMP/CTSM/issues/1078) + atm_timestep = nint(dt) + if (abs(atm_timestep - dt) > (1.e-5 * dt)) then + call wrf_error_fatal('ctsm_init: expect dt representable as integer') + end if + + call lilac_init1() + call lilac_init2( & + mpicom = comp_comm, & + atm_global_index = atm_global_index, & + atm_lons = atm_lons, & + atm_lats = atm_lats, & + atm_global_nx = atm_global_nx, & + atm_global_ny = atm_global_ny, & + atm_calendar = atm_calendar, & + atm_timestep = atm_timestep, & + atm_start_year = atm_start_year, & + atm_start_mon = atm_start_mon, & + atm_start_day = atm_start_day, & + atm_start_secs = atm_start_secs_since_midnight, & + starttype_in = atm_starttype, & + fields_needed_from_data = [ & + lilac_a2l_Faxa_bcphidry, lilac_a2l_Faxa_bcphodry, lilac_a2l_Faxa_bcphiwet, & + lilac_a2l_Faxa_ocphidry, lilac_a2l_Faxa_ocphodry, lilac_a2l_Faxa_ocphiwet, & + lilac_a2l_Faxa_dstwet1, lilac_a2l_Faxa_dstdry1, & + lilac_a2l_Faxa_dstwet2, lilac_a2l_Faxa_dstdry2, & + lilac_a2l_Faxa_dstwet3, lilac_a2l_Faxa_dstdry3, & + lilac_a2l_Faxa_dstwet4, lilac_a2l_Faxa_dstdry4]) + + end subroutine ctsm_init + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + subroutine ctsm_run( & + ! bounds + ids, ide, jds, jde, & + ims, ime, jms, jme, kms, & + its, ite, jts, jte, & + + ! restart flag + restart_flag, & + + ! general information + dt, xland, xice, xice_threshold, & + + ! atm -> lnd variables + dz8w, ht, u_phy, v_phy, p8w, t_phy, th_phy, & + qv_curr, rainbl, sr, & + glw, swvisdir, swvisdif, swnirdir, swnirdif, & + + ! lnd -> atm variables + tsk, t2, qsfc, albedo, & + ust, hfx, lh, qfx, emiss, z0, znt) + + use lilac_mod, only : lilac_run + use ctsm_LilacCouplingFieldIndices + + integer, intent(in) :: ids + integer, intent(in) :: ide + integer, intent(in) :: jds + integer, intent(in) :: jde + integer, intent(in) :: ims + integer, intent(in) :: ime + integer, intent(in) :: jms + integer, intent(in) :: jme + integer, intent(in) :: kms + integer, intent(in) :: its + integer, intent(in) :: ite + integer, intent(in) :: jts + integer, intent(in) :: jte + + + logical, intent(in) :: restart_flag ! restart flag sent from WRF to CTSM + + real , intent(in) :: dt ! timestep [s] + real , intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water) + real , intent(in) :: xice(ims: , jms: ) ! fraction of grid that is seaice + real , intent(in) :: xice_threshold ! fraction of grid determining seaice + + real, intent(in) :: dz8w(ims: , kms: , jms: ) ! thickness of atmo layers [m] + real, intent(in) :: ht(ims: , jms: ) ! terrain height [m] + real, intent(in) :: u_phy(ims: , kms: , jms: ) ! 3D U wind component [m/s] + real, intent(in) :: v_phy(ims: , kms: , jms: ) ! 3D V wind component [m/s] + real, intent(in) :: p8w(ims: , kms: , jms: ) ! 3D pressure, valid at interface [Pa] + real, intent(in) :: t_phy(ims: , kms: , jms: ) ! 3D atmospheric temperature valid at mid-levels [K] + real, intent(in) :: th_phy(ims: , kms: , jms: ) ! 3D atmospheric temperature valid at mid-levels [K] + real, intent(in) :: qv_curr(ims: , kms: , jms: ) ! 3D water vapor mixing ratio [kg/kg_dry] + real, intent(in) :: rainbl(ims: , jms: ) ! total input precipitation [mm] + real, intent(in) :: sr(ims: , jms: ) ! frozen precipitation ratio [-] + real, intent(in) :: glw(ims: , jms: ) ! longwave down at surface [W m-2] + real, intent(in) :: swvisdir(ims: , jms: ) ! vis direct beam solar rad onto surface [W m-2] + real, intent(in) :: swvisdif(ims: , jms: ) ! vis diffuse solar rad onto surface [W m-2] + real, intent(in) :: swnirdir(ims: , jms: ) ! nir direct beam solar rad onto surface [W m-2] + real, intent(in) :: swnirdif(ims: , jms: ) ! nir diffuse solar rad onto surface [W m-2] + + real, intent(inout) :: tsk(ims: , jms: ) ! surface temperature [K] + real, intent(inout) :: t2(ims: , jms: ) ! diagnostic 2-m temperature [K] + real, intent(inout) :: qsfc(ims: , jms: ) ! bulk surface specific humidity + real, intent(inout) :: albedo(ims: , jms: ) ! total grid albedo [-] + real, intent(inout) :: ust(ims: , jms: ) ! u* in similarity theory [m/s] + real, intent(inout) :: hfx(ims: , jms: ) ! sensible heat flux [W m-2] + real, intent(inout) :: lh(ims: , jms: ) ! latent heat flux [W m-2] + real, intent(inout) :: qfx(ims: , jms: ) ! latent heat flux [kg s-1 m-2] + real, intent(inout) :: emiss(ims: , jms: ) ! surface emissivity [between 0 and 1] + real, intent(inout) :: z0(ims: , jms: ) ! background roughness length [m] + real, intent(inout) :: znt(ims: , jms: ) ! thermal time-varying roughness length [m] + + logical, save :: first_call = .true. ! true if and only if this is the first time this subroutine has been called in this run + + integer :: i, j + + real :: landmask(ims:ime, jms:jme) ! 1 over land, 0 over non-land (ocean, sea ice) + real :: zlvl(ims:ime, jms:jme) ! mid-point of bottom atm layer [m] + real :: forc_u(ims:ime, jms:jme) ! u wind component, bottom atm layer [m/s] + real :: forc_v(ims:ime, jms:jme) ! v wind component, bottom atm layer [m/s] + real :: pbot(ims:ime, jms:jme) ! surface pressure [Pa] + real :: forc_th(ims:ime, jms:jme) ! potential temperature [K] + real :: forc_t(ims:ime, jms:jme) ! temperature [K] + real :: forc_q(ims:ime, jms:jme) ! specific humidity [kg/kg] + real :: precip_rate(ims:ime, jms:jme) ! total precipitation rate [mm/s] + real :: rain_convective_rate(ims:ime, jms:jme) ! rate of convective rain [mm/s] + real :: rain_largescale_rate(ims:ime, jms:jme) ! rate of large-scale rain [mm/s] + real :: snow_convective_rate(ims:ime, jms:jme) ! rate of convective snow [mm/s] + real :: snow_largescale_rate(ims:ime, jms:jme) ! rate of large-scale snow [mm/s] + + real :: qref(ims:ime, jms:jme) ! 2m surface specific humidity calculated by CTSM [kg/kg] + real :: albedo_visdir(ims:ime, jms:jme) ! albedo calculated by CTSM, visible direct [-] + real :: albedo_visdif(ims:ime, jms:jme) ! albedo calculated by CTSM, visible diffuse [-] + real :: albedo_nirdir(ims:ime, jms:jme) ! albedo calculated by CTSM, near IR direct [-] + real :: albedo_nirdif(ims:ime, jms:jme) ! albedo calculated by CTSM, near IR diffuse [-] + + logical :: no_negative = .true. ! flag for converting negative values to zero + + ! ------------------------------------------------------------------------ + ! Calculate landmask and export it to ctsm via lilac. + ! + ! For efficiency, only do this the first time this subroutine is called: we assume + ! landmask stays constant throughout the run (there can be transitions between open + ! ocean and sea ice, but the land itself stays fixed). + ! ------------------------------------------------------------------------ + + if (first_call) then + do j = jts, jte + do i = its, ite + if (is_land(xland(i,j), xice(i,j), xice_threshold)) then + landmask(i,j) = 1. + else + landmask(i,j) = 0. + end if + end do + end do + call export_to_lilac(lilac_a2l_Sa_landfrac, landmask) + end if + + ! ------------------------------------------------------------------------ + ! Calculate derived variables, and 2-d versions of 3-d fields + ! ------------------------------------------------------------------------ + + ! dz8w = thickness of full levels; we want the mid-point of the bottom-most level + zlvl(its:ite, jts:jte) = 0.5 * dz8w(its:ite, 1, jts:jte) + + forc_u(its:ite, jts:jte) = u_phy(its:ite, 1, jts:jte) + forc_v(its:ite, jts:jte) = v_phy(its:ite, 1, jts:jte) + + pbot(its:ite, jts:jte) = p8w(its:ite, 1, jts:jte) + + forc_th(its:ite, jts:jte) = th_phy(its:ite, 1, jts:jte) + + forc_t(its:ite, jts:jte) = t_phy(its:ite, 1, jts:jte) + + ! Convert from mixing ratio to specific humidity + forc_q(its:ite, jts:jte) = qv_curr(its:ite, 1, jts:jte)/(1.0 + qv_curr(its:ite, 1, jts:jte)) + + ! Separate total precip into rain and snow. Arbitrarily assign all precipitation to + ! convective (CTSM requires separate convective vs. large-scale precipitation, but + ! then just adds the two together). + precip_rate(its:ite, jts:jte) = rainbl(its:ite, jts:jte)/dt + snow_convective_rate(its:ite, jts:jte) = precip_rate(its:ite, jts:jte) * sr(its:ite, jts:jte) + snow_largescale_rate(its:ite, jts:jte) = 0. + rain_convective_rate(its:ite, jts:jte) = precip_rate(its:ite, jts:jte) * (1. - sr(its:ite, jts:jte)) + rain_largescale_rate(its:ite, jts:jte) = 0. + + ! ------------------------------------------------------------------------ + ! Export data to ctsm via lilac + ! ------------------------------------------------------------------------ + + call export_to_lilac(lilac_a2l_Sa_z, zlvl) + call export_to_lilac(lilac_a2l_Sa_topo, ht) + call export_to_lilac(lilac_a2l_Sa_u, forc_u) + call export_to_lilac(lilac_a2l_Sa_v, forc_v) + call export_to_lilac(lilac_a2l_Sa_ptem, forc_th) + call export_to_lilac(lilac_a2l_Sa_pbot, pbot) + call export_to_lilac(lilac_a2l_Sa_tbot, forc_t) + call export_to_lilac(lilac_a2l_Sa_shum, forc_q) + call export_to_lilac(lilac_a2l_Faxa_lwdn, glw) + call export_to_lilac(lilac_a2l_Faxa_rainc, rain_convective_rate) + call export_to_lilac(lilac_a2l_Faxa_rainl, rain_largescale_rate) + call export_to_lilac(lilac_a2l_Faxa_snowc, snow_convective_rate) + call export_to_lilac(lilac_a2l_Faxa_snowl, snow_largescale_rate) + call export_to_lilac(lilac_a2l_Faxa_swndr, swnirdir, no_negative) + call export_to_lilac(lilac_a2l_Faxa_swvdr, swvisdir, no_negative) + call export_to_lilac(lilac_a2l_Faxa_swndf, swnirdif, no_negative) + call export_to_lilac(lilac_a2l_Faxa_swvdf, swvisdif, no_negative) + + ! ------------------------------------------------------------------------ + ! Run ctsm via lilac + ! ------------------------------------------------------------------------ + + ! FIXME(wjs, 2020-01-01) Use correct values for restart and stop alarms + call lilac_run( & + write_restarts_now = restart_flag, & + stop_now = .false.) + + ! ------------------------------------------------------------------------ + ! Import data from ctsm via lilac + ! ------------------------------------------------------------------------ + + call import_from_lilac(lilac_l2a_Sl_t, tsk) + call import_from_lilac(lilac_l2a_Sl_tref, t2) + call import_from_lilac(lilac_l2a_Sl_qref, qref) + + call import_from_lilac(lilac_l2a_Sl_avsdr, albedo_visdir) + call import_from_lilac(lilac_l2a_Sl_anidr, albedo_nirdir) + call import_from_lilac(lilac_l2a_Sl_avsdf, albedo_visdif) + call import_from_lilac(lilac_l2a_Sl_anidf, albedo_nirdif) + call calculate_total_albedo( & + ims=ims, ime=ime, jms=jms, jme=jme, & + its=its, ite=ite, jts=jts, jte=jte, & + xland = xland, & + xice = xice, & + xice_threshold = xice_threshold, & + albedo_visdir = albedo_visdir, & + albedo_nirdir = albedo_nirdir, & + albedo_visdif = albedo_visdif, & + albedo_nirdif = albedo_nirdif, & + swvisdir = swvisdir, & + swnirdir = swnirdir, & + swvisdif = swvisdif, & + swnirdif = swnirdif, & + albedo = albedo) + + call import_from_lilac(lilac_l2a_Sl_fv, ust) + + call import_from_lilac(lilac_l2a_Fall_sen, hfx) + call import_from_lilac(lilac_l2a_Fall_lat, lh) + call import_from_lilac(lilac_l2a_Fall_evap, qfx) + + ! At least for now, use CTSM's Sl_z0m (momentum roughness length) for both background + ! roughness length and thermal time-varying roughness length, even though it's + ! possible that those two should differ in some way. + call import_from_lilac(lilac_l2a_Sl_z0m, z0) + + ! ------------------------------------------------------------------------ + ! Calculate derived variables + ! ------------------------------------------------------------------------ + + do j = jts, jte + do i = its, ite + if (is_land(xland(i,j), xice(i,j), xice_threshold)) then + + ! Convert from specific humidity to mixing ratio. Note that qref is specific + ! humidity at 2m, whereas qsfc is supposed to be specified right at the + ! surface. So there isn't a perfect correspondence between the two, but + ! given that qsfc is just being used as a diagnostic quantity when using + ! CTSM (for now), we won't worry about this. + qsfc(i,j) = qref(i,j) / (1. - qref(i,j)) + + ! CTSM assumes an emissivity of 1 + emiss(i,j) = 1. + + ! At least for now, use CTSM's Sl_z0m (momentum roughness length) for both + ! background roughness length and thermal time-varying roughness length, even + ! though it's possible that those two should differ in some way. + znt(i,j) = z0(i,j) + + ! Flip signs of heat and moisture fluxes: CTSM sends as positive downwards; + ! WRF expects them to be positive upwards + hfx(i,j) = -hfx(i,j) + lh(i,j) = -lh(i,j) + qfx(i,j) = -qfx(i,j) + end if + end do + end do + + first_call = .false. + + contains + + subroutine export_to_lilac(field_index, var_2d , no_negative) + ! Reshape var_2d to 1d for LILAC, then set the appropriate atm2lnd variable in LILAC + + use ctsm_LilacCouplingFields, only : lilac_atm2lnd + + integer, intent(in) :: field_index ! one of the field indices defined in ctsm_LilacCouplingFieldIndices + real, intent(in) :: var_2d(ims: , jms: ) + + real(r8), allocatable :: var_1d(:) + + logical , optional, intent(in) :: no_negative ! convert negative values to zero + + call convert_2d_to_1d( & + ide=ide, jde=jde, & + ims=ims, ime=ime, & + jms=jms, jme=jme, & + its=its, ite=ite, & + jts=jts, jte=jte, & + var_2d=var_2d, var_1d=var_1d) + + ! zero-ing out negative value when we no_negative flag. + if (present (no_negative)) then + if (no_negative) then + var_1d = 0.5* (var_1d +abs(var_1d)) + end if + end if + + call lilac_atm2lnd( & + field_index = field_index, & + data = var_1d) + + end subroutine export_to_lilac + + subroutine import_from_lilac(field_index, var_2d) + ! Get the appropriate lnd2atm variable from LILAC, then reshape 1d LILAC variable to + ! fill var_2d + + use ctsm_LilacCouplingFields, only : lilac_lnd2atm + + integer, intent(in) :: field_index ! one of the field indices defined in ctsm_LilacCouplingFieldIndices + real, intent(inout) :: var_2d(ims: , jms: ) + + real(r8), allocatable :: var_1d(:) + integer :: num_points + + call get_num_points( & + ide=ide, jde=jde, & + its=its, ite=ite, & + jts=jts, jte=jte, & + num_points = num_points) + + allocate(var_1d(num_points)) + + call lilac_lnd2atm( & + field_index = field_index, & + data = var_1d) + + call convert_1d_to_2d( & + ide=ide, jde=jde, & + ims=ims, ime=ime, & + jms=jms, jme=jme, & + its=its, ite=ite, & + jts=jts, jte=jte, & + xland=xland, xice=xice, xice_threshold=xice_threshold, & + var_1d=var_1d, var_2d=var_2d) + + end subroutine import_from_lilac + + end subroutine ctsm_run + !----------------------------------------------------------------------- + + !----------------------------------------------------------------------- + subroutine calculate_total_albedo( & + ims, ime, jms, jme, its, ite, jts, jte, & + xland, xice, xice_threshold, & + albedo_visdir, albedo_nirdir, albedo_visdif, albedo_nirdif, & + swvisdir, swnirdir, swvisdif, swnirdif, & + albedo) + ! Calculate total albedo from its 4 components, based on ratio of incoming SW in the + ! 4 components. + + integer, intent(in) :: ims + integer, intent(in) :: ime + integer, intent(in) :: jms + integer, intent(in) :: jme + integer, intent(in) :: its + integer, intent(in) :: ite + integer, intent(in) :: jts + integer, intent(in) :: jte + + real, intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water) + real, intent(in) :: xice(ims: , jms: ) ! fraction of grid that is seaice + real, intent(in) :: xice_threshold ! fraction of grid determining seaice + + real, intent(in) :: albedo_visdir(ims: , jms: ) ! albedo, visible direct [-] + real, intent(in) :: albedo_nirdir(ims: , jms: ) ! albedo, near IR direct [-] + real, intent(in) :: albedo_visdif(ims: , jms: ) ! albedo, visible diffuse [-] + real, intent(in) :: albedo_nirdif(ims: , jms: ) ! albedo, near IR diffuse [-] + real, intent(in) :: swvisdir(ims: , jms: ) ! solar rad onto surface, visible direct [W m-2] + real, intent(in) :: swnirdir(ims: , jms: ) ! solar rad onto surface, near IR direct [W m-2] + real, intent(in) :: swvisdif(ims: , jms: ) ! solar rad onto surface, visible diffuse [W m-2] + real, intent(in) :: swnirdif(ims: , jms: ) ! solar rad onto surface, near IR diffuse [W m-2] + + real, intent(inout) :: albedo(ims: , jms: ) ! total surface albedo [-] + + real, parameter :: albedo_default = 0.3 + + integer :: i, j + real :: sw_tot ! total solar rad onto surface [W m-2] + + do j = jms, jme + do i = ims, ime + if (is_land(xland(i,j), xice(i,j), xice_threshold)) then + sw_tot = swvisdir(i,j) + swnirdir(i,j) + swvisdif(i,j) + swnirdif(i,j) + if (sw_tot > 0.) then + albedo(i,j) = & + albedo_visdir(i,j) * (swvisdir(i,j) / sw_tot) + & + albedo_nirdir(i,j) * (swnirdir(i,j) / sw_tot) + & + albedo_visdif(i,j) * (swvisdif(i,j) / sw_tot) + & + albedo_nirdif(i,j) * (swnirdif(i,j) / sw_tot) + else + ! Night; albedo shouldn't matter; use coefficients from module_sf_clm + albedo(i,j) = & + albedo_visdir(i,j) * 0.35 + & + albedo_nirdir(i,j) * 0.35 + & + albedo_visdif(i,j) * 0.15 + & + albedo_nirdif(i,j) * 0.15 + end if + + if (abs(albedo(i,j) - 1.) < 1.e-5) then + ! CTSM gives albedo values of 1 at night. To avoid problems (in case CTSM + ! thinks it's night but the albedo is still needed in WRF), use some more + ! reasonable value in this case. + albedo(i,j) = albedo_default + end if + end if + end do + end do + end subroutine calculate_total_albedo + +#endif + ! endif WRF_USE_CTSM + +end module module_sf_ctsm diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index d082f547ad..1f3f74da07 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -94,6 +94,7 @@ SUBROUTINE surface_driver( & #if (EM_CORE==1) ! & ,lakemask, lakeflag & !lake & ,lakemask & !lake + , restart_flag & ! restart_flag #endif ! cyl ocean variable ,OM_TMP,OM_S,OM_U,OM_V,OM_DEPTH,OM_ML,OM_LON & @@ -311,6 +312,7 @@ SUBROUTINE surface_driver( & ,RUCLSMSCHEME & ,PXLSMSCHEME & ,CLMSCHEME & + ,CTSMSCHEME & ,SSIBSCHEME & !ssib ,MYNNSFCSCHEME & ,OMLSCHEME & @@ -354,6 +356,9 @@ SUBROUTINE surface_driver( & USE module_sf_noah_seaice_drv #ifdef WRF_USE_CLM USE module_sf_clm +#endif +#ifdef WRF_USE_CTSM + USE module_sf_ctsm, only : ctsm_run #endif USE module_sf_ssib ! ssib USE module_sf_ruclsm @@ -1350,7 +1355,8 @@ SUBROUTINE surface_driver( & tksatu3d real, dimension(ims:ime,jms:jme ),intent(in) :: lakedepth2d #if (EM_CORE==1) - real , dimension(ims:ime,jms:jme ) :: lakemask + real , dimension(ims:ime,jms:jme ) :: lakemask + logical , intent(in) :: restart_flag ! INTEGER :: lakeflag #endif ! logical, dimension(ims:ime,jms:jme ),intent(in) :: lake @@ -3861,6 +3867,184 @@ SUBROUTINE surface_driver( & ! end of CLM scehme ! ------------------------------------------------------------------- +#endif + +#ifdef WRF_USE_CTSM + CASE (CTSMSCHEME) + + IF (MYJ) call wrf_error_fatal('CTSM is not currently compatible with MYJ. Please pick a different PBL scheme,') + + IF (PRESENT(qv_curr) .AND. PRESENT(rainbl) .AND. & + ! PRESENT(emiss) .AND. PRESENT(t2) .AND. & + ! PRESENT(declin) .AND. PRESENT(coszen) .AND. & + ! PRESENT(hrang) .AND. PRESENT( xlat_urb2d) .AND. & + ! PRESENT(dzr) .AND. & + ! PRESENT( dzb) .AND. PRESENT(dzg) .AND. & + ! PRESENT(tr_urb2d) .AND. PRESENT(tb_urb2d) .AND. & + ! PRESENT(tg_urb2d) .AND. PRESENT(tc_urb2d) .AND. & + ! PRESENT(qc_urb2d) .AND. PRESENT(uc_urb2d) .AND. & + ! PRESENT(xxxr_urb2d) .AND. PRESENT(xxxb_urb2d) .AND. & + ! PRESENT(xxxg_urb2d) .AND. & + ! PRESENT(xxxc_urb2d) .AND. PRESENT(trl_urb3d) .AND. & + ! PRESENT(tbl_urb3d) .AND. PRESENT(tgl_urb3d) .AND. & + ! PRESENT(sh_urb2d) .AND. PRESENT(lh_urb2d) .AND. & + ! PRESENT(g_urb2d) .AND. PRESENT(rn_urb2d) .AND. & + ! PRESENT(ts_urb2d) .AND. & + ! PRESENT(frc_urb2d) .AND. PRESENT(utype_urb2d) .AND. & + .TRUE. ) THEN + + +!------------------------------------------------------------------ + IF ( FRACTIONAL_SEAICE == 1) THEN + ! The fields passed to LSM need to represent the full ice values, not + ! the fractional values. Convert ALBEDO and EMISS from the blended value + ! to a value representing only the sea-ice portion. Albedo over open + ! water is taken to be 0.08. Emissivity over open water is taken to be 0.98 + DO j = j_start(ij) , j_end(ij) + DO i = i_start(ij) , i_end(ij) + IF ( ( XICE(I,J) .GE. XICE_THRESHOLD ) .AND. ( XICE(i,j) .LE. 1 ) ) THEN + ALBEDO(I,J) = (ALBEDO(I,J)-(1.-XICE(I,J))*0.08)/XICE(I,J) + EMISS(I,J) = (EMISS(I,J)-(1.-XICE(I,J))*0.98)/XICE(I,J) + ENDIF + ENDDO + ENDDO + IF ( isisfc ) THEN + ! Use surface layer routine values from the ice portion of grid point + ELSE + ! + ! We don't have surface layer routine values at this time, so + ! just use what we have. Use ice component of TSK + ! + CALL get_local_ice_tsk( ims, ime, jms, jme, & + i_start(ij), i_end(ij), & + j_start(ij), j_end(ij), & + itimestep, .false., tice2tsk_if2cold, & + XICE, XICE_THRESHOLD, & + SST, TSK, TSK_SEA, TSK_LOCAL ) + + DO j = j_start(ij) , j_end(ij) + DO i = i_start(ij) , i_end(ij) + TSK(i,j) = TSK_LOCAL(i,j) + ENDDO + ENDDO + ENDIF + ENDIF + + + CALL ctsm_run( & + ! bounds + ids=ids, ide=ide, jds=jds, jde=jde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, & + its=i_start(ij), ite=i_end(ij), jts=j_start(ij), jte=j_end(ij), & + + ! restart flag + restart_flag=restart_flag,& + + ! general information + dt = dt, & + xland = xland, & + xice = xice, & + xice_threshold = xice_threshold, & + + ! atm -> lnd variables + dz8w = dz8w, & + ht = ht, & + u_phy = u_phy, & + v_phy = v_phy, & + p8w = p_phy, & + t_phy = t_phy, & + th_phy = th_phy, & + qv_curr = qv_curr, & + rainbl = rainbl, & + sr = sr, & + glw = glw, & + swvisdir = swvisdir, & + swvisdif = swvisdif, & + swnirdir = swnirdir, & + swnirdif = swnirdif, & + + ! lnd -> atm variables + tsk = tsk, & + t2 = t2, & + qsfc = qsfc, & + albedo = albedo, & + ust = ust, & + hfx = hfx, & + lh = lh, & + qfx = qfx, & + emiss = emiss, & + z0 = z0, & + znt = znt) + + call seaice_noah( SEAICE_ALBEDO_OPT, SEAICE_ALBEDO_DEFAULT, SEAICE_THICKNESS_OPT, & + & SEAICE_THICKNESS_DEFAULT, SEAICE_SNOWDEPTH_OPT, & + & SEAICE_SNOWDEPTH_MAX, SEAICE_SNOWDEPTH_MIN, & + & t_phy, qv_curr, p8w, dz8w, num_soil_layers, dt, frpcpn, sr, & + & glw, swdown, rainbl, snoalb, qgh, xice, xice_threshold, & + & albsi, icedepth, snowsi, & + & tslb, emiss, albedo, z0, tsk, snow, snowc, snowh, & + & chs, chs2, cqs2, & + & br, znt, lh, hfx, qfx, potevp, grdflx, qsfc, acsnow, & + & acsnom, snopcx, sfcrunoff, noahres, & + & sf_urban_physics, b_t_bep, b_q_bep, rho, & + & ids,ide, jds,jde, kds,kde, & + & ims,ime, jms,jme, kms,kme, & + & i_start(ij),i_end(ij), j_start(ij),j_end(ij), kts,kte ) + + IF ( FRACTIONAL_SEAICE == 1 ) THEN + ! LSM Returns full land/ice values, no fractional values. + ! We return to a fractional component here. SFLX currently hard-wires + ! emissivity over sea ice to 0.98, the same value as over open water, so + ! the fractional consideration doesn't have any effect for emissivity. + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + IF ( ( XICE(I,J) .GE. XICE_THRESHOLD) .AND. ( XICE(i,j) .LE. 1.0 ) ) THEN + albedo(i,j) = ( albedo(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * 0.08 ) + emiss(i,j) = ( emiss(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * 0.98 ) + ENDIF + ENDDO + ENDDO + + IF ( isisfc ) THEN + DO j=j_start(ij),j_end(ij) + DO i=i_start(ij),i_end(ij) + IF ( ( XICE(I,J) .GE. XICE_THRESHOLD) .AND. ( XICE(i,j) .LE. 1.0 ) ) THEN + ! Weighted average of fields between ice-cover values and open-water values. + flhc(i,j) = ( flhc(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * flhc_sea(i,j) ) + flqc(i,j) = ( flqc(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * flqc_sea(i,j) ) + cpm(i,j) = ( cpm(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * cpm_sea(i,j) ) + cqs2(i,j) = ( cqs2(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * cqs2_sea(i,j) ) + chs2(i,j) = ( chs2(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * chs2_sea(i,j) ) + chs(i,j) = ( chs(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * chs_sea(i,j) ) + qsfc(i,j) = ( qsfc(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * qsfc_sea(i,j) ) + qgh(i,j) = ( qgh(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * qgh_sea(i,j) ) + qz0(i,j) = ( qz0(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * qz0_sea(i,j) ) + ! print *,'hfx =',hfx_sea(170,20) + ! print *,'XICE =',XICE(170,20) + ! print *,'QSFC =',QSFC(170,20) + hfx(i,j) = ( hfx(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * hfx_sea(i,j) ) + qfx(i,j) = ( qfx(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * qfx_sea(i,j) ) + lh(i,j) = ( lh(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * lh_sea(i,j) ) +!save old tsk_ice + tsk_save(i,j) = tsk(i,j) + tsk(i,j) = ( tsk(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * tsk_sea(i,j) ) + ENDIF + ENDDO + ENDDO + ELSE + DO j = j_start(ij) , j_end(ij) + DO i = i_start(ij) , i_end(ij) + IF ( ( XICE(I,J) .GE. XICE_THRESHOLD ) .AND. ( XICE(i,j) .LE. 1.0 ) ) THEN + ! Compute TSK as the open-water and ice-cover average + tsk(i,j) = ( tsk(i,j) * XICE(i,j) ) + ( (1.0-XICE(i,j)) * tsk_sea(i,j) ) + ENDIF + ENDDO + ENDDO + ENDIF + ENDIF + + ENDIF + #endif CASE (SSIBSCHEME) diff --git a/share/input_wrf.F b/share/input_wrf.F index c162dfc798..13fd3d4886 100644 --- a/share/input_wrf.F +++ b/share/input_wrf.F @@ -740,6 +740,10 @@ SUBROUTINE input_wrf ( fid , grid , config_flags , switch , ierr ) ! All is well. Noah-MP and Noah have compatible wrfinput files. ELSE IF ( ( config_flags%sf_surface_physics == NOAHMPSCHEME ) .and. ( itmp == LSMSCHEME ) ) then ! All is well. Noah-MP and Noah have compatible wrfinput files. +#if ( NMM_CORE != 1 ) + ELSE IF ( ( config_flags%sf_surface_physics == CTSMSCHEME ) ) then + ! All is well. CTSM doesn't care about the wrfinput file. +#endif ELSE call wrf_message("----------------- ERROR -------------------") WRITE(wrf_err_message,'("namelist : sf_surface_physics = ",I10)') config_flags%sf_surface_physics diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index db28b315a3..23674025c1 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -2107,6 +2107,30 @@ END FUNCTION bep_bem_nbui_max END IF #endif +!----------------------------------------------------------------------- +! The CTSM scheme may not even be compiled, so make sure it is not allowed +! to be run if the code is not available. +!----------------------------------------------------------------------- + +#if !defined ( WRF_USE_CTSM ) + oops = 0 + DO i = 1, model_config_rec % max_dom + IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE + IF ( model_config_rec%sf_surface_physics(i) .EQ. CTSMSCHEME ) THEN + oops = oops + 1 + END IF + ENDDO ! Loop over domains + IF ( oops .GT. 0 ) THEN + wrf_err_message = '--- ERROR: The CTSM surface scheme was requested in the namelist.input file.' + CALL wrf_debug ( 0, TRIM(wrf_err_message) ) + wrf_err_message = '--- ERROR: However, the WRF CTSM scheme was not compiled in WRF.' + CALL wrf_debug ( 0, TRIM(wrf_err_message) ) + wrf_err_message = '--- ERROR: Please read doc/README.CTSM for how to compile WRF with CTSM.' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + count_fatal_error = count_fatal_error + 1 + END IF +#endif + !----------------------------------------------------------------------- ! grav_settling = 1 must be turned off for mp_physics=28. !----------------------------------------------------------------------- @@ -3029,6 +3053,11 @@ SUBROUTINE set_physics_rconfigs model_config_rec % num_soil_layers = 2 ELSE IF ( model_config_rec % sf_surface_physics(1) .EQ. CLMSCHEME ) THEN model_config_rec % num_soil_layers = 10 +#if (EM_CORE == 1) + ELSE IF ( model_config_rec % sf_surface_physics(1) .EQ. CTSMSCHEME ) THEN + ! Using 4 for the sake of the sea ice scheme + model_config_rec % num_soil_layers = 4 +#endif ELSE IF ( model_config_rec % sf_surface_physics(1) .EQ. SSIBSCHEME ) THEN model_config_rec % num_soil_layers = 3 #if (NMM_CORE == 1) diff --git a/wrftladj/solve_em_ad.F b/wrftladj/solve_em_ad.F index cb02ec72aa..1ab1b56328 100644 --- a/wrftladj/solve_em_ad.F +++ b/wrftladj/solve_em_ad.F @@ -149,6 +149,7 @@ SUBROUTINE solve_em_ad ( grid , config_flags & ! Flag for producing diagnostic fields (e.g., radar reflectivity) LOGICAL :: diag_flag + LOGICAL :: restart_flag ! tells if it is a restart timestep to write restart files #if (WRF_CHEM==1) ! Index cross-referencing array for tendency accumulation @@ -314,6 +315,13 @@ SUBROUTINE solve_em_ad ( grid , config_flags & ! End of adaptive time step modifications !----------------------------------------------------------------------------- ! +! Set restart flag value history output time +!----------------------------------------------------------------------------- + restart_flag = .false. + if ( Is_alarm_tstep(grid%domain_clock, grid%alarms(restart_alarm)) ) then + restart_flag = .true. + endif +! ! Set diagnostic flag value history output time !----------------------------------------------------------------------------- ! if ( Is_alarm_tstep(grid_clock, grid_alarms(HISTORY_ALARM)) ) then @@ -677,6 +685,7 @@ SUBROUTINE solve_em_ad ( grid , config_flags & , k_start , k_end & , f_flux & , aerocu & + , restart_flag & ) #ifdef DM_PARALLEL