diff --git a/Makefile b/Makefile index 3974c7e12..935bc7484 100644 --- a/Makefile +++ b/Makefile @@ -16,11 +16,11 @@ gnu: # BUILDTARGET GNU Fortran, C, and C++ compilers "CC_SERIAL = gcc" \ "CXX_SERIAL = g++" \ "FFLAGS_PROMOTION = -fdefault-real-8 -fdefault-double-8" \ - "FFLAGS_OPT = -std=f2008 -O3 -ffree-line-length-none -fconvert=big-endian -ffree-form" \ + "FFLAGS_OPT = -std=f2008 -O3 -fallow-argument-mismatch -ffree-line-length-none -fconvert=big-endian -ffree-form" \ "CFLAGS_OPT = -O3" \ "CXXFLAGS_OPT = -O3" \ "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -std=f2008 -g -ffree-line-length-none -fconvert=big-endian -ffree-form -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \ + "FFLAGS_DEBUG = -g -ffree-line-length-none -fallow-argument-mismatch -fconvert=big-endian -ffree-form -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \ "CFLAGS_DEBUG = -g" \ "CXXFLAGS_DEBUG = -g" \ "LDFLAGS_DEBUG = -g" \ @@ -154,7 +154,7 @@ nvhpc: # BUILDTARGET NVIDIA HPC SDK "FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \ "CFLAGS_DEBUG = -O0 -g -traceback" \ "CXXFLAGS_DEBUG = -O0 -g -traceback" \ - "LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \ + "LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \ "FFLAGS_OMP = -mp" \ "CFLAGS_OMP = -mp" \ "FFLAGS_ACC = -Mnofma -acc -gpu=cc70,cc80 -Minfo=accel" \ @@ -184,7 +184,7 @@ pgi: # BUILDTARGET PGI compiler suite "FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \ "CFLAGS_DEBUG = -O0 -g -traceback" \ "CXXFLAGS_DEBUG = -O0 -g -traceback" \ - "LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \ + "LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \ "FFLAGS_OMP = -mp" \ "CFLAGS_OMP = -mp" \ "FFLAGS_ACC = -Mnofma -acc -Minfo=accel" \ @@ -216,7 +216,7 @@ pgi-summit: # BUILDTARGET PGI compiler suite w/OpenACC options for ORNL Summit "FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \ "CFLAGS_DEBUG = -O0 -g -traceback" \ "CXXFLAGS_DEBUG = -O0 -g -traceback" \ - "LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \ + "LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \ "FFLAGS_OMP = -mp" \ "CFLAGS_OMP = -mp" \ "PICFLAG = -fpic" \ @@ -667,10 +667,10 @@ intel: # BUILDTARGET Intel oneAPI Fortran, C, and C++ compiler suite "CFLAGS_OPT = -O3" \ "CXXFLAGS_OPT = -O3" \ "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -g -convert big_endian -free -check bounds,pointers,arg_temp_created,format,shape,contiguous -fpe0 -traceback" \ + "FFLAGS_DEBUG = -g -convert big_endian -free -check all -fpe0 -traceback" \ "CFLAGS_DEBUG = -g -traceback" \ "CXXFLAGS_DEBUG = -g -traceback" \ - "LDFLAGS_DEBUG = -g -traceback" \ + "LDFLAGS_DEBUG = -g -fpe0 -traceback" \ "FFLAGS_OMP = -qopenmp" \ "CFLAGS_OMP = -qopenmp" \ "PICFLAG = -fpic" \ @@ -690,7 +690,7 @@ intel-xd2000: "CC_SERIAL = cc" \ "CXX_SERIAL = CC" \ "FFLAGS_PROMOTION = -real-size 64" \ - "FFLAGS_OPT = -O3 -convert big_endian -FR -march=core-avx2 -mtune=core-avx2" \ + "FFLAGS_OPT = -O3 -convert big_endian -FR" \ "CFLAGS_OPT = -O3 -std=gnu90" \ "CXXFLAGS_OPT = -O3" \ "LDFLAGS_OPT = -O3" \ @@ -715,7 +715,7 @@ intel2-xd2000: "CC_SERIAL = cc" \ "CXX_SERIAL = CC" \ "FFLAGS_PROMOTION = -real-size 64" \ - "FFLAGS_OPT = -O2 -convert big_endian -FR -march=core-avx2 -mtune=core-avx2" \ + "FFLAGS_OPT = -O2 -convert big_endian -FR" \ "CFLAGS_OPT = -O2 -std=gnu90" \ "CXXFLAGS_OPT = -O2" \ "LDFLAGS_OPT = -O2" \ @@ -741,11 +741,21 @@ gfortran-xd2000: # BUILDTARGET GNU Fortran, C, and C++ compilers "CC_SERIAL = cc" \ "CXX_SERIAL = CC" \ "FFLAGS_PROMOTION = -fdefault-real-8 -fdefault-double-8" \ - "FFLAGS_OPT = -O3 -ffree-line-length-none -fconvert=big-endian -ffree-form -fallow-argument-mismatch" \ - "CFLAGS_OPT = -O3" \ - "CXXFLAGS_OPT = -O3" \ - "LDFLAGS_OPT = -O3" \ - "FFLAGS_DEBUG = -g -ffree-line-length-none -fconvert=big-endian -free -ffree-form -fallow-argument-mismatch -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \ + "FFLAGS_OPT = -O3 -fopenmp \ + -march=native -mtune=native \ + -fno-fast-math \ + -ffp-contract=off \ + -ffree-line-length-none -fconvert=big-endian -ffree-form \ + -fallow-argument-mismatch" \ + "CFLAGS_OPT = -O3 -march=native -mtune=native" \ + "CXXFLAGS_OPT = -O3 -march=native -mtune=native" \ + "LDFLAGS_OPT = -O3" \ + "FFLAGS_DEBUG = -O0 -g \ + -fcheck=all -fbacktrace \ + -finit-real=nan -finit-integer=-999 \ + -ffpe-trap=invalid,zero,overflow \ + -ffree-line-length-none -fconvert=big-endian -ffree-form \ + -fallow-argument-mismatch" \ "CFLAGS_DEBUG = -g" \ "CXXFLAGS_DEBUG = -g" \ "LDFLAGS_DEBUG = -g" \ @@ -762,6 +772,65 @@ gfortran-xd2000: # BUILDTARGET GNU Fortran, C, and C++ compilers "OPENACC = $(OPENACC)" \ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" ) +cray-xd2000: + ( $(MAKE) all \ + "FC_PARALLEL = ftn" \ + "CC_PARALLEL = cc" \ + "CXX_PARALLEL = CC" \ + "FC_SERIAL = ftn" \ + "CC_SERIAL = cc" \ + "CXX_SERIAL = CC" \ + "FFLAGS_PROMOTION = -sreal64 " \ + "FFLAGS_OPT = -Ofast -hcpu=x86-turin -hipa1 -ffree" \ + "CFLAGS_OPT = -Ofast" \ + "CXXFLAGS_OPT = -Ofast" \ + "LDFLAGS_OPT = -hbyteswapio" \ + "FFLAGS_DEBUG = -eD -O0 -ffree" \ + "CFLAGS_DEBUG = -O0 -g -Weverything" \ + "CXXFLAGS_DEBUG = -O0 -g -Weverything" \ + "LDFLAGS_DEBUG = -eD -O0 -hbyteswapio" \ + "FFLAGS_OMP = -homp" \ + "CFLAGS_OMP = -fopenmp" \ + "FFLAGS_ACC =" \ + "CFLAGS_ACC =" \ + "BUILD_TARGET = $(@)" \ + "CORE = $(CORE)" \ + "DEBUG = $(DEBUG)" \ + "USE_PAPI = $(USE_PAPI)" \ + "OPENMP = $(OPENMP)" \ + "OPENACC = $(OPENACC)" \ + "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" ) + +nvhpc-xd2000: # BUILDTARGET NVIDIA HPC SDK + ( $(MAKE) all \ + "FC_PARALLEL = ftn" \ + "CC_PARALLEL = cc" \ + "CXX_PARALLEL = CC" \ + "FC_SERIAL = ftn" \ + "CC_SERIAL = cc" \ + "CXX_SERIAL = CC" \ + "FFLAGS_PROMOTION = -r8" \ + "FFLAGS_OPT = -gopt -O4 -byteswapio -Mfree" \ + "CFLAGS_OPT = -gopt -O3" \ + "CXXFLAGS_OPT = -gopt -O3" \ + "LDFLAGS_OPT = -gopt -O3" \ + "FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Mchkstk -byteswapio -Mfree -traceback" \ + "CFLAGS_DEBUG = -O0 -g -traceback" \ + "CXXFLAGS_DEBUG = -O0 -g -traceback" \ + "LDFLAGS_DEBUG = -O0 -g -traceback" \ + "FFLAGS_OMP = -mp" \ + "CFLAGS_OMP = -mp" \ + "FFLAGS_ACC = -Mnofma -acc -gpu=cc70,cc80 -Minfo=accel" \ + "CFLAGS_ACC =" \ + "PICFLAG = -fpic" \ + "BUILD_TARGET = $(@)" \ + "CORE = $(CORE)" \ + "DEBUG = $(DEBUG)" \ + "USE_PAPI = $(USE_PAPI)" \ + "OPENMP = $(OPENMP)" \ + "OPENACC = $(OPENACC)" \ + "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DCPRPGI" ) + CPPINCLUDES = FCINCLUDES = LIBS = @@ -821,10 +890,10 @@ else # Not using PIO, using SMIOL endif ifneq "$(NETCDF)" "" -ifneq ($(wildcard $(NETCDF)/lib/libnetcdf.*), ) +ifneq ($(wildcard $(NETCDF)/lib), ) NETCDFLIBLOC = lib endif -ifneq ($(wildcard $(NETCDF)/lib64/libnetcdf.*), ) +ifneq ($(wildcard $(NETCDF)/lib64), ) NETCDFLIBLOC = lib64 endif CPPINCLUDES += -I$(NETCDF)/include @@ -845,10 +914,10 @@ endif ifneq "$(PNETCDF)" "" -ifneq ($(wildcard $(PNETCDF)/lib/libpnetcdf.*), ) +ifneq ($(wildcard $(PNETCDF)/lib), ) PNETCDFLIBLOC = lib endif -ifneq ($(wildcard $(PNETCDF)/lib64/libpnetcdf.*), ) +ifneq ($(wildcard $(PNETCDF)/lib64), ) PNETCDFLIBLOC = lib64 endif CPPINCLUDES += -I$(PNETCDF)/include @@ -1008,11 +1077,11 @@ else # else ifdef $(TIMER_LIB) endif # endif ifdef $(TIMER_LIB) ifeq "$(TAU)" "true" - LINKER=tau_f90.sh + LNK=tau_f90.sh CPPINCLUDES += -DMPAS_TAU -DMPAS_TAU_TIMERS TAU_MESSAGE="TAU Hooks are on." else - LINKER=$(FC) + LNK=$(FC) TAU_MESSAGE="TAU Hooks are off." endif @@ -1042,6 +1111,28 @@ else OPENACC_MESSAGE="MPAS was built without OpenACC accelerator support." endif +ifneq ($(wildcard .mpas_core_*), ) # CHECK FOR BUILT CORE + +ifneq ($(wildcard .mpas_core_$(CORE)), ) # CHECK FOR SAME CORE AS ATTEMPTED BUILD. + override AUTOCLEAN=false + CONTINUE=true +else + LAST_CORE=`cat .mpas_core_*` + +ifeq "$(AUTOCLEAN)" "true" # CHECK FOR CLEAN PRIOR TO BUILD OF A NEW CORE. + CONTINUE=true + AUTOCLEAN_MESSAGE="Infrastructure was cleaned prior to building ." +else + CONTINUE=false +endif # END OF AUTOCLEAN CHECK + +endif # END OF CORE=LAST_CORE CHECK + +else + + override AUTOCLEAN=false + CONTINUE=true +endif # END IF BUILT CORE CHECK ifneq ($(wildcard namelist.$(NAMELIST_SUFFIX)), ) # Check for generated namelist file. NAMELIST_MESSAGE="A default namelist file (namelist.$(NAMELIST_SUFFIX).defaults) has been generated, but namelist.$(NAMELIST_SUFFIX) has not been modified." @@ -1098,119 +1189,12 @@ report_builds: @echo "CORE=$(CORE)" endif +ifeq "$(CONTINUE)" "true" all: mpas_main - +else +all: clean_core endif -# -# The rebuild_check target determines whether the shared framework or $(CORE) were -# previously compiled with incompatible options, and stops the build with an error -# message if so. -# -rebuild_check: - @# - @# Write current build options to a file .build_opts.tmp, to later be - @# compared with build options use for the shared framework or core. - @# Only build options that affect compatibility are written, while options - @# like $(RM), $(BUILD_TARGET), and $(CORE) are not. - @# - $(shell printf "FC=$(FC)\n$\ - CC=$(CC)\n$\ - CXX=$(CXX)\n$\ - SFC=$(SFC)\n$\ - SCC=$(SCC)\n$\ - CFLAGS=$(CFLAGS)\n$\ - CXXFLAGS=$(CXXFLAGS)\n$\ - FFLAGS=$(FFLAGS)\n$\ - LDFLAGS=$(LDFLAGS)\n$\ - CPPFLAGS=$(CPPFLAGS)\n$\ - LIBS=$(LIBS)\n$\ - CPPINCLUDES=$(CPPINCLUDES)\n$\ - OPENMP=$(OPENMP)\n$\ - OPENMP_OFFLOAD=$(OPENMP_OFFLOAD)\n$\ - OPENACC=$(OPENACC)\n$\ - TAU=$(TAU)\n$\ - PICFLAG=$(PICFLAG)\n$\ - TIMER_LIB=$(TIMER_LIB)\n$\ - GEN_F90=$(GEN_F90)\n" | sed 's/-DMPAS_EXE_NAME=[^[:space:]]*//' | sed 's/-DMPAS_NAMELIST_SUFFIX=[^[:space:]]*//' | sed 's/-DCORE_[^[:space:]]*//' | sed 's/-DMPAS_GIT_VERSION=[^[:space:]]*//' > .build_opts.tmp ) - - @# - @# PREV_BUILD is set to "OK" if the shared framework and core are either - @# clean or were previously compiled with compatible options. Otherwise, - @# PREV_BUILD is set to "shared framework" if the shared framework was - @# built with incompatible options, or "$(CORE) core" if the core was - @# built with incompatible options. - @# - $(eval PREV_BUILD := $(shell $\ - if [ -f ".build_opts.framework" ]; then $\ - cmp -s .build_opts.tmp .build_opts.framework; $\ - if [ $$? -eq 0 ]; then $\ - stat=0; $\ - else $\ - stat=1; $\ - x="shared framework"; $\ - if [ "$(AUTOCLEAN)" = "true" ]; then $\ - cp .build_opts.tmp .build_opts.framework; $\ - fi; $\ - fi $\ - else $\ - stat=0; $\ - cp .build_opts.tmp .build_opts.framework; $\ - fi; $\ - : ; $\ - : At this this point, stat is already set, and we should only ; $\ - : set it to 1 but never to 0, as that might mask an incompatibility ; $\ - : in the framework build. ; $\ - : ; $\ - if [ -f ".build_opts.$(CORE)" ]; then $\ - cmp -s .build_opts.tmp .build_opts.$(CORE); $\ - if [ $$? -ne 0 ]; then $\ - stat=1; $\ - if [ "$$x" = "" ]; then $\ - x="$(CORE) core"; $\ - else $\ - x="$$x and $(CORE) core"; $\ - fi; $\ - if [ "$(AUTOCLEAN)" = "true" ]; then $\ - cp .build_opts.tmp .build_opts.$(CORE); $\ - fi; $\ - fi; $\ - else $\ - if [ $$stat -eq 0 ]; then $\ - cp .build_opts.tmp .build_opts.$(CORE); $\ - fi; $\ - fi; $\ - rm -f .build_opts.tmp; $\ - if [ $$stat -eq 1 ]; then $\ - printf "$$x"; $\ - else $\ - printf "OK"; $\ - fi; $\ - )) - - $(if $(findstring and,$(PREV_BUILD)),$(eval VERB=were),$(eval VERB=was)) -ifeq "$(AUTOCLEAN)" "true" - $(if $(findstring framework,$(PREV_BUILD)),$(eval AUTOCLEAN_DEPS+=clean_shared)) - $(if $(findstring core,$(PREV_BUILD)),$(eval AUTOCLEAN_DEPS+=clean_core)) - $(if $(findstring OK,$(PREV_BUILD)), $(eval override AUTOCLEAN=false), ) - $(eval AUTOCLEAN_MESSAGE=The $(PREV_BUILD) $(VERB) cleaned and re-compiled.) -else - $(if $(findstring OK,$(PREV_BUILD)), \ - , \ - $(info ************************************************************************) \ - $(info The $(PREV_BUILD) $(VERB) previously compiled with ) \ - $(info incompatible options. Please do one of the following:) \ - $(info ) \ - $(info - Clean the $(CORE) core, which will also cause the shared) \ - $(info framework to be cleaned; then compile the $(CORE) core.) \ - $(info ) \ - $(info or)\ - $(info ) \ - $(info - Add AUTOCLEAN=true to the build command to automatically clean) \ - $(info and re-compile the $(PREV_BUILD).) \ - $(info ) \ - $(info ************************************************************************) \ - $(error )) endif @@ -1350,7 +1334,7 @@ ifeq "$(OPENACC)" "true" endif # OPENACC eq true -pio_test: openmp_test openacc_test pnetcdf_test +pio_test: openmp_test openacc_test @# @# PIO_VERS will be set to: @# 0 if no working PIO library was detected (and .piotest.log will contain error messages) @@ -1436,10 +1420,9 @@ mpi_f08_test: $(info Checking for mpi_f08 support...) $(eval MPAS_MPI_F08 := $(shell $\ printf "program main\n$\ - & use mpi_f08, only : MPI_Init, MPI_Comm, MPI_INTEGER, MPI_Datatype\n$\ + & use mpi_f08, only : MPI_Init, MPI_Comm\n$\ & integer :: ierr\n$\ & type (MPI_Comm) :: comm\n$\ - & type (MPI_Datatype), parameter :: MPI_INTEGERKIND = MPI_INTEGER\n$\ & call MPI_Init(ierr)\n$\ end program main\n" | sed 's/&/ /' > mpi_f08.f90; $\ $\ @@ -1458,63 +1441,26 @@ mpi_f08_test: $(if $(findstring 1,$(MPAS_MPI_F08)), $(eval MPI_F08_MESSAGE = "Using the mpi_f08 module."), ) $(if $(findstring 1,$(MPAS_MPI_F08)), $(info mpi_f08 module detected.)) - -pnetcdf_test: - @# - @# Create test C programs that look for PNetCDF header file and some symbols in it - @# -ifneq "$(PNETCDF)" "" - @echo "Checking for a working PnetCDF library..." - @printf "#include \"pnetcdf.h\"\n\ - &#include \"mpi.h\"\n\ - &int main(){\n\ - & int err, ncid;\n\ - & err = ncmpi_create(MPI_COMM_WORLD, \"foo.nc\", NC_NOCLOBBER, MPI_INFO_NULL, &ncid);\n\ - & return 0;\n\ - &}\n" | sed 's/&/ /' > pnetcdf.c - @( $(CC) pnetcdf.c $(CPPINCLUDES) $(CFLAGS) $(LDFLAGS) -L$(PNETCDF)/$(PNETCDFLIBLOC) -lpnetcdf -o pnetcdf.out > pnetcdf.log 2>&1; \ - if [ $$? -eq 0 ] ; then \ - echo "$(CC) can compile test PnetCDF C program."; \ - else \ - echo "*********************************************************"; \ - echo "ERROR: Test PnetCDF C program could not be compiled by $(CC)."; \ - echo "Please ensure you have a working PnetCDF library installed."; \ - echo ""; \ - echo "The following compilation command failed with errors:" ; \ - echo "$(CC) pnetcdf.c $(CPPINCLUDES) $(CFLAGS) $(LDFLAGS) -L$(PNETCDF)/$(PNETCDFLIBLOC) -lpnetcdf -o pnetcdf.out"; \ - echo ""; \ - echo "Test program pnetcdf.c and output pnetcdf.log have been left"; \ - echo "in the top-level MPAS directory for further debugging"; \ - echo "*********************************************************"; \ - rm -f pnetcdf.out; exit 1; \ - fi ) - - @rm -f pnetcdf.c pnetcdf.out pnetcdf.log -else - @echo "*********************************************************"; \ - echo "ERROR: The PNETCDF environment variable isn't set."; \ - echo "Please set this variable to where PnetCDF is installed."; \ - echo "*********************************************************"; \ - exit 1 -endif - - ifneq "$(PIO)" "" -MAIN_DEPS = rebuild_check openmp_test openacc_test pnetcdf_test pio_test mpi_f08_test +MAIN_DEPS = openmp_test openacc_test pio_test mpi_f08_test override CPPFLAGS += "-DMPAS_PIO_SUPPORT" else -MAIN_DEPS = rebuild_check openmp_test openacc_test pnetcdf_test mpi_f08_test +MAIN_DEPS = openmp_test openacc_test mpi_f08_test IO_MESSAGE = "Using the SMIOL library." override CPPFLAGS += "-DMPAS_SMIOL_SUPPORT" endif + mpas_main: $(MAIN_DEPS) +ifeq "$(AUTOCLEAN)" "true" + $(RM) .mpas_core_* +endif cd src; $(MAKE) FC="$(FC)" \ CC="$(CC)" \ CXX="$(CXX)" \ SFC="$(SFC)" \ SCC="$(SCC)" \ - LINKER="$(LINKER)" \ + LNK="$(LNK)" \ CFLAGS="$(CFLAGS)" \ CXXFLAGS="$(CXXFLAGS)" \ FFLAGS="$(FFLAGS)" \ @@ -1527,11 +1473,11 @@ mpas_main: $(MAIN_DEPS) FCINCLUDES="$(FCINCLUDES)" \ CORE="$(CORE)"\ AUTOCLEAN="$(AUTOCLEAN)" \ - AUTOCLEAN_DEPS="$(AUTOCLEAN_DEPS)" \ GEN_F90="$(GEN_F90)" \ NAMELIST_SUFFIX="$(NAMELIST_SUFFIX)" \ EXE_NAME="$(EXE_NAME)" + @echo "$(EXE_NAME)" > .mpas_core_$(CORE) if [ -e src/$(EXE_NAME) ]; then mv src/$(EXE_NAME) .; fi ( cd src/core_$(CORE); $(MAKE) ROOT_DIR="$(PWD)" post_build ) @echo "*******************************************************************************" @@ -1553,13 +1499,11 @@ endif @echo $(IO_MESSAGE) @echo "*******************************************************************************" clean: - cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)" AUTOCLEAN="$(AUTOCLEAN)" + cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)" + $(RM) .mpas_core_* $(RM) $(EXE_NAME) $(RM) namelist.$(NAMELIST_SUFFIX).defaults $(RM) streams.$(NAMELIST_SUFFIX).defaults - if [ -f .build_opts.framework ]; then $(RM) .build_opts.framework; fi - if [ -f .build_opts.$(CORE) ]; then $(RM) .build_opts.$(CORE); fi - core_error: @echo "" @echo "*******************************************************************************" @@ -1570,6 +1514,26 @@ core_error: exit 1 error: errmsg +clean_core: + @echo "" + @echo "*******************************************************************************" + @echo " The MPAS infrastructure is currently built for the $(LAST_CORE) core." + @echo " Before building the $(CORE) core, please do one of the following." + @echo "" + @echo "" + @echo " To remove the $(LAST_CORE)_model executable and clean the MPAS infrastructure, run:" + @echo " make clean CORE=$(LAST_CORE)" + @echo "" + @echo " To preserve all executables except $(CORE)_model and clean the MPAS infrastructure, run:" + @echo " make clean CORE=$(CORE)" + @echo "" + @echo " Alternatively, AUTOCLEAN=true can be appended to the make command to force a clean," + @echo " build a new $(CORE)_model executable, and preserve all other executables." + @echo "" + @echo "*******************************************************************************" + @echo "" + exit 1 + else # CORE IF all: error @@ -1597,7 +1561,7 @@ errmsg: @echo " DEBUG=true - builds debug version. Default is optimized version." @echo " USE_PAPI=true - builds version using PAPI for timers. Default is off." @echo " TAU=true - builds version using TAU hooks for profiling. Default is off." - @echo " AUTOCLEAN=true - Enables automatic cleaning and re-compilation of code as needed." + @echo " AUTOCLEAN=true - forces a clean of infrastructure prior to build new core." @echo " GEN_F90=true - Generates intermediate .f90 files through CPP, and builds with them." @echo " TIMER_LIB=opt - Selects the timer library interface to be used for profiling the model. Options are:" @echo " TIMER_LIB=native - Uses native built-in timers in MPAS" diff --git a/README.md b/README.md index 8a736f304..46a664adf 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# MONAN-v2.0.0-rc +# MONAN ## Model for Ocean-laNd-Atmosphere predictioN @@ -6,10 +6,10 @@ MONAN is a community model of the Unified Earth System that has as its main obje The MONAN Model is managed by a scientific committee appointed by INPE's director and has its initial version structure (0.1.0) based on the dynamic core of the MPAS 8.0.1 Model. part of the physics used by MONAN is obtained from the MPAS model and another part obtained from other sources or developed by the community. The MPAS model can be found at the link [GitHub - MPAS-Dev/MPAS-Model: Repository for MPAS models and shared framework releases.](https://github.com/MPAS-Dev/MPAS-Model) +## History -History -==== -- Version 2.0.0-rc (Release Candidate) - This GF version includes a new formulation for the gross entrainment rate that depends on the model grid spacing (Zhao et al 2024, GRL https://doi.org/10.1029/2024GL110735). This may help a smoother transition from non-resolved to resolved scales (grey-zone for deep convection). Also, this formulation, together with the cold-pool parameterization, significantly improves the organization of convection over tropical regions (Freitas, 2024 JAMES and 2026 WCO5). So, this version is the starting point to retuning MPAS/MONAN 2.0 for the next candidate for regional and global operation on a medium-range time scale. +- Version 2.0.0-rc (Release Candidate) - Includes: MPAS-Model 8.3.1 version, GF convection parametrization improvements, Check control numeric stability, Cumulative surfaces fluxes variables, Spatially heterogenious soil color maps for Noahmp, Soil layer bug fix, NVHPC compilation support. Now, it is default the new land surface (NOAH_MP) and Gravity Wave Drag (UGWP) for the suites 'mesoscale_reference_monan' and 'convection_permitting_monan'. +- Version 1.4.4 - Cold-Start Dynamics Substepping for numerical stabilization and conditional dynamic substepping during the first N full integration steps; option for compilation using a Cray compiler; fixed debug error for compilation using the Intel compiler and compilation flag adjustments for the GNU compiler. - Version 1.4.3-rc (Release Candidate) - This development contains further configuration for the pre-operational MONAN global 10km uniform resolution, adapted to run on the new CRAY supercomputer. It contains new tuning for the GF scheme and sets WSM6 as the default cloud microphysics parameterization. For this configuration, the namelist.atmosphere should be set to config_physics_suite = 'mesoscale_reference_monan'. - Version 1.4.2-rc (Release Candidate) - Speed up by about 7%. Cleanup to become the initial version of the C3P Community Cloud-Convection Parameterization. Joint development between INPE and NOAA/GSL. Removed files not needed anymore. Adding effects of PCW (Neelin et al. 2009) and vertical shear of horizontal wind on the entrainment rate. This should improve model simulations of MCSs. Additional trigger function based on Xie et al 2019. Additional comments and references. - Version 1.4.1-rc (Release Candidate) - Terrain height (ter, calculated previously in the pre processing with init_atmosphere_model) included in the Registry.xml's input section so that can be read from 'init' file and post processed. kubota relhum evalute modification on mpas_isobaric_diagnostics.F. @@ -29,6 +29,8 @@ History - Version 0.2.0 - Included variables and new isobaric levels. - Version 0.1.0 - Initial version structure (0.1.0) based on the dynamic core of the MPAS 8.0.1 Model. +## MPAS + The Model for Prediction Across Scales (MPAS) is a collaborative project for developing atmosphere, ocean, and other earth-system simulation components for use in climate, regional climate, and weather studies. The primary development @@ -44,8 +46,7 @@ high-level data types, communication routines, and I/O routines. By using MPAS, developers can leverage pre-existing code and focus more on development of their model. -BUILDING -======== +## BUILDING This README is provided as a brief introduction to the MPAS framework. It does not provide details about each specific model, nor does it provide building @@ -63,8 +64,7 @@ core's user's guide, which can be found at the following web sites: [MPAS-Seaice](http://mpas-dev.github.io/sea_ice/releases.html) -Code Layout ----------- +## Code Layout Within the MPAS repository, code is laid out as follows. Sub-directories are only described below the src directory. diff --git a/namelist_monan/namelist.atmosphere b/namelist_monan/namelist.atmosphere index 9b786b74d..24c4d7d7f 100644 --- a/namelist_monan/namelist.atmosphere +++ b/namelist_monan/namelist.atmosphere @@ -6,6 +6,8 @@ config_split_dynamics_transport = true config_number_of_sub_steps = 2 config_dynamics_split_steps = 3 + config_coldstart_substeps = 3 + config_coldstart_steps_to_substep = 30 config_h_mom_eddy_visc2 = 0.0 config_h_mom_eddy_visc4 = 0.0 config_v_mom_eddy_visc2 = 0.0 diff --git a/namelist_monan/namelist.atmosphere.TEMPLATE b/namelist_monan/namelist.atmosphere.TEMPLATE index 2a102cb94..c060a3e7b 100755 --- a/namelist_monan/namelist.atmosphere.TEMPLATE +++ b/namelist_monan/namelist.atmosphere.TEMPLATE @@ -6,6 +6,8 @@ config_split_dynamics_transport = true config_number_of_sub_steps = 2 config_dynamics_split_steps = 3 + config_coldstart_substeps = 3 + config_coldstart_steps_to_substep = 30 config_h_mom_eddy_visc2 = 0.0 config_h_mom_eddy_visc4 = 0.0 config_v_mom_eddy_visc2 = 0.0 diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 86e37f24a..e4388f2f6 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -111,6 +111,16 @@ description="When config_split_dynamics_transport = T, the number of RK steps per transport step" possible_values="Positive integer values"/> + + + + + + @@ -591,6 +603,8 @@ + + @@ -704,6 +718,12 @@ + + + + + + @@ -783,6 +803,10 @@ + + + + @@ -883,6 +907,8 @@ + + @@ -1011,6 +1037,7 @@ + @@ -1240,6 +1267,8 @@ + + @@ -2253,7 +2282,7 @@ + + + + + + + + @@ -3073,6 +3119,18 @@ description="latent heat flux at the surface from ocean spray" packages="bl_mynn_in;bl_ysu_in"/> + + + + + + - - - - - - - - - - - - - - - + + + + @@ -3563,6 +3627,12 @@ + + + + @@ -3863,7 +3933,10 @@ + description="dominant soil texture category"/> + + @@ -3871,6 +3944,9 @@ + + @@ -3898,6 +3974,12 @@ + + + + diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 7f9be1196..be881df39 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -409,6 +409,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) integer :: nCells, nEdges, nVertices, nVertLevels integer :: thread character(len=StrKIND), pointer :: mminlu + character(len=StrKIND), pointer :: mminsc integer, pointer :: nThreads integer, dimension(:), pointer :: cellThreadStart, cellThreadEnd @@ -562,6 +563,16 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) write(mminlu,'(a)') 'USGS' end if + ! Before calling physics_init, ensure that mminsc contains the name of the soil colour dataset + call mpas_pool_get_array(sfc_input, 'mminsc', mminsc) + if (len_trim(mminsc) == 0) then + call mpas_log_write('****************************************************************') + call mpas_log_write('No information on soil colour dataset is available.') + call mpas_log_write('Assume that we are using ''DEFAULT_RAD_NOAH''.') + call mpas_log_write('****************************************************************') + write(mminsc,'(a)') 'DEFAULT_RAD_NOAH' + end if + if (moist_physics) then !initialization of some input variables in registry: @@ -572,7 +583,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) call mpas_pool_get_subpool(block % structs, 'atm_input', atm_input) call mpas_pool_get_subpool(block % structs, 'output_noahmp', output_noahmp) call physics_tables_init(dminfo, block % configs) - call physics_registry_init(mesh, block % configs, sfc_input) + call physics_registry_init(mesh, block % dimensions, block % configs, sfc_input) call physics_run_init(block % configs, mesh, state, clock, stream_manager) !initialization of all physics: @@ -1000,6 +1011,17 @@ subroutine atm_do_timestep(domain, dt, itimestep) real (kind=RKIND) :: xtime_s integer :: ierr + ! Cold-start substepping — controlled via namelist (nhyd_model record): + ! config_coldstart_substeps : number of dynamics sub-steps per full step (1 = disabled) + ! config_coldstart_steps_to_substep : number of full steps in which substepping is active (0 = disabled) + ! Physics always uses the full dt; only dynamics is sub-stepped. + integer, pointer :: config_coldstart_substeps + integer, pointer :: config_coldstart_steps_to_substep + integer :: step, integration_steps, itimestep_local + real (kind=RKIND) :: integration_dt + type (mpas_pool_type), pointer :: state_substep + logical :: substep_active + clock => domain % clock mpas_log_info => domain % logInfo @@ -1019,7 +1041,40 @@ subroutine atm_do_timestep(domain, dt, itimestep) endif #endif - call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group) + ! --- Cold-start substepping logic (Skamarock/JCSDA, ported for MONAN 1.4.5-rc) --- + ! Determine whether this full timestep falls within the substepping window + call mpas_pool_get_config(domain % blocklist % configs, 'config_coldstart_substeps', config_coldstart_substeps) + call mpas_pool_get_config(domain % blocklist % configs, 'config_coldstart_steps_to_substep', config_coldstart_steps_to_substep) + + substep_active = .false. + integration_dt = dt + integration_steps = 1 + itimestep_local = itimestep + + if (itimestep .le. config_coldstart_steps_to_substep .and. config_coldstart_substeps .gt. 1) then + substep_active = .true. + integration_dt = dt / real(config_coldstart_substeps, kind=RKIND) + integration_steps = config_coldstart_substeps + itimestep_local = (itimestep - 1) * config_coldstart_substeps + 1 + end if + + do step = 1, integration_steps + + if (substep_active) & + call mpas_log_write(' [cold-start] substep $i of $i (effective dt = $r s)', & + intArgs=(/step, integration_steps/), realArgs=(/integration_dt/)) + + call atm_timestep(domain, integration_dt, currTime, itimestep_local, exchange_halo_group) + + ! Between sub-steps: promote time level 2 → time level 1 for state fields + if (step .lt. integration_steps) then + call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state_substep) + call mpas_pool_shift_time_levels(state_substep) + itimestep_local = itimestep_local + 1 + end if + + end do + ! --- End cold-start substepping --- end subroutine atm_do_timestep diff --git a/src/core_atmosphere/physics/Registry_noahmp.xml b/src/core_atmosphere/physics/Registry_noahmp.xml index 89d980f72..ec4c7b030 100644 --- a/src/core_atmosphere/physics/Registry_noahmp.xml +++ b/src/core_atmosphere/physics/Registry_noahmp.xml @@ -16,7 +16,7 @@ - diff --git a/src/core_atmosphere/physics/mpas_atmphys_control.F b/src/core_atmosphere/physics/mpas_atmphys_control.F index bcdacad74..96d225f3d 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_control.F +++ b/src/core_atmosphere/physics/mpas_atmphys_control.F @@ -157,24 +157,24 @@ subroutine physics_namelist_check(configs) if (trim(config_microp_scheme) == 'suite') config_microp_scheme = 'mp_wsm6' if (trim(config_convection_scheme) == 'suite') config_convection_scheme = 'cu_gf_monan' if (trim(config_pbl_scheme) == 'suite') config_pbl_scheme = 'bl_mynn' - if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ysu_gwdo' + if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ugwp_gwdo' if (trim(config_radt_lw_scheme) == 'suite') config_radt_lw_scheme = 'rrtmg_lw' if (trim(config_radt_sw_scheme) == 'suite') config_radt_sw_scheme = 'rrtmg_sw' if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction_monan' if (trim(config_sfclayer_scheme) == 'suite') config_sfclayer_scheme = 'sf_mynn' - if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'sf_noah' - + if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'sf_noahmp' + else if (trim(config_physics_suite) == 'convection_permitting_monan') then if (trim(config_microp_scheme) == 'suite') config_microp_scheme = 'mp_thompson' if (trim(config_convection_scheme) == 'suite') config_convection_scheme = 'cu_gf_monan' if (trim(config_pbl_scheme) == 'suite') config_pbl_scheme = 'bl_mynn' - if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ysu_gwdo' + if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ugwp_gwdo' if (trim(config_radt_lw_scheme) == 'suite') config_radt_lw_scheme = 'rrtmg_lw' if (trim(config_radt_sw_scheme) == 'suite') config_radt_sw_scheme = 'rrtmg_sw' if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction_monan' if (trim(config_sfclayer_scheme) == 'suite') config_sfclayer_scheme = 'sf_mynn' - if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'sf_noah' + if (trim(config_lsm_scheme) == 'suite') config_lsm_scheme = 'sf_noahmp' else if (trim(config_physics_suite) == 'none') then @@ -358,11 +358,12 @@ subroutine physics_namelist_check(configs) end subroutine physics_namelist_check !================================================================================================================= - subroutine physics_registry_init(mesh,configs,sfc_input) + subroutine physics_registry_init(mesh,dims,configs,sfc_input) !================================================================================================================= !input arguments: type(mpas_pool_type),intent(in):: mesh + type(mpas_pool_type),intent(in):: dims type(mpas_pool_type),intent(in):: configs !inout arguments: @@ -372,40 +373,52 @@ subroutine physics_registry_init(mesh,configs,sfc_input) logical,pointer:: config_do_restart character(len=StrKIND),pointer:: config_lsm_scheme integer,pointer:: nCells + integer, pointer :: nSoilLevels integer,dimension(:),pointer:: landmask - real(kind=RKIND),dimension(:,:),pointer:: dzs + real(kind=RKIND),dimension(:,:), pointer :: dzs + real(kind=RKIND),dimension(:) , pointer :: dzstop, kzs !local variables: integer:: iCell + integer:: iSoil !----------------------------------------------------------------------------------------------------------------- call mpas_pool_get_config(configs,'config_do_restart',config_do_restart) call mpas_pool_get_config(configs,'config_lsm_scheme',config_lsm_scheme) + call mpas_pool_get_dimension(dims, 'nSoilLevels', nSoilLevels) call mpas_pool_get_dimension(mesh,'nCells',nCells) call mpas_pool_get_array(sfc_input,'landmask',landmask) call mpas_pool_get_array(sfc_input,'dzs' , dzs ) + call mpas_pool_get_array(sfc_input,'dzstop' , dzstop ) + call mpas_pool_get_array(sfc_input,'kzs' , kzs ) -!initialization of input variables, if needed: + + +!initialization of input variables, if needed: if(.not. config_do_restart) then lsm_select: select case(trim(config_lsm_scheme)) + case("sf_noah","sf_noahmp") + ! Report soil layers as defined during the initialisation step. + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Soil depth layers (physics_registry_init)') + call mpas_log_write('---~---') + iCell = min(2,nCells) + do iSoil=1, nSoilLevels + call mpas_log_write(' Layer $i - Thickness = $r' & + ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell)/)) + end do + call mpas_log_write('---~---') + call mpas_log_write(' ') + + case default - case("sf_noah","sf_noahmp") - !initialize the thickness of the soil layers for the Noah scheme: - do iCell = 1, nCells - dzs(1,iCell) = 0.10_RKIND - dzs(2,iCell) = 0.30_RKIND - dzs(3,iCell) = 0.60_RKIND - dzs(4,iCell) = 1.00_RKIND - enddo - - case default - end select lsm_select endif diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver.F b/src/core_atmosphere/physics/mpas_atmphys_driver.F index 26eff1738..ef8f34796 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver.F @@ -90,7 +90,8 @@ end subroutine halo_exchange_routine ! driver_sfclayer : driver for surface layer scheme. ! update_convection_step1 : updates lifetime of deep convective clouds in Kain-Fritsch scheme. ! update_convection_step2 : updates accumulated precipitation output from convection schemes. -! update_radiation_diagnostics: updates accumualted radiation diagnostics from radiation schemes. +! update_radiation_diagnostics: updates accumulated radiation diagnostics from radiation schemes. +! update_sfcflux_diagnostics : updates accumulated surface flux diagnostics from LSM/Sfc layer schemes. ! ! add-ons and modifications to sourcecode: ! ---------------------------------------- @@ -165,6 +166,7 @@ subroutine physics_driver(domain,itimestep,xtime_s,exchange_halo_group) logical, pointer:: config_oml1d real(kind=RKIND),pointer:: config_bucket_radt + real(kind=RKIND),pointer:: config_bucket_sflx !local variables: type(block_type),pointer:: block @@ -189,6 +191,7 @@ subroutine physics_driver(domain,itimestep,xtime_s,exchange_halo_group) call mpas_pool_get_config(domain%configs,'config_radt_sw_scheme' ,config_radt_sw_scheme ) call mpas_pool_get_config(domain%configs,'config_sfclayer_scheme' ,config_sfclayer_scheme ) call mpas_pool_get_config(domain%configs,'config_bucket_radt' ,config_bucket_radt ) + call mpas_pool_get_config(domain%configs,'config_bucket_sflx' ,config_bucket_sflx ) call mpas_pool_get_config(domain%configs,'config_bucket_update' ,config_bucket_update ) call mpas_pool_get_config(domain%configs,'config_frac_seaice' ,config_frac_seaice ) call mpas_pool_get_config(domain%configs,'config_oml1d' ,config_oml1d ) @@ -338,6 +341,18 @@ subroutine physics_driver(domain,itimestep,xtime_s,exchange_halo_group) !$OMP END PARALLEL DO endif + + !call to accumulate surface flux diagnostics if needed: + if(config_bucket_update /= 'none' .or. config_bucket_sflx > 0._RKIND) then +!$OMP PARALLEL DO + do thread=1,nThreads + call update_sfcflux_diagnostics(block%configs,mesh,diag_physics, & + cellSolveThreadStart(thread), cellSolveThreadEnd(thread)) + end do +!$OMP END PARALLEL DO + endif + + !call to pbl schemes: if(config_pbl_scheme .ne. 'off' .and. config_sfclayer_scheme .ne. 'off') then call allocate_pbl(block%configs) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F index 0116dcf56..ba5faf286 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F @@ -139,6 +139,7 @@ subroutine allocate_lsm if(.not.allocated(cpm_p) ) allocate(cpm_p(ims:ime,jms:jme) ) if(.not.allocated(cqs2_p) ) allocate(cqs2_p(ims:ime,jms:jme) ) if(.not.allocated(isltyp_p) ) allocate(isltyp_p(ims:ime,jms:jme) ) + if(.not.allocated(isctyp_p) ) allocate(isctyp_p(ims:ime,jms:jme) ) if(.not.allocated(ivgtyp_p) ) allocate(ivgtyp_p(ims:ime,jms:jme) ) if(.not.allocated(glw_p) ) allocate(glw_p(ims:ime,jms:jme) ) if(.not.allocated(grdflx_p) ) allocate(grdflx_p(ims:ime,jms:jme) ) @@ -216,6 +217,7 @@ subroutine deallocate_lsm if(allocated(gsw_p) ) deallocate(gsw_p ) if(allocated(hfx_p) ) deallocate(hfx_p ) if(allocated(isltyp_p) ) deallocate(isltyp_p ) + if(allocated(isctyp_p) ) deallocate(isctyp_p ) if(allocated(ivgtyp_p) ) deallocate(ivgtyp_p ) if(allocated(lai_p) ) deallocate(lai_p ) if(allocated(lh_p) ) deallocate(lh_p ) @@ -279,7 +281,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) character(len=StrKIND),pointer:: config_microp_scheme, & config_convection_scheme - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, & grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc, & @@ -336,6 +338,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) call mpas_pool_get_array(diag_physics,'znt' ,znt ) call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isctyp' ,isctyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'shdmin' ,shdmin ) call mpas_pool_get_array(sfc_input,'shdmax' ,shdmax ) @@ -409,6 +412,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) znt_p(i,j) = znt(i) isltyp_p(i,j) = isltyp(i) + isctyp_p(i,j) = isctyp(i) ivgtyp_p(i,j) = ivgtyp(i) shdmin_p(i,j) = shdmin(i) shdmax_p(i,j) = shdmax(i) @@ -490,7 +494,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) !local pointers: character(len=StrKIND),pointer:: config_microp_scheme - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, & grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc, & @@ -544,6 +548,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) call mpas_pool_get_array(diag_physics,'znt' ,znt ) call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isctyp' ,isctyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'shdmin' ,shdmin ) call mpas_pool_get_array(sfc_input,'shdmax' ,shdmax ) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F index 7b93e7cf6..11b8c4d06 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F @@ -58,7 +58,7 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm integer:: i,its,ite integer:: n,ns,nsoil,nsnow,nzsnow - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer:: latCell,lonCell real(kind=RKIND),dimension(:),pointer:: shdmax,shdmin,vegfra,tmn,xice,xland @@ -842,14 +842,10 @@ subroutine lsm_noahmp_toMPAS(diag_physics,diag_physics_noahmp,output_noahmp,sfc_ snicexy(ns,i) = mpas_noahmp%snicexy(i,n) snliqxy(ns,i) = mpas_noahmp%snliqxy(i,n) enddo - do ns = 1,nsnow + do ns = 1,nzsnow n = ns - nsnow zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n) enddo - do ns = nsnow+1,nzsnow - n = ns - nsoil + 1 - zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n) - enddo enddo diff --git a/src/core_atmosphere/physics/mpas_atmphys_init.F b/src/core_atmosphere/physics/mpas_atmphys_init.F index 4cf54fefb..1aa2569d0 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_init.F +++ b/src/core_atmosphere/physics/mpas_atmphys_init.F @@ -125,11 +125,13 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state i_acswupb,i_acswupbc,i_acswupt,i_acswuptc, & i_aclwdnb,i_aclwdnbc,i_aclwdnt,i_aclwdntc, & i_aclwupb,i_aclwupbc,i_aclwupt,i_aclwuptc + integer,dimension(:),pointer:: i_achfx,i_acqfx,i_aclh,i_acgrdflx,i_acnoahres real(kind=RKIND),dimension(:),pointer:: acswdnb,acswdnbc,acswdnt,acswdntc, & acswupb,acswupbc,acswupt,acswuptc, & aclwdnb,aclwdnbc,aclwdnt,aclwdntc, & aclwupb,aclwupbc,aclwupt,aclwuptc + real(kind=RKIND),dimension(:),pointer :: achfx,acqfx,aclh,acgrdflx,acnoahres real(kind=RKIND),dimension(:),pointer:: nsteps_accum,ndays_accum,tday_accum, & tyear_accum,tyear_mean real(kind=RKIND),dimension(:),pointer:: sst,sstsk,tmn,xice,xicem @@ -191,6 +193,12 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state call mpas_pool_get_array(diag_physics,'i_aclwupt' ,i_aclwupt ) call mpas_pool_get_array(diag_physics,'i_aclwuptc' ,i_aclwuptc ) + call mpas_pool_get_array(diag_physics,'i_achfx' ,i_achfx ) + call mpas_pool_get_array(diag_physics,'i_acqfx' ,i_acqfx ) + call mpas_pool_get_array(diag_physics,'i_aclh' ,i_aclh ) + call mpas_pool_get_array(diag_physics,'i_acgrdflx' ,i_acgrdflx ) + call mpas_pool_get_array(diag_physics,'i_acnoahres' ,i_acnoahres ) + call mpas_pool_get_array(diag_physics,'acswdnb' ,acswdnb ) call mpas_pool_get_array(diag_physics,'acswdnbc' ,acswdnbc ) call mpas_pool_get_array(diag_physics,'acswdnt' ,acswdnt ) @@ -208,6 +216,12 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state call mpas_pool_get_array(diag_physics,'aclwupt' ,aclwupt ) call mpas_pool_get_array(diag_physics,'aclwuptc' ,aclwuptc ) + call mpas_pool_get_array(diag_physics,'achfx' ,achfx ) + call mpas_pool_get_array(diag_physics,'acqfx' ,acqfx ) + call mpas_pool_get_array(diag_physics,'aclh' ,aclh ) + call mpas_pool_get_array(diag_physics,'acgrdflx' ,acgrdflx ) + call mpas_pool_get_array(diag_physics,'acnoahres' ,acnoahres ) + call mpas_pool_get_array(diag_physics,'nsteps_accum',nsteps_accum) call mpas_pool_get_array(diag_physics,'ndays_accum' ,ndays_accum ) call mpas_pool_get_array(diag_physics,'tday_accum' ,tday_accum ) @@ -307,6 +321,25 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state enddo endif +! Initialisation of accumulated surface fluxes and corresponding counters (i_ac*). Counters i_ac* track +! the number of times the absolute value of the accumulated surface fluxes exceed their prescribed +! threshold values. + if(.not. config_do_restart .and. config_lsm_scheme.ne.'off') then + do iCell = 1, nCellsSolve + i_achfx (iCell) = 0 + i_acqfx (iCell) = 0 + i_aclh (iCell) = 0 + i_acgrdflx (iCell) = 0 + i_acnoahres(iCell) = 0 + + achfx (iCell) = 0._RKIND + acqfx (iCell) = 0._RKIND + aclh (iCell) = 0._RKIND + acgrdflx (iCell) = 0._RKIND + acnoahres (iCell) = 0._RKIND + enddo + endif + !initialization of xicem: if(.not.config_do_restart) then ! call mpas_log_write('--- initialization of xicem:') diff --git a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F index d6dc1bc0c..c7a943f14 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F +++ b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F @@ -285,9 +285,14 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs) integer :: iCell,iSoil integer, pointer :: nCellsSolve,nSoilLevels,nFGSoilLevels integer, pointer :: config_nsoillevels + real(kind=RKIND), pointer :: config_dzstop + real(kind=RKIND), pointer :: config_kzs real(kind=RKIND),dimension(:,:),pointer:: dzs_fg,zs_fg real(kind=RKIND),dimension(:,:),pointer:: dzs,zs + real(kind=RKIND),dimension(:),pointer:: dzstop,kzs + logical :: config_fine + logical :: default_layers !----------------------------------------------------------------------------------------------------------------- !call mpas_log_write('') @@ -301,11 +306,39 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs) call mpas_pool_get_array(fg, 'dzs_fg', dzs_fg) call mpas_pool_get_array(fg, 'zs', zs) call mpas_pool_get_array(fg, 'dzs', dzs) + call mpas_pool_get_array(fg, 'dzstop', dzstop) + call mpas_pool_get_array(fg, 'kzs', kzs) call mpas_pool_get_config(configs, 'config_nsoillevels', config_nsoillevels) - if(config_nsoillevels .ne. 4) & - call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.') + + call mpas_pool_get_config(configs, 'config_dzstop', config_dzstop) + call mpas_pool_get_config(configs, 'config_kzs' , config_kzs) + + ! This will check that there are enough soil layers, the thickness of the top soil + ! layer is reasonable and that the exponential factor is positively defined. The only + ! exception to the latter rule is when the number of soil layers is fixed at four. In + ! this case, setting both config_dzstop and config_kzs to zero will enable the default + ! soil layers by Noah and Noah-MP. + config_fine = ( ( config_nsoillevels >= 4 ) .and. ( config_dzstop >= 0.001_RKIND ) & + .and. ( config_kzs > 0.0_RKIND ) ) + default_layers = ( ( config_nsoillevels == 4 ) .and. ( config_dzstop == 0.0_RKIND ) & + .and. ( config_kzs == 0.0_RKIND ) ) + + if (.not. ( config_fine .or. default_layers ) ) then + call mpas_log_write( ' ' ) + call mpas_log_write( '---~---' ) + call mpas_log_write( ' Invalid soil depth settings.' ) + call mpas_log_write( '---~---' ) + call mpas_log_write( ' config_nsoillevels = $i (It must be at least 4).' & + , intArgs=(/config_nsoillevels/)) + call mpas_log_write( ' config_dztop = $r (It must be >= 0.001).' & + , realArgs=(/config_dzstop/)) + call mpas_log_write( ' config_kzs = $r (It must be >= 0.0).' & + , realArgs=(/ config_kzs /)) + call mpas_log_write(' ') + call physics_error_fatal(' At least one of the settings for soil depth is incorrect.') + end if do iCell = 1, nCellsSolve iSoil = 1 @@ -318,10 +351,21 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs) enddo do iCell = 1, nCellsSolve - dzs(1,iCell) = 0.10_RKIND - dzs(2,iCell) = 0.30_RKIND - dzs(3,iCell) = 0.60_RKIND - dzs(4,iCell) = 1.00_RKIND + dzstop(iCell) = config_dzstop + kzs(iCell) = config_kzs + + if (default_layers) then + ! This is the original set up for layers, which will be used by default. + dzs(1,iCell) = 0.10_RKIND + dzs(2,iCell) = 0.30_RKIND + dzs(3,iCell) = 0.60_RKIND + dzs(4,iCell) = 1.00_RKIND + else + ! User-defined soil layers. + do iSoil = 1, nSoilLevels + dzs(iSoil,iCell) = dzstop(iCell) * exp( kzs(iCell) * real(iSoil-1,kind=RKIND) ) + end do + end if iSoil = 1 zs(iSoil,iCell) = 0.5_RKIND * dzs(iSoil,iCell) @@ -333,6 +377,18 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs) enddo + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Soil depth layers (init_soil_layers_depth)') + call mpas_log_write('---~---') + iCell = min(2,nCellsSolve) + do iSoil=1, nSoilLevels + call mpas_log_write(' Layer $i -- Thickness = $r -- MidDepth = $r ' & + ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell)/)) + end do + call mpas_log_write('---~---') + call mpas_log_write(' ') + end subroutine init_soil_layers_depth !================================================================================================================= @@ -351,7 +407,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) type (mpas_pool_type), intent(inout) :: fg !local variables: - integer:: iCell,ifgSoil,iSoil + integer:: iCell,ifgSoil,iSoil,iCellShow integer, pointer:: nCellsSolve,nFGSoilLevels,nSoilLevels integer:: num_sm,num_st integer,dimension(:),pointer:: landmask @@ -409,8 +465,8 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) call mpas_log_write('Error in interpolation of sm_fg to MPAS grid: num_sm = $i', messageType=MPAS_LOG_CRIT, intArgs=(/num_sm/)) endif - if(config_nsoillevels .ne. 4) & - call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.') + if(config_nsoillevels < 4) & + call physics_error_fatal('NOAH and Noah-MP require at least 4 soil layers. Correct config_nsoillevels.') if(.not.allocated(zhave) ) allocate(zhave(nFGSoilLevels+2,nCellsSolve) ) if(.not.allocated(st_input)) allocate(st_input(nFGSoilLevels+2,nCellsSolve)) @@ -429,13 +485,13 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) sm_input(ifgSoil+1,iCell) = sm_fg(ifgSoil,iCell) enddo - zhave(nFGSoilLevels+2,iCell) = 300._RKIND/100._RKIND + zhave(nFGSoilLevels+2,iCell) = max( zs(nSoilLevels,iCell) + dzs(nSoilLevels,iCell), 300._RKIND/100._RKIND) st_input(nFGSoilLevels+2,iCell) = tmn(iCell) sm_input(nFGSoilLevels+2,iCell) = sm_input(nFGSoilLevels,iCell) if(iCell .eq. 1) then do ifgSoil = 1,nFGSoilLevels+2 - call mpas_log_write('$i $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/)) + call mpas_log_write(' ifgSoil = $i -- zhave = $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/)) enddo endif @@ -443,14 +499,15 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) !... interpolate the soil temperature, soil moisture, and soil liquid temperature to the four ! layers used in the NOAH land surface scheme: - + iCellShow = 0 do iCell = 1, nCellsSolve if(landmask(iCell) .eq. 1) then + if (iCellShow == 0) iCellShow = iCell noah: do iSoil = 1 , nSoilLevels input: do ifgSoil = 1 , nFGSoilLevels+2-1 - if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', & + if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', & intArgs=(/iSoil,ifgSoil/), & realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/)) @@ -461,7 +518,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) (st_input(ifgSoil,iCell) * (zhave(ifgSoil+1,iCell)-zs(iSoil,iCell)) & + st_input(ifgSoil+1,iCell) * (zs(iSoil,iCell)-zhave(ifgSoil,iCell))) & / (zhave(ifgSoil+1,iCell)-zhave(ifgSoil,iCell)) - if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', & + if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', & intArgs=(/iSoil,ifgSoil/), & realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/)) @@ -506,6 +563,18 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs) enddo + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Soil initial conditions (init_soil_layers_properties)') + call mpas_log_write('---~---') + iCell = max(1,iCellShow) + do iSoil=1, nSoilLevels + call mpas_log_write(' Layer $i -- Thickness = $r -- midDepth = $r -- Temperature = $r -- Moisture = $r' & + ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell),tslb(iSoil,iCell),smois(iSoil,iCell)/)) + end do + call mpas_log_write('---~---') + call mpas_log_write(' ') + if(allocated(zhave) ) deallocate(zhave ) if(allocated(st_input)) deallocate(st_input) if(allocated(sm_input)) deallocate(sm_input) diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F index b4d7c38ab..f7c28a3c5 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F @@ -91,7 +91,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) character(len=StrKIND),pointer:: mminlu,mminsl integer,pointer:: nCells,nSoilLevels - integer,dimension(:),pointer:: ivgtyp,isltyp + integer,dimension(:),pointer:: ivgtyp,isltyp,isctyp real(kind=RKIND),dimension(:),pointer:: snoalb,snow,snowh real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o @@ -119,6 +119,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) call mpas_pool_get_dimension(mesh,'nSoilLevels',nSoilLevels) call mpas_pool_get_array(sfc_input,'isltyp', isltyp) + call mpas_pool_get_array(sfc_input,'isctyp', isctyp) call mpas_pool_get_array(sfc_input,'ivgtyp', ivgtyp) call mpas_pool_get_array(sfc_input,'sh2o' , sh2o ) call mpas_pool_get_array(sfc_input,'smois' , smois ) @@ -150,6 +151,20 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// & "of ISLTYP. Is this field in the input?" ) + ! Make sure all cells have reasonable soil colour indices. + errflag = 0 + do iCell = 1, nCells + if(isctyp(iCell) < 1) then + errflag = 1 + write(err_message,*) "module_sf_noahlsm.F: lsminit: out of range ISCTYP ", & + iCell,isctyp(iCell),isltyp(iCell),ivgtyp(iCell) + call physics_message(err_message) + endif + end do + if (errflag == 1) & + call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// & + "of ISCTYP. Is this field in the input?" ) + !initializes soil liquid water content SH2O: do iCell = 1, nCells bx = bb(isltyp(iCell)) diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F index 19ad287d0..05c6c8b65 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F @@ -46,6 +46,7 @@ subroutine init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,o !--- local variables and arrays: character(len=StrKIND),pointer:: mminlu + character(len=StrKIND),pointer:: mminsc integer:: ns @@ -73,6 +74,10 @@ subroutine init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,o call mpas_pool_get_array(sfc_input,'mminlu',mminlu) mpas_noahmp%llanduse = mminlu +!--- read NoahmpTable.TBL: + call mpas_pool_get_array(sfc_input,'mminsc',mminsc) + mpas_noahmp%lsoilcol = mminsc + !call mpas_log_write(' ') !call mpas_log_write('--- enter subroutine NoahmpReadTable:') call NoahmpReadTable(mpas_noahmp) @@ -234,7 +239,7 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu logical,pointer:: urban_physics integer,pointer:: nsoilcomps - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp integer,dimension(:),pointer:: isnowxy integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi @@ -322,11 +327,14 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu !--- initialization of time-invariant surface variables needed in subroutine NoahmpInitMain: call mpas_pool_get_array(sfc_input,'dzs' ,dzs ) call mpas_pool_get_array(sfc_input,'isltyp',isltyp) + call mpas_pool_get_array(sfc_input,'isctyp',isctyp) call mpas_pool_get_array(sfc_input,'ivgtyp',ivgtyp) do i = its, ite mpas_noahmp%isltyp(i) = isltyp(i) + mpas_noahmp%isctyp(i) = isctyp(i) mpas_noahmp%ivgtyp(i) = ivgtyp(i) + enddo do ns = 1, nsoil mpas_noahmp%dzs(ns) = dzs(ns,its) diff --git a/src/core_atmosphere/physics/mpas_atmphys_manager.F b/src/core_atmosphere/physics/mpas_atmphys_manager.F index 205714872..05415a573 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_manager.F +++ b/src/core_atmosphere/physics/mpas_atmphys_manager.F @@ -63,6 +63,11 @@ module mpas_atmphys_manager character(len=*), parameter:: acradtAlarmID = 'acradt' type(MPAS_TimeInterval_Type):: acradtTimeStep +!defines alarm to check if the accumulated surface flux diagnostics are greater than +!their maximum allowed value: + character(len=*), parameter:: acsflxAlarmID = 'acsflx' + type(MPAS_TimeInterval_Type):: acsflxTimeStep + !defines alarm to compute some physics diagnostics, such as radar reflectivity: character(len=*), parameter:: diagAlarmID = 'diag' @@ -156,7 +161,8 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) character(len=StrKIND),pointer:: config_convection_scheme, & config_radt_lw_scheme, & - config_radt_sw_scheme + config_radt_sw_scheme, & + config_lsm_scheme character(len=StrKIND),pointer:: config_conv_interval, & config_radtlw_interval, & @@ -187,6 +193,7 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) call mpas_pool_get_config(domain%blocklist%configs,'config_convection_scheme',config_convection_scheme) call mpas_pool_get_config(domain%blocklist%configs,'config_radt_lw_scheme' ,config_radt_lw_scheme ) call mpas_pool_get_config(domain%blocklist%configs,'config_radt_sw_scheme' ,config_radt_sw_scheme ) + call mpas_pool_get_config(domain%blocklist%configs,'config_lsm_scheme' ,config_lsm_scheme ) call mpas_pool_get_config(domain%blocklist%configs,'config_conv_interval' ,config_conv_interval ) call mpas_pool_get_config(domain%blocklist%configs,'config_radtlw_interval',config_radtlw_interval) @@ -361,6 +368,16 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) call mpas_log_write('--- time to apply limit to accumulated radiation diags. L_ACRADT = $l',logicArgs=(/l_acradt/)) endif +!check to see if it is time to apply limit to the surface flux diagnostics: + if(trim(config_lsm_scheme) /= "off") then + l_acsflx = .false. + if(mpas_is_alarm_ringing(clock,acsflxAlarmID,acsflxTimeStep,ierr=ierr)) then + call mpas_reset_clock_alarm(clock,acsflxAlarmID,acsflxTimeStep,ierr=ierr) + l_acsflx = .true. + endif + call mpas_log_write('--- time to apply limit to accumulated surface flux diags. L_ACSFLX =$l',logicArgs=(/l_acsflx/)) + endif + !check to see if it is time to calculate additional physics diagnostics: l_diags = .false. if (mpas_is_alarm_ringing(clock,diagAlarmID,ierr=ierr)) then @@ -612,6 +629,17 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call physics_error_fatal('subroutine physics_init: error creating alarm radiation limit') endif +!set alarm to check if the accumulated surface flux diagnostics is greater than its maximum +!allowed value: + if(config_bucket_update /= "none") then + call mpas_set_timeInterval(acsflxTimeStep,dt=dt,ierr=ierr) + call mpas_set_timeInterval(alarmTimeStep,timeString=config_bucket_update,ierr=ierr) + alarmStartTime = startTime + alarmTimeStep + call mpas_add_clock_alarm(clock,acsflxAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) + if(ierr /= 0) & + call physics_error_fatal('subroutine physics_init: error creating alarm surface-flux limit') + endif + !set alarm to calculate physics diagnostics on IO outpt only: call MPAS_stream_mgr_get_property(stream_manager, 'output', MPAS_STREAM_PROPERTY_RECORD_INTV, stream_interval, & direction=MPAS_STREAM_OUTPUT, ierr=ierr) @@ -708,6 +736,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) l_camlw = .false. l_acrain = .false. l_acradt = .false. + l_acsflx = .false. !initialization for CAM radiation schemes only: if(trim(config_radt_lw_scheme) .eq. "cam_lw" .or. & diff --git a/src/core_atmosphere/physics/mpas_atmphys_update.F b/src/core_atmosphere/physics/mpas_atmphys_update.F index 421a90e28..f05f6cca4 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_update.F +++ b/src/core_atmosphere/physics/mpas_atmphys_update.F @@ -16,7 +16,8 @@ module mpas_atmphys_update implicit none private public:: physics_update, & - update_radiation_diagnostics + update_radiation_diagnostics, & + update_sfcflux_diagnostics !Update diagnostics. @@ -27,6 +28,7 @@ module mpas_atmphys_update ! ----------------------------------- ! physics_update : not used. ! update_radiation_diagnostics: update accumulated radiation diagnostics. +! update_sfcflux_diagnostics: update accumulated surface flux diagnostics. ! ! add-ons and modifications to sourcecode: ! ---------------------------------------- @@ -253,6 +255,99 @@ subroutine update_radiation_diagnostics(configs,mesh,diag_physics,its,ite) end subroutine update_radiation_diagnostics +!================================================================================================================= + subroutine update_sfcflux_diagnostics(configs,mesh,diag_physics,its,ite) +!================================================================================================================= + +!input arguments: + type(mpas_pool_type),intent(in):: configs + type(mpas_pool_type),intent(in):: mesh + integer,intent(in):: its,ite + +!inout arguments: + type(mpas_pool_type),intent(inout):: diag_physics + +!local pointers: + integer,pointer:: nCellsSolve + integer,dimension(:),pointer :: i_achfx, i_acqfx, i_aclh, i_acgrdflx, i_acnoahres + + real(kind=RKIND),pointer:: bucket_sflx + real(kind=RKIND),dimension(:),pointer:: hfx , qfx , lh , grdflx , noahres + real(kind=RKIND),dimension(:),pointer:: achfx, acqfx, aclh, acgrdflx, acnoahres + +!local variables and arrays: + integer:: iCell + +!----------------------------------------------------------------------------------------------------------------- + + call mpas_pool_get_config(configs,'config_bucket_sflx',bucket_sflx) + + call mpas_pool_get_dimension(mesh,'nCellsSolve',nCellsSolve) + + call mpas_pool_get_array(diag_physics,'i_achfx' , i_achfx ) + call mpas_pool_get_array(diag_physics,'i_acqfx' , i_acqfx ) + call mpas_pool_get_array(diag_physics,'i_aclh' , i_aclh ) + call mpas_pool_get_array(diag_physics,'i_acgrdflx' , i_acgrdflx ) + call mpas_pool_get_array(diag_physics,'i_acnoahres' , i_acnoahres ) + + call mpas_pool_get_array(diag_physics,'achfx' , achfx ) + call mpas_pool_get_array(diag_physics,'acqfx' , acqfx ) + call mpas_pool_get_array(diag_physics,'aclh' , aclh ) + call mpas_pool_get_array(diag_physics,'acgrdflx' , acgrdflx ) + call mpas_pool_get_array(diag_physics,'acnoahres' , acnoahres ) + + call mpas_pool_get_array(diag_physics,'hfx' , hfx ) + call mpas_pool_get_array(diag_physics,'qfx' , qfx ) + call mpas_pool_get_array(diag_physics,'lh' , lh ) + call mpas_pool_get_array(diag_physics,'grdflx' , grdflx ) + call mpas_pool_get_array(diag_physics,'noahres' , noahres ) + + do iCell = its, ite + achfx (iCell) = achfx (iCell) + hfx (ICell) * dt_dyn + acqfx (iCell) = acqfx (iCell) + qfx (ICell) * dt_dyn + aclh (iCell) = aclh (iCell) + lh (ICell) * dt_dyn + acgrdflx (iCell) = acgrdflx (iCell) + grdflx (ICell) * dt_dyn + acnoahres(iCell) = acnoahres(iCell) + noahres(ICell) * dt_dyn + end do + + if (l_acsflx .and. bucket_sflx > 0._RKIND) then + do iCell = its, ite + + ! For the time being we apply the bucket logic to the surface fluxes. However, + ! these fluxes can be either positive or negative, so we must apply a two-sided + ! bucket to avoid positive and negative overflow. + if(abs(achfx(iCell)) > bucket_sflx) then + i_achfx(iCell) = i_achfx(iCell) + 1 + achfx (iCell) = sign(achfx(iCell),1.) * (abs(achfx(iCell) - bucket_sflx)) + end if + + if(abs(acqfx(iCell)) > bucket_sflx) then + i_acqfx(iCell) = i_acqfx(iCell) + 1 + acqfx (iCell) = sign(acqfx(iCell),1.) * (abs(acqfx(iCell) - bucket_sflx)) + end if + + if(abs(aclh (iCell)) > bucket_sflx) then + i_aclh(iCell) = i_aclh(iCell) + 1 + aclh (iCell) = sign(aclh(iCell),1.) * (abs(aclh(iCell) - bucket_sflx)) + end if + + if(abs(acgrdflx(iCell)) > bucket_sflx) then + i_acgrdflx(iCell) = i_acgrdflx(iCell) + 1 + acgrdflx (iCell) = & + sign(acgrdflx(iCell),1.) * (abs(acgrdflx(iCell) - bucket_sflx)) + end if + + if(abs(acnoahres(iCell)) > bucket_sflx) then + i_acnoahres(iCell) = i_acnoahres(iCell) + 1 + acnoahres (iCell) = & + sign(acnoahres(iCell),1.) * (abs(acnoahres(iCell) - bucket_sflx)) + end if + end do + + end if + + end subroutine update_sfcflux_diagnostics + !================================================================================================================= end module mpas_atmphys_update !================================================================================================================= diff --git a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F index 6e2057d5c..4d6a59e2d 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F +++ b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F @@ -151,7 +151,7 @@ subroutine physics_update_sst(dminfo,config_frac_seaice,mesh,sfc_input,diag_phys call mpas_pool_get_array(sfc_input,'isice' ,isice ) call mpas_pool_get_array(sfc_input,'iswater' ,iswater ) - call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'landmask' ,landmask ) call mpas_pool_get_array(sfc_input,'vegfra' ,vegfra ) diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index 3b997ccc8..fcba4ef15 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -148,6 +148,7 @@ module mpas_atmphys_vars logical:: l_diags !controls when to calculate physics diagnostics. logical:: l_acrain !when .true., limit to accumulated rain is applied. logical:: l_acradt !when .true., limit to lw and sw radiation is applied. + logical:: l_acsflx !when .true., limit to surface fluxes is applied. logical:: l_mp_tables !when .true., read look-up tables for Thompson cloud microphysics scheme. integer,public:: ids,ide,jds,jde,kds,kde @@ -894,6 +895,7 @@ module mpas_atmphys_vars integer,dimension(:,:),allocatable:: & isltyp_p, &!dominant soil type category [-] + isctyp_p, &!dominant soil colour category [-] ivgtyp_p !dominant vegetation category [-] real(kind=RKIND),dimension(:),allocatable:: & diff --git a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md index ae45c39f5..6427e6b8f 100644 --- a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md +++ b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md @@ -400,7 +400,7 @@ - but in the file listed in the namelist as: HRLDAS_SETUP_FILE = " - The initialization fields are: SNOW,CANWAT,TSK,TSLB,SMOIS - - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,DZS,ZS + - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,ISCTYP,DZS,ZS - This file can also contains some optional fields: LAI - NOTE: a WRF input file can be used as a HRLDAS_SETUP_FILE diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 index 2de35ed9c..2acfd6e01 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 @@ -62,7 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) ! config domain variable noahmp%config%domain%SurfaceType = 1 noahmp%config%domain%NumSwRadBand = 2 - noahmp%config%domain%SoilColor = 4 + noahmp%config%domain%SoilColor = NoahmpIO%ISCTYP(I) noahmp%config%domain%NumCropGrowStage = 8 noahmp%config%domain%FlagSoilProcess = NoahmpIO%calculate_soil noahmp%config%domain%NumSoilTimeStep = NoahmpIO%soil_update_steps @@ -75,6 +75,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) noahmp%config%domain%SoilTimeStep = NoahmpIO%DTBL * NoahmpIO%soil_update_steps noahmp%config%domain%GridSize = NoahmpIO%DX noahmp%config%domain%LandUseDataName = NoahmpIO%LLANDUSE + noahmp%config%domain%SoilColDataName = NoahmpIO%LSOILCOL noahmp%config%domain%VegType = NoahmpIO%IVGTYP(I) noahmp%config%domain%CropType = NoahmpIO%CROPCAT(I) noahmp%config%domain%IndicatorIceSfc = NoahmpIO%ICE @@ -114,10 +115,10 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) if ( noahmp%config%nmlist%OptSoilProperty == 1 ) then noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! soil type same in all layers elseif ( noahmp%config%nmlist%OptSoilProperty == 2 ) then - noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1 - noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2 - noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3 - noahmp%config%domain%SoilType(4) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layer4 + noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1 + noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2 + noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3 + noahmp%config%domain%SoilType(4:NumSoilLayer) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layers 4 and below elseif ( noahmp%config%nmlist%OptSoilProperty == 3 ) then noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! to initialize with default endif diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 index 66458d638..49a709f12 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 @@ -37,6 +37,7 @@ subroutine NoahmpIOVarFinalizeDefault(NoahmpIO) if ( allocated (NoahmpIO%zsoil) ) deallocate ( NoahmpIO%zsoil ) ! depth to soil interfaces [m] if ( allocated (NoahmpIO%ivgtyp) ) deallocate ( NoahmpIO%ivgtyp ) ! vegetation type if ( allocated (NoahmpIO%isltyp) ) deallocate ( NoahmpIO%isltyp ) ! soil type + if ( allocated (NoahmpIO%isctyp) ) deallocate ( NoahmpIO%isctyp ) ! soil colour class if ( allocated (NoahmpIO%vegfra) ) deallocate ( NoahmpIO%vegfra ) ! vegetation fraction [] if ( allocated (NoahmpIO%tmn) ) deallocate ( NoahmpIO%tmn ) ! deep soil temperature [K] if ( allocated (NoahmpIO%xland) ) deallocate ( NoahmpIO%xland ) ! =2 ocean; =1 land/seaice diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 index 4f3c3f4f2..cfd91bb91 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 @@ -41,6 +41,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) if ( .not. allocated (NoahmpIO%zsoil) ) allocate ( NoahmpIO%zsoil (1:nsoil ) ) ! depth to soil interfaces [m] if ( .not. allocated (NoahmpIO%ivgtyp) ) allocate ( NoahmpIO%ivgtyp (its:ite ) ) ! vegetation type if ( .not. allocated (NoahmpIO%isltyp) ) allocate ( NoahmpIO%isltyp (its:ite ) ) ! soil type + if ( .not. allocated (NoahmpIO%isctyp) ) allocate ( NoahmpIO%isctyp (its:ite ) ) ! soil colour class if ( .not. allocated (NoahmpIO%vegfra) ) allocate ( NoahmpIO%vegfra (its:ite ) ) ! vegetation fraction [] if ( .not. allocated (NoahmpIO%tmn) ) allocate ( NoahmpIO%tmn (its:ite ) ) ! deep soil temperature [K] if ( .not. allocated (NoahmpIO%xland) ) allocate ( NoahmpIO%xland (its:ite ) ) ! =2 ocean; =1 land/seaice @@ -469,6 +470,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) NoahmpIO%ice = undefined_int NoahmpIO%ivgtyp = undefined_int NoahmpIO%isltyp = undefined_int + NoahmpIO%isctyp = undefined_int NoahmpIO%isnowxy = undefined_int NoahmpIO%coszen = undefined_real NoahmpIO%xlat = undefined_real diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 index 05a29d703..2abf99670 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 @@ -68,6 +68,7 @@ module NoahmpIOVarType integer :: soil_update_steps ! number of model time steps to update soil process integer, allocatable, dimension(:) :: ivgtyp ! vegetation type integer, allocatable, dimension(:) :: isltyp ! soil type + integer, allocatable, dimension(:) :: isctyp ! soil colour class real(kind=kind_noahmp), allocatable, dimension(:) :: coszen ! cosine zenith angle real(kind=kind_noahmp), allocatable, dimension(:) :: xlat ! latitude [rad] real(kind=kind_noahmp), allocatable, dimension(:,:) :: dz8w ! thickness of atmo layers [m] @@ -516,6 +517,7 @@ module NoahmpIOVarType CHARACTER(LEN=256) :: mminsl = 'STAS' ! soil classification CHARACTER(LEN=256) :: llanduse ! (=USGS, using USGS landuse classification) + CHARACTER(LEN=256) :: lsoilcol ! (=Noah, using Noah's native soil colour classification) !------------------------------------------------------------------------ ! Timing: diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 index eb01ceb2f..949ea8c90 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 @@ -28,7 +28,7 @@ subroutine NoahmpReadTable(NoahmpIO) !------------------------------------------------------- integer, parameter :: MVT = 27 ! number of vegetation types integer, parameter :: MBAND = 2 ! number of radiation bands - integer, parameter :: MSC = 8 ! number of soil texture + integer, parameter :: MSC = 21 ! number of soil colour types integer, parameter :: MAX_SOILTYP = 30 ! max number of soil types integer, parameter :: NCROP = 5 ! number of crop types integer, parameter :: NSTAGE = 8 ! number of crop growth stages @@ -98,12 +98,18 @@ subroutine NoahmpReadTable(NoahmpIO) CZIL_DATA ! radiation parameters + character(len=256) :: RAD_DATASET_DESCRIPTION + integer :: NSC real(kind=kind_noahmp) :: BETADS, BETAIS, EICE real(kind=kind_noahmp), dimension(MBAND) :: ALBICE, ALBLAK, OMEGAS real(kind=kind_noahmp), dimension(2) :: EG real(kind=kind_noahmp), dimension(MSC) :: ALBSAT_VIS, ALBSAT_NIR, ALBDRY_VIS, ALBDRY_NIR + namelist / noahmp_rad_categories / RAD_DATASET_DESCRIPTION, NSC namelist / noahmp_rad_parameters / ALBSAT_VIS, ALBSAT_NIR, ALBDRY_VIS, ALBDRY_NIR, ALBICE, ALBLAK, OMEGAS, & BETADS, BETAIS, EG, EICE + namelist / noahmp_clm_rad_categories / RAD_DATASET_DESCRIPTION, NSC + namelist / noahmp_clm_rad_parameters / ALBSAT_VIS, ALBSAT_NIR, ALBDRY_VIS, ALBDRY_NIR, ALBICE, ALBLAK, OMEGAS, & + BETADS, BETAIS, EG, EICE ! global parameters real(kind=kind_noahmp) :: CO2, O2, TIMEAN, FSATMX, Z0SNO, SSI, SNOW_RET_FAC ,SNOW_EMIS, SWEMX, TAU0, & @@ -847,21 +853,54 @@ subroutine NoahmpReadTable(NoahmpIO) if (ierr /= 0) then write(*,'("WARNING: Cannot find file NoahmpTable.TBL")') endif - read(15,noahmp_rad_parameters) + + + + DATASET_IDENTIFIER = NoahmpIO%LSOILCOL + + inquire( file='NoahmpTable.TBL', exist=file_named ) + if ( file_named ) then + open(15, file="NoahmpTable.TBL", status='old', form='formatted', action='read', iostat=ierr) + else + open(15, status='old', form='formatted', action='read', iostat=ierr) + end if + if ( ierr /= 0 ) then + write(*,'("WARNING: Cannot find file NoahmpTable.TBL")') + endif + + select case (trim(DATASET_IDENTIFIER)) + case ("DEFAULT_RAD_NOAH") + read(15,noahmp_rad_categories,iostat=ierr) + if ( ierr /= 0 ) then + write(*,'("WARNING: Namelist ''noahmp_rad_categories'' not found; using the default values.")') + write(*,'("WARNING: Update your NoahmpTable.TBL. This will eventually become an error.")') + + RAD_DATASET_DESCRIPTION = "DEFAULT_RAD_NOAH" ! radiation (soil colour) type dataset + NSC = 8 ! total number of soil colour categories in Noah + rewind(15) + end if + read(15,noahmp_rad_parameters) + case ("MODIFIED_RAD_CLM_NOAH") + read(15,noahmp_clm_rad_categories) + read(15,noahmp_clm_rad_parameters) + case default + write(*,'("WARNING: Unrecognised DATASET_IDENTIFIER in subroutine ReadNoahmpTable")') + write(*,'("WARNING: DATASET_IDENTIFIER = ''", A, "''")') trim(DATASET_IDENTIFIER) + end select close(15) ! assign values - NoahmpIO%ALBSAT_TABLE(:,1) = ALBSAT_VIS ! saturated soil albedos: 1=vis, 2=nir - NoahmpIO%ALBSAT_TABLE(:,2) = ALBSAT_NIR ! saturated soil albedos: 1=vis, 2=nir - NoahmpIO%ALBDRY_TABLE(:,1) = ALBDRY_VIS ! dry soil albedos: 1=vis, 2=nir - NoahmpIO%ALBDRY_TABLE(:,2) = ALBDRY_NIR ! dry soil albedos: 1=vis, 2=nir - NoahmpIO%ALBICE_TABLE = ALBICE - NoahmpIO%ALBLAK_TABLE = ALBLAK - NoahmpIO%OMEGAS_TABLE = OMEGAS - NoahmpIO%BETADS_TABLE = BETADS - NoahmpIO%BETAIS_TABLE = BETAIS - NoahmpIO%EG_TABLE = EG - NoahmpIO%EICE_TABLE = EICE + NoahmpIO%ALBSAT_TABLE(1:NSC,1) = ALBSAT_VIS(1:NSC) ! saturated soil albedos: 1=vis, 2=nir + NoahmpIO%ALBSAT_TABLE(1:NSC,2) = ALBSAT_NIR(1:NSC) ! saturated soil albedos: 1=vis, 2=nir + NoahmpIO%ALBDRY_TABLE(1:NSC,1) = ALBDRY_VIS(1:NSC) ! dry soil albedos: 1=vis, 2=nir + NoahmpIO%ALBDRY_TABLE(1:NSC,2) = ALBDRY_NIR(1:NSC) ! dry soil albedos: 1=vis, 2=nir + NoahmpIO%ALBICE_TABLE = ALBICE + NoahmpIO%ALBLAK_TABLE = ALBLAK + NoahmpIO%OMEGAS_TABLE = OMEGAS + NoahmpIO%BETADS_TABLE = BETADS + NoahmpIO%BETAIS_TABLE = BETAIS + NoahmpIO%EG_TABLE = EG + NoahmpIO%EICE_TABLE = EICE !---------------- NoahmpTable.TBL global parameters inquire( file='NoahmpTable.TBL', exist=file_named ) diff --git a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL index c9d37c5b4..e25c233ff 100644 --- a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL +++ b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL @@ -404,6 +404,11 @@ LAI_DEC = 4.0, 4.5, 0.0, 0.0, 2.0, 0.0, 0.0, 0.2, 0.3, 0.4, 0.2, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 1.0, 0.6, 0.0, / +&noahmp_rad_categories + RAD_DATASET_DESCRIPTION = "DEFAULT_RAD_NOAH" ! radiation (soil colour) type dataset + NSC = 8 ! total number of soil colour categories in Noah +/ + &noahmp_rad_parameters !------------------------------------------------------------------------------ ! soil color: 1 2 3 4 5 6 7 8 soil color index for soil albedo @@ -421,6 +426,34 @@ EICE = 0.98 ! emissivity ice surface / +&noahmp_clm_rad_categories + RAD_DATASET_DESCRIPTION = "MODIFIED_RAD_CLM_NOAH" ! radiation (soil colour) type dataset + NSC = 21 ! total number of soil colour categories in CLM +/ + +&noahmp_clm_rad_parameters + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + ! soil color index for soil albedo + ! soil color: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + ! saturated soil albedo at visible band + ALBSAT_VIS = 0.25, 0.23, 0.21, 0.20, 0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.10, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.00 + ! saturated soil albedo at NIR band + ALBSAT_NIR = 0.50, 0.46, 0.42, 0.40, 0.38, 0.36, 0.34, 0.32, 0.30, 0.28, 0.26, 0.24, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.08, 0.00 + ! dry soil albedo at visible band + ALBDRY_VIS = 0.36, 0.34, 0.32, 0.31, 0.30, 0.29, 0.28, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.08, 0.00 + ! dry soil albedo at NIR band + ALBDRY_NIR = 0.61, 0.57, 0.53, 0.51, 0.49, 0.48, 0.45, 0.43, 0.41, 0.39, 0.37, 0.35, 0.33, 0.31, 0.29, 0.27, 0.25, 0.23, 0.21, 0.16, 0.00 + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + ALBICE = 0.80, 0.55 ! albedo land ice: 1=vis, 2=nir + ALBLAK = 0.60, 0.40 ! albedo frozen lakes: 1=vis, 2=nir + OMEGAS = 0.8 , 0.4 ! two-stream parameter omega for snow + BETADS = 0.5 ! two-stream parameter betad for snow + BETAIS = 0.5 ! two-stream parameter betaI for snow + EG = 0.97, 0.98 ! emissivity soil surface 1-soil;2-lake + EICE = 0.98 ! emissivity ice surface +/ + &noahmp_global_parameters ! atmospheric constituants CO2 = 395.0e-06 ! CO2 partial pressure diff --git a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90 index 7b5e83913..a082d9834 100644 --- a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90 @@ -5,6 +5,8 @@ module BalanceErrorCheckGlacierMod use Machine use NoahmpVarType use ConstantDefineMod + use NoahmpFatalErrorMod,only: Noahmp_error_fatal + use mpas_log implicit none @@ -77,19 +79,22 @@ subroutine BalanceWaterCheckGlacier(noahmp) #ifndef WRF_HYDRO if ( abs(WaterBalanceError) > 0.1 ) then - if ( WaterBalanceError > 0) then - write(*,*) "The model is gaining water (WaterBalanceError is positive)" - else - write(*,*) "The model is losing water (WaterBalanceError is negative)" - endif - write(*,*) "WaterBalanceError = ",WaterBalanceError, "kg m{-2} timestep{-1}" - write(*, & - '(" GridIndexI GridIndexJ WaterStorageTotEnd WaterStorageTotBeg PrecipTotRefHeight & - EvapGroundNet RunoffSurface RunoffSubsurface")') - write(*,'(i6,1x,i6,1x,2f15.3,9f11.5)') GridIndexI, GridIndexJ, WaterStorageTotEnd, WaterStorageTotBeg, & - PrecipTotRefHeight*MainTimeStep, EvapGroundNet*MainTimeStep, & - RunoffSurface*MainTimeStep, RunoffSubsurface*MainTimeStep - stop "Error: Water budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP water budget conservation error (glacier):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' WaterStorageTotBeg = $r ', realArgs = (/ WaterStorageTotBeg /) ) + call mpas_log_write(' WaterStorageTotEnd = $r ', realArgs = (/ WaterStorageTotEnd /) ) + call mpas_log_write(' WaterBalanceError = $r ', realArgs = (/ WaterBalanceError /) ) + call mpas_log_write(' (positive value above means water gain).' ) + call mpas_log_write(' PrecipTotRefHeight = $r ', realArgs = (/ PrecipTotRefHeight*MainTimeStep /) ) + call mpas_log_write(' EvapGroundNet = $r ', realArgs = (/ EvapGroundNet*MainTimeStep /) ) + call mpas_log_write(' RunoffSurface = $r ', realArgs = (/ RunoffSurface*MainTimeStep /) ) + call mpas_log_write(' RunoffSubsurface = $r ', realArgs = (/ RunoffSubsurface*MainTimeStep /) ) + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Water budget problem in NoahMP LSM (glacier)") endif #endif @@ -133,27 +138,45 @@ subroutine BalanceEnergyCheckGlacier(noahmp) RadSwBalanceError = RadSwDownRefHeight - (RadSwAbsSfc + RadSwReflSfc) ! print out diagnostics when error is large if ( abs(RadSwBalanceError) > 0.01 ) then - write(*,*) "GridIndexI, GridIndexJ = ", GridIndexI, GridIndexJ - write(*,*) "RadSwBalanceError = ", RadSwBalanceError - write(*,*) "RadSwDownRefHeight = ", RadSwDownRefHeight - write(*,*) "RadSwReflSfc = ", RadSwReflSfc - write(*,*) "RadSwAbsGrd = ", RadSwAbsGrd - write(*,*) "RadSwAbsSfc = ", RadSwAbsSfc - stop "Error: Solar radiation budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP solar radiation budget conservation error (glacier):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) ) + call mpas_log_write(' (positive value above means energy gain).' ) + call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) ) + call mpas_log_write(' RadSwDownRefHeight = $r ', realArgs = (/ RadSwDownRefHeight /) ) + call mpas_log_write(' RadSwReflSfc = $r ', realArgs = (/ RadSwReflSfc /) ) + call mpas_log_write(' RadSwAbsGrd = $r ', realArgs = (/ RadSwAbsGrd /) ) + call mpas_log_write(' RadSwAbsSfc = $r ', realArgs = (/ RadSwAbsSfc /) ) + call mpas_log_write('---~---') + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Solar radiation budget problem in NoahMP LSM (glacier)") endif ! error in surface energy balance should be <0.01 W/m2 EnergyBalanceError = RadSwAbsGrd + HeatPrecipAdvSfc - (RadLwNetSfc + HeatSensibleSfc + HeatLatentGrd + HeatGroundTot) ! print out diagnostics when error is large if ( abs(EnergyBalanceError) > 0.01 ) then - write(*,*) 'EnergyBalanceError = ', EnergyBalanceError, ' at GridIndexI,GridIndexJ: ', GridIndexI, GridIndexJ - write(*,'(a17,F10.4)' ) "Net longwave: ", RadLwNetSfc - write(*,'(a17,F10.4)' ) "Total sensible: ", HeatSensibleSfc - write(*,'(a17,F10.4)' ) "Ground evap: ", HeatLatentGrd - write(*,'(a17,F10.4)' ) "Total ground: ", HeatGroundTot - write(*,'(a17,4F10.4)') "Precip advected: ", HeatPrecipAdvSfc - write(*,'(a17,F10.4)' ) "absorbed shortwave: ", RadSwAbsGrd - stop "Error: Surface energy budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP energy budget conservation error (glacier):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' EnergyBalanceError = $r ', realArgs = (/ EnergyBalanceError /) ) + call mpas_log_write(' (positive value above energy gain).' ) + call mpas_log_write(' Net longwave = $r ', realArgs = (/ RadLwNetSfc /) ) + call mpas_log_write(' Total sensible = $r ', realArgs = (/ HeatSensibleSfc /) ) + call mpas_log_write(' Ground evap = $r ', realArgs = (/ HeatLatentGrd /) ) + call mpas_log_write(' Total ground = $r ', realArgs = (/ HeatGroundTot /) ) + call mpas_log_write(' Precip advected = $r ', realArgs = (/ HeatPrecipAdvSfc /) ) + call mpas_log_write(' Absorbed shortwave = $r ', realArgs = (/ RadSwAbsGrd /) ) + call mpas_log_write('---~---') + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Energy budget problem in NoahMP LSM (glacier)") endif end associate diff --git a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90 index f076e2a5e..ce58f61d4 100644 --- a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90 @@ -5,6 +5,8 @@ module BalanceErrorCheckMod use Machine use NoahmpVarType use ConstantDefineMod + use NoahmpFatalErrorMod,only: Noahmp_error_fatal + use mpas_log implicit none @@ -68,6 +70,7 @@ subroutine BalanceWaterCheck(noahmp) ! local variable integer :: LoopInd ! loop index + real(kind=kind_noahmp) :: SoilWaterStorage ! -------------------------------------------------------------------- associate( & @@ -116,10 +119,12 @@ subroutine BalanceWaterCheck(noahmp) ! only water balance check for every soil timestep ! Error in water balance should be < 0.1 mm if ( SurfaceType == 1 ) then ! soil - WaterStorageTotEnd = CanopyLiqWater + CanopyIce + SnowWaterEquiv + WaterStorageAquifer + SoilWaterStorage = 0.0 do LoopInd = 1, NumSoilLayer - WaterStorageTotEnd = WaterStorageTotEnd + SoilMoisture(LoopInd) * ThicknessSnowSoilLayer(LoopInd) * 1000.0 + SoilWaterStorage = SoilWaterStorage + SoilMoisture(LoopInd) * ThicknessSnowSoilLayer(LoopInd) * 1000.0 enddo + WaterStorageTotEnd = SoilWaterStorage + CanopyLiqWater + CanopyIce + SnowWaterEquiv + WaterStorageAquifer + ! accumualted water change (only for canopy and snow during non-soil timestep) SfcWaterTotChgAcc = SfcWaterTotChgAcc + (WaterStorageTotEnd - WaterStorageTotBeg) ! snow, canopy, and soil water change PrecipTotAcc = PrecipTotAcc + PrecipTotRefHeight * MainTimeStep ! accumulated precip @@ -134,21 +139,44 @@ subroutine BalanceWaterCheck(noahmp) TileDrain) #ifndef WRF_HYDRO if ( abs(WaterBalanceError) > 0.1 ) then - if ( WaterBalanceError > 0 ) then - write(*,*) "The model is gaining water (WaterBalanceError is positive)" - else - write(*,*) "The model is losing water (WaterBalanceError is negative)" - endif - write(*,*) 'WaterBalanceError = ',WaterBalanceError, "kg m{-2} timestep{-1}" - write(*, & - '(" GridIndexI GridIndexJ SfcWaterTotChgAcc PrecipTotRefHeightAcc IrrigationRateMicro & - IrrigationRateFlood EvapCanopyNetAcc EvapGroundNetAcc TranspirationAcc RunoffSurface & - RunoffSubsurface WaterTableDepth TileDrain")') - write(*,'(i6,i6,f10.3,10f10.5)') GridIndexI, GridIndexJ, SfcWaterTotChgAcc, PrecipTotAcc, & - IrrigationRateMicro*1000.0, IrrigationRateFlood*1000.0, & - EvapCanopyNetAcc, EvapGroundNetAcc, TranspirationAcc, RunoffSurface, & - RunoffSubsurface, WaterTableDepth, TileDrain - stop "Error: Water budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP water budget conservation error (land):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' WaterStorageTotBeg = $r ', realArgs = (/ WaterStorageTotBeg /) ) + call mpas_log_write(' WaterStorageTotEnd = $r ', realArgs = (/ WaterStorageTotEnd /) ) + call mpas_log_write(' WaterBalanceError = $r ', realArgs = (/ WaterBalanceError /) ) + call mpas_log_write(' (positive value above means water gain).' ) + call mpas_log_write('--- State:') + call mpas_log_write(' SoilWaterStorage = $r ', realArgs = (/ SoilWaterStorage /) ) + call mpas_log_write(' CanopyLiqWater = $r ', realArgs = (/ CanopyLiqWater /) ) + call mpas_log_write(' CanopyIce = $r ', realArgs = (/ CanopyIce /) ) + call mpas_log_write(' SnowWaterEquiv = $r ', realArgs = (/ SnowWaterEquiv /) ) + call mpas_log_write(' WaterStorageAquifer = $r ', realArgs = (/ WaterStorageAquifer /) ) + call mpas_log_write('--- Fluxes:') + call mpas_log_write(' SfcWaterTotChgAcc = $r ', realArgs = (/ SfcWaterTotChgAcc /) ) + call mpas_log_write(' PrecipTotAcc = $r ', realArgs = (/ PrecipTotAcc /) ) + call mpas_log_write(' IrrigationRateMicro*1000.0 = $r ', realArgs = (/ IrrigationRateMicro*1000.0 /) ) + call mpas_log_write(' IrrigationRateFlood*1000.0 = $r ', realArgs = (/ IrrigationRateFlood*1000.0 /) ) + call mpas_log_write(' EvapCanopyNetAcc = $r ', realArgs = (/ EvapCanopyNetAcc /) ) + call mpas_log_write(' EvapGroundNetAcc = $r ', realArgs = (/ EvapGroundNetAcc /) ) + call mpas_log_write(' TranspirationAcc = $r ', realArgs = (/ TranspirationAcc /) ) + call mpas_log_write(' RunoffSurface = $r ', realArgs = (/ RunoffSurface /) ) + call mpas_log_write(' RunoffSubsurface = $r ', realArgs = (/ RunoffSubsurface /) ) + call mpas_log_write(' WaterTableDepth = $r ', realArgs = (/ WaterTableDepth /) ) + call mpas_log_write(' TileDrain = $r ', realArgs = (/ TileDrain /) ) + call mpas_log_write('--- Soil state:') + call mpas_log_write(' NumSoilLayer = $i ',intArgs=(/ NumSoilLayer /)) + do LoopInd=1,NumSoilLayer + call mpas_log_write(' Layer $i - Thickness = $r ; SoilMoisture = $r' & + , intArgs = (/ LoopInd /) & + , realArgs = (/ ThicknessSnowSoilLayer(LoopInd), SoilMoisture(LoopInd) /) ) + end do + call mpas_log_write('---~---') + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Water budget problem in NoahMP LSM (land)") endif #endif endif ! FlagSoilProcess @@ -208,23 +236,31 @@ subroutine BalanceEnergyCheck(noahmp) RadSwBalanceError = RadSwDownRefHeight - (RadSwAbsSfc + RadSwReflSfc) ! print out diagnostics when error is large if ( abs(RadSwBalanceError) > 0.01 ) then - write(*,*) "GridIndexI, GridIndexJ = ", GridIndexI, GridIndexJ - write(*,*) "RadSwBalanceError = ", RadSwBalanceError - write(*,*) "VEGETATION ---------" - write(*,*) "RadSwDownRefHeight * VegFrac = ", RadSwDownRefHeight*VegFrac - write(*,*) "VegFrac*RadSwAbsVeg + RadSwAbsGrd = ", VegFrac*RadSwAbsVeg+RadSwAbsGrd - write(*,*) "VegFrac*RadSwReflVeg + RadSwReflGrd = ", VegFrac*RadSwReflVeg+RadSwReflGrd - write(*,*) "GROUND -------" - write(*,*) "(1 - VegFrac) * RadSwDownRefHeight = ", (1.0-VegFrac)*RadSwDownRefHeight - write(*,*) "(1 - VegFrac) * RadSwAbsGrd = ", (1.0-VegFrac)*RadSwAbsGrd - write(*,*) "(1 - VegFrac) * RadSwReflGrd = ", (1.0-VegFrac)*RadSwReflGrd - write(*,*) "RadSwReflVeg = ", RadSwReflVeg - write(*,*) "RadSwReflGrd = ", RadSwReflGrd - write(*,*) "RadSwReflSfc = ", RadSwReflSfc - write(*,*) "RadSwAbsVeg = ", RadSwAbsVeg - write(*,*) "RadSwAbsGrd = ", RadSwAbsGrd - write(*,*) "RadSwAbsSfc = ", RadSwAbsSfc - stop "Error: Solar radiation budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP solar radiation budget conservation error (land):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) ) + call mpas_log_write(' (positive value above means energy gain).' ) + call mpas_log_write('--- Vegetation ') + call mpas_log_write(' RadSwDownRefHeight * VegFrac = $r ', realArgs = (/ RadSwDownRefHeight*VegFrac /) ) + call mpas_log_write(' VegFrac*RadSwAbsVeg + RadSwAbsGrd = $r ', realArgs = (/ VegFrac*RadSwAbsVeg+RadSwAbsGrd /) ) + call mpas_log_write(' VegFrac*RadSwReflVeg + RadSwReflGrd = $r ', realArgs = (/ VegFrac*RadSwReflVeg+RadSwReflGrd /) ) + call mpas_log_write('--- Ground ') + call mpas_log_write(' (1 - VegFrac) * RadSwDownRefHeight = $r ', realArgs = (/ (1.0-VegFrac)*RadSwDownRefHeight /) ) + call mpas_log_write(' (1 - VegFrac) * RadSwAbsGrd = $r ', realArgs = (/ (1.0-VegFrac)*RadSwAbsGrd /) ) + call mpas_log_write(' (1 - VegFrac) * RadSwReflGrd = $r ', realArgs = (/ (1.0-VegFrac)*RadSwReflGrd /) ) + call mpas_log_write(' RadSwReflVeg = $r ', realArgs = (/ RadSwReflVeg /) ) + call mpas_log_write(' RadSwReflGrd = $r ', realArgs = (/ RadSwReflGrd /) ) + call mpas_log_write(' RadSwReflSfc = $r ', realArgs = (/ RadSwReflSfc /) ) + call mpas_log_write(' RadSwAbsVeg = $r ', realArgs = (/ RadSwAbsVeg /) ) + call mpas_log_write(' RadSwAbsGrd = $r ', realArgs = (/ RadSwAbsGrd /) ) + call mpas_log_write(' RadSwAbsSfc = $r ', realArgs = (/ RadSwAbsSfc /) ) + call mpas_log_write('---~---') + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Solar radiation budget problem in NoahMP LSM (land)") endif ! error in surface energy balance should be <0.01 W/m2 @@ -233,19 +269,32 @@ subroutine BalanceEnergyCheck(noahmp) HeatLatentTransp + HeatGroundTot + HeatLatentIrriEvap + HeatCanStorageChg) ! print out diagnostics when error is large if ( abs(EnergyBalanceError) > 0.01 ) then - write(*,*) 'EnergyBalanceError = ', EnergyBalanceError, ' at GridIndexI,GridIndexJ: ', GridIndexI, GridIndexJ - write(*,'(a17,F10.4)' ) "Net solar: ", RadSwAbsSfc - write(*,'(a17,F10.4)' ) "Net longwave: ", RadLwNetSfc - write(*,'(a17,F10.4)' ) "Total sensible: ", HeatSensibleSfc - write(*,'(a17,F10.4)' ) "Canopy evap: ", HeatLatentCanopy - write(*,'(a17,F10.4)' ) "Ground evap: ", HeatLatentGrd - write(*,'(a17,F10.4)' ) "Transpiration: ", HeatLatentTransp - write(*,'(a17,F10.4)' ) "Total ground: ", HeatGroundTot - write(*,'(a17,F10.4)' ) "Sprinkler: ", HeatLatentIrriEvap - write(*,'(a17,F10.4)' ) "Canopy heat storage change: ", HeatCanStorageChg - write(*,'(a17,4F10.4)') "Precip advected: ", HeatPrecipAdvSfc,HeatPrecipAdvCanopy,HeatPrecipAdvVegGrd,HeatPrecipAdvBareGrd - write(*,'(a17,F10.4)' ) "Veg fraction: ", VegFrac - stop "Error: Energy budget problem in NoahMP LSM" + call mpas_log_write(' ') + call mpas_log_write('---~---') + call mpas_log_write(' Noah-MP energy budget conservation error (land):') + call mpas_log_write('---~---') + call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) ) + call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) ) + call mpas_log_write(' EnergyBalanceError = $r ', realArgs = (/ EnergyBalanceError /) ) + call mpas_log_write(' (positive value above energy gain).' ) + call mpas_log_write(' Net solar = $r ', realArgs = (/ RadSwAbsSfc /) ) + call mpas_log_write(' Net longwave = $r ', realArgs = (/ RadLwNetSfc /) ) + call mpas_log_write(' Total sensible = $r ', realArgs = (/ HeatSensibleSfc /) ) + call mpas_log_write(' Canopy evap = $r ', realArgs = (/ HeatLatentCanopy /) ) + call mpas_log_write(' Ground evap = $r ', realArgs = (/ HeatLatentGrd /) ) + call mpas_log_write(' Transpiration = $r ', realArgs = (/ HeatLatentTransp /) ) + call mpas_log_write(' Total ground = $r ', realArgs = (/ HeatGroundTot /) ) + call mpas_log_write(' Sprinkler = $r ', realArgs = (/ HeatLatentIrriEvap /) ) + call mpas_log_write(' Canopy heat storage change = $r ', realArgs = (/ HeatCanStorageChg /) ) + + call mpas_log_write(' Precip advected (surface) = $r ', realArgs = (/ HeatPrecipAdvSfc /) ) + call mpas_log_write(' Precip advected (canopy) = $r ', realArgs = (/ HeatPrecipAdvCanopy /) ) + call mpas_log_write(' Precip advected (veg. grnd) = $r ', realArgs = (/ HeatPrecipAdvVegGrd /) ) + call mpas_log_write(' Precip advected (bare grnd) = $r ', realArgs = (/ HeatPrecipAdvBareGrd /) ) + call mpas_log_write(' Veg fraction = $r ', realArgs = (/ VegFrac /) ) + call mpas_log_write('---~---') + call mpas_log_write(' ') + call Noahmp_error_fatal("Error: Energy budget problem in NoahMP LSM (land)") endif end associate diff --git a/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarInitMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarInitMod.F90 index 5c8af537b..2f9c5e9e2 100644 --- a/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarInitMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarInitMod.F90 @@ -49,6 +49,7 @@ subroutine ConfigVarInitDefault(noahmp) ! config domain variable noahmp%config%domain%LandUseDataName = "MODIFIED_IGBP_MODIS_NOAH" + noahmp%config%domain%SoilColDataName = "DEFAULT_RAD_NOAH" noahmp%config%domain%FlagUrban = .false. noahmp%config%domain%FlagCropland = .false. noahmp%config%domain%FlagDynamicCrop = .false. diff --git a/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarType.F90 b/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarType.F90 index dc7979f3c..1ea5013d6 100644 --- a/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarType.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/src/ConfigVarType.F90 @@ -128,6 +128,7 @@ module ConfigVarType type :: domain_type character(len=256) :: LandUseDataName ! landuse dataset name (USGS or MODIFIED_IGBP_MODIS_NOAH) + character(len=256) :: SoilColDataName ! soil colour dataset name (DEFAULT_RAD_NOAH or MODIFIED_RAD_CLM_NOAH) logical :: FlagUrban ! flag for urban grid logical :: FlagCropland ! flag to identify croplands logical :: FlagDynamicCrop ! flag to activate dynamic crop model diff --git a/src/core_atmosphere/physics/physics_noahmp/src/Makefile b/src/core_atmosphere/physics/physics_noahmp/src/Makefile index 675bdf9df..f4e2b2bb3 100644 --- a/src/core_atmosphere/physics/physics_noahmp/src/Makefile +++ b/src/core_atmosphere/physics/physics_noahmp/src/Makefile @@ -31,6 +31,7 @@ OBJS = ConstantDefineMod.o \ IrrigationInfilPhilipMod.o \ IrrigationMicroMod.o \ MatrixSolverTriDiagonalMod.o \ + NoahmpFatalErrorMod.o \ RunoffSubSurfaceDrainageMod.o \ RunoffSubSurfaceEquiWaterTableMod.o \ RunoffSubSurfaceGroundWaterMod.o \ @@ -168,6 +169,7 @@ IrrigationInfilPhilipMod.o: ../utility/Machine.o NoahmpVarType.o Const IrrigationMicroMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ IrrigationInfilPhilipMod.o MatrixSolverTriDiagonalMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o +NoahmpFatalErrorMod.o: ../utility/Machine.o RunoffSubSurfaceDrainageMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o RunoffSubSurfaceEquiWaterTableMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ WaterTableEquilibriumMod.o @@ -301,10 +303,10 @@ BiochemNatureVegMainMod.o: ../utility/Machine.o NoahmpVarType.o Co BiochemCropMainMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o CarbonFluxCropMod.o \ CropGrowDegreeDayMod.o CropPhotosynthesisMod.o IrrigationPrepareMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o IrrigationTriggerMod.o -BalanceErrorCheckMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o +BalanceErrorCheckMod.o: ../utility/Machine.o NoahmpFatalErrorMod.o NoahmpVarType.o ConstantDefineMod.o GeneralInitMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o GroundWaterMmfMod.o: ../utility/Machine.o NoahmpVarType.o ../drivers/hrldas/NoahmpIOVarType.o -BalanceErrorCheckGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o +BalanceErrorCheckGlacierMod.o: ../utility/Machine.o NoahmpFatalErrorMod.o NoahmpVarType.o ConstantDefineMod.o EnergyMainGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowCoverGlacierMod.o \ GroundRoughnessPropertyGlacierMod.o GroundThermalPropertyGlacierMod.o \ SurfaceAlbedoGlacierMod.o SurfaceRadiationGlacierMod.o SurfaceEmissivityGlacierMod.o \ @@ -347,5 +349,5 @@ clean: $(RM) *.i .F90.o: - $(FC) $(CPPFLAGS) $(COREDEF) $(FFLAGS_NONSTD) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I../utility -I../../../../framework + $(FC) $(CPPFLAGS) $(COREDEF) $(FFLAGS_NONSTD) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I../utility -I../../../../framework -I../../../../external/esmf_time_f90 diff --git a/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90 new file mode 100644 index 000000000..b42c4be30 --- /dev/null +++ b/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90 @@ -0,0 +1,39 @@ +module NoahmpFatalErrorMod + ! This module defines an error message for Noah-MP that allows for gracefully failing. + implicit none + private + public :: NoahMP_error_fatal + + ! I would rather call mpas_derived_types from within this module, but I tried so many + ! things and I still cannot figure out a way of loading it (somehow mpas_log seems to + ! be accessible, but not mpas_derived_types). For the time being I will duplicate the + ! values, fully aware that this is not the correct thing to do, but please, if you + ! know how to make this neat, let me know... + + integer, parameter :: NOAHMP_LOG_OUT = 1 !< code for message type "output" + integer, parameter :: NOAHMP_LOG_WARN = 2 !< code for message type "warning" + integer, parameter :: NOAHMP_LOG_ERR = 3 !< code for message type "error" + integer, parameter :: NOAHMP_LOG_CRIT = 4 !< code for message type "critical error" + + contains + + !---~--- + ! Print error message then gracefully exit. + !---~--- + subroutine Noahmp_error_fatal(str) + + use mpas_log , only : mpas_log_write + + ! input arguments: + character(len=*),intent(in):: str + + + call mpas_log_write(' ' , messageType=NOAHMP_LOG_ERR) + call mpas_log_write('---~---' , messageType=NOAHMP_LOG_ERR) + call mpas_log_write(' Noah-MP FATAL ERROR', messageType=NOAHMP_LOG_ERR) + call mpas_log_write('---~---' , messageType=NOAHMP_LOG_ERR) + call mpas_log_write(trim(str) , messageType=NOAHMP_LOG_ERR) + call mpas_log_write('Noah-MP abort' , messageType=NOAHMP_LOG_CRIT) + end subroutine Noahmp_error_fatal + !---~--- +end module NoahmpFatalErrorMod diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 7b8905fff..831976fbf 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -171,6 +171,11 @@ description="The soil category classification to use" possible_values="`STATSGO' or `BNU'"/> + + + + + + + + + + + @@ -543,10 +563,12 @@ + + @@ -598,6 +620,8 @@ + + @@ -844,9 +868,15 @@ + + + + @@ -1218,6 +1248,14 @@ + + + + * isice = 24 !> * isurban = 1 !> * isoilwater = 14 + !> * islcolwater = 9 ! !----------------------------------------------------------------------- function mpas_geotile_mgr_init(mgr, path) result(ierr) @@ -119,7 +120,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) integer, pointer :: tile_z_start, tile_z_end integer, pointer :: signed integer, pointer :: tile_bdr - integer, pointer :: iswater, islake, isice, isurban, isoilwater + integer, pointer :: iswater, islake, isice, isurban, isoilwater, islcolwater integer, pointer :: category_min, category_max integer :: err_level real (kind=RKIND), pointer :: dx ! Grid spacing in the x-direction @@ -188,7 +189,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) ! ! If this is a categorical field, then check to see if it has category_max and category_min, - ! and then set the defaults of iswater, islake, isice, isurban and isoilwater + ! and then set the defaults of iswater, islake, isice, isurban, isoilwater and islcolwater ! call mpas_pool_get_config(mgr % pool, 'type', fieldType) if (fieldType == 'categorical') then @@ -219,12 +220,14 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) isice => null() isurban => null() isoilwater => null() + islcolwater => null() call mpas_pool_get_config(mgr % pool, 'iswater', iswater) call mpas_pool_get_config(mgr % pool, 'islake', islake) call mpas_pool_get_config(mgr % pool, 'isice', isice) call mpas_pool_get_config(mgr % pool, 'isurban', isurban) call mpas_pool_get_config(mgr % pool, 'isoilwater', isoilwater) + call mpas_pool_get_config(mgr % pool, 'islcolwater', islcolwater) if (.not. associated(iswater)) then call mpas_pool_add_config(mgr % pool, 'iswater', 16) @@ -245,6 +248,10 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) if (.not. associated(isoilwater)) then call mpas_pool_add_config(mgr % pool, 'isoilwater', 14) endif + + if (.not. associated(islcolwater)) then + call mpas_pool_add_config(mgr % pool, 'islcolwater', 1) + end if endif ! diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index 44e58bc78..9a3a3845e 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -74,6 +74,7 @@ subroutine init_atm_setup_case(domain, stream_manager) character (len=StrKIND), pointer :: config_specified_zeta_levels character(len=StrKIND), pointer :: mminlu + character(len=StrKIND), pointer :: mminsc character(len=StrKIND), pointer :: xtime real (kind=RKIND) :: dt real (kind=RKIND), pointer :: Time @@ -257,6 +258,20 @@ subroutine init_atm_setup_case(domain, stream_manager) write(mminlu,'(a)') 'USGS' end if + ! + ! If at this point the mminsc variable is blank, we assume that the static interp step was + ! not run, and that we are working with a static file created before there was a choice + ! of soil colour datasets; in this case, the dataset was almost necessarily DEFAULT_RAD_NOAH + ! + call mpas_pool_get_array(mesh, 'mminsc', mminsc) + if (len_trim(mminsc) == 0) then + call mpas_log_write('****************************************************************') + call mpas_log_write('No information on soil colour dataset is available.') + call mpas_log_write('Assume that we are using ''DEFAULT_RAD_NOAH''.') + call mpas_log_write('****************************************************************') + write(mminsc,'(a)') 'DEFAULT_RAD_NOAH' + end if + call init_atm_case_gfs(block_ptr, mesh, nCells, nEdges, nVertLevels, fg, state, & diag, diag_physics, block_ptr % dimensions, block_ptr % configs) diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 0c540094d..2720643f4 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -93,6 +93,7 @@ end subroutine interp_accumulation_function real (kind=RKIND) :: soilcomp_msgval = 255.0_RKIND ! Modified later based on index file for soilcomp integer, dimension(:), pointer :: lu_index integer, dimension(:), pointer :: soilcat_top + integer, dimension(:), pointer :: soilcol_idx integer, dimension(:), pointer :: nhs integer, dimension(:,:), allocatable:: ncat ! Landmask is used by the accumulation function for maxsnoalb and soilcomp, @@ -122,6 +123,7 @@ subroutine init_atm_static(mesh, dims, configs) character(kind=c_char), dimension(StrKIND+1) :: c_fname character(len=StrKIND), pointer :: config_geog_data_path character(len=StrKIND), pointer :: config_landuse_data + character(len=StrKIND), pointer :: config_soilcol_data character(len=StrKIND), pointer :: config_topo_data character(len=StrKIND), pointer :: config_vegfrac_data character(len=StrKIND), pointer :: config_albedo_data @@ -155,6 +157,8 @@ subroutine init_atm_static(mesh, dims, configs) integer, pointer :: isice_lu, iswater_lu integer :: iswater_soil + integer :: iswater_slcol + integer :: isdefault_slcol integer, pointer :: nCells, nCellsSolve, nEdges, nVertices, maxEdges logical, pointer :: on_a_sphere real (kind=RKIND), pointer :: sphere_radius @@ -185,9 +189,11 @@ subroutine init_atm_static(mesh, dims, configs) real (kind=RKIND), pointer :: missing_value integer, dimension(:), pointer :: lu_index integer, dimension(:), pointer :: soilcat_top + integer, dimension(:), pointer :: soilcol_idx integer, dimension(:), pointer :: landmask integer, dimension(:), pointer :: bdyMaskCell character(len=StrKIND), pointer :: mminlu + character(len=StrKIND), pointer :: mminsc real (kind=RKIND) :: xPixel, yPixel, zPixel @@ -217,6 +223,7 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_pool_get_config(configs, 'config_geog_data_path', config_geog_data_path) call mpas_pool_get_config(configs, 'config_landuse_data', config_landuse_data) call mpas_pool_get_config(configs, 'config_soilcat_data', config_soilcat_data) + call mpas_pool_get_config(configs, 'config_soilcol_data', config_soilcol_data) call mpas_pool_get_config(configs, 'config_topo_data', config_topo_data) call mpas_pool_get_config(configs, 'config_vegfrac_data', config_vegfrac_data) call mpas_pool_get_config(configs, 'config_albedo_data', config_albedo_data) @@ -267,9 +274,11 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_pool_get_array(mesh, 'lu_index', lu_index) call mpas_pool_get_array(mesh, 'mminlu', mminlu) + call mpas_pool_get_array(mesh, 'mminsc', mminsc) call mpas_pool_get_array(mesh, 'isice_lu', isice_lu) call mpas_pool_get_array(mesh, 'iswater_lu', iswater_lu) call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top) + call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx) call mpas_pool_get_array(mesh, 'landmask', landmask) call mpas_pool_get_array(mesh, 'snoalb', snoalb) call mpas_pool_get_array(mesh, 'greenfrac', greenfrac) @@ -366,6 +375,23 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT) end select surface_input_select0 +! +! Set soil colour parameters for soil albedo +! + slcol_input_select0: select case(trim(config_soilcol_data)) + case('DEFAULT_RAD_NOAH') + write(mminsc,'(a)') 'DEFAULT_RAD_NOAH' + case('MODIFIED_RAD_CLM_NOAH') + write(mminsc,'(a)') 'MODIFIED_RAD_CLM_NOAH' + case default + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Invalid soil colour dataset '''//trim(config_soilcol_data) & + //''' selected for config_soilcol_data', messageType=MPAS_LOG_ERR) + call mpas_log_write(' Possible options are: ''DEFAULT_RAD_NOAH'', ''MODIFIED_RAD_CLM_NOAH''', messageType=MPAS_LOG_ERR) + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT) + end select slcol_input_select0 + ! ! Interpolate HGT ! @@ -443,17 +469,53 @@ subroutine init_atm_static(mesh, dims, configs) supersample_fac=supersample_fac_30s) call mpas_log_write('--- end interpolate SOILCAT_TOP') +! +! Interpolate SOILCOL_IDX +! + slcol_input_select1: select case(trim(config_soilcol_data)) + case('DEFAULT_RAD_NOAH') + call mpas_log_write('Using 9-class default NOAH soil colour classes') + geog_sub_path = 'soilcolour_30s/' + ! Set default category and default water category for NOAH soil colour. + ! These will be overwritten if soil colour data sets exist. + isdefault_slcol = 4 + iswater_slcol = 1 + + case('MODIFIED_RAD_CLM_NOAH') + call mpas_log_write('Using 21-class MODIS 30-arc-second land cover dataset') + geog_sub_path = 'clm_soilcolour_21class_30s/' + + ! Set default category and default water category for NOAH soil colour + ! These will be overwritten if soil colour data sets exist. + isdefault_slcol = 14 + iswater_slcol = 21 + case default + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Invalid soil colour dataset '''//trim(config_soilcol_data) & + //''' selected for config_landuse_data', messageType=MPAS_LOG_ERR) + call mpas_log_write(' Possible options are: ''DEFAULT_RAD_NOAH'', ''MODIFIED_RAD_CLM_NOAH''', messageType=MPAS_LOG_ERR) + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT) + end select slcol_input_select1 + + ! Call interpolation routine for soil colour. + call mpas_log_write('--- start interpolate SOILCOL_IDX') + call interp_soilcol(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), & + isdefault_slcol, iswater_slcol, supersample_fac=supersample_fac_30s) + call mpas_log_write('--- end interpolate SOILCOL_IDX') + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! where (lu_index == isice_lu) soilcat_top = 16 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! CORRECT INCONSISTENT SOIL AND LAND USE DATA +! CORRECT INCONSISTENT SOIL TYPE, SOIL COLOUR AND LAND USE DATA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do iCell = 1,nCells - if (lu_index(iCell) == iswater_lu .or. & - soilcat_top(iCell) == iswater_soil) then + if (lu_index(iCell) == iswater_lu .or. & + soilcat_top(iCell) == iswater_soil .or. & + soilcol_idx(iCell) == iswater_slcol ) then if (lu_index(iCell) /= iswater_lu) then call mpas_log_write('Turning lu_index into water at $i', intArgs=(/iCell/)) lu_index(iCell) = iswater_lu @@ -462,6 +524,10 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_log_write('Turning soilcat_top into water at $i', intArgs=(/iCell/)) soilcat_top(iCell) = iswater_soil end if + if (soilcol_idx(iCell) /= iswater_slcol) then + call mpas_log_write('Turning soilcol_idx into water at $i', intArgs=(/iCell/)) + soilcol_idx(iCell) = iswater_slcol + end if end if end do @@ -1936,6 +2002,104 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl end subroutine interp_soilcat + !*********************************************************************** + ! + ! routine interp_soilcol + ! + !> \brief Interpolate soil colour indices + !> \author Marcos Longo + !> \date 15 October 2025 + !> \details + !> This sub-routine is based on interp_soilcat. The procedure interpolates soil + !> colour indices by using the init_atm_map_static_data routine and then by + !> accumulating the pixel values into each cell using + !> category_interp_accumulation. + !> + !> Variable mesh should be an mpas_pool that contains nCells and soilcol_idx; + !> kdtree should be an initialised mpas_kd_type tree with dimensions + !> (xCell, yCell, zCell); and geog_data_path should be the path to the + !> soil colour index data set. The values used by this data set to flag + !> pixels with no data (e.g., pixels fully covered with water or ice). Some + !> data sets may not have the no data flag, but we harmonise this data set + !> with the land use class mapping. + !> + !> In case the path is not found, instead of issuing an error, the routine will + !> simply assume constant soil colour indices across the entire domain. This + !> is done because data sets may or may not be available, so we add the fall + !> back is implemented for back compatibility. + !----------------------------------------------------------------------- + subroutine interp_soilcol(mesh, kdtree, geog_data_path, isdefault_slcol, iswater_slcol, & + supersample_fac) + + implicit none + + ! Input variables + type (mpas_pool_type), intent(inout) :: mesh + type (mpas_kd_type), pointer, intent(in) :: kdtree + character (len=*), intent(in) :: geog_data_path + integer, intent(in) :: isdefault_slcol + integer, intent(inout) :: iswater_slcol + integer, intent(in), optional :: supersample_fac + + ! Local variables + type (mpas_geotile_mgr_type) :: mgr + integer, pointer :: nCells + integer, pointer :: iswater_slcol_ptr + + real (kind=RKIND), pointer :: scalefactor + + integer :: iCell + integer :: ierr + + ierr = mgr % init(trim(geog_data_path)) + if (ierr /= 0) then + call mpas_log_write("--> Warning: error occured initalising interpolation for "//trim(geog_data_path)//".") + call mpas_log_write(" Assuming constant soil colour for the entire domain.") + + call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx) + + ! Define the entire domain with the default soil colour index + soilcol_idx(1:nCells) = isdefault_slcol + + return + end if + + call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx) + call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor ) + call mpas_pool_get_config(mgr % pool, 'category_min', category_min ) + call mpas_pool_get_config(mgr % pool, 'category_max', category_max ) + call mpas_pool_get_config(mgr % pool, 'islcolwater' , iswater_slcol_ptr) + + iswater_slcol = iswater_slcol_ptr + + allocate(ncat(category_min:category_max, nCells)) + ncat(:,:) = 0 + + call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, & + supersample_fac=supersample_fac) + + do iCell = 1, nCells + ! Variable ncat will tally the classes occurring in the grid cell. We select the class that has the + ! highest count as the value representing the grid cell. We use maxloc for this, however, maxloc returns + ! the index corresponding to the maximum value starting from 1. We shift the result to make sure the + ! category value correctly maps between category_min and category_max. + soilcol_idx(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min + end do + deallocate(ncat) + + ierr = mgr % finalize() + if (ierr /= 0) then + call mpas_log_write("Error occured finalising interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT) + return + end if + + nullify(category_min) + nullify(category_max) + + end subroutine interp_soilcol + !*********************************************************************** ! ! routine derive_landmask diff --git a/src/core_init_atmosphere/mpas_parse_geoindex.F b/src/core_init_atmosphere/mpas_parse_geoindex.F index 753ed4ee8..0be2c8098 100644 --- a/src/core_init_atmosphere/mpas_parse_geoindex.F +++ b/src/core_init_atmosphere/mpas_parse_geoindex.F @@ -131,7 +131,8 @@ function mpas_parse_index(path, geo_pool) result(ierr) .or. trim(lhs) == 'description' & .or. trim(lhs) == 'row_order' & .or. trim(lhs) == 'endian' & - .or. trim(lhs) == 'mminlu' ) then + .or. trim(lhs) == 'mminlu' & + .or. trim(lhs) == 'mminsc' ) then char_t = rhs call mpas_pool_add_config(geo_pool, trim(lhs), char_t) @@ -171,6 +172,7 @@ function mpas_parse_index(path, geo_pool) result(ierr) .or. trim(lhs) == 'isice' & .or. trim(lhs) == 'isurban' & .or. trim(lhs) == 'isoilwater' & + .or. trim(lhs) == 'islcolwater' & .or. trim(lhs) == 'filename_digits' ) then ! Because each compiler handles reporting type errors when transferring diff --git a/src/framework/mpas_io.F b/src/framework/mpas_io.F index 09514a366..bf5a89804 100644 --- a/src/framework/mpas_io.F +++ b/src/framework/mpas_io.F @@ -868,7 +868,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('PIO_SUPPORT: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if !call mpas_log_write('Inquired about variable ID $i.', intArgs=(/new_fieldlist_node % fieldhandle % fieldid/) ) @@ -907,7 +907,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('SMIOL_SUPPORT_A: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if new_fieldlist_node % fieldhandle % field_type = smiol_type @@ -955,7 +955,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('SMIOL_SUPPORT_B: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if new_fieldlist_node% fieldhandle % ndims = smiol_ndims