diff --git a/.circleci/config.yml b/.circleci/config.yml index 6b7ad7dddf85..5b107ce77773 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -7,7 +7,7 @@ jobs: - checkout - run: bash .circleci/setup.sh - run: bash .circleci/checkout.sh - - run: wget https://dabdceba-6d04-11e5-ba46-22000b92c6ec.e.globus.org/containers/public/e3sm.sif + - run: bash .circleci/container.sh - run: command: singularity exec --hostname singularity e3sm.sif .circleci/run.sh no_output_timeout: 60m diff --git a/.circleci/container.sh b/.circleci/container.sh new file mode 100755 index 000000000000..b2f1d7b90d35 --- /dev/null +++ b/.circleci/container.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +wget -t 3 http://portal.nersc.gov/project/e3sm/lukasz/e3sm.sif +if [ $? -eq 0 ]; then + exit 0; +fi + +wget -t 3 https://dabdceba-6d04-11e5-ba46-22000b92c6ec.e.globus.org/containers/public/e3sm.sif +if [ $? -eq 0 ]; then + exit 0; +fi + +exit -1 diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index df1e63022b34..000000000000 --- a/.travis.yml +++ /dev/null @@ -1,46 +0,0 @@ -os: linux - -dist: bionic - -git: - submodules: false - -python: - - "3.7" - -addons: - apt: - packages: - - flawfinder - - squashfs-tools - - build-essential - - uuid-dev - - libuuid1 - - libffi-dev - - libssl-dev - - libssl1.0.0 - - libarchive-dev - - libgpgme11-dev - - libseccomp-dev - - pkg-config - - cryptsetup-bin - homebrew: - packages: - - squashfs - update: true - -sudo: required - -before_install: - # Replace all ssh URLs to submodules with HTTP URLs - - sed -i 's/git@github.com:/https:\/\/github.com\//' .gitmodules - - git submodule update --init - # Download the E3SM Singularity container - - wget https://dabdceba-6d04-11e5-ba46-22000b92c6ec.e.globus.org/containers/public/e3sm.sif - # Install Singularity - - sudo chmod u+x .travis/*.sh - - .travis/setup.sh - -script: - # Build a case - - travis_wait 30 singularity exec --hostname singularity e3sm.sif .travis/run.sh diff --git a/.travis/run.sh b/.travis/run.sh deleted file mode 100755 index 744d53d798ab..000000000000 --- a/.travis/run.sh +++ /dev/null @@ -1,51 +0,0 @@ -#!/bin/bash - -case="master.A_WCYCL1850.ne4_oQU240.baseline" -compset="A_WCYCL1850" -res="ne4_oQU240" - -bldlog_dir="${HOME}/projects/e3sm/scratch/$case/bld" - -shopt -s extglob -shopt -s nullglob - -function command_yellow -{ - echo -e "\e[33m\$ $1\e[0m" - eval $1 - return $? -} - -function command_yellow_rc -{ - command_yellow "$1" - rc=$? - if [ $rc -ne 0 ]; then - exit $rc - fi -} - -function print_bldlog -{ - bldlog=$1 - bldlog_path=(${bldlog_dir}/$bldlog.bldlog.+([[:digit:]])-+([[:digit:]])) - if [ -n "${bldlog_path}" ]; then - command_yellow "cat ${bldlog_path}" - fi -} - -command_yellow_rc "mkdir -p $HOME/projects/e3sm/cesm-inputdata" -command_yellow_rc "cd cime/scripts" -command_yellow_rc "./create_newcase --case $case --compset $compset --res $res" -command_yellow_rc "cd $case" -command_yellow_rc "./case.setup" -command_yellow "./case.build" -rc=$? -if [ $rc -ne 0 ]; then - print_bldlog "csm_share" - print_bldlog "pio" - print_bldlog "mct" - print_bldlog "gptl" - print_bldlog "e3sm" - exit $rc -fi diff --git a/.travis/setup.sh b/.travis/setup.sh deleted file mode 100755 index 6023ec7c3be6..000000000000 --- a/.travis/setup.sh +++ /dev/null @@ -1,34 +0,0 @@ -#!/bin/bash - -# Remove the pre-installed Go 1.11 from the environment -go version -for p in $(xargs -n1 -d: <<< $PATH); do - echo $p | grep '/go' > /dev/null - if [ $? -ne 0 ]; then - TMPPATH="$p:$TMPPATH" - fi -done -export PATH="$TMPPATH" -unset GOPATH -unset GOROOT - -# Install Go 1.13.5 -export VERSION=1.13.5 OS=linux ARCH=amd64 && \ - wget https://dl.google.com/go/go$VERSION.$OS-$ARCH.tar.gz && \ - sudo tar -C /usr/local -xzf go$VERSION.$OS-$ARCH.tar.gz && \ - rm go$VERSION.$OS-$ARCH.tar.gz - -export GOPATH="${HOME}/go" -export PATH="/usr/local/go/bin:${PATH}:${GOPATH}/bin" -go version - -# Install Singularity -export VERSION=3.5.3 && # adjust this as necessary \ - wget https://github.com/sylabs/singularity/releases/download/v${VERSION}/singularity-${VERSION}.tar.gz && \ - tar -xzf singularity-${VERSION}.tar.gz && \ - cd singularity -./mconfig && \ - make -C ./builddir && \ - sudo make -C ./builddir install - -singularity --version diff --git a/cime_config/allactive/config_pesall.xml b/cime_config/allactive/config_pesall.xml index 4b3c52c6179c..96c5d9139e01 100644 --- a/cime_config/allactive/config_pesall.xml +++ b/cime_config/allactive/config_pesall.xml @@ -8387,12 +8387,12 @@ - -compset A_WCYCL* -res ne30pg*EC30to60* on 14 nodes pure-MPI, ~6 sypd + -compset A_WCYCL* -res ne30pg*EC30to60* on 14 nodes pure-MPI, ~8.1 sypd 675 - 128 - 128 - 576 + 64 + 64 + 640 192 704 @@ -8406,21 +8406,21 @@ 0 - 576 - 576 + 640 + 640 0 704 0 - -compset A_WCYCL* -res ne30pg*EC30to60* on 26 nodes pure-MPI, ~10 sypd + -compset A_WCYCL* -res ne30pg*EC30to60* on 28 nodes pure-MPI, ~15.3 sypd 1350 - 128 - 128 - 1280 - 256 + 64 + 64 + 1344 + 384 1408 @@ -8433,21 +8433,21 @@ 0 - 1280 - 1280 + 1344 + 1344 0 1408 0 - - -compset A_WCYCL* -res ne30pg*EC30to60* on 51 nodes pure-MPI, ~14 sypd + + -compset A_WCYCL* -res ne30pg*EC30to60* on 53 nodes pure-MPI, ~26.4 sypd - 2700 + 2752 192 192 2560 - 512 + 640 2752 @@ -8467,15 +8467,15 @@ 0 - - -compset A_WCYCL* -res ne30pg*EC30to60* on 60 nodes pure-MPI, ~24 sypd + + -compset A_WCYCL* -res ne30pg*EC30to60* on 71 nodes pure-MPI, ~31.5 sypd - 3200 - 640 - 640 - 2560 - 640 - 3200 + 3648 + 192 + 192 + 3456 + 896 + 3648 1 @@ -8487,22 +8487,22 @@ 0 - 2560 - 2560 + 3456 + 3456 0 - 3200 + 3648 0 - - -compset A_WCYCL* -res ne30pg*EC30to60* on 80 nodes pure-MPI, ~30 sypd + + -compset A_WCYCL* -res ne30pg*EC30to60* on 105 nodes pure-MPI, ~42.3 sypd - 4160 + 5400 320 320 - 3840 - 960 - 4160 + 5120 + 1280 + 5440 1 @@ -8514,22 +8514,53 @@ 0 - 3840 - 3840 + 5120 + 5120 + 0 + 5440 + 0 + + + + + + + + -compset A_WCYCL* -res ne30pg*EC30to60* on 8 debug Q nodes threaded, 0.8 sypd + + 384 + 64 + 64 + 320 + 128 + 384 + + + 2 + 2 + 2 + 2 + 2 + 2 + + + 0 + 320 + 320 0 - 4160 + 384 0 - - -compset A_WCYCL* -res ne30pg*EC30to60* on 105 nodes pure-MPI, ~38 sypd + + -compset A_WCYCL* -res ne30pg*EC30to60* on 128 default Q nodes pure-MPI, 3.8 sypd 5400 - 320 - 320 - 5120 + 512 + 512 + 6400 1280 - 5440 + 6912 1 @@ -8541,10 +8572,10 @@ 0 - 5120 - 5120 + 6400 + 6400 0 - 5440 + 6912 0 diff --git a/cime_config/machines/config_batch.xml b/cime_config/machines/config_batch.xml index d2e45aac9cb3..de69e9e0f8ee 100644 --- a/cime_config/machines/config_batch.xml +++ b/cime_config/machines/config_batch.xml @@ -72,7 +72,7 @@ - + @@ -271,6 +271,15 @@ + + + --gres=gpu:1 + + + gpu + + + acme-small diff --git a/cime_config/machines/config_compilers.xml b/cime_config/machines/config_compilers.xml index a571095b9696..52c9f32ab0d9 100644 --- a/cime_config/machines/config_compilers.xml +++ b/cime_config/machines/config_compilers.xml @@ -886,6 +886,38 @@ flags should be captured within MPAS CMake files. $ENV{PNETCDF_PATH} + + + -DGPU + -DHAVE_NANOTIME -DBIT64 -DHAVE_SLASHPROC -DHAVE_GETTIMEOFDAY + + gpfs + + -O2 -Mvect=nosimd + + + -O2 -Mvect=nosimd -DSUMMITDEV_PGI -DHAVE_IEEE_ARITHMETIC -DNO_R16 + + + -DSUMMITDEV_PGI -DHAVE_IEEE_ARITHMETIC -DNO_R16 + + + -Minline -ta=tesla:ccall,fastmath,loadcache:L1,unroll,fma,managed,deepcopy,nonvvm -Mcuda -Minfo=accel + + + $SHELL{$ENV{NETCDF_FORTRAN_PATH}/bin/nf-config --flibs} -llapack -lblas + $SHELL{$ENV{NETCDF_C_PATH}/bin/nc-config --libs} + + TRUE + FORTRAN + + -lstdc++ + + $ENV{NETCDF_C_PATH} + $ENV{NETCDF_FORTRAN_PATH} + $ENV{PNETCDF_PATH} + + -fopenmp @@ -2018,6 +2050,9 @@ flags should be captured within MPAS CMake files. + + --host=cray + -DARCH_MIC_KNL diff --git a/cime_config/machines/config_machines.xml b/cime_config/machines/config_machines.xml index a867a3fb6ce6..b9ff7819f4e2 100644 --- a/cime_config/machines/config_machines.xml +++ b/cime_config/machines/config_machines.xml @@ -1323,6 +1323,75 @@ + + ANL/LCRC Linux Cluster: 6x 128c EPYC nodes with 8x A100 GPUs + gpulogin.* + LINUX + pgigpu + openmpi + e3sm + /lcrc/group/e3sm + .* + /lcrc/group/e3sm/$USER/scratch/swing + /lcrc/group/e3sm/data/inputdata + /lcrc/group/e3sm/data/inputdata/atm/datm7 + /lcrc/group/e3sm/$USER/archive/$CASE + /lcrc/group/e3sm/baselines/swing/$COMPILER + /lcrc/group/e3sm/soft/tools/cprnc/cprnc + 8 + e3sm_gpu + 4 + slurm + E3SM + 128 + 128 + TRUE + + srun + + -l -n {{ total_tasks }} -N {{ num_nodes }} -K + $SHELL{if [ 128 -ge `./xmlquery --value MAX_MPITASKS_PER_NODE` ]; then echo "--cpu_bind=cores"; else echo "--cpu_bind=threads";fi;} + -c $SHELL{echo 256/ {{ tasks_per_node }} |bc} + -m plane={{ tasks_per_node }} + + + + /gpfs/fs1/soft/swing/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.3.0/lmod-8.3-5tuyfdb/lmod/lmod/init/sh + /gpfs/fs1/soft/swing/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.3.0/lmod-8.3-5tuyfdb/lmod/lmod/init/csh + /gpfs/fs1/soft/swing/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.3.0/lmod-8.3-5tuyfdb/lmod/lmod/init/env_modules_python.py + /gpfs/fs1/soft/swing/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.3.0/lmod-8.3-5tuyfdb/lmod/lmod/libexec/lmod python + module + module + + + cmake/3.21.1-e5i6eks + + + nvhpc/20.9-37zsymt + cuda/11.1.1-nkh7mm7 + openmpi/4.1.1-r6ebr2e + netcdf-c/4.7.4-zppo53l + netcdf-cxx/4.2-wjm7fye + netcdf-fortran/4.5.3-srsajjs + parallel-netcdf/1.12.1-75szceu + + + $CIME_OUTPUT_ROOT/$CASE/run + $CIME_OUTPUT_ROOT/$CASE/bld + 0.1 + 1000 + + $SHELL{dirname $(dirname $(which nc-config))} + $SHELL{dirname $(dirname $(which nf-config))} + $SHELL{dirname $(dirname $(which pnetcdf_version))} + /lcrc/group/e3sm/soft/perl/5.26.0/bin:$ENV{PATH} + + + 64M + cores + + + ANL/LCRC Cluster, Cray CS400, 352-nodes Xeon Phi 7230 KNLs 64C/1.3GHz + 672-nodes Xeon E5-2695v4 Broadwells 36C/2.10GHz, Intel Omni-Path network, SLURM batch system, Lmod module environment. beboplogin.* @@ -1614,7 +1683,7 @@ 0.1 1000 - /projects/ccsm/e3sm/libs/perl5 + /projects/ccsm/e3sm/soft/perl/5.26.0/bin:$ENV{PATH} /projects/ccsm/e3sm/tools/mpas 1 -e PMI_LABEL_ERROUT=1 diff --git a/cime_config/tests.py b/cime_config/tests.py index bbaffaf8faef..7edc16c2b569 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -35,6 +35,7 @@ "ERS.f09_g16.IELMBC", "SMS.r05_r05.I1850ELMCN.elm-qian_1948", "SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-lulcc_sville", + "SMS_Ly2_P1x1.1x1_smallvilleIA.IELMCNCROP.elm-per_crop", "ERS.r05_r05.RMOSGPCC.mosart-gpcc_1972", "ERS.MOS_USRDAT.RMOSGPCC.mosart-mos_usrdat", "SMS.MOS_USRDAT.RMOSGPCC.mosart-unstructure", @@ -166,11 +167,11 @@ "SMS_Ld2.ne30_oECv3.BGCEXP_CNTL_CNPRDCTC_1850.elm-bgcexp", "SMS_D_Ld1.T62_oEC60to30v3.DTESTM", "SMS_D_Ld1.ne30pg2_r05_EC30to60E2r2.A_WCYCL1850S_CMIP6", - "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMFXX", - "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_use_VT", + "ERP_Ln9.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_fixed_subcycle", + "ERS_Ln9.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_use_VT", "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_use_ESMT", "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMFOMP", - "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMF1-RCEMIP", + "ERS_Ln9.ne4pg2_ne4pg2.F-MMF1-RCEMIP", ) }, @@ -206,12 +207,12 @@ "time" : "02:00:00", "tests" : ( # MMF tests - "SMS_D_Ln3_P96x1.ne4pg2_ne4pg2.F-MMF1", - "SMS_Ln3_P96x1.ne4pg2_ne4pg2.F-MMFXX-AQP1", - "SMS_Ln3_P96x1.ne4pg2_ne4pg2.F-MMFXX-RCEMIP", - "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMF1.eam-mmf_crmout", - "ERS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_crmout", - "SMS_Ln9_P96x1.ne4pg2_ne4pg2.F-MMF2-ECPP", + "SMS_D_Ln3.ne4pg2_ne4pg2.F-MMF1", + "SMS_Ln3.ne4pg2_ne4pg2.F-MMFXX-AQP1", + "SMS_Ln3.ne4pg2_ne4pg2.F-MMFXX-RCEMIP", + "ERS_Ln9.ne4pg2_ne4pg2.F-MMF1.eam-mmf_crmout", + "ERS_Ln9.ne4pg2_ne4pg2.F-MMFXX.eam-mmf_crmout", + "SMS_Ln9.ne4pg2_ne4pg2.F-MMF2-ECPP", # non-MMF tests with RRTMGP "ERP_Ln9.ne4pg2_ne4pg2.FC5AV1C-L.eam-rrtmgp", ) diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index 808b0536bfc4..38a8c23dab4c 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -3981,18 +3981,6 @@ Initialize with observed aerosols from IOP file. Default: FALSE - -Assume night time (turn off SW rad). -Default: FALSE - - - -Turn off LW radiation computation. -Default: FALSE - - Turn off microphysics computation. diff --git a/components/eam/cime_config/testdefs/testmods_dirs/eam/mmf_fixed_subcycle/shell_commands b/components/eam/cime_config/testdefs/testmods_dirs/eam/mmf_fixed_subcycle/shell_commands new file mode 100644 index 000000000000..6d08a87ad1cc --- /dev/null +++ b/components/eam/cime_config/testdefs/testmods_dirs/eam/mmf_fixed_subcycle/shell_commands @@ -0,0 +1 @@ +./xmlchange --append -id CAM_CONFIG_OPTS -val " -cppdefs ' -DMMF_DIR_NS -DMMF_FIXED_SUBCYCLE ' " diff --git a/components/eam/src/control/runtime_opts.F90 b/components/eam/src/control/runtime_opts.F90 index 21e317a06d4d..0004ea82a003 100644 --- a/components/eam/src/control/runtime_opts.F90 +++ b/components/eam/src/control/runtime_opts.F90 @@ -173,8 +173,6 @@ module runtime_opts logical :: scm_diurnal_avg logical :: scm_crm_mode logical :: scm_observed_aero -logical :: swrad_off -logical :: lwrad_off logical :: precip_off contains @@ -326,7 +324,7 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in ) namelist /cam_inparm/ iopfile,scm_iop_srf_prop,scm_relaxation, & scm_relaxation_low, scm_relaxation_high, & scm_diurnal_avg,scm_crm_mode,scm_clubb_iop_name, & - scm_observed_aero,swrad_off,lwrad_off, precip_off + scm_observed_aero, precip_off !----------------------------------------------------------------------- @@ -370,8 +368,6 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in ) scm_diurnal_avg_out=scm_diurnal_avg, & scm_crm_mode_out=scm_crm_mode, & scm_observed_aero_out=scm_observed_aero, & - swrad_off_out=swrad_off, & - lwrad_off_out=lwrad_off, & precip_off_out=precip_off, & scm_clubb_iop_name_out=scm_clubb_iop_name) end if @@ -449,8 +445,6 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in ) scm_diurnal_avg_in=scm_diurnal_avg, & scm_crm_mode_in=scm_crm_mode, & scm_observed_aero_in=scm_observed_aero, & - swrad_off_in=swrad_off, & - lwrad_off_in=lwrad_off, & precip_off_in=precip_off, & scm_clubb_iop_name_in=scm_clubb_iop_name) end if diff --git a/components/eam/src/control/scamMod.F90 b/components/eam/src/control/scamMod.F90 index c554b50506d4..1c2c35004e9c 100644 --- a/components/eam/src/control/scamMod.F90 +++ b/components/eam/src/control/scamMod.F90 @@ -183,8 +183,6 @@ module scamMod logical*4, public :: scm_iop_srf_prop ! use the specified surface properties logical*4, public :: scm_relaxation! use relaxation logical*4, public :: scm_observed_aero ! use observed aerosols in SCM file - logical*4, public :: swrad_off ! turn off SW radiation (assume night) - logical*4, public :: lwrad_off ! turn off LW radiation logical*4, public :: precip_off ! turn off precipitation processes logical*4, public :: use_replay ! use e3sm generated forcing logical*4, public :: use_3dfrc ! use 3d forcing @@ -202,10 +200,10 @@ module scamMod subroutine scam_default_opts( scmlat_out,scmlon_out,iopfile_out, & - single_column_out,scm_iop_srf_prop_out, scm_relaxation_out, & - scm_relaxation_low_out, scm_relaxation_high_out, & + single_column_out,scm_iop_srf_prop_out, scm_relaxation_out, & + scm_relaxation_low_out, scm_relaxation_high_out, & scm_diurnal_avg_out, scm_crm_mode_out, scm_observed_aero_out, & - swrad_off_out, lwrad_off_out, precip_off_out, scm_clubb_iop_name_out) + precip_off_out, scm_clubb_iop_name_out) !----------------------------------------------------------------------- real(r8), intent(out), optional :: scmlat_out,scmlon_out character*(max_path_len), intent(out), optional :: iopfile_out @@ -215,8 +213,6 @@ subroutine scam_default_opts( scmlat_out,scmlon_out,iopfile_out, & logical, intent(out), optional :: scm_diurnal_avg_out logical, intent(out), optional :: scm_crm_mode_out logical, intent(out), optional :: scm_observed_aero_out - logical, intent(out), optional :: swrad_off_out - logical, intent(out), optional :: lwrad_off_out logical, intent(out), optional :: precip_off_out real(r8), intent(out), optional :: scm_relaxation_low_out real(r8), intent(out), optional :: scm_relaxation_high_out @@ -233,8 +229,6 @@ subroutine scam_default_opts( scmlat_out,scmlon_out,iopfile_out, & if ( present(scm_diurnal_avg_out) ) scm_diurnal_avg_out = .false. if ( present(scm_crm_mode_out) ) scm_crm_mode_out = .false. if ( present(scm_observed_aero_out)) scm_observed_aero_out = .false. - if ( present(swrad_off_out)) swrad_off_out = .false. - if ( present(lwrad_off_out)) lwrad_off_out = .false. if ( present(precip_off_out)) precip_off_out = .false. if ( present(scm_clubb_iop_name_out) ) scm_clubb_iop_name_out = ' ' @@ -242,9 +236,9 @@ end subroutine scam_default_opts subroutine scam_setopts( scmlat_in, scmlon_in,iopfile_in,single_column_in, & scm_iop_srf_prop_in, scm_relaxation_in, & - scm_relaxation_low_in, scm_relaxation_high_in, & + scm_relaxation_low_in, scm_relaxation_high_in, & scm_diurnal_avg_in, scm_crm_mode_in, scm_observed_aero_in, & - swrad_off_in, lwrad_off_in, precip_off_in, scm_clubb_iop_name_in) + precip_off_in, scm_clubb_iop_name_in) !----------------------------------------------------------------------- real(r8), intent(in), optional :: scmlon_in, scmlat_in character*(max_path_len), intent(in), optional :: iopfile_in @@ -254,8 +248,6 @@ subroutine scam_setopts( scmlat_in, scmlon_in,iopfile_in,single_column_in, & logical, intent(in), optional :: scm_diurnal_avg_in logical, intent(in), optional :: scm_crm_mode_in logical, intent(in), optional :: scm_observed_aero_in - logical, intent(in), optional :: swrad_off_in - logical, intent(in), optional :: lwrad_off_in logical, intent(in), optional :: precip_off_in character(len=*), intent(in), optional :: scm_clubb_iop_name_in real(r8), intent(in), optional :: scm_relaxation_low_in @@ -295,14 +287,6 @@ subroutine scam_setopts( scmlat_in, scmlon_in,iopfile_in,single_column_in, & if (present (scm_observed_aero_in)) then scm_observed_aero=scm_observed_aero_in endif - - if (present (swrad_off_in)) then - swrad_off=swrad_off_in - endif - - if (present (lwrad_off_in)) then - lwrad_off=lwrad_off_in - endif if (present (precip_off_in)) then precip_off=precip_off_in diff --git a/components/eam/src/dynamics/se/stepon.F90 b/components/eam/src/dynamics/se/stepon.F90 index ef83b6550955..44c3b7071566 100644 --- a/components/eam/src/dynamics/se/stepon.F90 +++ b/components/eam/src/dynamics/se/stepon.F90 @@ -14,7 +14,7 @@ module stepon use constituents, only: pcnst, cnst_name, cnst_longname use cam_abortutils, only: endrun use ppgrid, only: begchunk, endchunk - use physconst, only: zvir, cappa + use physconst, only: zvir, cappa, gravit use physics_types, only: physics_state, physics_tend use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_startf, t_stopf, t_barrierf @@ -28,6 +28,7 @@ module stepon use scamMod, only: use_iop, doiopupdate, single_column, & setiopupdate, readiopdata use element_mod, only: element_t + use element_ops, only: get_field, get_field_i use shr_const_mod, only: SHR_CONST_PI implicit none @@ -148,6 +149,10 @@ subroutine stepon_init(dyn_in, dyn_out ) call addfld('DYN_V' ,(/ 'lev' /), 'A', 'm/s', 'Meridional Velocity', gridname='GLL') call addfld('DYN_OMEGA',(/ 'lev' /), 'A', 'Pa/s', 'Vertical Velocity', gridname='GLL' ) call addfld('DYN_PS' ,horiz_only, 'A', 'Pa', 'Surface pressure', gridname='GLL') + call addfld('DYN_PNH' ,(/ 'lev' /), 'A', 'Pa', 'Nonhydrostatic pressure',gridname='GLL') + call addfld('DYN_W' ,(/ 'ilev' /),'A', 'm/s', 'Vertical velocity', gridname='GLL') + call addfld('DYN_Z3' ,(/ 'ilev' /),'A', 'm', 'Geopotential Height (above sea level)', gridname='GLL') + call addfld('DYN_MU' ,(/ 'ilev' /),'A', 'Pa/Pa','dPNH/dPH', gridname='GLL') end subroutine stepon_init @@ -238,6 +243,7 @@ subroutine stepon_run2(phys_state, phys_tend, dyn_in, dyn_out ) real(r8) :: rec2dt real(r8) :: dp(np,np,nlev),fq,fq0,qn0, ftmp(npsq,nlev,2) real(r8) :: tmp_dyn(np,np,nlev,nelemd) + real(r8) :: tmp_dyn_i(np,np,nlevp) real(r8) :: fmtmp(np,np,nlev) real(r8) :: p_m(np,np,nlev) ! temporary midpoint pressure for DYN_OMEGA output real(r8) :: p_i(np,np,nlevp) ! temporary interface pressure for DYN_OMEGA output @@ -412,6 +418,33 @@ subroutine stepon_run2(phys_state, phys_tend, dyn_in, dyn_out ) call outfld('DIV',tmp_dyn(1,1,1,ie),npsq,ie) enddo endif + + if (hist_fld_active('DYN_PNH')) then + do ie=1,nelemd + ! time level ntQ is not used + call get_field(dyn_in%elem(ie),'pnh',tmp_dyn(:,:,:,ie),hvcoord,tl_f,1) + call outfld('DYN_PNH',tmp_dyn(1,1,1,ie),npsq,ie) + enddo + endif + if (hist_fld_active('DYN_Z3')) then + do ie=1,nelemd + call get_field_i(dyn_in%elem(ie),'geo_i',tmp_dyn_i(:,:,:),hvcoord,tl_f) + call outfld('DYN_Z3',tmp_dyn_i(:,:,:) / gravit ,npsq,ie) + enddo + endif + if (hist_fld_active('DYN_W')) then + do ie=1,nelemd + call get_field_i(dyn_in%elem(ie),'w_i',tmp_dyn_i(:,:,:),hvcoord,tl_f) + call outfld('DYN_W',tmp_dyn_i(:,:,:),npsq,ie) + enddo + endif + if (hist_fld_active('DYN_MU')) then + do ie=1,nelemd + call get_field_i(dyn_in%elem(ie),'mu_i',tmp_dyn_i(:,:,:),hvcoord,tl_f) + call outfld('DYN_MU',tmp_dyn_i(:,:,:),npsq,ie) + enddo + endif + if (hist_fld_active('FU') .or. hist_fld_active('FV') ) then do ie=1,nelemd do j=1,np diff --git a/components/eam/src/physics/crm/rrtmg/radiation.F90 b/components/eam/src/physics/crm/rrtmg/radiation.F90 index ece6fe4f382d..adf5c382ba52 100644 --- a/components/eam/src/physics/crm/rrtmg/radiation.F90 +++ b/components/eam/src/physics/crm/rrtmg/radiation.F90 @@ -34,8 +34,7 @@ module radiation use error_messages, only: handle_err use cam_control_mod, only: lambm0, obliqr, mvelpp, eccen use scamMod, only: scm_crm_mode, single_column,have_cld,cldobs,& - have_clwp,clwpobs,have_tg,tground,swrad_off,& - lwrad_off + have_clwp,clwpobs,have_tg,tground use perf_mod, only: t_startf, t_stopf use cam_logfile, only: iulog @@ -324,20 +323,24 @@ function radiation_do(op, timestep) end if select case (op) - - case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always - - case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always - - case default - call endrun('radiation_do: unknown operation:'//op) - + case ('sw') ! do a shortwave heating calc this timestep? + if (iradsw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradsw == 1 & + .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if + case ('lw') ! do a longwave heating calc this timestep? + if (iradlw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradlw == 1 & + .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if + case default + call endrun('radiation_do: unknown operation:'//op) end select end function radiation_do @@ -364,18 +367,17 @@ real(r8) function radiation_nextsw_cday() nstep = get_nstep() dtime = get_step_size() offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') + if (iradsw/=0) then + do while (.not. dosw) + nstep = nstep + 1 + offset = offset + dtime + if (radiation_do('sw', nstep)) then + radiation_nextsw_cday = get_curr_calday(offset=offset) + dosw = .true. + end if + end do end if - + end function radiation_nextsw_cday !================================================================================================ @@ -1330,14 +1332,6 @@ subroutine radiation_tend( state, ptend,pbuf, & calday = get_curr_calday() call zenith (calday, clat, clon, coszrs, ncol, dt_avg) - ! We can bypass the shortwave calculation by setting the cosine of the solar - ! zenith angle to zero for all columns, because the shortwave code collapses - ! the inputs to daytime-only arrays. In case the swrad_off flag is set then, - ! we force shortwave to not be calculated by setting coszrs = 0. - if (swrad_off) then - coszrs(:)=0._r8 - endif - ! Output cosine solar zenith angle call outfld('COSZRS', coszrs(1:ncol), ncol, lchnk) @@ -2106,21 +2100,6 @@ subroutine radiation_tend( state, ptend,pbuf, & clm_seed, lu, ld ) call t_stopf ('rad_rrtmg_lw') - if (lwrad_off) then - qrl(:,:) = 0._r8 - qrlc(:,:) = 0._r8 - flns(:) = 0._r8 - flnt(:) = 0._r8 - flnsc(:) = 0._r8 - flntc(:) = 0._r8 - cam_out%flwds(:) = 0._r8 - flut(:) = 0._r8 - flutc(:) = 0._r8 - fnl(:,:) = 0._r8 - fcnl(:,:) = 0._r8 - fldsc(:) = 0._r8 - end if !lwrad_off - do i=1,ncol lwcf(i)=flutc(i) - flut(i) end do diff --git a/components/eam/src/physics/crm/rrtmgp/radiation.F90 b/components/eam/src/physics/crm/rrtmgp/radiation.F90 index b7f9623f8af8..214a72f10385 100644 --- a/components/eam/src/physics/crm/rrtmgp/radiation.F90 +++ b/components/eam/src/physics/crm/rrtmgp/radiation.F90 @@ -14,7 +14,7 @@ module radiation use shr_kind_mod, only: r8=>shr_kind_r8, cl=>shr_kind_cl use ppgrid, only: pcols, pver, pverp, begchunk, endchunk use cam_abortutils, only: endrun - use scamMod, only: scm_crm_mode, single_column, swrad_off + use scamMod, only: scm_crm_mode, single_column use rad_constituents, only: N_DIAG use radconstants, only: & nswbands, nlwbands, & @@ -364,13 +364,21 @@ function radiation_do(op, timestep) ! and the longwave. In practice, these are probably usually the same. select case (op) case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + if (iradsw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradsw == 1 & + .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + if (iradlw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradlw == 1 & + .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if case default call endrun('radiation_do: unknown operation:'//op) end select @@ -401,18 +409,17 @@ real(r8) function radiation_nextsw_cday() nstep = get_nstep() dtime = get_step_size() offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') + if (iradsw/=0) then + do while (.not. dosw) + nstep = nstep + 1 + offset = offset + dtime + if (radiation_do('sw', nstep)) then + radiation_nextsw_cday = get_curr_calday(offset=offset) + dosw = .true. + end if + end do end if - + end function radiation_nextsw_cday !================================================================================================ @@ -1447,11 +1454,6 @@ subroutine radiation_tend(state_in,ptend, pbuf, cam_out, cam_in, & call set_cosine_solar_zenith_angle(state, dt_avg, coszrs(1:ncol)) call outfld('COSZRS', coszrs(1:ncol), ncol, state%lchnk) - ! If the swrad_off flag is set, meaning we should not do SW radiation, then - ! we just set coszrs to zero everywhere. TODO: why not just set dosw false - ! and skip the loop? - if (swrad_off) coszrs(:) = 0._r8 - ! Gather night/day column indices for subsetting SW inputs; we only want to ! do the shortwave radiative transfer during the daytime to save ! computational cost (and because RRTMGP will fail for cosine solar zenith diff --git a/components/eam/src/physics/crm/sam/MICRO_SAM1MOM/microphysics.F90 b/components/eam/src/physics/crm/sam/MICRO_SAM1MOM/microphysics.F90 index db7d0ebe05bb..2c029ef2d390 100644 --- a/components/eam/src/physics/crm/sam/MICRO_SAM1MOM/microphysics.F90 +++ b/components/eam/src/physics/crm/sam/MICRO_SAM1MOM/microphysics.F90 @@ -507,6 +507,10 @@ subroutine precip_fall(ncrms,hydro_type, omega, ind) nprec = 1 endif +#ifdef MMF_FIXED_SUBCYCLE + nprec = 4 +#endif + ! loop over iterations do iprec = 1,nprec !$acc parallel loop gang vector collapse(4) async(asyncid) diff --git a/components/eam/src/physics/crm/sam/kurant.F90 b/components/eam/src/physics/crm/sam/kurant.F90 index c50ff4a9126b..8dcf6449ffa0 100644 --- a/components/eam/src/physics/crm/sam/kurant.F90 +++ b/components/eam/src/physics/crm/sam/kurant.F90 @@ -63,6 +63,10 @@ subroutine kurant(ncrms) !$acc wait(asyncid) ncycle = max(ncycle,max(1,ceiling(cfl/0.7D0))) +#ifdef MMF_FIXED_SUBCYCLE + ncycle = max_ncycle +#endif + if(ncycle.gt.max_ncycle) then if(masterproc) print *,'kurant() - the number of cycles exceeded max_ncycle = ',max_ncycle do icrm = 1 , ncrms diff --git a/components/eam/src/physics/crm/samomp/MICRO_SAM1MOM/microphysics.F90 b/components/eam/src/physics/crm/samomp/MICRO_SAM1MOM/microphysics.F90 index 7fcf5b7bfb44..9ee18d3cbfdd 100644 --- a/components/eam/src/physics/crm/samomp/MICRO_SAM1MOM/microphysics.F90 +++ b/components/eam/src/physics/crm/samomp/MICRO_SAM1MOM/microphysics.F90 @@ -547,6 +547,10 @@ subroutine precip_fall(ncrms,hydro_type, omega, ind) nprec = 1 endif +#ifdef MMF_FIXED_SUBCYCLE + nprec = 4 +#endif + ! loop over iterations do iprec = 1,nprec !$omp target teams distribute parallel do collapse(4) diff --git a/components/eam/src/physics/crm/samomp/kurant.F90 b/components/eam/src/physics/crm/samomp/kurant.F90 index 7aa59c512b06..2c631af7f1e8 100644 --- a/components/eam/src/physics/crm/samomp/kurant.F90 +++ b/components/eam/src/physics/crm/samomp/kurant.F90 @@ -62,6 +62,10 @@ subroutine kurant(ncrms) !$omp taskwait ncycle = max(ncycle,max(1,ceiling(cfl/0.7D0))) +#ifdef MMF_FIXED_SUBCYCLE + ncycle = max_ncycle +#endif + if(ncycle.gt.max_ncycle) then !$omp target update from(wm) !$omp target update from(uhm) diff --git a/components/eam/src/physics/crm/samxx/kurant.cpp b/components/eam/src/physics/crm/samxx/kurant.cpp index ae0de556b8f9..a3b2ccea99b2 100644 --- a/components/eam/src/physics/crm/samxx/kurant.cpp +++ b/components/eam/src/physics/crm/samxx/kurant.cpp @@ -68,6 +68,10 @@ void kurant () { ncycle = max(ncycle,max(1,static_cast(ceil(cfl/0.7)))); +#ifdef MMF_FIXED_SUBCYCLE + ncycle = max_ncycle; +#endif + if(ncycle > max_ncycle) { std::cout << "\nkurant() - the number of cycles exceeded max_ncycle = "<< max_ncycle << std::endl; exit(-1); diff --git a/components/eam/src/physics/crm/samxx/microphysics.cpp b/components/eam/src/physics/crm/samxx/microphysics.cpp index 53f9b3bbeeb5..8a68fda2fd46 100644 --- a/components/eam/src/physics/crm/samxx/microphysics.cpp +++ b/components/eam/src/physics/crm/samxx/microphysics.cpp @@ -108,6 +108,10 @@ void precip_fall(int hydro_type, real4d &omega) { } //if(nprec > 1){std::cout << nprec << std::endl;} + +#ifdef MMF_FIXED_SUBCYCLE + nprec = 4; +#endif for(int iprec = 1; iprec<=nprec; iprec++) { diff --git a/components/eam/src/physics/rrtmg/radiation.F90 b/components/eam/src/physics/rrtmg/radiation.F90 index f78000ca3ba8..fe5cfb8ce699 100644 --- a/components/eam/src/physics/rrtmg/radiation.F90 +++ b/components/eam/src/physics/rrtmg/radiation.F90 @@ -24,8 +24,7 @@ module radiation use error_messages, only: handle_err use cam_control_mod, only: lambm0, obliqr, mvelpp, eccen use scamMod, only: scm_crm_mode, single_column,have_cld,cldobs,& - have_clwp,clwpobs,have_tg,tground,swrad_off,& - lwrad_off + have_clwp,clwpobs,have_tg,tground use perf_mod, only: t_startf, t_stopf use cam_logfile, only: iulog @@ -304,20 +303,24 @@ function radiation_do(op, timestep) end if select case (op) - - case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always - - case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always - - case default - call endrun('radiation_do: unknown operation:'//op) - + case ('sw') ! do a shortwave heating calc this timestep? + if (iradsw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradsw == 1 & + .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if + case ('lw') ! do a longwave heating calc this timestep? + if (iradlw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradlw == 1 & + .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if + case default + call endrun('radiation_do: unknown operation:'//op) end select end function radiation_do @@ -344,18 +347,17 @@ real(r8) function radiation_nextsw_cday() nstep = get_nstep() dtime = get_step_size() offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') + if (iradsw/=0) then + do while (.not. dosw) + nstep = nstep + 1 + offset = offset + dtime + if (radiation_do('sw', nstep)) then + radiation_nextsw_cday = get_curr_calday(offset=offset) + dosw = .true. + end if + end do end if - + end function radiation_nextsw_cday !================================================================================================ @@ -1111,10 +1113,6 @@ subroutine radiation_tend(state,ptend, pbuf, & call get_rlat_all_p(lchnk, ncol, clat) call get_rlon_all_p(lchnk, ncol, clon) call zenith (calday, clat, clon, coszrs, ncol, dt_avg) - - if (swrad_off) then - coszrs(:)=0._r8 ! coszrs is only output for zenith - endif call output_rad_data( pbuf, state, cam_in, landm, coszrs ) @@ -1464,22 +1462,7 @@ subroutine radiation_tend(state,ptend, pbuf, & clm_seed, lu, ld ) call t_stopf ('rad_rrtmg_lw') - if (lwrad_off) then - qrl(:,:) = 0._r8 - qrlc(:,:) = 0._r8 - flns(:) = 0._r8 - flnt(:) = 0._r8 - flnsc(:) = 0._r8 - flntc(:) = 0._r8 - cam_out%flwds(:) = 0._r8 - flut(:) = 0._r8 - flutc(:) = 0._r8 - fnl(:,:) = 0._r8 - fcnl(:,:) = 0._r8 - fldsc(:) = 0._r8 - end if !lwrad_off - - do i=1,ncol + do i=1,ncol lwcf(i)=flutc(i) - flut(i) end do diff --git a/components/eam/src/physics/rrtmgp/radiation.F90 b/components/eam/src/physics/rrtmgp/radiation.F90 index be5f36952242..f078c8324117 100644 --- a/components/eam/src/physics/rrtmgp/radiation.F90 +++ b/components/eam/src/physics/rrtmgp/radiation.F90 @@ -15,7 +15,7 @@ module radiation use ppgrid, only: pcols, pver, pverp, begchunk, endchunk use cam_abortutils, only: endrun use cam_history_support, only: add_hist_coord - use scamMod, only: scm_crm_mode, single_column, swrad_off + use scamMod, only: scm_crm_mode, single_column use rad_constituents, only: N_DIAG use radconstants, only: nswbands, nlwbands, & get_band_index_sw, get_band_index_lw, test_get_band_index, & @@ -358,13 +358,21 @@ function radiation_do(op, timestep) ! and the longwave. In practice, these are probably usually the same. select case (op) case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + if (iradsw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradsw == 1 & + .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + if (iradlw==0) then + radiation_do = .false. + else + radiation_do = nstep == 0 .or. iradlw == 1 & + .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + end if case default call endrun('radiation_do: unknown operation:'//op) end select @@ -395,18 +403,17 @@ real(r8) function radiation_nextsw_cday() nstep = get_nstep() dtime = get_step_size() offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') + if (iradsw/=0) then + do while (.not. dosw) + nstep = nstep + 1 + offset = offset + dtime + if (radiation_do('sw', nstep)) then + radiation_nextsw_cday = get_curr_calday(offset=offset) + dosw = .true. + end if + end do end if - + end function radiation_nextsw_cday !================================================================================================ @@ -1281,10 +1288,6 @@ subroutine radiation_tend(state, ptend, pbuf, cam_out, cam_in, & ! Get cosine solar zenith angle for current time step. call set_cosine_solar_zenith_angle(state, dt_avg, coszrs(1:ncol)) call outfld('COSZRS', coszrs(1:ncol), ncol, state%lchnk) - ! If the swrad_off flag is set, meaning we should not do SW radiation, then - ! we just set coszrs to zero everywhere. TODO: why not just set dosw false - ! and skip the loop? - if (swrad_off) coszrs(:) = 0._r8 ! Do shortwave cloud optics calculations call t_startf('rad_cld_optics_sw') diff --git a/components/elm/bld/ELMBuildNamelist.pm b/components/elm/bld/ELMBuildNamelist.pm index 19cd56e724fd..545f90d15c35 100755 --- a/components/elm/bld/ELMBuildNamelist.pm +++ b/components/elm/bld/ELMBuildNamelist.pm @@ -1350,7 +1350,7 @@ sub setup_cmdl_maxpft { my $val; my $var = "maxpft"; my %maxpatchpft; - $maxpatchpft{'.true.'} = 25; + $maxpatchpft{'.true.'} = 51; $maxpatchpft{'.false.'} = 17; if ( $opts->{$var} ne "default") { $val = $opts->{$var}; diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml index 2d1bbc17a37b..4f3552de6963 100644 --- a/components/elm/bld/namelist_files/namelist_defaults.xml +++ b/components/elm/bld/namelist_files/namelist_defaults.xml @@ -110,7 +110,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case). for the CLM2 data in the CESM distribution --> lnd/clm2/paramdata/clm_params_c180301.nc -lnd/clm2/paramdata/clm_params_c180301.nc +lnd/clm2/paramdata/clm_params_c210219.nc lnd/clm2/paramdata/clm_params.cbgc.c07292018.nc lnd/clm2/paramdata/clm_params_c180524.nc @@ -296,7 +296,7 @@ lnd/clm2/surfdata_map/surfdata_10x15_mp24_simyr2000_c130927.nc lnd/clm2/surfdata_map/surfdata_1x1_numaIA_mp24_simyr2000_c130927.nc -lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp24_simyr2000_c130927.nc +lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp50_simyr2000_c210617.nc lnd/clm2/surfdata_map/surfdata_5x5_amazon_simyr2000_c130927.nc diff --git a/components/elm/cime_config/config_compsets.xml b/components/elm/cime_config/config_compsets.xml index 90008e45c782..2f7bcca3d645 100644 --- a/components/elm/cime_config/config_compsets.xml +++ b/components/elm/cime_config/config_compsets.xml @@ -94,6 +94,11 @@ 2000_DATM%QIA_ELM%CN-CROP_SICE_SOCN_MOSART_SGLC_SWAV + + ICBELMCNCROP + 1850_SATM_ELM%CN-CROP_SICE_SOCN_MOSART_SGLC_SWAV + + IRCP45ELMBGC RCP4_DATM%QIA_ELM%BGC_SICE_SOCN_MOSART_SGLC_SWAV diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm index 97834bbba99d..a9090338e383 100644 --- a/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm +++ b/components/elm/cime_config/testdefs/testmods_dirs/elm/lulcc_sville/user_nl_elm @@ -1,4 +1,4 @@ -flanduse_timeseries = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/landuse.timeseries_1x1_smallvilleIA_mp24_hist_simyr1850-1855_c141229.nc' -fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp24_simyr2000_c130927.nc' +flanduse_timeseries = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/landuse.timeseries_1x1_smallvilleIA_mp50_hist_simyr1850-1855_c210617.nc' +fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp50_simyr2000_c210617.nc' check_dynpft_consistency = .false. diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm new file mode 100644 index 000000000000..9188c64325a4 --- /dev/null +++ b/components/elm/cime_config/testdefs/testmods_dirs/elm/per_crop/user_nl_elm @@ -0,0 +1,2 @@ +fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1x1_smallvilleIA_mp50_simyr2000_miscanthus_c200606.nc' +do_budgets = .false. diff --git a/components/elm/src/biogeochem/AllocationMod.F90 b/components/elm/src/biogeochem/AllocationMod.F90 index 12064748aaee..6f83850cc7b4 100644 --- a/components/elm/src/biogeochem/AllocationMod.F90 +++ b/components/elm/src/biogeochem/AllocationMod.F90 @@ -321,6 +321,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp use elm_varctl , only : carbonphosphorus_only! use pftvarcon , only: npcropmin, declfact, bfact, aleaff, arootf, astemf, noveg use pftvarcon , only: arooti, fleafi, allconsl, allconss, grperc, grpnow, nsoybean + use pftvarcon , only: percrop use elm_varpar , only: nlevdecomp use elm_varcon , only: nitrif_n2o_loss_frac, secspday ! @@ -612,7 +613,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp if (ivt(p) >= npcropmin) then ! skip 2 generic crops - if (croplive(p)) then + if (croplive(p) .and. percrop(ivt(p)) == 0.0_r8 ) then ! same phases appear in subroutine CropPhenology ! Phase 1 completed: @@ -725,6 +726,22 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp f5 = arepr(p) / aleaf(p) g1 = 0.25_r8 + else if (croplive(p) .and. percrop(ivt(p)) == 1.0_r8) then + arepr(p) = 0._r8 + aroot(p) = max(0._r8, min(1._r8, arooti(ivt(p)) - & + (arooti(ivt(p)) - arootf(ivt(p))) * & + min(1._r8, hui(p)/gddmaturity(p)))) + fleaf = fleafi(ivt(p)) * (exp(-bfact(ivt(p))) - & + exp(-bfact(ivt(p))*hui(p)/gddmaturity(p))) / & ! replacing huigrain with gddmaturity since huigrain does not exist for perennial crops + (exp(-bfact(ivt(p)))-1) ! fraction alloc to leaf (from J Norman alloc curve) + aleaf(p) = max(1.e-5_r8, (1._r8 - aroot(p)) * fleaf) + astem(p) = 1._r8 - arepr(p) - aleaf(p) - aroot(p) + + f1 = aroot(p) / aleaf(p) + f3 = astem(p) / aleaf(p) + f5 = arepr(p) / aleaf(p) + g1 = 0.25_r8 + else ! .not croplive f1 = 0._r8 f3 = 0._r8 diff --git a/components/elm/src/biogeochem/EcosystemDynMod.F90 b/components/elm/src/biogeochem/EcosystemDynMod.F90 index d9f7797c300b..bdcf93e64637 100644 --- a/components/elm/src/biogeochem/EcosystemDynMod.F90 +++ b/components/elm/src/biogeochem/EcosystemDynMod.F90 @@ -465,6 +465,7 @@ end subroutine EcosystemDynNoLeaching1 subroutine EcosystemDynNoLeaching2(bounds, & num_soilc, filter_soilc, & num_soilp, filter_soilp, num_pcropp, filter_pcropp, doalb, & + num_ppercropp, filter_ppercropp, & cnstate_vars, & atm2lnd_vars, & canopystate_vars, soilstate_vars, crop_vars, ch4_vars, & @@ -515,6 +516,8 @@ subroutine EcosystemDynNoLeaching2(bounds, & integer , intent(in) :: filter_soilp(:) ! filter for soil patches integer , intent(in) :: num_pcropp ! number of prog. crop patches in filter integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches + integer , intent(in) :: num_ppercropp ! number of prog perennial crop patches in filter + integer , intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches logical , intent(in) :: doalb ! true = surface albedo calculation time step type(cnstate_type) , intent(inout) :: cnstate_vars type(atm2lnd_type) , intent(in) :: atm2lnd_vars @@ -579,7 +582,7 @@ subroutine EcosystemDynNoLeaching2(bounds, & event = 'Phenology' call t_start_lnd(event) call Phenology(num_soilc, filter_soilc, num_soilp, filter_soilp, & - num_pcropp, filter_pcropp, doalb, atm2lnd_vars, & + num_pcropp, filter_pcropp, num_ppercropp, filter_ppercropp, doalb, atm2lnd_vars, & crop_vars, canopystate_vars, soilstate_vars, & cnstate_vars ) call t_stop_lnd(event) diff --git a/components/elm/src/biogeochem/PhenologyMod.F90 b/components/elm/src/biogeochem/PhenologyMod.F90 index f331414e2061..bf46863084ab 100644 --- a/components/elm/src/biogeochem/PhenologyMod.F90 +++ b/components/elm/src/biogeochem/PhenologyMod.F90 @@ -32,6 +32,12 @@ module PhenologyMod use VegetationType , only : veg_pp use VegetationDataType , only : veg_es, veg_ef, veg_cs, veg_cf, veg_ns, veg_nf use VegetationDataType , only : veg_ps, veg_pf + use CNCarbonStateType , only : carbonstate_type + use CNCarbonFluxType , only : carbonflux_type + use CNNitrogenStateType , only : nitrogenstate_type + use CNNitrogenFluxType , only : nitrogenflux_type + use PhosphorusStateType , only : phosphorusstate_type + use PhosphorusFluxType , only : phosphorusflux_type !!!Added for gpu timing info use timeinfoMod @@ -252,7 +258,7 @@ end subroutine readPhenolParams !----------------------------------------------------------------------- subroutine Phenology (num_soilc, filter_soilc, num_soilp, filter_soilp, & - num_pcropp, filter_pcropp, doalb, atm2lnd_vars, & + num_pcropp, filter_pcropp, num_ppercropp, filter_ppercropp, doalb, atm2lnd_vars, & crop_vars, canopystate_vars, soilstate_vars, & cnstate_vars) ! @@ -268,6 +274,8 @@ subroutine Phenology (num_soilc, filter_soilc, num_soilp, filter_soilp, & integer , intent(in) :: filter_soilp(:) ! filter for soil patches integer , intent(in) :: num_pcropp ! number of prog. crop patches in filter integer , intent(in) :: filter_pcropp(:)! filter for prognostic crop patches + integer , intent(in) :: num_ppercropp ! number of prog perennial crop patches in filter + integer , intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches logical , intent(in) :: doalb ! true if time for sfc albedo calc type(atm2lnd_type) , intent(in) :: atm2lnd_vars type(crop_type) , intent(inout) :: crop_vars @@ -302,6 +310,11 @@ subroutine Phenology (num_soilc, filter_soilc, num_soilp, filter_soilp, & ) end if + if (doalb .and. num_ppercropp > 0 ) then + call PerennialCropPhenology(num_ppercropp, filter_ppercropp, & + crop_vars, cnstate_vars) + end if + ! the same onset and offset routines are called regardless of ! phenology type - they depend only on onset_flag, offset_flag, bglfr, and bgtr @@ -314,6 +327,11 @@ subroutine Phenology (num_soilc, filter_soilc, num_soilp, filter_soilp, & cnstate_vars) end if + if (num_ppercropp > 0 ) then + call CNPerennialCropHarvest(num_ppercropp, filter_ppercropp, & + num_soilc, filter_soilc, crop_vars, cnstate_vars) + end if + call CNOffsetLitterfall(num_soilp, filter_soilp, & cnstate_vars) @@ -1712,7 +1730,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & harvdate(p) = NOT_Harvested if (ivt(p)==nsoybean .or. ivt(p) == nsoybeanirrig) gddmaturity(p)=min(gdd1020(p),hybgdd(ivt(p))) - if (ivt(p)==ncorn .or. ivt(p)==ncornirrig) gddmaturity(p)=max(950._r8, min(gdd820(p)*0.85_r8, hybgdd(ivt(p)))) + if (ivt(p)==ncorn .or. ivt(p)==ncornirrig) then + gddmaturity(p)=max(950._r8, min(gdd820(p)*0.85_r8, hybgdd(ivt(p)))) + end if if (ivt(p)==nscereal .or. ivt(p) == nscerealirrig) gddmaturity(p)=min(gdd020(p),hybgdd(ivt(p))) leafc_xfer(p) = 1._r8 ! initial seed at planting to appear @@ -1910,6 +1930,215 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & end subroutine CropPhenology +!----------------------------------------------------------------------- + subroutine PerennialCropPhenology(num_ppercropp, filter_ppercropp, & + crop_vars, cnstate_vars) + + ! !DESCRIPTION: + ! Code based on ORCHIDEE-MICT-BIOENERGY model (Li et al., 2018) + ! !USES: + use shr_const_mod , only : SHR_CONST_TKFRZ + use clm_time_manager , only : get_curr_calday, get_days_per_year + use pftvarcon , only : gddmin, hybgdd + use pftvarcon , only : minplanttemp, planttemp, senestemp, min_days_senes + use elm_varcon , only : spval, secspday + ! + ! !ARGUMENTS: + integer , intent(in) :: num_ppercropp ! number of prog perennial crop patches in filter + integer , intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches + type(crop_type) , intent(inout) :: crop_vars + type(cnstate_type) , intent(inout) :: cnstate_vars + ! + ! LOCAL VARAIBLES: + integer jday ! julian day of the year + integer fp,p ! patch indices + integer c ! column indices + integer g ! gridcell indices + integer h ! hemisphere indices + integer t ! topographic indices + real(r8) dayspyr ! days per year + real(r8) ndays_on ! number of days to fertilize + real(r8):: soilt + real(r8):: dt + !------------------------------------------------------------------------ + + associate( & + ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type + + t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) + + leaf_long => veg_vp%leaf_long , & ! Input: [real(r8) (:) ] leaf longevity (yrs) + froot_long => veg_vp%froot_long , & ! Input: [real(r8) (:) ] fine root longevity (yrs) + leafcn => veg_vp%leafcn , & ! Input: [real(r8) (:) ] leaf C:N (gC/gN) + leafcp => veg_vp%leafcp , & ! Input: [real(r8) (:) ] leaf C:P (gC/gP) + fertnitro => veg_vp%fertnitro , & ! Input: [real(r8) (:) ] max fertilizer to be applied in total (kgN/m2) + + t10 => veg_es%t_a10 , & ! Input: [real(r8) (:) ] 10-day running mean of the 2 m temperature (K) + a10tmin => veg_es%t_a10min , & ! Input: [real(r8) (:) ] 10-day running mean of min 2-m temperature + + harvdate => crop_vars%harvdate_patch , & ! Output: [integer (:) ] harvest date + harvday => crop_vars%harvday_patch , & ! Ouptut: [real(r8) ):) ] harvest day + croplive => crop_vars%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested + crpyld => crop_vars%crpyld_patch , & ! Output: [real(r8) ):) ] harvested crop (bu/acre) + dmyield => crop_vars%dmyield_patch , & ! Output: [real(r8) ):) ] dry matter harvested crop (t/ha) + + bglfr_leaf => cnstate_vars%bglfr_leaf_patch , & ! Output: [real(r8) (:) ] background leaf litterfall rate (1/s) + bglfr_froot => cnstate_vars%bglfr_froot_patch , & ! Output: [real(r8) (:) ] background fine root litterfall rate (1/s) + bgtr => cnstate_vars%bgtr_patch , & ! Output: [real(r8) (:) ] background transfer growth rate (1/s) + lgsf => cnstate_vars%lgsf_patch , & ! Output: [real(r8) (:) ] long growing season factor [0-1] + onset_flag => cnstate_vars%onset_flag_patch , & ! Output: [real(r8) (:) ] onset flag + onset_counter => cnstate_vars%onset_counter_patch , & ! Output: [real(r8) (:) ] onset counter + onset_gddflag => cnstate_vars%onset_gddflag_patch , & ! Output: [real(r8) (:) ] onset growing degree days flag + onset_gdd => cnstate_vars%onset_gdd_patch , & ! Output: [real(r8) (:) ] onset growing degree days + offset_flag => cnstate_vars%offset_flag_patch , & ! Output: [real(r8) (:) ] offset flag + offset_counter => cnstate_vars%offset_counter_patch , & ! Output: [real(r8) (:) ] offset counter + gddmaturity => cnstate_vars%gddmaturity_patch , & ! Output: [real(r8) (:) ] gdd needed to harvest + + leafc_xfer => veg_cs%leafc_xfer , & ! Output: [real(r8) (:) ] (gC/m2) leaf C transfer + + crop_seedc_to_leaf => veg_cf%crop_seedc_to_leaf , & ! Output: [real(r8) (:) ] (gC/m2/s) seed source to PFT-level + + fert => veg_nf%fert , & ! Output: [real(r8) (:) ] (gN/m2/s) fertilizer applied each timestep + fert_counter => veg_nf%fert_counter , & ! Output: [real(r8) (:) ] >0 fertilize; <=0 not (seconds) + + leafn_xfer => veg_ns%leafn_xfer , & ! Output: [real(r8) (:) ] (gN/m2) leaf N transfer + + crop_seedn_to_leaf => veg_nf%crop_seedn_to_leaf , & ! Output: [real(r8) (:) ] (gN/m2/s) seed source to PFT-level + + leafp_xfer => veg_ps%leafp_xfer , & ! Output: [real(r8) (:) ] (gP/m2) leaf P transfer + + crop_seedp_to_leaf => veg_pf%crop_seedp_to_leaf & ! Output [real(r8) (:) ] (gP/m2/s) seed source to PFT-level + ) + + ! get time info + dt = dtime_mod + dayspyr = get_days_per_year() + jday = get_curr_calday() + + ndays_on = 20._r8 ! number of days to fertilize + + do fp = 1, num_ppercropp ! prognostic perennial crops loop + p = filter_ppercropp(fp) + c = veg_pp%column(p) + g = veg_pp%gridcell(p) + t = veg_pp%topounit(p) + h = inhemi(p) + + ! background litterfall and transfer rates; long growing season factor + bglfr_leaf(p) = 0._r8 ! this value changes later in a crop's life cycle + bglfr_froot(p) = 0._r8 ! this value changes later in a crop's life cycle + bgtr(p) = 0._r8 + lgsf(p) = 0._r8 + + ! B.Drewniak - zero our yield calculator + crpyld(p) = 0._r8 + dmyield(p) = 0._r8 + + if (.not. croplive(p)) then + + ! Perennial crop is planted when temperature conditions are met + if ( t10(p) /= spval .and. a10tmin(p) /= spval .and. & + t10(p) > planttemp(ivt(p)) .and. & + a10tmin(p) > minplanttemp(ivt(p)) .and. & + jday >= minplantjday(ivt(p),h) .and. & + jday <= maxplantjday(ivt(p),h) ) then + croplive(p) = .true. + gddmaturity(p) = hybgdd(ivt(p)) + offset_flag(p) = 1._r8 + offset_counter(p) = dt + harvdate(p) = NOT_Harvested + harvday(p) = NOT_Harvested + + leafc_xfer(p) = 1._r8 ! initial seed at planting to appear + leafn_xfer(p) = leafc_xfer(p) / leafcn(ivt(p)) + leafp_xfer(p) = leafc_xfer(p) / leafcp(ivt(p)) + + crop_seedc_to_leaf(p) = leafc_xfer(p)/dt + crop_seedn_to_leaf(p) = leafn_xfer(p)/dt + crop_seedp_to_leaf(p) = leafp_xfer(p)/dt + end if + end if ! crop not live + + if (croplive(p)) then + + ! update onset_counter and test for the end of the onset period + if (onset_flag(p) == 1.0_r8) then + ! decrement counter for onset period + onset_counter(p) = onset_counter(p) - dt + + ! fertilizer counter + if (fert_counter(p) <= 0._r8) then + fert(p) = 0._r8 + else ! continue same fert application every timestep + fert_counter(p) = fert_counter(p) - dt + end if + + ! if end of the onset period, reset phenology flags and indices + if (onset_counter(p) <= 0.0_r8 .and. & + (t10(p) /= spval .and. t10(p) < senestemp(ivt(p)))) then + + onset_flag(p) = 0.0_r8 + onset_counter(p) = 0.0_r8 + offset_flag(p) = 1._r8 + offset_counter(p) = dt + if (harvdate(p) >= NOT_Harvested) harvdate(p) = jday + if (harvday(p) >= NOT_Harvested) harvday(p) = jday + croplive(p) = .false. + bglfr_leaf(p) = 1._r8/(leaf_long(ivt(p)) * dayspyr * secspday) + bglfr_froot(p) = 1._r8/(froot_long(ivt(p)) * dayspyr * secspday) + + end if + end if ! onset flag + + ! test for switching from senescence period to growth period + if (offset_flag(p) == 1.0_r8) then + ! Test to turn on growing degree-day sum, if off + if (onset_gddflag(p) == 0._r8) then + onset_gddflag(p) = 1._r8 + onset_gdd(p) = 0._r8 + end if + + ! if the gdd flag is set, and if the soil is above freezing + ! then accumulate growing degree days for onset trigger + soilt = t_soisno(c,3) + if (onset_gddflag(p) == 1.0_r8 .and. soilt > SHR_CONST_TKFRZ) then + onset_gdd(p) = onset_gdd(p) + (soilt-SHR_CONST_TKFRZ)*fracday + end if + + ! if accumulated gdd has exceeded gdd required for leaf onset + ! then set onset flag + ! for perennial crops gddmin tracks gdd required for leaf onset + if (onset_gdd(p) > gddmin(ivt(p))) then + + ! Initialize onset counter + onset_flag(p) = 1.0_r8 + offset_flag(p) = 0._r8 + harvdate(p) = NOT_Harvested + harvday(p) = NOT_Harvested + onset_gddflag(p) = 0.0_r8 + onset_gdd(p) = 0.0_r8 + onset_counter(p) = min_days_senes(ivt(p)) * secspday + fert_counter(p) = ndays_on * secspday + fert(p) = fertnitro(ivt(p)) * 1000._r8 / fert_counter(p) + end if + end if ! offset flag + else ! crop not live + ! next 2 lines conserve mass if leaf*_xfer > 0 due to interpinic + crop_seedc_to_leaf(p) = crop_seedc_to_leaf(p) - leafc_xfer(p)/dt + crop_seedn_to_leaf(p) = crop_seedn_to_leaf(p) - leafn_xfer(p)/dt + crop_seedp_to_leaf(p) = crop_seedp_to_leaf(p) - leafp_xfer(p)/dt + onset_counter(p) = 0._r8 + leafc_xfer(p) = 0._r8 + leafn_xfer(p) = leafc_xfer(p) / leafcn(ivt(p)) + leafp_xfer(p) = leafc_xfer(p) / leafcp(ivt(p)) + end if ! croplive + + end do ! prognostic perennial crops loop + + end associate + + end subroutine PerennialCropPhenology + !----------------------------------------------------------------------- subroutine CropPhenologyInit(bounds) ! @@ -1918,7 +2147,7 @@ subroutine CropPhenologyInit(bounds) ! initialized, and after ecophyscon file is read in. ! ! !USES: - use pftvarcon , only: npcropmin, npcropmax, mnNHplantdate + use pftvarcon , only: npcropmin, npcropmax, nppercropmin, nppercropmax, mnNHplantdate use pftvarcon , only: mnSHplantdate, mxNHplantdate use pftvarcon , only: mxSHplantdate use clm_time_manager, only: get_calday @@ -1953,6 +2182,15 @@ subroutine CropPhenologyInit(bounds) maxplantjday(n,inSH) = int( get_calday( mxSHplantdate(n), 0 ) ) end do + do n = nppercropmin, nppercropmax + minplantjday(n,inNH) = int( get_calday( mnNHplantdate(n), 0 ) ) + maxplantjday(n,inNH) = int( get_calday( mxNHplantdate(n), 0 ) ) + end do + do n = nppercropmin, nppercropmax + minplantjday(n,inSH) = int( get_calday( mnSHplantdate(n), 0 ) ) + maxplantjday(n,inSH) = int( get_calday( mxSHplantdate(n), 0 ) ) + end do + ! Figure out what hemisphere each PFT is in do p = bounds%begp, bounds%endp g = veg_pp%gridcell(p) @@ -2321,6 +2559,8 @@ subroutine CNOnsetGrowth (num_soilp, filter_soilp, & ! Determines the flux of stored C and N from transfer pools to display ! pools during the phenological onset period. ! add flux for phosphorus - X.YANG + ! !USES: + use pftvarcon , only : percrop ! ! !ARGUMENTS: !$acc routine seq @@ -2401,7 +2641,7 @@ subroutine CNOnsetGrowth (num_soilp, filter_soilp, & ! The transfer rate is a linearly decreasing function of time, ! going to zero on the last timestep of the onset period - if (onset_counter(p) == dt) then + if (onset_counter(p) == dt .or. percrop(ivt(p)) == 1.0_r8) then t1 = 1.0_r8 / dt else t1 = 2.0_r8 / (onset_counter(p)) @@ -2575,6 +2815,102 @@ subroutine CNCropHarvest (num_pcropp, filter_pcropp,& end associate end subroutine CNCropHarvest +!---------------------------------------------------------------------- + subroutine CNPerennialCropHarvest (num_ppercropp, filter_ppercropp, num_soilc, filter_soilc, & + crop_vars, cnstate_vars) + ! + ! !DESCRIPTION: + ! This routine handles harvest for agriculture vegetation types, such as + ! corn, soybean, and wheat. This routine allows harvest to be calculated + ! instead of in the OffsetLitterfall subroutine. The harvest index is + ! determined based on the LPJ model. + ! + ! !ARGUMENTS: + integer, intent(in) :: num_ppercropp ! number of prog perennial crop patches in filter + integer, intent(in) :: filter_ppercropp(:) ! filter for prognostic perennial crop patches + integer, intent(in) :: num_soilc ! number of soil columns in filter + integer, intent(in) :: filter_soilc(:) ! soil column filter + + type(crop_type) , intent(inout) :: crop_vars + type(cnstate_type) , intent(inout) :: cnstate_vars + ! + ! !LOCAL VARIABLES: + ! local pointers to implicit in scalars + integer :: p ! indices + integer :: fp ! lake filter pft index + real(r8):: t1 ! temporary variable + real(r8):: cgrain ! amount of carbon in the grain + real(r8):: dt + !------------------------------------------------------------------------- + associate(& + ivt => veg_pp%itype , & ! Input: [integer (:)] pft vegetation type + offset_flag => cnstate_vars%offset_flag_patch , & ! Input: [real(r8) (:)] offset flag + offset_counter => cnstate_vars%offset_counter_patch, & ! Input: [real(r8) (:)] offset days counter + + presharv => veg_vp%presharv , & ! Input: [real(r8) (:)] proportion of residue harvested + fyield => veg_vp%fyield , & ! Input: [real(r8) (:)] fraction of grain actually harvested + convfact => veg_vp%convfact , & ! Input: [real(r8) (:)] conversation factor from gC/m2 to bu/acre + + leafc => veg_cs%leafc , & ! Input: [real(r8) (:)] leaf C (gC/m2) + livestemc => veg_cs%livestemc , & ! Input: [real(r8) (:)] livestem C (gC/m2) + leafn => veg_ns%leafn , & ! Input: [real(r8) (:)] leaf N (gN/m2) + livestemn => veg_ns%livestemn , & ! Input: [real(r8) (:)] livestem N (gN/m2) + leafp => veg_ps%leafp , & ! Input: [real(r8) (:)] leaf P (gP/m2) + livestemp => veg_ps%livestemp , & ! Input: [real(r8) (:)] livestem P (gP/m2) + cpool_to_livestemc => veg_cf%cpool_to_livestemc , & ! Input: [real(r8) (:)] allocation to live stem C (gC/m2/s) + cpool_to_leafc => veg_cf%cpool_to_leafc , & ! Input: [real(r8) (:)] allocation to leaf C (gC/m2/s) + npool_to_leafn => veg_nf%npool_to_leafn , & ! Input: [real(r8) (:)] allocation to leaf N (gN/m2/s) + npool_to_livestemn => veg_nf%npool_to_livestemn , & ! Input: [real(r8) (:)] allocation to live stem N (gN/m2/s) + ppool_to_leafp => veg_pf%ppool_to_leafp , & ! Input: [real(r8) (:)] allocation to leaf P (gP/m2/s) + ppool_to_livestemp => veg_pf%ppool_to_livestemp , & ! Input: [real(r8) (:)] allocation to live stem P (gP/m2/s) + hrv_leafc_to_prod1c => veg_cf%hrv_leafc_to_prod1c , & ! Input: [real(r8) (:)] crop leaf C harvested + hrv_livestemc_to_prod1c => veg_cf%hrv_livestemc_to_prod1c , & ! Input: [real(r8) (:)] crop live stem C harvested + hrv_leafn_to_prod1n => veg_nf%hrv_leafn_to_prod1n , & ! Input: [real(r8) (:)] crop leaf N harvested + hrv_livestemn_to_prod1n => veg_nf%hrv_livestemn_to_prod1n , & ! Input: [real(r8) (:)] crop live stem N harvested + hrv_leafp_to_prod1p => veg_pf%hrv_leafp_to_prod1p , & ! Input: [real(r8) (:)] crop leaf P harvested + hrv_livestemp_to_prod1p => veg_pf%hrv_livestemp_to_prod1p , & ! Input: [real(r8) (:)] crop live stem P harvested + crpyld => crop_vars%crpyld_patch , & ! InOut: [real(r8) ):)] harvested crop (bu/acre) + dmyield => crop_vars%dmyield_patch & ! InOut: [real(r8) ):)] dry matter harvested crop (t/ha) + ) + + cgrain = 0.50_r8 + do fp = 1, num_ppercropp + p = filter_ppercropp(fp) + ! only calculate during the offset period + if (offset_flag(p) == 1._r8) then + + if (offset_counter(p) == dt) then + t1 = 1._r8 / dt + ! calculate yield (crpyld = bu/acre and dmyield = t/ha) + ! for perennial bioenergy grass yield comes from leaf and stem harvest + crpyld(p) = presharv(ivt(p)) * (leafc(p) + cpool_to_leafc(p)*dt + livestemc(p) + cpool_to_livestemc(p)*dt) * fyield(ivt(p)) * convfact(ivt(p)) / (cgrain * 1000) + dmyield(p) = presharv(ivt(p)) * (leafc(p) + cpool_to_leafc(p)*dt + livestemc(p) + cpool_to_livestemc(p)*dt) * fyield(ivt(p)) * 0.01 / cgrain + + !calculate harvested carbon and nitrogen; remaining goes into litterpool + hrv_leafc_to_prod1c(p) = presharv(ivt(p)) * ((t1 * leafc(p)) + cpool_to_leafc(p)) + hrv_livestemc_to_prod1c(p) = presharv(ivt(p)) * ((t1 * livestemc(p)) + cpool_to_livestemc(p)) + + ! Do the same for Nitrogen + hrv_leafn_to_prod1n(p) = presharv(ivt(p)) * ((t1 * leafn(p)) + npool_to_leafn(p)) + hrv_livestemn_to_prod1n(p) = presharv(ivt(p)) * ((t1 * livestemn(p)) + npool_to_livestemn(p)) + + ! Do the same for Phosphorus + hrv_leafp_to_prod1p(p) = presharv(ivt(p)) * ((t1 * leafp(p)) + ppool_to_leafp(p)) + hrv_livestemp_to_prod1p(p) = presharv(ivt(p)) * ((t1 * livestemp(p)) + ppool_to_livestemp(p)) + + end if ! offseddt_counter + + end if ! offset_flag + end do + + ! gather all pft-level fluxes from harvest to the column + ! for C and N inputs + + call CNCropHarvestPftToColumn(num_soilc, filter_soilc, cnstate_vars) + + end associate + end subroutine CNPerennialCropHarvest + !----------------------------------------------------------------------- subroutine CNOffsetLitterfall (num_soilp, filter_soilp, & cnstate_vars) diff --git a/components/elm/src/biogeochem/VegStructUpdateMod.F90 b/components/elm/src/biogeochem/VegStructUpdateMod.F90 index 3f2ea15476d8..706e89daab1b 100644 --- a/components/elm/src/biogeochem/VegStructUpdateMod.F90 +++ b/components/elm/src/biogeochem/VegStructUpdateMod.F90 @@ -39,6 +39,7 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & ! !USES: use pftvarcon , only : noveg, nc3crop, nc3irrig, nbrdlf_evr_shrub, nbrdlf_dcd_brl_shrub use pftvarcon , only : ncorn, ncornirrig, npcropmin, ztopmx, laimx + use pftvarcon , only : nmiscanthus, nmiscanthusirrig, nswitchgrass, nswitchgrassirrig use clm_time_manager , only : get_rad_step_size use elm_varctl , only : spinup_state, spinup_mortality_factor ! @@ -197,7 +198,9 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & if (tlai(p) >= laimx(ivt(p))) peaklai(p) = 1 ! used in CNAllocation - if (ivt(p) == ncorn .or. ivt(p) == ncornirrig) then + if (ivt(p) == ncorn .or. ivt(p) == ncornirrig .or. & + ivt(p) == nmiscanthus .or. ivt(p) == nmiscanthusirrig .or. & + ivt(p) == nswitchgrass .or. ivt(p) == nswitchgrassirrig) then tsai(p) = 0.1_r8 * tlai(p) else tsai(p) = 0.2_r8 * tlai(p) diff --git a/components/elm/src/external_models/emi/src/elm_stub/utils/clm_varpar.F90 b/components/elm/src/external_models/emi/src/elm_stub/utils/clm_varpar.F90 index e4220fa94fc2..bc39d27b24bf 100644 --- a/components/elm/src/external_models/emi/src/elm_stub/utils/clm_varpar.F90 +++ b/components/elm/src/external_models/emi/src/elm_stub/utils/clm_varpar.F90 @@ -45,17 +45,19 @@ module clm_varpar integer, parameter :: ndst = 4 ! number of dust size classes (BGC only) integer, parameter :: dst_src_nbr = 3 ! number of size distns in src soil (BGC only) integer, parameter :: sz_nbr = 200 ! number of sub-grid bins in large bin of dust size distribution (BGC only) - integer, parameter :: mxpft = 24 ! maximum number of PFT's for any mode; + integer, parameter :: mxpft = 50 ! maximum number of PFT's for any mode; ! FIX(RF,032414) might we set some of these automatically from reading pft-physiology? integer, parameter :: numveg = 16 ! number of veg types (without specific crop) integer, parameter :: nlayer = 3 ! number of VIC soil layer --Added by AWang integer :: nlayert ! number of VIC soil layer + 3 lower thermal layers integer :: numpft = mxpft ! actual # of patches (without bare) - integer :: numcft = 10 ! actual # of crops + integer :: numcft = 36 ! actual # of crops logical :: crop_prog = .true. ! If prognostic crops is turned on integer :: maxpatch_urb= 5 ! max number of urban patches (columns) in urban landunit + integer :: mxpft_nc ! maximum number of PFT's when use_crop=False; + integer :: maxpatch_pft ! max number of plant functional types in naturally vegetated landunit (namelist setting) integer, parameter :: nsoilorder = 15 ! number of soil orders @@ -109,11 +111,12 @@ subroutine clm_varpar_init() if (use_crop) then numpft = mxpft ! actual # of patches (without bare) - numcft = 10 ! actual # of crops + numcft = 36 ! actual # of crops crop_prog = .true. ! If prognostic crops is turned on else numpft = numveg ! actual # of patches (without bare) numcft = 2 ! actual # of crops + mxpft_nc = 24 ! maximum number of PFT's when use_crop=False; crop_prog = .false. ! If prognostic crops is turned on end if diff --git a/components/elm/src/external_models/fates b/components/elm/src/external_models/fates index 36e60bb31b33..6c3e30ebe633 160000 --- a/components/elm/src/external_models/fates +++ b/components/elm/src/external_models/fates @@ -1 +1 @@ -Subproject commit 36e60bb31b330ab7d2d75066174f4a140a39cfec +Subproject commit 6c3e30ebe6335dd5ac891dfb567d5749f17c7cdd diff --git a/components/elm/src/main/elm_driver.F90 b/components/elm/src/main/elm_driver.F90 index f8e3de0a647e..5bc3f5d05b00 100644 --- a/components/elm/src/main/elm_driver.F90 +++ b/components/elm/src/main/elm_driver.F90 @@ -1078,11 +1078,11 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) end if !if (use_elm_interface) !-------------------------------------------------------------------------------- - call EcosystemDynNoLeaching2(bounds_clump, & filter(nc)%num_soilc, filter(nc)%soilc, & filter(nc)%num_soilp, filter(nc)%soilp, & filter(nc)%num_pcropp, filter(nc)%pcropp, doalb, & + filter(nc)%num_ppercropp, filter(nc)%ppercropp, & cnstate_vars, atm2lnd_vars, & canopystate_vars, soilstate_vars, crop_vars, ch4_vars, & photosyns_vars, soilhydrology_vars, energyflux_vars, & diff --git a/components/elm/src/main/elm_varpar.F90 b/components/elm/src/main/elm_varpar.F90 index 8fb069d0de8f..0b1a1974691c 100644 --- a/components/elm/src/main/elm_varpar.F90 +++ b/components/elm/src/main/elm_varpar.F90 @@ -47,17 +47,19 @@ module elm_varpar integer, parameter :: ndst = 4 ! number of dust size classes (BGC only) integer, parameter :: dst_src_nbr = 3 ! number of size distns in src soil (BGC only) integer, parameter :: sz_nbr = 200 ! number of sub-grid bins in large bin of dust size distribution (BGC only) - integer, parameter :: mxpft = 24 ! maximum number of PFT's for any mode; + integer, parameter :: mxpft = 50 ! maximum number of PFT's for any mode; ! FIX(RF,032414) might we set some of these automatically from reading pft-physiology? integer, parameter :: numveg = 16 ! number of veg types (without specific crop) integer, parameter :: nlayer = 3 ! number of VIC soil layer --Added by AWang integer :: nlayert ! number of VIC soil layer + 3 lower thermal layers integer :: numpft = mxpft ! actual # of patches (without bare) - integer :: numcft = 10 ! actual # of crops + integer :: numcft = 36 ! actual # of crops logical :: crop_prog = .true. ! If prognostic crops is turned on integer :: maxpatch_urb= 5 ! max number of urban patches (columns) in urban landunit + integer :: mxpft_nc ! maximum number of PFT's when use_crop=False; + integer :: maxpatch_pft ! max number of plant functional types in naturally vegetated landunit (namelist setting) integer, parameter :: nsoilorder = 15 ! number of soil orders @@ -113,11 +115,12 @@ subroutine elm_varpar_init() if (use_crop) then numpft = mxpft ! actual # of patches (without bare) - numcft = 10 ! actual # of crops + numcft = 36 ! actual # of crops crop_prog = .true. ! If prognostic crops is turned on else numpft = numveg ! actual # of patches (without bare) numcft = 2 ! actual # of crops + mxpft_nc = 24 ! maximum number of PFT's when use_crop=False; crop_prog = .false. ! If prognostic crops is turned on end if diff --git a/components/elm/src/main/filterMod.F90 b/components/elm/src/main/filterMod.F90 index 318bc40928ae..6d6c2bca1177 100644 --- a/components/elm/src/main/filterMod.F90 +++ b/components/elm/src/main/filterMod.F90 @@ -30,6 +30,8 @@ module filterMod integer, pointer :: pcropp(:) ! prognostic crop filter (pfts) integer :: num_pcropp ! number of pfts in prognostic crop filter + integer, pointer :: ppercropp(:) ! prognostic perennial crop filter (pfts) + integer :: num_ppercropp ! number of pfts in prognostic perennial crop filter integer, pointer :: soilnopcropp(:) ! soil w/o prog. crops (pfts) integer :: num_soilnopcropp ! number of pfts in soil w/o prog crops @@ -210,6 +212,7 @@ subroutine allocFiltersOneGroup(this_filter) allocate(this_filter(nc)%nourbanl(bounds%endl-bounds%begl+1)) allocate(this_filter(nc)%pcropp(bounds%endp-bounds%begp+1)) + allocate(this_filter(nc)%ppercropp(bounds%endp-bounds%begp+1)) allocate(this_filter(nc)%soilnopcropp(bounds%endp-bounds%begp+1)) allocate(this_filter(nc)%icemecc(bounds%endc-bounds%begc+1)) @@ -268,7 +271,7 @@ subroutine setFiltersOneGroup(bounds, this_filter, include_inactive, icemask_grc ! ! !USES: use decompMod , only : BOUNDS_LEVEL_CLUMP - use pftvarcon , only : npcropmin + use pftvarcon , only : npcropmin, nppercropmin use landunit_varcon, only : istsoil, istcrop, istice_mec use column_varcon, only : icol_road_perv ! @@ -284,6 +287,8 @@ subroutine setFiltersOneGroup(bounds, this_filter, include_inactive, icemask_grc integer :: fl ! lake filter index integer :: fnl,fnlu ! non-lake filter index integer :: fs ! soil filter index + integer :: fc, fpc ! crop and perennial crop filter index + integer :: fnc ! non-crop filter index integer :: f, fn ! general indices integer :: g !gridcell index !------------------------------------------------------------------------ @@ -391,24 +396,31 @@ subroutine setFiltersOneGroup(bounds, this_filter, include_inactive, icemask_grc ! Create prognostic crop and soil w/o prog. crop filters at pft-level ! according to where the crop model should be used - fl = 0 - fnl = 0 + fc = 0 + fpc = 0 + fnc = 0 do p = bounds%begp,bounds%endp if (veg_pp%active(p) .or. include_inactive) then - if (veg_pp%itype(p) >= npcropmin) then !skips 2 generic crop types - fl = fl + 1 - this_filter(nc)%pcropp(fl) = p - else + if (veg_pp%itype(p) < npcropmin) then l =veg_pp%landunit(p) if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then - fnl = fnl + 1 - this_filter(nc)%soilnopcropp(fnl) = p + fnc = fnc + 1 + this_filter(nc)%soilnopcropp(fnc) = p + end if + else + if (veg_pp%itype(p) < nppercropmin) then + fc = fc + 1 + this_filter(nc)%pcropp(fc) = p + else if (veg_pp%itype(p) >= nppercropmin) then + fpc = fpc + 1 + this_filter(nc)%ppercropp(fpc) = p end if end if end if end do - this_filter(nc)%num_pcropp = fl - this_filter(nc)%num_soilnopcropp = fnl ! This wasn't being set before... + this_filter(nc)%num_pcropp = fc + this_filter(nc)%num_ppercropp = fpc + this_filter(nc)%num_soilnopcropp = fnc ! This wasn't being set before... ! Create landunit-level urban and non-urban filters diff --git a/components/elm/src/main/pftvarcon.F90 b/components/elm/src/main/pftvarcon.F90 index e86c99d7408e..071bf28d25c1 100644 --- a/components/elm/src/main/pftvarcon.F90 +++ b/components/elm/src/main/pftvarcon.F90 @@ -10,6 +10,7 @@ module pftvarcon use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use elm_varpar , only : mxpft, numrad, ivis, inir, cft_lb, cft_ub + use elm_varpar , only: mxpft_nc use elm_varctl , only : iulog, use_vertsoilc use elm_varpar , only : nlevdecomp_full, nsoilorder use elm_varctl , only : nu_com @@ -37,6 +38,7 @@ module pftvarcon integer :: nc3_nonarctic_grass !value for C3 non-arctic grass integer :: nc4_grass !value for C4 grass integer :: npcropmin !value for first crop + integer :: nppercropmin !value for first perennial crop integer :: ncorn !value for corn, rain fed (rf) integer :: ncornirrig !value for corn, irrigated (ir) integer :: nscereal !value for spring temperate cereal (rf) @@ -46,9 +48,36 @@ module pftvarcon integer :: nsoybean !value for soybean (rf) integer :: nsoybeanirrig !value for soybean (ir) integer :: npcropmax !value for last prognostic crop in list + integer :: nppercropmax !value for last prognostic perennial crop in list integer :: nc3crop !value for generic crop (rf) integer :: nc3irrig !value for irrigated generic crop (ir) - + integer :: ncassava !value for cassava, rain fed (rf) + integer :: ncassavairrig !value for cassava, irrigated (ir) + integer :: ncotton !value for cotton, rain fed (rf) + integer :: ncottonirrig !value for cotton, irrigated (ir) + integer :: nfoddergrass !value for foddergrass, rain fed (rf) + integer :: nfoddergrassirrig !value for foddergrass, irrigated (ir) + integer :: noilpalm !value for oilpalm, rain fed (rf) + integer :: noilpalmirrig !value for oilpalm, irrigated (ir) + integer :: nograins !value for other grains, rain fed (rf) + integer :: nograinsirrig !value for other grains, irrigated (ir) + integer :: nrapeseed !value for rapeseed, rain fed (rf) + integer :: nrapeseedirrig !value for rapeseed, irrigated (ir) + integer :: nrice !value for rice, rain fed (rf) + integer :: nriceirrig !value for rice, irrigated (ir) + integer :: nrtubers !value for root tubers, rain fed (rf) + integer :: nrtubersirrig !value for root tubers, irrigated (ir) + integer :: nsugarcane !value for sugarcane, rain fed (rf) + integer :: nsugarcaneirrig !value for sugarcane, irrigated (ir) + integer :: nmiscanthus !value for miscanthus, rain fed (rf) + integer :: nmiscanthusirrig !value for miscanthus, irrigated (ir) + integer :: nswitchgrass !value for switchgrass, rain fed (rf) + integer :: nswitchgrassirrig !value for switchgrass, irrigated (ir) + integer :: npoplar !value for poplar, rain fed (rf) + integer :: npoplarirrig !value for poplar, irrigated (ir) + integer :: nwillow !value for willow, rain fed (rf) + integer :: nwillowirrig !value for willow, irrigated (ir) + ! Number of crop functional types actually used in the model. This includes each CFT for ! which is_pft_known_to_model is true. Note that this includes irrigated crops even if ! irrigation is turned off in this run: it just excludes crop types that aren't handled @@ -67,6 +96,7 @@ module pftvarcon real(r8), allocatable :: roota_par(:) !CLM rooting distribution parameter [1/m] real(r8), allocatable :: rootb_par(:) !CLM rooting distribution parameter [1/m] real(r8), allocatable :: crop(:) !crop pft: 0. = not crop, 1. = crop pft + real(r8), allocatable :: percrop(:) !perennial crop pft: 0. = not annual crop, 1. = perennial crop pft real(r8), allocatable :: irrigated(:) !irrigated pft: 0. = not, 1. = irrigated real(r8), allocatable :: smpso(:) !soil water potential at full stomatal opening (mm) real(r8), allocatable :: smpsc(:) !soil water potential at full stomatal closure (mm) @@ -130,6 +160,8 @@ module pftvarcon integer , allocatable :: mxSHplantdate(:)!maximum planting date for SouthHemisphere (YYYYMMDD) real(r8), allocatable :: planttemp(:) !planting temperature used in Phenology (K) real(r8), allocatable :: minplanttemp(:) !mininum planting temperature used in Phenology (K) + real(r8), allocatable :: senestemp(:) !senescence temperature for perennial crops used in Phenology (K) + real(r8), allocatable :: min_days_senes(:) !minimum leaf age to allow for leaf senescence real(r8), allocatable :: froot_leaf(:) !allocation parameter: new fine root C per new leaf C (gC/gC) real(r8), allocatable :: stem_leaf(:) !allocation parameter: new stem c per new leaf C (gC/gC) real(r8), allocatable :: croot_stem(:) !allocation parameter: new coarse root C per new stem C (gC/gC) @@ -175,7 +207,7 @@ module pftvarcon real(r8), allocatable :: root_dmx(:) !maximum root depth integer, parameter :: pftname_len = 40 ! max length of pftname - character(len=pftname_len) :: pftname(0:mxpft) !PFT description + character(len=:), allocatable :: pftname(:) !PFT description real(r8), parameter :: reinickerp = 1.6_r8 !parameter in allometric equation real(r8), parameter :: dwood = 2.5e5_r8 !cn wood density (gC/m3); lpj:2.0e5 @@ -346,6 +378,32 @@ subroutine pftconrd expected_pftnames(22) = 'irrigated_winter_temperate_cereal ' expected_pftnames(23) = 'soybean ' expected_pftnames(24) = 'irrigated_soybean ' + expected_pftnames(25) = 'cassava ' + expected_pftnames(26) = 'irrigated_cassava ' + expected_pftnames(27) = 'cotton ' + expected_pftnames(28) = 'irrigated_cotton ' + expected_pftnames(29) = 'foddergrass ' + expected_pftnames(30) = 'irrigated_foddergrass ' + expected_pftnames(31) = 'oilpalm ' + expected_pftnames(32) = 'irrigated_oilpalm ' + expected_pftnames(33) = 'other_grains ' + expected_pftnames(34) = 'irrigated_other_grains ' + expected_pftnames(35) = 'rapeseed ' + expected_pftnames(36) = 'irrigated_rapeseed ' + expected_pftnames(37) = 'rice ' + expected_pftnames(38) = 'irrigated_rice ' + expected_pftnames(39) = 'root_tubers ' + expected_pftnames(40) = 'irrigated_root_tubers ' + expected_pftnames(41) = 'sugarcane ' + expected_pftnames(42) = 'irrigated_sugarcane ' + expected_pftnames(43) = 'miscanthus ' + expected_pftnames(44) = 'irrigated_miscanthus ' + expected_pftnames(45) = 'switchgrass ' + expected_pftnames(46) = 'irrigated_switchgrass ' + expected_pftnames(47) = 'poplar ' + expected_pftnames(48) = 'irrigated_poplar ' + expected_pftnames(49) = 'willow ' + expected_pftnames(50) = 'irrigated_willow ' allocate( dleaf (0:mxpft) ) allocate( c3psn (0:mxpft) ) @@ -359,13 +417,14 @@ subroutine pftconrd allocate( roota_par (0:mxpft) ) allocate( rootb_par (0:mxpft) ) allocate( crop (0:mxpft) ) + allocate( percrop (0:mxpft) ) allocate( irrigated (0:mxpft) ) allocate( smpso (0:mxpft) ) allocate( smpsc (0:mxpft) ) allocate( fnitr (0:mxpft) ) allocate( slatop (0:mxpft) ) allocate( dsladlai (0:mxpft) ) - allocate( leafcn (0:mxpft) ) + allocate( leafcn (0:mxpft) ) allocate( flnr (0:mxpft) ) allocate( woody (0:mxpft) ) allocate( lflitcn (0:mxpft) ) @@ -413,6 +472,8 @@ subroutine pftconrd allocate( mxSHplantdate (0:mxpft) ) allocate( planttemp (0:mxpft) ) allocate( minplanttemp (0:mxpft) ) + allocate( senestemp (0:mxpft) ) + allocate( min_days_senes (0:mxpft) ) allocate( froot_leaf (0:mxpft) ) allocate( stem_leaf (0:mxpft) ) allocate( croot_stem (0:mxpft) ) @@ -455,6 +516,12 @@ subroutine pftconrd allocate( fyield (0:mxpft) ) allocate( root_dmx (0:mxpft) ) + if (use_crop) then + allocate(character(pftname_len) :: pftname(0:mxpft)) + else + allocate(character(pftname_len) :: pftname(0:mxpft_nc)) + end if + allocate( VMAX_PLANT_NH4(0:mxpft) ) allocate( VMAX_PLANT_NO3(0:mxpft) ) allocate( VMAX_PLANT_P(0:mxpft) ) @@ -653,7 +720,13 @@ subroutine pftconrd if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) call ncd_io('fyield',fyield, 'read', ncid, readvar=readv, posNOTonfile=.true.) if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) - endif + call ncd_io('percrop', percrop, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) + call ncd_io('senescence_temp', senestemp, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) + call ncd_io('min_days_senescence', min_days_senes, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) + end if if(use_dynroot)then call ncd_io('root_dmx',root_dmx, 'read', ncid, readvar=readv, posNOTonfile=.true.) if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(__FILE__, __LINE__)) @@ -936,6 +1009,7 @@ subroutine pftconrd ! or other modules that co-exist while FATES is on, we may want to preserve these pft definitions ! on non-fates columns. For now, they are incompatible, and this check is warranted (rgk 04-2017) if(.not. use_fates)then + if(.not. use_crop .and. i > mxpft_nc) EXIT ! exit the do loop if ( trim(adjustl(pftname(i))) /= trim(expected_pftnames(i)) )then write(iulog,*)'pftconrd: pftname is NOT what is expected, name = ', & trim(pftname(i)), ', expected name = ', trim(expected_pftnames(i)) @@ -968,30 +1042,103 @@ subroutine pftconrd if ( trim(pftname(i)) == 'irrigated_winter_temperate_cereal' ) nwcerealirrig = i if ( trim(pftname(i)) == 'soybean' ) nsoybean = i if ( trim(pftname(i)) == 'irrigated_soybean' ) nsoybeanirrig = i + if ( trim(pftname(i)) == 'cassava' ) ncassava = i + if ( trim(pftname(i)) == 'irrigated_cassava' ) ncassavairrig = i + if ( trim(pftname(i)) == 'cotton' ) ncotton = i + if ( trim(pftname(i)) == 'irrigated_cotton' ) ncottonirrig = i + if ( trim(pftname(i)) == 'foddergrass' ) nfoddergrass = i + if ( trim(pftname(i)) == 'irrigated_foddergrass' ) nfoddergrassirrig = i + if ( trim(pftname(i)) == 'oilpalm' ) noilpalm = i + if ( trim(pftname(i)) == 'irrigated_oilpalm' ) noilpalmirrig = i + if ( trim(pftname(i)) == 'other_grains' ) nograins = i + if ( trim(pftname(i)) == 'irrigated_other_grains' ) nograinsirrig = i + if ( trim(pftname(i)) == 'rapeseed' ) nrapeseed = i + if ( trim(pftname(i)) == 'irrigated_rapeseed' ) nrapeseedirrig = i + if ( trim(pftname(i)) == 'rice' ) nrice = i + if ( trim(pftname(i)) == 'irrigated_rice' ) nriceirrig = i + if ( trim(pftname(i)) == 'root_tubers' ) nrtubers = i + if ( trim(pftname(i)) == 'irrigated_root_tubers' ) nrtubersirrig = i + if ( trim(pftname(i)) == 'sugarcane' ) nsugarcane = i + if ( trim(pftname(i)) == 'irrigated_sugarcane' ) nsugarcaneirrig = i + if ( trim(pftname(i)) == 'miscanthus' ) nmiscanthus = i + if ( trim(pftname(i)) == 'irrigated_miscanthus' ) nmiscanthusirrig = i + if ( trim(pftname(i)) == 'switchgrass' ) nswitchgrass = i + if ( trim(pftname(i)) == 'irrigated_switchgrass' ) nswitchgrassirrig = i + if ( trim(pftname(i)) == 'poplar' ) npoplar = i + if ( trim(pftname(i)) == 'irrigated_poplar' ) npoplarirrig = i + if ( trim(pftname(i)) == 'willow' ) nwillow = i + if ( trim(pftname(i)) == 'irrigated_willow' ) nwillowirrig = i end do ntree = nbrdlf_dcd_brl_tree ! value for last type of tree npcropmin = ncorn ! first prognostic crop - npcropmax = nsoybeanirrig ! last prognostic crop in list + if( .not. use_crop) then + npcropmax = nsoybeanirrig ! last prognostic crop in list + else + npcropmax = nsugarcaneirrig ! last prognostic crop in list + nppercropmin = nmiscanthus ! first prognostic perennial crop + nppercropmax = nwillowirrig ! last prognostic perennial crop in list + end if call set_is_pft_known_to_model() call set_num_cfts_known_to_model() if( .not. use_fates ) then - if ( npcropmax /= mxpft )then - call endrun(msg=' ERROR: npcropmax is NOT the last value'//errMsg(__FILE__, __LINE__)) + if( .not. use_crop) then + if ( npcropmax /= mxpft_nc )then + call endrun(msg=' ERROR: npcropmax is NOT the last value'//errMsg(__FILE__, __LINE__)) + end if + else + if ( nppercropmax /= mxpft )then + call endrun(msg=' ERROR: nppercropmax is NOT the last value'//errMsg(__FILE__, __LINE__)) + end if end if do i = 0, mxpft - if ( irrigated(i) == 1.0_r8 .and. (i == nc3irrig .or. & - i == ncornirrig .or. & - i == nscerealirrig .or. & - i == nwcerealirrig .or. & - i == nsoybeanirrig) )then - ! correct - else if ( irrigated(i) == 0.0_r8 )then - ! correct + if(.not. use_crop .and. i > mxpft_nc) EXIT ! exit the do loop + if( .not. use_crop) then + if ( irrigated(i) == 1.0_r8 .and. (i == nc3irrig .or. & + i == ncornirrig .or. & + i == nscerealirrig .or. & + i == nwcerealirrig .or. & + i == nsoybeanirrig ) ) then + ! correct + else if ( irrigated(i) == 0.0_r8 )then + ! correct + else + call endrun(msg=' ERROR: irrigated has wrong values'//errMsg(__FILE__, __LINE__)) + end if else - call endrun(msg=' ERROR: irrigated has wrong values'//errMsg(__FILE__, __LINE__)) + if ( irrigated(i) == 1.0_r8 .and. (i == nc3irrig .or. & + i == ncornirrig .or. & + i == nscerealirrig .or. & + i == nwcerealirrig .or. & + i == nsoybeanirrig .or. & + i == ncassavairrig .or. & + i == ncottonirrig .or. & + i == nfoddergrassirrig .or. & + i == noilpalmirrig .or. & + i == nograinsirrig .or. & + i == nrapeseedirrig .or. & + i == nriceirrig .or. & + i == nrtubersirrig .or. & + i == nsugarcaneirrig .or. & + i == nmiscanthusirrig .or. & + i == nswitchgrassirrig .or. & + i == npoplarirrig .or. & + i == nwillowirrig ) )then + ! correct + else if ( irrigated(i) == 0.0_r8 )then + ! correct + else + call endrun(msg=' ERROR: irrigated has wrong values'//errMsg(__FILE__, __LINE__)) + end if + if ( percrop(i) == 1.0_r8 .and. (i >= nmiscanthus .and. i <= nppercropmax))then + ! correct + else if ( percrop(i) == 0.0_r8 )then + ! correct + else + call endrun(msg=' ERROR: perennial crop has wrong values'//errMsg(__FILE__, __LINE__)) + end if end if if ( crop(i) == 1.0_r8 .and. (i >= nc3crop .and. i <= npcropmax) )then ! correct diff --git a/components/elm/src/main/restFileMod.F90 b/components/elm/src/main/restFileMod.F90 index bb13d52fcbf9..279aefc34fa9 100644 --- a/components/elm/src/main/restFileMod.F90 +++ b/components/elm/src/main/restFileMod.F90 @@ -1116,8 +1116,9 @@ subroutine restFile_add_ipft_metadata(ncid) ! Add global metadata defining pft types ! ! !USES: - use elm_varpar, only : natpft_lb, mxpft, cft_lb, cft_ub + use elm_varpar, only : natpft_lb, mxpft, cft_lb, cft_ub, mxpft_nc use pftvarcon , only : pftname_len, pftname + use elm_varctl, only : use_crop ! ! !ARGUMENTS: type(file_desc_t), intent(inout) :: ncid ! local file id @@ -1131,6 +1132,7 @@ subroutine restFile_add_ipft_metadata(ncid) !----------------------------------------------------------------------- do ptype = natpft_lb, mxpft + if(.not. use_crop .and. ptype > mxpft_nc) EXIT ! exit the do loop attname = att_prefix // pftname(ptype) call ncd_putatt(ncid, ncd_global, attname, ptype) end do diff --git a/components/homme/src/share/control_mod.F90 b/components/homme/src/share/control_mod.F90 index fc33ac15c4f6..bd535a5c0aad 100644 --- a/components/homme/src/share/control_mod.F90 +++ b/components/homme/src/share/control_mod.F90 @@ -258,6 +258,24 @@ module control_mod real (kind=real_kind), public :: dcmip16_mu_q = -1 ! additional uniform viscosity (scalar tracers); -1 implies it defaults to dcmip16_mu_s value real (kind=real_kind), public :: interp_lon0 = 0.0d0 +!PLANAR + real (kind=real_kind), private, parameter :: tol_zero=1e-10 !tolerance to determine if lx,ly,sx,sy are set + + real (kind=real_kind), public :: bubble_T0 = 270.0 !bubble ref state + real (kind=real_kind), public :: bubble_dT = 0.5 !bubble dTheta + real (kind=real_kind), public :: bubble_xycenter = 0.0 !bubble xy position + real (kind=real_kind), public :: bubble_zcenter = 3000.0 !bubble z position + real (kind=real_kind), public :: bubble_ztop = 10000.0 !bubble z top + real (kind=real_kind), public :: bubble_xyradius = 2000.0!bubble radius along x or y axis + real (kind=real_kind), public :: bubble_zradius = 1500.0 !bubble radius along z axis + logical, public :: bubble_cosine = .TRUE. !bubble uniform or cosine + logical, public :: bubble_moist = .FALSE. ! + real (kind=real_kind), public :: bubble_moist_dq = 0.0 !bubble dQ parameter + integer, public :: bubble_prec_type = 0 !0 kessler, 1 rj + logical, protected :: case_planar_bubble = .FALSE. + + public :: set_planar_defaults + contains function timestep_make_parameters_consistent(par, rsplit, qsplit, & @@ -622,4 +640,79 @@ subroutine test_timestep_make_parameters_consistent(par, nerr) write(iulog,'(a,i2)') 'test_timestep_make_parameters_consistent nerr', nerr end subroutine test_timestep_make_parameters_consistent + +subroutine set_planar_defaults() + +use physical_constants, only: Lx, Ly, Sx, Sy + +!since defaults here depend on test, they cannot be set before ctl_nl is read, unlike some other parameters, bubble_*, etc. +!if true, most likely lx,ly,sx,sy weren't set in ctl_nl + if ( abs(lx).le.tol_zero .and. abs(ly).le.tol_zero & + .and. abs(sx).le.tol_zero .and. abs(sy).le.tol_zero )then + if (test_case == "planar_dbl_vrtx") then + Lx = 5000.0D0 * 1000.0D0 + Ly = 5000.0D0 * 1000.0D0 + Sx = 0.0D0 + Sy = 0.0D0 + else if (test_case == "planar_hydro_gravity_wave") then + Lx = 6000.0D0 * 1000.0D0 + Ly = 6000.0D0 * 1000.0D0 + Sx = -3000.0D0 * 1000.0D0 + Sy = -3000.0D0 * 1000.0D0 + else if (test_case == "planar_nonhydro_gravity_wave") then + Lx = 300.0D0 * 1000.0D0 + Ly = 300.0D0 * 1000.0D0 + Sx = -150.0D0 * 1000.0D0 + Sy = -150.0D0 * 1000.0D0 + else if (test_case == "planar_hydro_mtn_wave") then + Lx = 240.0D0 * 1000.0D0 + Ly = 240.0D0 * 1000.0D0 + Sx = 0.0D0 * 1000.0D0 + Sy = 0.0D0 * 1000.0D0 + else if (test_case == "planar_nonhydro_mtn_wave") then + Lx = 144.0D0 * 1000.0D0 + Ly = 144.0D0 * 1000.0D0 + Sx = 0.0D0 * 1000.0D0 + Sy = 0.0D0 * 1000.0D0 + else if (test_case == "planar_schar_mtn_wave") then + Lx = 100.0D0 * 1000.0D0 + Ly = 100.0D0 * 1000.0D0 + Sx = 0.0D0 * 1000.0D0 + Sy = 0.0D0 * 1000.0D0 + else if (test_case == "planar_density_current" .OR. test_case == "planar_moist_density_current") then + Lx = 51.2D0 * 1000.0D0 + Ly = 51.2D0 * 1000.0D0 + Sx = -25.6D0 * 1000.0D0 + Sy = -25.6D0 * 1000.0D0 + else if (test_case == "planar_rising_bubble" ) then + Lx = 2.0D0 * 10000.0D0 + Ly = 2.0D0 * 10000.0D0 + Sx = -10000.0D0 + Sy = -10000.0D0 +! THESE ARE WRONG AND NEED TO BE FIXED WHEN THESE CASES ARE ACTUALLY IMPLEMENTED.... +!else if (test_case == "planar_baroclinic_instab" .OR. test_case == "planar_moist_baroclinic_instab") then +! Lx = 5000.0D0 * 1000.0D0 +! Ly = 5000.0D0 * 1000.0D0 +! Sx = 0.0D0 +! Sy = 0.0D0 +! else if (test_case == "planar_tropical_cyclone") then +! Lx = 5000.0D0 * 1000.0D0 +! Ly = 5000.0D0 * 1000.0D0 +! Sx = 0.0D0 +! Sy = 0.0D0 +! else if (test_case == "planar_supercell") then +! Lx = 5000.0D0 * 1000.0D0 +! Ly = 5000.0D0 * 1000.0D0 +! Sx = 0.0D0 +! Sy = 0.0D0 + + endif + endif !if lx,ly,sx,sy are not set in nl + + if (test_case == "planar_rising_bubble" ) then + case_planar_bubble = .TRUE. + end if + +end subroutine set_planar_defaults + end module control_mod diff --git a/components/homme/src/share/namelist_mod.F90 b/components/homme/src/share/namelist_mod.F90 index acaf89a89620..c4dd7787886e 100644 --- a/components/homme/src/share/namelist_mod.F90 +++ b/components/homme/src/share/namelist_mod.F90 @@ -15,6 +15,7 @@ module namelist_mod use arkode_mod, only: rel_tol, abs_tol, calc_nonlinear_stats, use_column_solver #endif use physical_constants, only : rearth, rrearth, DD_PI + use physical_constants, only : scale_factor, scale_factor_inv, domain_size, laplacian_rigid_factor use physical_constants, only : Sx, Sy, Lx, Ly, dx, dy, dx_ref, dy_ref @@ -103,6 +104,26 @@ module namelist_mod se_fv_phys_remap_alg, & timestep_make_subcycle_parameters_consistent + +!PLANAR setup +#ifndef CAM + use control_mod, only: & + set_planar_defaults,& + bubble_T0, & + bubble_dT, & + bubble_xycenter, & + bubble_zcenter, & + bubble_ztop, & + bubble_xyradius, & + bubble_zradius, & + bubble_cosine, & + bubble_moist, & + bubble_moist_dq, & + bubble_prec_type, & + case_planar_bubble +#endif + + #ifndef CAM use control_mod, only: & pertlim, & @@ -318,6 +339,21 @@ subroutine readnl(par) dcmip2_x_xi, & !dcmip2-x mountain wavelength (m) dcmip4_moist, & !dcmip4 moist, 0 or 1 dcmip4_X !dcmip4 scaling factor, nondim +!PLANAR + namelist /ctl_nl/ & + lx, ly, & + sx, sy, & + bubble_T0, & + bubble_dT, & + bubble_xycenter, & + bubble_zcenter, & + bubble_ztop, & + bubble_xyradius, & + bubble_zradius, & + bubble_cosine, & + bubble_moist, & + bubble_moist_dq, & + bubble_prec_type namelist /vert_nl/ & vform, & vfile_mid, & @@ -422,7 +458,6 @@ subroutine readnl(par) ! write(iulog,*) "masterproc=",par%masterproc if (par%masterproc) then - write(iulog,*)"reading ctl namelist..." #if defined(CAM) unitn=getunit() @@ -491,11 +526,26 @@ subroutine readnl(par) #ifndef CAM - +!checks before the next NL +!check on ne if (topology == "plane" .and. ne /= 0) then call abortmp('cannot set ne for planar topology, use ne_x and ne_y instead') end if +!setting default PLANAR values if they are not set in ctl_nl +call set_planar_defaults() + +!checks on planar tests +if (test_case(1:7)=="planar_") then + +!check on Lx, Ly +if(lx .le. 0.d0 .or. ly .le. 0.d0) then +call abortmp('for planar tests, planar_lx and planar_ly should be >0') +endif + +endif + + #ifdef _PRIM write(iulog,*) "reading physics namelist..." if (test_case == "held_suarez0" .or. & @@ -766,8 +816,32 @@ subroutine readnl(par) #ifndef HOMME_WITHOUT_PIOLIBRARY call MPI_bcast(mesh_file,MAX_FILE_LEN,MPIChar_t ,par%root,par%comm,ierr) #endif - call MPI_bcast(theta_hydrostatic_mode ,1,MPIlogical_t,par%root,par%comm,ierr) + +!PLANAR + call MPI_bcast(lx ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(ly ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(sx ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(sy ,1,MPIreal_t ,par%root,par%comm,ierr) call MPI_bcast(planar_slice ,1,MPIlogical_t,par%root,par%comm,ierr) + +!PLANAR +#ifndef CAM + call MPI_bcast(bubble_T0 ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_dT ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_xycenter,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_zcenter ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_ztop ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_xyradius ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_zradius ,1,MPIreal_t ,par%root,par%comm,ierr) + call MPI_bcast(bubble_cosine ,1,MPIlogical_t,par%root,par%comm,ierr) + call MPI_bcast(bubble_moist ,1,MPIlogical_t,par%root,par%comm,ierr) + call MPI_bcast(bubble_moist_dq ,1,MPIreal_t,par%root,par%comm,ierr) + call MPI_bcast(bubble_prec_type, 1, MPIinteger_t, par%root,par%comm,ierr) + + call MPI_bcast(case_planar_bubble,1,MPIlogical_t,par%root,par%comm,ierr) +#endif + + call MPI_bcast(theta_hydrostatic_mode ,1,MPIlogical_t,par%root,par%comm,ierr) call MPI_bcast(transport_alg ,1,MPIinteger_t,par%root,par%comm,ierr) call MPI_bcast(semi_lagrange_cdr_alg ,1,MPIinteger_t,par%root,par%comm,ierr) call MPI_bcast(semi_lagrange_cdr_check ,1,MPIlogical_t,par%root,par%comm,ierr) @@ -834,6 +908,7 @@ subroutine readnl(par) use_moisture = ( moisture /= "dry") if (qsize<1) use_moisture = .false. + if(par%masterproc) print *, "use moisture:", use_moisture ! use maximum available: @@ -923,84 +998,27 @@ subroutine readnl(par) scale_factor = 1.0D0 scale_factor_inv = 1.0D0 laplacian_rigid_factor = 0.0D0 !this eliminates the correction to ensure the Laplacian doesn't damp rigid motion - if (test_case == "planar_dbl_vrtx") then - Lx = 5000.0D0 * 1000.0D0 - Ly = 5000.0D0 * 1000.0D0 - Sx = 0.0D0 - Sy = 0.0D0 - else if (test_case == "planar_hydro_gravity_wave") then - Lx = 6000.0D0 * 1000.0D0 - Ly = 6000.0D0 * 1000.0D0 - Sx = -3000.0D0 * 1000.0D0 - Sy = -3000.0D0 * 1000.0D0 - else if (test_case == "planar_nonhydro_gravity_wave") then - Lx = 300.0D0 * 1000.0D0 - Ly = 300.0D0 * 1000.0D0 - Sx = -150.0D0 * 1000.0D0 - Sy = -150.0D0 * 1000.0D0 - else if (test_case == "planar_hydro_mtn_wave") then - Lx = 240.0D0 * 1000.0D0 - Ly = 240.0D0 * 1000.0D0 - Sx = 0.0D0 * 1000.0D0 - Sy = 0.0D0 * 1000.0D0 - else if (test_case == "planar_nonhydro_mtn_wave") then - Lx = 144.0D0 * 1000.0D0 - Ly = 144.0D0 * 1000.0D0 - Sx = 0.0D0 * 1000.0D0 - Sy = 0.0D0 * 1000.0D0 - else if (test_case == "planar_schar_mtn_wave") then - Lx = 100.0D0 * 1000.0D0 - Ly = 100.0D0 * 1000.0D0 - Sx = 0.0D0 * 1000.0D0 - Sy = 0.0D0 * 1000.0D0 - else if (test_case == "planar_density_current" .OR. test_case == "planar_moist_density_current") then - Lx = 51.2D0 * 1000.0D0 - Ly = 51.2D0 * 1000.0D0 - Sx = -25.6D0 * 1000.0D0 - Sy = -25.6D0 * 1000.0D0 - else if (test_case == "planar_rising_bubble" .OR. test_case == "planar_moist_rising_bubble") then - Lx = 1.0D0 * 1000.0D0 - Ly = 1.0D0 * 1000.0D0 - Sx = 0.0D0 - Sy = 0.0D0 -! THESE ARE WRONG AND NEED TO BE FIXED WHEN THESE CASES ARE ACTUALLY IMPLEMENTED.... -!else if (test_case == "planar_baroclinic_instab" .OR. test_case == "planar_moist_baroclinic_instab") then -! Lx = 5000.0D0 * 1000.0D0 -! Ly = 5000.0D0 * 1000.0D0 -! Sx = 0.0D0 -! Sy = 0.0D0 -! else if (test_case == "planar_tropical_cyclone") then -! Lx = 5000.0D0 * 1000.0D0 -! Ly = 5000.0D0 * 1000.0D0 -! Sx = 0.0D0 -! Sy = 0.0D0 -! else if (test_case == "planar_supercell") then -! Lx = 5000.0D0 * 1000.0D0 -! Ly = 5000.0D0 * 1000.0D0 -! Sx = 0.0D0 -! Sy = 0.0D0 - end if ! makes the y-direction cells identical in size to the x-dir cells ! this is important for hyperviscosity, etc. ! Also adjusts Sy so y-dir domain is centered at 0 - if (planar_slice .eqv. .true.) then - Ly = (Lx / ne_x) * ne_y - Sy = -Ly/2.0D0 - end if + if (planar_slice .eqv. .true.) then + Ly = (Lx / ne_x) * ne_y + Sy = -Ly/2.0D0 + end if - domain_size = Lx * Ly - dx = Lx/ne_x - dy = Ly/ne_y - dx_ref = 1.0D0/ne_x - dy_ref = 1.0D0/ne_y + domain_size = Lx * Ly + dx = Lx/ne_x + dy = Ly/ne_y + dx_ref = 1.0D0/ne_x + dy_ref = 1.0D0/ne_y - else if (geometry == "sphere") then + else if (geometry == "sphere") then scale_factor = rearth scale_factor_inv = rrearth domain_size = 4.0D0*DD_PI laplacian_rigid_factor = rrearth - end if + end if ! if plane !logic around different hyperviscosity options if (hypervis_power /= 0) then diff --git a/components/homme/src/share/planar_mod.F90 b/components/homme/src/share/planar_mod.F90 index b741c00cbab4..64bf02a7c0fa 100644 --- a/components/homme/src/share/planar_mod.F90 +++ b/components/homme/src/share/planar_mod.F90 @@ -20,6 +20,8 @@ module planar_mod use control_mod, only : hypervis_scaling, cubed_sphere_map use spacecurve_mod, only : GilbertCurve + use physical_constants, only: Lx, Ly, Sx, Sy, dx, dy, dx_ref, dy_ref + implicit none private @@ -262,7 +264,6 @@ end subroutine plane_init_atomic subroutine coordinates_atomic(elem,gll_points) use element_mod, only : element_t - use physical_constants, only: dx, dy, dx_ref, dy_ref, Sx, Sy, Lx, Ly type (element_t) :: elem real (kind=longdouble_kind) :: gll_points(np) @@ -310,7 +311,6 @@ end subroutine coordinates_atomic subroutine Dmap(D, a,b, corners3D, ref_map, cartp, facenum) - use physical_constants, only: dx, dy real (kind=real_kind), intent(out) :: D(2,2) real (kind=real_kind), intent(in) :: a,b type (cartesian3D_t) :: corners3D(4) !x,y,z coords of element corners @@ -363,7 +363,6 @@ end subroutine Dmap subroutine metric_atomic(elem,gll_points,alpha) use element_mod, only : element_t - use physical_constants, only: Lx, Ly, dx, dy type (element_t) :: elem real(kind=real_kind) :: alpha @@ -597,7 +596,6 @@ end subroutine coriolis_init_atomic subroutine plane_set_corner_coordinates(elem) use element_mod, only : element_t - use physical_constants, only : dx_ref, dy_ref type (element_t) :: elem ! Local variables diff --git a/components/homme/src/test_mod.F90 b/components/homme/src/test_mod.F90 index 82042c45dad3..eddc4bbe01c5 100644 --- a/components/homme/src/test_mod.F90 +++ b/components/homme/src/test_mod.F90 @@ -4,7 +4,7 @@ module test_mod -use control_mod, only: test_case, sub_case, dt_remap_factor, runtype +use control_mod, only: test_case, sub_case, dt_remap_factor, runtype, bubble_moist use dimensions_mod, only: np, nlev, nlevp, qsize use derivative_mod, only: derivative_t, gradient_sphere use element_mod, only: element_t @@ -88,6 +88,7 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n case('planar_nonhydro_mtn_wave'); case('planar_schar_mtn_wave'); case('planar_rising_bubble'); + if (bubble_moist) call dcmip2016_init(); case('planar_density_current'); case('planar_baroclinic_instab'); case('planar_moist_rising_bubble'); @@ -235,6 +236,9 @@ subroutine compute_test_forcing(elem,hybrid,hvcoord,nt,ntQ,dt,nets,nete,tl) case('dcmip2016_test2'); call dcmip2016_test2_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl,2) case('dcmip2016_test3'); call dcmip2016_test3_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) + case('planar_rising_bubble'); + if (bubble_moist) call dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) + case('held_suarez0'); do ie=nets,nete call hs_forcing(elem(ie),hvcoord,nt,ntQ,dt) @@ -242,21 +246,7 @@ subroutine compute_test_forcing(elem,hybrid,hvcoord,nt,ntQ,dt,nets,nete,tl) endselect -!for ftype3 we scale tendencies by dp - if(ftype == 3) then - !initialize dp3d from ps - do ie=nets,nete - do k=1,nlev - dp(:,:)= ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & - ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*elem(ie)%state%ps_v(:,:,nt) - elem(ie)%derived%FT(:,:,k) = elem(ie)%derived%FT(:,:,k) * dp(:,:) - elem(ie)%derived%FM(:,:,1,k) = elem(ie)%derived%FM(:,:,1,k) * dp(:,:) - elem(ie)%derived%FM(:,:,2,k) = elem(ie)%derived%FM(:,:,2,k) * dp(:,:) - enddo - enddo - endif - -end subroutine +end subroutine compute_test_forcing !_____________________________________________________________________ diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90 index 3fed7d4c90cc..d616ff63af42 100644 --- a/components/homme/src/test_src/dcmip16_wrapper.F90 +++ b/components/homme/src/test_src/dcmip16_wrapper.F90 @@ -11,7 +11,7 @@ module dcmip16_wrapper use dcmip12_wrapper, only: pressure_thickness, set_tracers, get_evenly_spaced_z, set_hybrid_coefficients use control_mod, only: test_case, dcmip16_pbl_type, dcmip16_prec_type, use_moisture, theta_hydrostatic_mode,& - sub_case + sub_case, case_planar_bubble, bubble_prec_type use baroclinic_wave, only: baroclinic_wave_test use supercell, only: supercell_init, supercell_test, supercell_z use tropical_cyclone, only: tropical_cyclone_test @@ -40,7 +40,7 @@ module dcmip16_wrapper real(rl), parameter :: rh2o = 461.5d0, & ! Gas constant for water vapor (J/kg/K) Mvap = (Rwater_vapor/Rgas) - 1.d0 ! Constant for virtual temp. calc. (~0.608) -real(rl) :: sample_period = 60.0_rl +real(rl) :: sample_period = 2.0_rl real(rl) :: rad2dg = 180.0_rl/pi type :: PhysgridData_t @@ -355,7 +355,7 @@ subroutine dcmip2016_test3(elem,hybrid,hvcoord,nets,nete) enddo - sample_period = 60.0 ! sec + sample_period = 1 ! 60 orig sec end subroutine !_______________________________________________________________________ @@ -399,7 +399,9 @@ subroutine dcmip2016_append_measurements(max_w,max_precl,min_ps,tl,hybrid) pmin_ps = parallelMin(min_ps, hybrid) if (hybrid%masterthread) then - print *,"time=",time_at(tl%nstep)," pmax_w (m/s)=",pmax_w," pmax_precl (mm/day)=",pmax_precl*(1000.0)*(24.0*3600)," pmin_ps (Pa)=",pmin_ps + print *,"time=",time_at(tl%nstep)," pmax_w (m/s)=",pmax_w + print *,"time=",time_at(tl%nstep)," pmax_precl (mm/day)=",pmax_precl*(1000.0)*(24.0*3600) + print *,"time=",time_at(tl%nstep)," pmin_ps (Pa)=",pmin_ps open(unit=10,file=w_filename,form="formatted",position="append") write(10,'(99E24.15)') pmax_w @@ -444,8 +446,17 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) integer :: pbl_type, prec_type, qi integer, parameter :: test = 1 + logical :: toy_chemistry_on + + if (case_planar_bubble) then + toy_chemistry_on = .false. + prec_type = bubble_prec_type + if (qsize .ne. 3) call abortmp('ERROR: moist bubble test requires qsize=3') + else + toy_chemistry_on = .true. + prec_type = dcmip16_prec_type + endif - prec_type = dcmip16_prec_type pbl_type = dcmip16_pbl_type max_w = -huge(rl) @@ -468,8 +479,10 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) qi=1; qv = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp qi=2; qc = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp qi=3; qr = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp - qi=4; cl = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp - qi=5; cl2 = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp + if (toy_chemistry_on) then + qi=4; cl = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp + qi=5; cl2 = elem(ie)%state%Qdp(:,:,:,qi,ntQ)/dp + endif ! ensure positivity where(qv<0); qv=0; endwhere @@ -503,6 +516,7 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) th_c = theta_kess(i,j,nlev:1:-1) ! get forced versions of u,v,p,qv,qc,qr. rho is constant + lat=0.0 call DCMIP2016_PHYSICS(test, u_c, v_c, p_c, th_c, qv_c, qc_c, qr_c, rho_c, dt, z_c, zi_c, lat, nlev, & precl(i,j,ie), pbl_type, prec_type) @@ -514,12 +528,15 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) qr(i,j,:) = qr_c(nlev:1:-1) theta_kess(i,j,:) = th_c(nlev:1:-1) - lon = elem(ie)%spherep(i,j)%lon - lat = elem(ie)%spherep(i,j)%lat + if (toy_chemistry_on) then + lon = elem(ie)%spherep(i,j)%lon + lat = elem(ie)%spherep(i,j)%lat + + do k=1,nlev + call tendency_terminator( lat*rad2dg, lon*rad2dg, cl(i,j,k), cl2(i,j,k), dt, ddt_cl(i,j,k), ddt_cl2(i,j,k)) + enddo + endif - do k=1,nlev - call tendency_terminator( lat*rad2dg, lon*rad2dg, cl(i,j,k), cl2(i,j,k), dt, ddt_cl(i,j,k), ddt_cl2(i,j,k)) - enddo enddo; enddo; @@ -544,8 +561,12 @@ subroutine dcmip2016_test1_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl) elem(ie)%derived%FQ(:,:,:,2) = (rho_dry/rho)*dp*(qc-qc0)/dt elem(ie)%derived%FQ(:,:,:,3) = (rho_dry/rho)*dp*(qr-qr0)/dt + + if (toy_chemistry_on) then qi=4; elem(ie)%derived%FQ(:,:,:,qi) = dp*ddt_cl qi=5; elem(ie)%derived%FQ(:,:,:,qi) = dp*ddt_cl2 + endif + ! perform measurements of max w, and max prect max_w = max( max_w , maxval(w ) ) @@ -889,6 +910,7 @@ subroutine dcmip2016_test2_forcing(elem,hybrid,hvcoord,nets,nete,nt,ntQ,dt,tl, t zi_c = zi (i,j,nlevp:1:-1) th_c = theta_kess(i,j,nlev:1:-1) + lat=0.0 ! get forced versions of u,v,p,qv,qc,qr. rho is constant call DCMIP2016_PHYSICS(test, u_c, v_c, p_c, th_c, qv_c, qc_c, qr_c, rho_c, dt, z_c, zi_c, lat, nlev, & precl(i,j,ie), pbl_type, prec_type) diff --git a/components/homme/src/test_src/dry_planar.F90 b/components/homme/src/test_src/dry_planar.F90 index ba067dc74793..ffc9a5784ca1 100644 --- a/components/homme/src/test_src/dry_planar.F90 +++ b/components/homme/src/test_src/dry_planar.F90 @@ -10,7 +10,7 @@ module dry_planar_tests use parallel_mod, only: abortmp use element_ops, only: set_state, set_state_i, tests_finalize use kinds, only: rl=>real_kind, iulog - use physical_constants, only : dd_pi + use physical_constants, only : dd_pi, rgas, rwater_vapor use dcmip12_wrapper, only : set_tracers, get_evenly_spaced_z, get_evenly_spaced_p, set_hybrid_coefficients, pressure_thickness use dimensions_mod, only: np, nlev, nlevp, qsize, qsize_d, nelemd use element_state, only: nt=>timelevels @@ -18,12 +18,17 @@ module dry_planar_tests implicit none +!OG take from physical const mod instead + real(rl), parameter :: Rd = 287.0d0, & ! Ideal gas const dry air (J kg^-1 K^1) g = 9.80616d0, & ! Gravity (m s^2) cp = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1) p0 = 100000.d0, &! reference pressure (Pa) kappa = Rd/cp + real(rl), parameter :: bubble_const1=3.8, bubble_const2=17.27, bubble_const3=273.0, bubble_const4=36.0 + + real(rl):: zi(nlevp), zm(nlev) ! z coordinates real(rl):: ddn_hyai(nlevp), ddn_hybi(nlevp) ! vertical derivativess of hybrid coefficients real(rl):: ztop @@ -170,11 +175,228 @@ subroutine planar_rising_bubble_init(elem,hybrid,hvcoord,nets,nete) type(hvcoord_t), intent(inout) :: hvcoord ! hybrid vertical coordinates integer, intent(in) :: nets,nete ! start, end element index - call abortmp('planar rising bubble not yet implemented') + !last input is coriolis parameter + call bubble_init(elem,hybrid,hvcoord,nets,nete,0.d0) end subroutine planar_rising_bubble_init +! was not tested with preqx + +! A note on relation T = \Pi * theta where \Pi = (p/p0)^dry_kappa that is used below: +! depending on assumptions for Gibbs potential, definition of \Pi with dry kapa +! may not hold. It is correct for theta-l and is used (though may be +! inconsistent) in preqx in other places. +subroutine bubble_init(elem,hybrid,hvcoord,nets,nete,f) + + use control_mod, only: bubble_T0, bubble_dT, bubble_xycenter, bubble_zcenter, bubble_ztop, & + bubble_xyradius,bubble_zradius, bubble_cosine, & + bubble_moist, bubble_moist_dq, bubble_prec_type + use physical_constants, only: Lx, Ly, Sx, Sy + use element_ops, only: set_elem_state + + type(element_t), intent(inout), target :: elem(:) ! element array + type(hybrid_t), intent(in) :: hybrid ! hybrid parallel structure + type(hvcoord_t), intent(inout) :: hvcoord ! hybrid vertical coordinates + integer, intent(in) :: nets,nete ! start, end element index + real(rl), intent(in) :: f ! (const) Coriolis force + + integer :: i,j,k,ie,ii + real(rl):: x,y,offset + real(rl):: pi(nlevp), pm(nlev), dpm(nlev), th0(nlevp), th0m(nlev), ai(nlevp), bi(nlevp), rr, & + qi_s(nlevp), qm_s(nlev), Ti(nlevp), Tm(nlev), qi(nlevp) + + real(rl):: zero_mid_init(np,np,nlev),zero_int_init(np,np,nlevp), & + dp_init(np,np,nlev),ps_init(np,np), & + phis_init(np,np,nlevp),t_init(np,np,nlev),p_init(np,np,nlev), & + zi_init(np,np,nlevp), zm_init(np,np,nlev) + + if (qsize < 3 .and. bubble_moist) then + call abortmp('planar moist bubble requires at least 3 tracers') + endif + + if (hybrid%masterthread) then + write(iulog,*) 'initializing hot bubble with' + print *, 'Lx, Ly =', Lx, Ly + print *, 'Sx, Sy =', Sx, Sy + print *, 'bubble_T0', bubble_T0 + print *, 'bubble_dT', bubble_dT + print *, 'bubble_xycenter', bubble_xycenter + print *, 'bubble_zcenter', bubble_zcenter + print *, 'bubble_ztop', bubble_ztop + print *, 'bubble_xyradius', bubble_xyradius + print *, 'bubble_zradius', bubble_zradius + print *, 'bubble_cosine', bubble_cosine + print *, 'bubble_moist', bubble_moist + print *, 'bubble_moist_dq', bubble_moist_dq + print *, 'bubble_prec_type (0 is Kessler (default), 1 is RJ)', bubble_prec_type + endif + + !for the background state + !evenly spaced with reversed indexing, just like in homme eta coord + call get_evenly_spaced_z(zi,zm,0.0_rl,bubble_ztop) + + !for the background state + do k=1,nlevp + Ti(k) = bubble_T0 - zi(k)*g/cp + pi(k) = p0*( Ti(k)/bubble_T0 )**(1.0/kappa) + enddo + do k=1,nlev + Tm(k) = bubble_T0 - zm(k)*g/cp + pm(k) = p0*( Tm(k)/bubble_T0 )**(1.0/kappa) + dpm(k)=(pi(k+1)-pi(k)) + enddo + ! this T, z, pressure, and q1 (set below) together do not obey + ! homme's EOS, meaning if one to compute phi from EOS, it won't match z + ! one could have zi replaced with output from phi_from_eos(), + ! but it won't satisfy model-independent setup. + ! instead, we use uniform zi as above. + + !create hybrid coords now from pi + !note that this code should not depend on partitioning + !it depends on ref pressure profile, whcih is init-ed in the same way for all elements/gll +#if 0 + !almost sigma + ai(:) = 0.0; bi(:) = 0.0 + ai(1) = pi(1)/p0; bi(nlevp) = one + + do k=2,nlevp + bi(k) = pi(k)/pi(nlevp) + enddo +#else + !old version, hybrid + ai(:) = 0.0; bi(:) = 0.0 + ai(1) = pi(1)/p0; bi(nlevp) = 1.0 + + do k=2,nlev + bi(k) = 1.0 - zi(k)/zi(1) + !restore ai frop given pressure + ai(k)=(pi(k)-bi(k)*pi(nlevp))/p0 + enddo +#endif + + hvcoord%hyai = ai; hvcoord%hybi = bi + + !set : hyam hybm + hvcoord%hyam = 0.5 *(ai(2:nlev+1) + ai(1:nlev)) + hvcoord%hybm = 0.5 *(bi(2:nlev+1) + bi(1:nlev)) + + !call set_layer_locations: sets etam, etai, dp0, checks that Am=ai/2+ai/2 + call set_layer_locations(hvcoord, .true., hybrid%masterthread) + + !for the background state + if (bubble_moist) then + !build specific humidity at saturation qs as in dcmip2016-kessler + do k=1,nlevp + qi_s(k) = bubble_const1 / pi(k) * exp( bubble_const2 * (Ti(k) - bubble_const3) / ( Ti(k) - bubble_const4 ) ) + enddo + else + qi_s(1:nlevp) = 0.0 + endif + + do k=1,nlev + qm_s(k)=(qi_s(k+1)+qi_s(k)) / 2.0 + enddo + + !set some of element-size arrays + zero_mid_init(:,:,:) = 0.0 + zero_int_init(:,:,:) = 0.0 + do j=1,np; do i=1,np; + zi_init(i,j,1:nlevp) = zi(1:nlevp) + zm_init(i,j,1:nlev) = zm(1:nlev) + ps_init(i,j) = pi(nlevp) + p_init(i,j,1:nlev) = pm(1:nlev) + dp_init(i,j,1:nlev) = dpm(1:nlev) + enddo; enddo + + ! set the rest of conditions for each element + do ie = nets,nete + do j=1,np; do i=1,np + + ! get horizontal coordinates at column i,j + x = elem(ie)%spherep(i,j)%lon; y = elem(ie)%spherep(i,j)%lat + + do k=1,nlevp + + !intermediate value + rr=(x-bubble_xycenter) *(x-bubble_xycenter) / bubble_xyradius / bubble_xyradius + & + (zi(k)-bubble_zcenter)*(zi(k)-bubble_zcenter) / bubble_zradius / bubble_zradius + + if (planar_slice .eqv. .true.) then + !no y dependence + rr=sqrt(rr) + else + rr =sqrt( rr + & + (y-bubble_xycenter) * (y-bubble_xycenter) / bubble_xyradius / bubble_xyradius ) + endif + + !set pot. temperature on interfaces + if ( rr < 1.0 ) then + + if (bubble_cosine) then + offset = cos(rr*dd_pi / 2.0 ) + th0(k) = bubble_T0 + bubble_dT * offset + qi(k) = qi_s(k) + bubble_moist_dq * offset + else + !0/1 nonsmooth function + th0(k) = bubble_T0 + bubble_dT + qi(k) = qi_s(k) + bubble_moist_dq + endif + + else + + !set to reference profile + th0(k) = bubble_T0 + qi(k) = qi_s(k) + + endif + + enddo ! k loop + + !set theta on midlevels and then T from theta, exner + th0m(1:nlev) = (th0(1:nlev) + th0(2:nlevp) ) / 2.0 + t_init(i,j,1:nlev) = th0m(1:nlev) * ( pm(1:nlev)/p0 )**kappa + + !set Q before set_elem_state + if (bubble_moist) then + do k=1,nlev + elem(ie)%state%Q(i,j,k,1) = ( qi(k) + qi(k+1) ) / 2.0 + enddo + else + elem(ie)%state%Q(i,j,:,1) = 0.0 + end if + + !call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,zm,zi,g,elem,n0,n1,ntQ) + call set_elem_state(zero_mid_init,zero_mid_init,zero_mid_init, & + zero_int_init,t_init,ps_init,zero_mid_init(:,:,1), & + p_init,dp_init,zm_init,zi_init,g,elem(ie),1,3,-1) + + enddo; enddo !i,j loop + enddo !ie loop + + !indexing of Q, Qdp + !Q (np,np,nlev,qsize_d) + !Qdp (np,np,nlev,qsize_d,2) + if (bubble_moist) then + ii=2 + else + ii=1 + endif + + do ie = nets,nete + elem(ie)%fcor(:,:) = f + + !all but vapor + elem(ie)%state%Q(:,:,:,ii:qsize) = 0.0 + + !sets hydro phi from (perturbed) theta and pressure, checks for hydrostatic balance after that, saves a state + !call tests_finalize(elem(ie),hvcoord) + enddo + +end subroutine bubble_init + + + ! planar baroclinic instability subroutine planar_baroclinic_instab_init(elem,hybrid,hvcoord,nets,nete) @@ -196,17 +418,16 @@ SUBROUTINE gravity_wave (x,y,p,z,zcoords,hybrid_eta,hyam,hybm,d,u,v,w,t,t_mean,p !----------------------------------------------------------------------- ! input/output params parameters at given location !----------------------------------------------------------------------- - - real(rl), intent(in) :: x ! x (m) + real(rl), intent(in) :: x ! x (m) real(rl), intent(in) :: y ! y (m) real(rl), intent(inout) :: z ! Height (m) real(rl), intent(in) :: hyam ! A coefficient for hybrid-eta coordinate, at model level midpoint real(rl), intent(in) :: hybm ! B coefficient for hybrid-eta coordinate, at model level midpoint logical, intent(in) :: hybrid_eta ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used - real(rl), intent(inout) :: p ! Pressure (Pa) - integer, intent(in) :: zcoords ! 0 or 1 see below + real(rl), intent(inout) :: p ! Pressure (Pa) + integer, intent(in) :: zcoords ! 0 or 1 see below real(rl), intent(in) :: d ! Width for Pert - real(rl), intent(out) :: u ! x-dir wind (m s^-1) + real(rl), intent(out) :: u ! x-dir wind (m s^-1) real(rl), intent(out) :: v ! y-dir wind (m s^-1) real(rl), intent(out) :: w ! Vertical Velocity (m s^-1) real(rl), intent(out) :: t ! Temperature (K) @@ -217,86 +438,80 @@ SUBROUTINE gravity_wave (x,y,p,z,zcoords,hybrid_eta,hyam,hybm,d,u,v,w,t,t_mean,p real(rl), intent(out) :: rho_mean ! density (kg m^-3) real(rl), intent(out) :: q ! Specific Humidity (kg/kg) - ! if zcoords = 1, then we use z and output z - ! if zcoords = 0, then we use p +! if zcoords = 1, then we use z and output z +! if zcoords = 0, then we use p !----------------------------------------------------------------------- ! test case parameters !----------------------------------------------------------------------- - real(rl), parameter :: & - u0 = 20.d0, & ! Reference Velocity - Tref = 300.d0, & ! Reference Temperature + Potential Temperature - Pref = 100000.d0, & ! Reference PS - ztop = 10000.d0, & ! Model Top + real(rl), parameter :: & + u0 = 20.d0, &! Reference Velocity + Tref = 300.d0, &! Reference Temperature + Potential Temperature + Pref = 100000.d0, &! Reference PS + ztop = 10000.d0, &! Model Top yc = 0.d0, & ! y-loc of Pert Center xc = 0.d0, & ! x-loc of Pert Center delta_theta = 1.d0, & ! Max Amplitude of Pert - Lz = 20000.d0, & ! Vertical Wavelength of Pert + Lz = 20000.d0, & ! Vertical Wavelength of Pert N = 0.01d0, & ! Brunt-Vaisala frequency - N2 = N*N, & ! Brunt-Vaisala frequency Squared + N2 = N*N, &! Brunt-Vaisala frequency Squared bigG = (g*g)/(N2*cp) ! Constant real(rl) :: height ! Model level height - real(rl) :: r2, s ! Shape of perturbation - real(rl) :: Ts ! Surface temperature + real(rl) :: r2, s ! Shape of perturbation + real(rl) :: Ts ! Surface temperature real(rl) :: t_pert ! temperature perturbation real(rl) :: theta_pert ! Pot-temp perturbation !----------------------------------------------------------------------- ! THE VELOCITIES !----------------------------------------------------------------------- + ! Zonal Velocity + u = u0 - ! Zonal Velocity - - u = u0 - - ! Meridional Velocity + ! Meridional Velocity + v = 0.d0 - v = 0.d0 - - ! Vertical Velocity = Vertical Pressure Velocity = 0 - - w = 0.d0 + ! Vertical Velocity = Vertical Pressure Velocity = 0 + w = 0.d0 !----------------------------------------------------------------------- ! PHIS (surface geopotential) !----------------------------------------------------------------------- - - phis = 0.d0 + phis = 0.d0 !----------------------------------------------------------------------- ! SURFACE TEMPERATURE !----------------------------------------------------------------------- - Ts = Tref + Ts = Tref !----------------------------------------------------------------------- ! PS (surface pressure) !----------------------------------------------------------------------- - ps = Pref + ps = Pref !----------------------------------------------------------------------- ! HEIGHT AND PRESSURE AND MEAN TEMPERATURE !----------------------------------------------------------------------- - if (zcoords .eq. 1) then + if (zcoords .eq. 1) then - height = z - p = ps*( (bigG/Ts)*exp(-N2*height/g)+1.d0 - (bigG/Ts) )**(cp/Rd) + height = z + p = ps*( (bigG/Ts)*exp(-N2*height/g)+1.d0 - (bigG/Ts) )**(cp/Rd) - else + else if (hybrid_eta) p = hyam*p0 + hybm*ps - height = (-g/N2)*log( (Ts/bigG)*( (p/ps)**(Rd/cp) - 1.d0 ) + 1.d0 ) - z = height + height = (-g/N2)*log( (Ts/bigG)*( (p/ps)**(Rd/cp) - 1.d0 ) + 1.d0 ) + z = height - endif + endif - t_mean = bigG*(1.d0 - exp(N2*height/g))+ Ts*exp(N2*height/g) + t_mean = bigG*(1.d0 - exp(N2*height/g))+ Ts*exp(N2*height/g) !----------------------------------------------------------------------- ! rho (density), unperturbed using the background temperature t_mean !----------------------------------------------------------------------- - - rho_mean = p/(Rd*t_mean) + rho_mean = p/(Rd*t_mean) !----------------------------------------------------------------------- ! POTENTIAL TEMPERATURE PERTURBATION, @@ -309,15 +524,15 @@ SUBROUTINE gravity_wave (x,y,p,z,zcoords,hybrid_eta,hyam,hybm,d,u,v,w,t,t_mean,p if (planar_slice .eqv. .true.) then r2 = (x-xc)**2 else - r2 = (x-xc)**2 + (y-yc)**2 -end if + r2 = (x-xc)**2 + (y-yc)**2 + end if - s = (d**2)/(d**2 + r2) + s = (d**2)/(d**2 + r2) - theta_pert = delta_theta*s*sin(2.d0*DD_PI*height/Lz) + theta_pert = delta_theta*s*sin(2.d0*DD_PI*height/Lz) - t_pert = theta_pert*(p/p0)**(Rd/cp) - t = t_mean + t_pert + t_pert = theta_pert*(p/p0)**(Rd/cp) + t = t_mean + t_pert rho = p/(Rd*t) @@ -325,7 +540,7 @@ SUBROUTINE gravity_wave (x,y,p,z,zcoords,hybrid_eta,hyam,hybm,d,u,v,w,t,t_mean,p ! initialize Q, set to zero !----------------------------------------------------------------------- - q = 0.d0 + q = 0.d0 END SUBROUTINE gravity_wave diff --git a/components/homme/src/test_src/moist_planar.F90 b/components/homme/src/test_src/moist_planar.F90 index 143dbf096a3a..3ee4b760abe5 100644 --- a/components/homme/src/test_src/moist_planar.F90 +++ b/components/homme/src/test_src/moist_planar.F90 @@ -34,7 +34,7 @@ subroutine planar_moist_rising_bubble_init(elem,hybrid,hvcoord,nets,nete) type(hvcoord_t), intent(inout) :: hvcoord ! hybrid vertical coordinates integer, intent(in) :: nets,nete ! start, end element index - call abortmp('planar moist rising bubble not yet implemented') + call abortmp('call planar_rising_bubble with moist options instead') end subroutine planar_moist_rising_bubble_init diff --git a/components/homme/src/theta-l/share/element_ops.F90 b/components/homme/src/theta-l/share/element_ops.F90 index 2cab9a762d52..7d75e7fc3bca 100644 --- a/components/homme/src/theta-l/share/element_ops.F90 +++ b/components/homme/src/theta-l/share/element_ops.F90 @@ -154,8 +154,11 @@ subroutine get_field_i(elem,name,field,hvcoord,nt) field = elem%state%w_i(:,:,1:nlevp,nt) endif case('geo_i'); - !is phinh_i set for HY runs or should here be an abort? - field = elem%state%phinh_i(:,:,1:nlevp,nt) + if(theta_hydrostatic_mode) then + call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,nt),elem%state%dp3d(:,:,:,nt),field) + else + field = elem%state%phinh_i(:,:,1:nlevp,nt) + endif case('mu_i'); call get_dpnh_dp_i(elem,field,hvcoord,nt) case('pnh_i'); diff --git a/components/homme/test/bubble/README b/components/homme/test/bubble/README new file mode 100644 index 000000000000..ce14e59bb40a --- /dev/null +++ b/components/homme/test/bubble/README @@ -0,0 +1,17 @@ + +committed namelists for moist runs 'rain' for nlev=60. precl is very sensitive to anything, including nlev. + + +run ncl like this: + +;;; SET fili, tt, shortc, ztop, fff +;;; as file to read, time frame, shortc name for the plot suffix +;;; ztop is that same as in the run, fff as field + +ncl plot-in-height.ncl fili=\"planar_rising_bubble1.nc\" tt=45 ztop=1500 shortc=\"dry\" fff=\"Th\" +ncl plot-in-height.ncl fili=\"planar_rising_bubble1.nc\" tt=70 ztop=20000 shortc=\"kessler\" fff=\"Th\" +ncl plot-in-height.ncl fili=\"planar_rising_bubble1.nc\" tt=70 ztop=20000 shortc=\"rj\" fff=\"Th\" + + + + diff --git a/components/homme/test/bubble/nh-dry.nl b/components/homme/test/bubble/nh-dry.nl new file mode 100644 index 000000000000..3a33e82c8c35 --- /dev/null +++ b/components/homme/test/bubble/nh-dry.nl @@ -0,0 +1,76 @@ +! +! namelist for planar nonhydrostatic gravity waves +!_______________________________________________________________________ +&ctl_nl + nthreads = 1 + partmethod = 4 ! mesh parition method: 4 = space filling curve + topology = "plane" ! mesh type: planar + geometry = "plane" ! mesh type: planar + test_case = "planar_rising_bubble" ! test identifier + theta_hydrostatic_mode = .false. +! theta_hydrostatic_mode = .true. + transport_alg = 0 + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'dry' + ne_x = 68 ! number of elements in x-dir + ne_y = 4 ! number of elements in y-dir + qsize = 0 ! num tracer fields + nmax = 45000 ! total number of steps: 600s / tstep + statefreq = 200 ! number of steps between screen dumps + restartfreq = -1 ! don't write restart files if < 0 + runtype = 0 ! 0 = new run + qsplit = -1 ! timesteps set via se_tstep, dt_remap_factor, dt_tracer_factor + rsplit = -1 + integration = 'explicit' ! explicit time integration + tstep_type = 9 ! 1 => default method +! tstep_type = 5 ! 1 => default method + tstep = 0.02 !! 20x20x60lev 0.1 ! dynamics timestep + dt_remap_factor = 20 ! remap every 1 time steps + dt_tracer_factor = 1 ! tracers run at dynamics time step + hypervis_order = 2 + hypervis_scaling = 3.2 +! hypervis_scaling = 0.0 + nu = 0.01 !0.216784 +!!!0.02 !!!!default EAM is 3.4e-8 for earth circ=40000 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hyperviz subcycling + hypervis_subcycle_tom = 1 + hypervis_subcycle_q = 1 + planar_slice=.true. + lx = 1000.0 + ly = 1000.0 + sx = -500.0 + sy = -500.0 + bubble_zcenter=200.0 + bubble_ztop=1500.0 + bubble_xyradius=250.0 + bubble_zradius =175.0 + bubble_cosine=.true. + bubble_T0=300.0 + bubble_dT=0.5 +/ +&vert_nl + vform = "ccm" ! vertical coordinate type "ccm"=hybrid pressure/terrain + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 1000, ! steps + output_varnames1 ='T','Th','ps','geo' ! variables to write to file + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ + + diff --git a/components/homme/test/bubble/nh-moist-kessler.nl b/components/homme/test/bubble/nh-moist-kessler.nl new file mode 100644 index 000000000000..13d4bc916ef1 --- /dev/null +++ b/components/homme/test/bubble/nh-moist-kessler.nl @@ -0,0 +1,82 @@ +! +! namelist for planar nonhydrostatic gravity waves +!_______________________________________________________________________ +&ctl_nl + nthreads = 1 + partmethod = 4 ! mesh parition method: 4 = space filling curve + topology = "plane" ! mesh type: planar + geometry = "plane" ! mesh type: planar + test_case = "planar_rising_bubble" ! test identifier + theta_hydrostatic_mode = .false. +! theta_hydrostatic_mode = .true. + transport_alg = 0 + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'moist' + ne_x = 40 !68 ! number of elements in x-dir + ne_y = 4 ! number of elements in y-dir + qsize = 3 ! num tracer fields + nmax = 14000 ! total number of steps: 600s / tstep + statefreq = 100 ! number of steps between screen dumps + restartfreq = -1 ! don't write restart files if < 0 + runtype = 0 ! 0 = new run + qsplit = -1 ! timesteps set via se_tstep, dt_remap_factor, dt_tracer_factor + rsplit = -1 + integration = 'explicit' ! explicit time integration + tstep_type = 9 ! 1 => default method +! tstep_type = 5 ! 1 => default method + tstep = 0.02 !! 20x20x60lev 0.1 ! dynamics timestep + dt_remap_factor = 1 ! remap every 1 time steps + dt_tracer_factor = 1 ! tracers run at dynamics time step + hypervis_order = 2 + hypervis_scaling = 3.0 +! hypervis_scaling = 0.0 + nu = 0.216784 !dry ran with 0.01 +!!!0.02 !!!!default EAM is 3.4e-8 for earth circ=40000 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hyperviz subcycling + hypervis_subcycle_tom = 1 + hypervis_subcycle_q = 1 + se_ftype=2 + limiter_option = 9 + planar_slice=.true. + lx = 20000.0 + ly = 1000.0 + sx = -10000.0 + sy = -500.0 + bubble_zcenter=2000.0 + bubble_ztop=20000.0 + bubble_xyradius=3000.0 + bubble_zradius =1000.0 + bubble_cosine=.true. +! bubble_cosine=.false. +! bubble_moist=.false. + bubble_moist=.true. + bubble_T0=300.0 + bubble_dT=2.0 + bubble_moist_dq=0.2 +/ +&vert_nl + vform = "ccm" ! vertical coordinate type "ccm"=hybrid pressure/terrain + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 200, ! steps + output_varnames1 ='T','Th','ps','geo','Q1','Q2','Q3' ! variables to write to file + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ + + diff --git a/components/homme/test/bubble/nh-moist-rj.nl b/components/homme/test/bubble/nh-moist-rj.nl new file mode 100644 index 000000000000..66d0438640c9 --- /dev/null +++ b/components/homme/test/bubble/nh-moist-rj.nl @@ -0,0 +1,83 @@ +! +! namelist for planar nonhydrostatic gravity waves +!_______________________________________________________________________ +&ctl_nl + nthreads = 1 + partmethod = 4 ! mesh parition method: 4 = space filling curve + topology = "plane" ! mesh type: planar + geometry = "plane" ! mesh type: planar + test_case = "planar_rising_bubble" ! test identifier + theta_hydrostatic_mode = .false. +! theta_hydrostatic_mode = .true. + transport_alg = 0 + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'moist' + ne_x = 68 ! number of elements in x-dir + ne_y = 4 ! number of elements in y-dir + qsize = 3 ! num tracer fields + nmax = 14000 ! total number of steps: 600s / tstep + statefreq = 100 ! number of steps between screen dumps + restartfreq = -1 ! don't write restart files if < 0 + runtype = 0 ! 0 = new run + qsplit = -1 ! timesteps set via se_tstep, dt_remap_factor, dt_tracer_factor + rsplit = -1 + integration = 'explicit' ! explicit time integration + tstep_type = 9 ! 1 => default method +! tstep_type = 5 ! 1 => default method + tstep = 0.016 !! 20x20x60lev 0.1 ! dynamics timestep + dt_remap_factor = 20 ! remap every 1 time steps + dt_tracer_factor = 1 ! tracers run at dynamics time step + hypervis_order = 2 + hypervis_scaling = 3.2 +! hypervis_scaling = 0.0 + nu = 0.01 ! 0.216784 !dry ran with 0.01 +!!!0.02 !!!!default EAM is 3.4e-8 for earth circ=40000 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hyperviz subcycling + hypervis_subcycle_tom = 1 + hypervis_subcycle_q = 1 + se_ftype=2 + limiter_option = 9 + planar_slice=.true. + lx = 20000.0 + ly = 1000.0 + sx = -10000.0 + sy = -500.0 + bubble_zcenter=2000.0 + bubble_ztop=20000.0 + bubble_xyradius=3000.0 + bubble_zradius =1000.0 + bubble_cosine=.true. +! bubble_cosine=.false. +! bubble_moist=.false. + bubble_moist=.true. + bubble_T0=300.0 + bubble_dT=2.0 + bubble_moist_dq=0.2 + bubble_prec_type=1 +/ +&vert_nl + vform = "ccm" ! vertical coordinate type "ccm"=hybrid pressure/terrain + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 200, ! steps + output_varnames1 ='T','Th','ps','geo','Q1','Q2','Q3' ! variables to write to file + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ + + diff --git a/components/homme/test/bubble/plot-in-height.ncl b/components/homme/test/bubble/plot-in-height.ncl new file mode 100644 index 000000000000..2f447e9e2276 --- /dev/null +++ b/components/homme/test/bubble/plot-in-height.ncl @@ -0,0 +1,148 @@ +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" + + +;run like this +; ncl plot-in-height.ncl fili=\"nh-dry-new.nc\" tt=45 ztop=1500 shortc=\"dry-new\" fff=\"Th\" + +begin + +;;; SET fili, tt, shortc, ztop, fff +;;; as file to read, time frame, shortc name for the plot suffix +;;; ztop is that same as in the run, fff as field + +print("file is "+fili) +print("short cut is "+shortc) +print("frame is "+tt) + +gg=9.80616 +P0=100000.0 +;ztop=20000.0 +numz=200 +;fff="Th" + + +;-- define file name + diri = "./" +; fili = "planar_nonhydro_gravity_wave1.nc" +; fili = "planar_rising_bubble1.nc" +; fili = "hydro.nc" + +;-- open file and read variable + f = addfile(diri+fili, "r") + +; var = f->Th(0,:,{0},{-150000:147487}) ;-- first time step, latitude=40N, longitude=0-60E. + lon_t = f->lon ;-- longitude=0-60E + lat_t = f->lat + lev_t = f->lev ;-- currently 17 levels + + hyam = f->hyam + hybm = f->hybm + P0mb = P0*0.01 + ps = f->ps + + lev_p = (/100, 200, 300, 400, 500, 600,700,800,900,1000/) + lev_p@units = "hPa" ; required for vinth2p + +nlev = dimsizes(lev_t) +nlat = dimsizes(lat_t) +nlon = dimsizes(lon_t) + +print(lev_t) + +;tt=19 +field=f->$fff$(tt,:,:) +phi=f->geo(tt,:,:) + +; vp = vinth2p (field(:,:), hyam,hybm, lev_p ,ps(tt,:), 2, P0mb, 1, False ) + + +selectedcols=ind(lat_t.eq.0) + +print(selectedcols) + +print(lat_t(selectedcols)) +print(lon_t(selectedcols)) + +nscols=dimsizes(selectedcols) +nplev=dimsizes(lev_p) + +print("num select cols="+nscols) + +vals_z=new((/numz,nscols/),double) +zvals=new((/numz/),double) + +;vals=new((/nplev,nscols/),double) +;psdummy=new((/nlat,1/),double) +;psdummy(:,0)=ps(tt,:) +;valsdummy=new((/nlev,nlat,1/),double) +;valsdummy(:,:,0)=field(:,:) +;vp = vinth2p (valsdummy, hyam,hybm, lev_p ,psdummy(:,:), 2, P0mb, 1, False ) + +;make z array +dz=(ztop - 0.0)/(numz-1) + +do ii=0,numz-1 +zvals(ii) = (ii)*dz +end do + +print("zvals="+zvals) + +do ii=0,nscols-1 +;print("ii="+ii+", aa="+aa(ii)) +;print("vals(0:nlev-1,ii)"+vals(0:nlev-1,ii)) +;print("f->Th(tt,0:nlev-1,aa(ii))"+f->Th(tt,0:nlev-1,aa(ii))) +;vals(0:nlev-1,ii)=f->Th(tt,0:nlev-1,aa(ii)) +;sample=field(:,aa(i)) +;vals(:,ii)=interpolated(:) + +zvalues=phi(::-1,selectedcols(ii))/gg +values=field(::-1,selectedcols(ii)) + +remapped_z = ftcurv(zvalues,values,zvals) +vals_z(:,ii) = remapped_z + +;if ( ii .ge. 27 .and. ii .le. 31 ) then +;print ("ii="+ii) +;print ("original values "+values) +;print ("original geo "+zvalues/9.81) +;print ("remapped"+remapped_z) +;end if + + +end do + + + +;-- define workstation + wks = gsn_open_wks("png",shortc+"-height-field-"+fff+"-tt-"+tostring(tt)) + +;-- set resources + res = True + res@tiMainString = "frame="+tostring(tt)+" "+fili ;-- title string + +; cmap = read_colormap_file("MPL_GnBu") ;-- read the color map + cmap = read_colormap_file("cmocean_haline") ;-- read the color map + + res@cnFillOn = True ;-- turn on color fill +; res@cnFillPalette = cmap(::-1,:) + res@cnFillPalette = cmap(::1,:) + res@cnLineLabelsOn = False ;-- turns off contour line labels + res@cnInfoLabelOn = False ;-- turns off contour info label + res@lbOrientation = "vertical" ;-- vertical label bar +; res@tiYAxisString = var@long_name+" [hPa]" + res@cnLinesOn = False + + res@cnMaxLevelCount = 30 + res@cnLevelSelectionMode = "EqualSpacedLevels" ; equally spaced contour levels + + res@sfXArray = lon_t(selectedcols) ;-- uses lon_t as plot x-axis +; res@sfYArray = lev_t/100 ;-- uses lev_t in hPa as plot y-axis + res@sfYArray = zvals ;-- uses lev_t in hPa as plot y-axis + +; res@trYReverse = True ;-- reverses y-axis + +;-- generate the plot + plot = gsn_csm_contour(wks,vals_z,res) + +end diff --git a/components/homme/test/reg_test/namelists/thetanh-dry-bubble.nl b/components/homme/test/reg_test/namelists/thetanh-dry-bubble.nl new file mode 100644 index 000000000000..65e598451975 --- /dev/null +++ b/components/homme/test/reg_test/namelists/thetanh-dry-bubble.nl @@ -0,0 +1,76 @@ +! +! namelist for planar nonhydrostatic gravity waves +!_______________________________________________________________________ +&ctl_nl + nthreads = 1 + partmethod = 4 ! mesh parition method: 4 = space filling curve + topology = "plane" ! mesh type: planar + geometry = "plane" ! mesh type: planar + test_case = "planar_rising_bubble" ! test identifier + theta_hydrostatic_mode = .false. +! theta_hydrostatic_mode = .true. + transport_alg = 0 + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'dry' + ne_x = 6 ! number of elements in x-dir + ne_y = 4 !2 does nto work? ! number of elements in y-dir + qsize = 0 ! num tracer fields + nmax = 1000 ! total number of steps: 600s / tstep + statefreq = 200 ! number of steps between screen dumps + restartfreq = -1 ! don't write restart files if < 0 + runtype = 0 ! 0 = new run + qsplit = -1 ! timesteps set via se_tstep, dt_remap_factor, dt_tracer_factor + rsplit = -1 + integration = 'explicit' ! explicit time integration + tstep_type = 9 ! 1 => default method +! tstep_type = 5 ! 1 => default method + tstep = 0.01 !! 20x20x60lev 0.1 ! dynamics timestep + dt_remap_factor = 10 ! remap every 1 time steps + dt_tracer_factor = 1 ! tracers run at dynamics time step + hypervis_order = 2 + hypervis_scaling = 3.2 +! hypervis_scaling = 0.0 + nu = 0.01 !0.216784 +!!!0.02 !!!!default EAM is 3.4e-8 for earth circ=40000 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hyperviz subcycling + hypervis_subcycle_tom = 1 + hypervis_subcycle_q = 1 + planar_slice=.true. + lx = 1000.0 + ly = 1000.0 + sx = -500.0 + sy = -500.0 + bubble_zcenter=200.0 + bubble_ztop=1500.0 + bubble_xyradius=250.0 + bubble_zradius =175.0 + bubble_cosine=.true. + bubble_T0=300.0 + bubble_dT=0.5 +/ +&vert_nl + vform = "ccm" ! vertical coordinate type "ccm"=hybrid pressure/terrain + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 1000, ! steps + output_varnames1 ='T','Th','ps','geo' ! variables to write to file + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ + + diff --git a/components/homme/test/reg_test/namelists/thetanh-moist-bubble.nl b/components/homme/test/reg_test/namelists/thetanh-moist-bubble.nl new file mode 100644 index 000000000000..6583961704cc --- /dev/null +++ b/components/homme/test/reg_test/namelists/thetanh-moist-bubble.nl @@ -0,0 +1,83 @@ +! +! namelist for planar nonhydrostatic gravity waves +!_______________________________________________________________________ +&ctl_nl + nthreads = 1 + partmethod = 4 ! mesh parition method: 4 = space filling curve + topology = "plane" ! mesh type: planar + geometry = "plane" ! mesh type: planar + test_case = "planar_rising_bubble" ! test identifier + theta_hydrostatic_mode = .false. +! theta_hydrostatic_mode = .true. + transport_alg = 0 + theta_advect_form = 1 + vert_remap_q_alg = 10 + moisture = 'moist' + ne_x = 6 ! number of elements in x-dir + ne_y = 4 ! number of elements in y-dir + qsize = 3 ! num tracer fields + nmax = 1000 ! total number of steps: 600s / tstep + statefreq = 200 ! number of steps between screen dumps + restartfreq = -1 ! don't write restart files if < 0 + runtype = 0 ! 0 = new run + qsplit = -1 ! timesteps set via se_tstep, dt_remap_factor, dt_tracer_factor + rsplit = -1 + integration = 'explicit' ! explicit time integration + tstep_type = 9 ! 1 => default method +! tstep_type = 5 ! 1 => default method + tstep = 0.016 !! 20x20x60lev 0.1 ! dynamics timestep + dt_remap_factor = 20 ! remap every 1 time steps + dt_tracer_factor = 1 ! tracers run at dynamics time step + hypervis_order = 2 + hypervis_scaling = 3.2 +! hypervis_scaling = 0.0 + nu = 0.01 ! 0.216784 !dry ran with 0.01 +!!!0.02 !!!!default EAM is 3.4e-8 for earth circ=40000 + nu_top = 0.0 + hypervis_order = 2 ! 2 = hyperviscosity + hypervis_subcycle = 1 ! 1 = no hyperviz subcycling + hypervis_subcycle_tom = 1 + hypervis_subcycle_q = 1 + se_ftype=2 + limiter_option = 9 + planar_slice=.true. + lx = 20000.0 + ly = 1000.0 + sx = -10000.0 + sy = -500.0 + bubble_zcenter=2000.0 + bubble_ztop=20000.0 + bubble_xyradius=3000.0 + bubble_zradius =1000.0 + bubble_cosine=.true. +! bubble_cosine=.false. +! bubble_moist=.false. + bubble_moist=.true. + bubble_T0=300.0 + bubble_dT=2.0 + bubble_moist_dq=0.2 + bubble_prec_type=0 +/ +&vert_nl + vform = "ccm" ! vertical coordinate type "ccm"=hybrid pressure/terrain + vanalytic = 1 ! set vcoords in initialization routine + vtop = 2.73919e-1 ! vertical coordinate at top of atm (z=10000m) +/ +&analysis_nl + output_dir = "./movies/" ! destination dir for netcdf file + output_timeunits = 0, ! 1=days, 2=hours, 0=timesteps + output_frequency = 1000, ! steps + output_varnames1 ='T','Th','ps','geo','Q1','Q2','Q3' ! variables to write to file + interp_type = 0 ! 0=native grid, 1=bilinear + output_type ='netcdf' ! netcdf or pnetcdf + num_io_procs = 1 + interp_nlat = 128 + interp_nlon = 256 + interp_gridtype = 2 ! gauss grid +/ +&prof_inparm + profile_outpe_num = 100 + profile_single_file = .true. +/ + + diff --git a/components/homme/test/reg_test/run_tests/test-list.cmake b/components/homme/test/reg_test/run_tests/test-list.cmake index 5a665486b8b3..961acfda81f3 100644 --- a/components/homme/test/reg_test/run_tests/test-list.cmake +++ b/components/homme/test/reg_test/run_tests/test-list.cmake @@ -38,6 +38,8 @@ SET(HOMME_TESTS thetanh-nhgw-slice.cmake preqx-nhgw-slice.cmake sweqx-dbvor.cmake + thetanh-moist-bubble.cmake + thetanh-dry-bubble.cmake ) IF (${HOMME_ENABLE_COMPOSE}) diff --git a/components/homme/test/reg_test/run_tests/thetanh-dry-bubble.cmake b/components/homme/test/reg_test/run_tests/thetanh-dry-bubble.cmake new file mode 100644 index 000000000000..63bc07196170 --- /dev/null +++ b/components/homme/test/reg_test/run_tests/thetanh-dry-bubble.cmake @@ -0,0 +1,19 @@ +############################################################### +# +# 200 steps of the low-resolution (33x33 elements) planar nonhydrostatic gravity wave test case +# +############################################################### + +# The name of this test (should be the basename of this file) +SET(TEST_NAME thetanh-dry-bubble) +# The specifically compiled executable that this test uses +SET(EXEC_NAME theta-l-nlev20-native) + +SET(NUM_CPUS 4) + +SET(NAMELIST_FILES ${HOMME_ROOT}/test/reg_test/namelists/thetanh-dry-bubble.nl) + +# compare all of these files against baselines: +SET(NC_OUTPUT_FILES planar_rising_bubble1.nc) + + diff --git a/components/homme/test/reg_test/run_tests/thetanh-moist-bubble.cmake b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble.cmake new file mode 100644 index 000000000000..33cfc855bbe7 --- /dev/null +++ b/components/homme/test/reg_test/run_tests/thetanh-moist-bubble.cmake @@ -0,0 +1,19 @@ +############################################################### +# +# 200 steps of the low-resolution (33x33 elements) planar nonhydrostatic gravity wave test case +# +############################################################### + +# The name of this test (should be the basename of this file) +SET(TEST_NAME thetanh-moist-bubble) +# The specifically compiled executable that this test uses +SET(EXEC_NAME theta-l-nlev20-native) + +SET(NUM_CPUS 4) + +SET(NAMELIST_FILES ${HOMME_ROOT}/test/reg_test/namelists/thetanh-moist-bubble.nl) + +# compare all of these files against baselines: +SET(NC_OUTPUT_FILES planar_rising_bubble1.nc) + + diff --git a/components/homme/test_execs/CMakeLists.txt b/components/homme/test_execs/CMakeLists.txt index 0f22e9de851a..b18dc4b56a27 100644 --- a/components/homme/test_execs/CMakeLists.txt +++ b/components/homme/test_execs/CMakeLists.txt @@ -209,6 +209,7 @@ IF(${BUILD_HOMME_THETA}) ADD_SUBDIRECTORY(theta-l-nlev60) ADD_SUBDIRECTORY(theta-l-nlev10-native) ADD_SUBDIRECTORY(theta-l-nlev20-native) +# ADD_SUBDIRECTORY(theta-l-nlev60-native) # ADD_SUBDIRECTORY(theta-l-nlev30-native) # ADD_SUBDIRECTORY(theta-l-nlev64-native) # ADD_SUBDIRECTORY(theta-l-nlev100-native) diff --git a/components/homme/test_execs/theta-l-nlev10-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev10-native/CMakeLists.txt index 0a58bc313090..153a550a8466 100644 --- a/components/homme/test_execs/theta-l-nlev10-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev10-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev10-native theta-l 4 4 10 TRUE TRUE 1 ) +createTestExec(theta-l-nlev10-native theta-l 4 4 10 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev100-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev100-native/CMakeLists.txt index ea3f7b4ab0e1..c56f7874e2e2 100644 --- a/components/homme/test_execs/theta-l-nlev100-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev100-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev100-native theta-l 4 4 100 TRUE TRUE 1 ) +createTestExec(theta-l-nlev100-native theta-l 4 4 100 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev128-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev128-native/CMakeLists.txt index b620797d7fe5..7def5635dd30 100644 --- a/components/homme/test_execs/theta-l-nlev128-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev128-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev128-native theta-l 4 4 128 TRUE TRUE 1 ) +createTestExec(theta-l-nlev128-native theta-l 4 4 128 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev140-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev140-native/CMakeLists.txt index f2eae2bc74b3..e143ff73ad24 100644 --- a/components/homme/test_execs/theta-l-nlev140-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev140-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev140-native theta-l 4 4 140 TRUE TRUE 1 ) +createTestExec(theta-l-nlev140-native theta-l 4 4 140 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev150-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev150-native/CMakeLists.txt index f2d63c6bac69..3c8dc15f97a1 100644 --- a/components/homme/test_execs/theta-l-nlev150-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev150-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev150-native theta-l 4 4 150 TRUE TRUE 1 ) +createTestExec(theta-l-nlev150-native theta-l 4 4 150 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev20-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev20-native/CMakeLists.txt index ad91de5ec627..711ca47ca921 100644 --- a/components/homme/test_execs/theta-l-nlev20-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev20-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev20-native theta-l 4 4 20 TRUE TRUE 1 ) +createTestExec(theta-l-nlev20-native theta-l 4 4 20 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev200-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev200-native/CMakeLists.txt index c7610c243350..6364c3a9394f 100644 --- a/components/homme/test_execs/theta-l-nlev200-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev200-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev200-native theta-l 4 4 200 TRUE TRUE 1 ) +createTestExec(theta-l-nlev200-native theta-l 4 4 200 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev256-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev256-native/CMakeLists.txt index d1c48fe5e5eb..ac73d8012886 100644 --- a/components/homme/test_execs/theta-l-nlev256-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev256-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev256-native theta-l 4 4 256 TRUE TRUE 1 ) +createTestExec(theta-l-nlev256-native theta-l 4 4 256 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev30-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev30-native/CMakeLists.txt index ef15fa1d61fc..5ad227ab856f 100644 --- a/components/homme/test_execs/theta-l-nlev30-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev30-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev30-native theta-l 4 4 30 TRUE TRUE 1 ) +createTestExec(theta-l-nlev30-native theta-l 4 4 30 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev300-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev300-native/CMakeLists.txt index 5a8693c2db3e..7de3eabbac8a 100644 --- a/components/homme/test_execs/theta-l-nlev300-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev300-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev300-native theta-l 4 4 300 TRUE TRUE 1 ) +createTestExec(theta-l-nlev300-native theta-l 4 4 300 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev60-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev60-native/CMakeLists.txt new file mode 100644 index 000000000000..8efc354f9a59 --- /dev/null +++ b/components/homme/test_execs/theta-l-nlev60-native/CMakeLists.txt @@ -0,0 +1,6 @@ +#_______________________________________________________________________ +# create executable needed by planar hgw (native output) + +thetal_setup() +# name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D +createTestExec(theta-l-nlev60-native theta-l 4 4 60 TRUE TRUE 3 ) diff --git a/components/homme/test_execs/theta-l-nlev64-native/CMakeLists.txt b/components/homme/test_execs/theta-l-nlev64-native/CMakeLists.txt index 7cbb8ec40671..2b71f17a9c3b 100644 --- a/components/homme/test_execs/theta-l-nlev64-native/CMakeLists.txt +++ b/components/homme/test_execs/theta-l-nlev64-native/CMakeLists.txt @@ -3,4 +3,4 @@ thetal_setup() # name target NP NC PLEV USE_PIO WITH_ENERGY QSIZE_D -createTestExec(theta-l-nlev64-native theta-l 4 4 64 TRUE TRUE 1 ) +createTestExec(theta-l-nlev64-native theta-l 4 4 64 TRUE TRUE 3 ) diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 1536d4b501af..6a21a74edf85 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -942,6 +942,7 @@ add_default($nl, 'config_use_activeTracers_ttd_forcing'); add_default($nl, 'config_use_debugTracers'); add_default($nl, 'config_reset_debugTracers_near_surface'); add_default($nl, 'config_reset_debugTracers_top_nLayers'); +add_default($nl, 'config_debugTracers_limit_advection'); add_default($nl, 'config_use_debugTracers_surface_bulk_forcing'); add_default($nl, 'config_use_debugTracers_surface_restoring'); add_default($nl, 'config_use_debugTracers_interior_restoring'); @@ -1048,6 +1049,38 @@ add_default($nl, 'config_use_MacroMoleculesTracers_ttd_forcing'); add_default($nl, 'config_use_MacroMoleculesTracers_surface_value'); add_default($nl, 'config_use_MacroMoleculesTracers_sea_ice_coupling'); +############################################# +# Namelist group: tracer_forcing_IFVTracers # +############################################# + +add_default($nl, 'config_use_IFVTracers'); +add_default($nl, 'config_IFVTracers_limit_advection'); +add_default($nl, 'config_IFVTracers_normalize_by_layerThickness'); +add_default($nl, 'config_use_IFVTracers_surface_bulk_forcing'); +add_default($nl, 'config_use_IFVTracers_surface_restoring'); +add_default($nl, 'config_use_IFVTracers_interior_restoring'); +add_default($nl, 'config_use_IFVTracers_exponential_decay'); +add_default($nl, 'config_use_IFVTracers_idealAge_forcing'); +add_default($nl, 'config_use_IFVTracers_ttd_forcing'); +add_default($nl, 'config_use_IFVTracers_surface_value'); +add_default($nl, 'config_use_IFVTracers_sea_ice_coupling'); + +############################################# +# Namelist group: tracer_forcing_IFHTracers # +############################################# + +add_default($nl, 'config_use_IFHTracers'); +add_default($nl, 'config_IFHTracers_limit_advection'); +add_default($nl, 'config_IFHTracers_normalize_by_layerThickness'); +add_default($nl, 'config_use_IFHTracers_surface_bulk_forcing'); +add_default($nl, 'config_use_IFHTracers_surface_restoring'); +add_default($nl, 'config_use_IFHTracers_interior_restoring'); +add_default($nl, 'config_use_IFHTracers_exponential_decay'); +add_default($nl, 'config_use_IFHTracers_idealAge_forcing'); +add_default($nl, 'config_use_IFHTracers_ttd_forcing'); +add_default($nl, 'config_use_IFHTracers_surface_value'); +add_default($nl, 'config_use_IFHTracers_sea_ice_coupling'); + ################################## # Namelist group: AM_globalStats # ################################## @@ -1637,6 +1670,8 @@ my @groups = qw(run_modes tracer_forcing_ecosystracers tracer_forcing_dmstracers tracer_forcing_macromoleculestracers + tracer_forcing_ifvtracers + tracer_forcing_ifhtracers am_globalstats am_surfaceareaweightedaverages am_watermasscensus diff --git a/components/mpas-ocean/bld/build-namelist-group-list b/components/mpas-ocean/bld/build-namelist-group-list index bc2657a2fb58..9754b982c53d 100644 --- a/components/mpas-ocean/bld/build-namelist-group-list +++ b/components/mpas-ocean/bld/build-namelist-group-list @@ -39,6 +39,8 @@ my @groups = qw(run_modes tracer_forcing_ecosystracers tracer_forcing_dmstracers tracer_forcing_macromoleculestracers + tracer_forcing_ifvtracers + tracer_forcing_ifhtracers am_globalstats am_surfaceareaweightedaverages am_watermasscensus diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index d361aa081137..1592b0d9c0bb 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -467,6 +467,7 @@ add_default($nl, 'config_salinity_restoring_under_sea_ice'); add_default($nl, 'config_use_debugTracers'); add_default($nl, 'config_reset_debugTracers_near_surface'); add_default($nl, 'config_reset_debugTracers_top_nLayers'); +add_default($nl, 'config_debugTracers_limit_advection'); add_default($nl, 'config_use_debugTracers_surface_bulk_forcing'); add_default($nl, 'config_use_debugTracers_surface_restoring'); add_default($nl, 'config_use_debugTracers_interior_restoring'); @@ -525,6 +526,38 @@ add_default($nl, 'config_use_MacroMoleculesTracers_ttd_forcing'); add_default($nl, 'config_use_MacroMoleculesTracers_surface_value'); add_default($nl, 'config_use_MacroMoleculesTracers_sea_ice_coupling'); +############################################# +# Namelist group: tracer_forcing_IFVTracers # +############################################# + +add_default($nl, 'config_use_IFVTracers'); +add_default($nl, 'config_IFVTracers_limit_advection'); +add_default($nl, 'config_IFVTracers_normalize_by_layerThickness'); +add_default($nl, 'config_use_IFVTracers_surface_bulk_forcing'); +add_default($nl, 'config_use_IFVTracers_surface_restoring'); +add_default($nl, 'config_use_IFVTracers_interior_restoring'); +add_default($nl, 'config_use_IFVTracers_exponential_decay'); +add_default($nl, 'config_use_IFVTracers_idealAge_forcing'); +add_default($nl, 'config_use_IFVTracers_ttd_forcing'); +add_default($nl, 'config_use_IFVTracers_surface_value'); +add_default($nl, 'config_use_IFVTracers_sea_ice_coupling'); + +############################################# +# Namelist group: tracer_forcing_IFHTracers # +############################################# + +add_default($nl, 'config_use_IFHTracers'); +add_default($nl, 'config_IFHTracers_limit_advection'); +add_default($nl, 'config_IFHTracers_normalize_by_layerThickness'); +add_default($nl, 'config_use_IFHTracers_surface_bulk_forcing'); +add_default($nl, 'config_use_IFHTracers_surface_restoring'); +add_default($nl, 'config_use_IFHTracers_interior_restoring'); +add_default($nl, 'config_use_IFHTracers_exponential_decay'); +add_default($nl, 'config_use_IFHTracers_idealAge_forcing'); +add_default($nl, 'config_use_IFHTracers_ttd_forcing'); +add_default($nl, 'config_use_IFHTracers_surface_value'); +add_default($nl, 'config_use_IFHTracers_sea_ice_coupling'); + ################################## # Namelist group: AM_globalStats # ################################## diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index a97c0f205336..ccef8c1e7d39 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -513,9 +513,10 @@ .false. -.false. +.true. .false. 20 +.false. .false. .false. .false. @@ -565,6 +566,32 @@ .false. .false. + +.true. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. + + +.true. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. +.false. + .true. 'output_interval' diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index a2a93dbfe61d..bc1974776c04 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -2355,6 +2355,14 @@ Valid values: Any positive integer value greater than or equal to 0. Default: Defined in namelist_defaults.xml + +if true, the debugTracers limiter for advection is on + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + if true, surface bulk forcing from coupler is added to surfaceTracerFlux in 'debugTracers' category @@ -2701,6 +2709,188 @@ Default: Defined in namelist_defaults.xml + + + +if true, the 'IFVGRP' category is enabled for the run + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, the IFVTracers limiter for advection is on + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, the IFVTracers Tend variables are divided by layerThickness + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface bulk forcing from coupler is added to surfaceTracerFlux in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface restoring source is applied to tracers in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, interior restoring source is applied to tracers in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, exponential decay source is applied to tracers in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, idealAge forcing source is applied to tracers in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, transit time distribution forcing source is applied to tracers in 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface value is computed for 'IFVGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, couple IFV fields with sea ice + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + + + + +if true, the 'IFHGRP' category is enabled for the run + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, the IFHTracers limiter for advection is on + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, the IFHTracers Tend variables are divided by layerThickness + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface bulk forcing from coupler is added to surfaceTracerFlux in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface restoring source is applied to tracers in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, interior restoring source is applied to tracers in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, exponential decay source is applied to tracers in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, idealAge forcing source is applied to tracers in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, transit time distribution forcing source is applied to tracers in 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, surface value is computed for 'IFHGRP' category + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +if true, couple IFH fields with sea ice + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + ') lines.append(' ') lines.append(' ') + + # These Impulse Function variables only appear if IFH or IFV are on. + # MRP note before merge: the next two variables are not needed if they are the sum of the others. + lines.append(' ') + lines.append(' ') + lines.append(' ') + lines.append(' ') + lines.append(' ') + lines.append(' ') + if ocn_iceberg == 'true': lines.append(' ') lines.append(' ') diff --git a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F index 5e11072a49cd..486edfd223cc 100644 --- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F +++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F @@ -118,6 +118,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: windStressBulkPKGActive logical, pointer :: variableBottomDragPKGActive logical, pointer :: tracerBudgetActive + logical, pointer :: IRFTracerBudgetPKGActive logical, pointer :: landIcePressurePKGActive logical, pointer :: landIceFluxesPKGActive logical, pointer :: landIceCouplingPKGActive @@ -174,6 +175,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_bulk_wind_stress logical, pointer :: config_use_bulk_thickness_flux logical, pointer :: config_compute_active_tracer_budgets + logical, pointer :: config_use_IFHTracers, config_use_IFVTracers character (len=StrKIND), pointer :: config_land_ice_flux_mode type (mpas_pool_iterator_type) :: groupItr @@ -265,6 +267,16 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ tracerBudgetActive = .true. end if + ! + ! test for tracer budget + ! + call mpas_pool_get_package(packagePool, 'IRFTracerBudgetPKGActive', IRFTracerBudgetPKGActive) + call mpas_pool_get_config(configPool, 'config_use_IFHTracers', config_use_IFHTracers) + call mpas_pool_get_config(configPool, 'config_use_IFVTracers', config_use_IFVTracers) + if ( config_use_IFHTracers.or.config_use_IFVTracers) then + IRFTracerBudgetPKGActive = .true. + end if + ! ! Test if chlorophyll, solar zenith angle, and clear sky radiation should be used ! diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 1bc75b28ffab..09cbf1ba7cae 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -1711,6 +1711,31 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp end parallel endif + if ( config_use_IFHTracers .and. config_IFHTracers_normalize_by_layerThickness) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k= 1, maxLevelCell(iCell) + IFHTracersHorAdvTend(:,k,iCell) = IFHTracersHorAdvTend(:,k,iCell) / layerThicknessNew(k,iCell) + IFHTracersHorMixTend(:,k,iCell) = IFHTracersHorMixTend(:,k,iCell) / layerThicknessNew(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + endif + if ( config_use_IFVTracers .and. config_IFVTracers_normalize_by_layerThickness) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k= 1, maxLevelCell(iCell) + IFVTracersVertAdvTend(:,k,iCell) = IFVTracersVertAdvTend(:,k,iCell) / layerThicknessNew(k,iCell) + IFVTracersVertMixTend(:,k,iCell) = IFVTracersVertMixTend(:,k,iCell) / layerThicknessNew(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + endif + call mpas_pool_begin_iteration(tracersPool) do while (mpas_pool_get_next_member(tracersPool, & groupItr) ) @@ -1951,6 +1976,29 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp end parallel end if + if ( config_use_IFHTracers.or.config_use_IFVTracers) then + call mpas_pool_begin_iteration(tracersPool) + do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) + if ( groupItr % memberType == MPAS_POOL_FIELD ) then + ! Reset Impulse Response Functions to original value + if ( trim(groupItr % memberName) == 'IFHTracers'.or. & + trim(groupItr % memberName) == 'IFVTracers' ) then + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroupCur, 1) + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroupNew, 2) + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + tracersGroupNew(:,k,iCell) = tracersGroupCur(:,k,iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + end if + end do + end if + call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, & scratchPool, tracersPool, 2) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index c724f7a9004d..2d1adf5086d4 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -148,7 +148,12 @@ module ocn_diagnostics_variables activeTracerVerticalAdvectionTendency, & activeTracerHorizontalAdvectionEdgeFlux, & activeTracerVerticalAdvectionTopFlux, & - activeTracerVertMixTendency + activeTracerVertMixTendency, & + IFHTracersHorAdvTend, & + IFVTracersVertAdvTend, & + IFHTracersHorMixTend, & + IFVTracersVertMixTend + real (kind=RKIND), pointer :: daysSinceStartOfSim character (len=StrKIND), pointer :: xtime, simulationStartTime @@ -298,6 +303,15 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e salinitySurfaceRestoringTendency) end if + if ( config_use_IFHTracers ) then + call mpas_pool_get_array(diagnosticsPool,'IFHTracersHorAdvTend',IFHTracersHorAdvTend) + call mpas_pool_get_array(diagnosticsPool,'IFHTracersHorMixTend',IFHTracersHorMixTend) + endif + if ( config_use_IFVTracers ) then + call mpas_pool_get_array(diagnosticsPool,'IFVTracersVertAdvTend',IFVTracersVertAdvTend) + call mpas_pool_get_array(diagnosticsPool,'IFVTracersVertMixTend',IFVTracersVertMixTend) + endif + ! Retrieve pointers for configurations where land_ice_flux_mode is enabled if ( trim(config_land_ice_flux_mode) /= 'off' ) then call mpas_pool_get_array(diagnosticsPool, 'topDrag', topDrag) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F index 2746a754d078..0372da8e543f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F @@ -108,7 +108,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:,:), pointer :: normalVelocity - real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroup + real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroup, tracersGroupNew integer, pointer :: nCells, nEdges, nVertices, nVertLevels @@ -172,6 +172,20 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ tracersGroup(:, maxLevelCell(iCell)+1:nVertLevels,iCell) = -1.0e34_RKIND tracersGroup(:, 1:minLevelCell(iCell)-1,iCell) = -1.0e34_RKIND end do + ! Copy Impulse Response Functions from time index 1 to time 2 + ! That is because we copy it from index 2 to 1 at the beginning + ! of each time step to reset the IRFs. + if ( trim(groupItr % memberName) == 'IFVTracers'.or. & + trim(groupItr % memberName) == 'IFHTracers' ) then + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroupNew, 2) + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, nVertLevels + tracersGroupNew(:,k,iCell) = tracersGroup(:,k,iCell) + end do + end do + !$omp end do + end if end if end if end do diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index 2b84140b2485..5fec0a775217 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -1057,6 +1057,20 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & !$omp end parallel endif ! compute budgets + if (trim(groupItr % memberName)=='IFHTracers') then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCellsAll + do k=1,nVertLevels + do n=1,nTracersGroup + IFHTracersHorMixTend(n,k,iCell) = tracerGroupTend(n,k,iCell) + enddo + enddo + enddo + !$omp end do + !$omp end parallel + endif + call ocn_tracer_hmix_tend(meshPool, layerThickEdge, & layerThickness, zMid, tracerGroup, & RediKappa, slopeTriadUp, & @@ -1081,6 +1095,20 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & !$omp end parallel endif ! compute budgets + if (trim(groupItr % memberName)=='IFHTracers') then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCellsAll + do k=1,nVertLevels + do n=1,nTracersGroup + IFHTracersHorMixTend(n,k,iCell) = tracerGroupTend(n,k,iCell) - IFHTracersHorMixTend(n,k,iCell) + enddo + enddo + enddo + !$omp end do + !$omp end parallel + endif + ! ! convert the surface tracer flux into a tracer tendency by ! distributing the flux across some number of surface layers diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F index 899e5876a521..ab751af3946c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F @@ -96,7 +96,8 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & !*** Local variables logical :: & - computeBudgets ! flag to compute active tracer budget + computeBudgets, &! flag to compute active tracer budget + limiterOn ! flag to limit tracer group real (kind=RKIND), dimension(:,:), pointer, contiguous :: & advCoefs, advCoefs3rd ! advection coefficients @@ -136,8 +137,17 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & ! determine whether active tracer budgets should be computed computeBudgets = (budgetDiagsOn .and. tracerGroupName == 'activeTracers') + if (tracerGroupName == 'activeTracers'.or. & + config_IFVTracers_limit_advection.and.(tracerGroupName == 'IFVTracers').or. & + config_IFHTracers_limit_advection.and.(tracerGroupName == 'IFHTracers').or. & + config_debugTracers_limit_advection.and.(tracerGroupName == 'debugTracers')) then + limiterOn = .true. + else + limiterOn = .false. + endif ! call specific advection routine based on choice of monotonicity + ! MRP note, for IRFs limiting is controlled within the mono subroutine if (monotonicOn) then call ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & normalThicknessFlux, w, dt, & @@ -147,13 +157,14 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & minLevelEdgeBot, maxLevelEdgeTop, & highOrderAdvectionMask, & edgeSignOnCell, meshPool, & - computeBudgets) + computeBudgets, tracerGroupName, & + limiterOn) else call ocn_tracer_advection_std_tend(tracers, advCoefs, advCoefs3rd, & nAdvCellsForEdge, advCellsForEdge, normalThicknessFlux, w, layerThickness, & layerThickness, dt, meshPool, tend, minLevelCell, maxLevelCell, & - minLevelEdgeBot, maxLevelEdgeTop, highOrderAdvectionMask, edgeSignOnCell) + minLevelEdgeBot, maxLevelEdgeTop, highOrderAdvectionMask, edgeSignOnCell, tracerGroupName) endif call mpas_timer_stop("tracer adv") diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_mono.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_mono.F index 4eac840d4df1..9cb3010c5cd6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_mono.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_mono.F @@ -98,7 +98,8 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & minLevelEdgeBot, maxLevelEdgeTop, & highOrderAdvectionMask, & edgeSignOnCell, meshPool, & - computeBudgets)!{{{ + computeBudgets, tracerGroupName, & + limiterOn)!{{{ !*** Input/Output parameters @@ -141,6 +142,10 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & meshPool !< [in] Mesh information logical, intent(in) :: & computeBudgets !< [in] Flag to compute active tracer budgets + character (len=*), intent(in) :: & + tracerGroupName !< [in] Tracer name to check if active tracers + logical, intent(in) :: & + limiterOn !< [in] Flag to limit this tracer group !*** Local variables @@ -456,6 +461,12 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end do !$omp end parallel + ! If the limiter is off for this tracer, set the scale factor to 1 for the + ! horizontal so no limiting occurs. + if (.not.limiterOn) then + flxIn(:,:) = 1.0_RKIND + flxOut(:,:) = 1.0_RKIND + endif #ifdef _ADV_TIMERS call mpas_timer_stop('scale factor build') call mpas_timer_start('rescale horiz fluxes') @@ -513,6 +524,18 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end do !$omp end parallel + if (tracerGroupName == 'IFHTracers') then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + IFHTracersHorAdvTend(iTracer,k,iCell) = workTend(k,iCell) + end do + end do ! iCell loop + !$omp end do + !$omp end parallel + end if + #ifdef _ADV_TIMERS call mpas_timer_stop('flux accumulate') call mpas_timer_start('advect diags') @@ -808,6 +831,13 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & #endif nCells = nCellsArray( 1 ) + ! If the limiter is off for this tracer, set the scale factor to 1 for the + ! vertical so no limiting occurs. + if (.not.limiterOn) then + flxIn(:,:) = 1.0_RKIND + flxOut(:,:) = 1.0_RKIND + endif + ! Accumulate the scaled high order vertical tendencies ! and the upwind tendencies !$omp parallel @@ -839,6 +869,18 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end do !$omp end parallel + if (tracerGroupName == 'IFVTracers') then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + IFVTracersVertAdvTend(iTracer,k,iCell) = workTend(k,iCell) + end do + end do ! iCell loop + !$omp end do + !$omp end parallel + end if + ! Compute advection diagnostics and monotonicity checks if requested #ifdef _ADV_TIMERS call mpas_timer_stop('flux accumulate') diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_std.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_std.F index c6aee55e4c94..a546c9c5259a 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_std.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_std.F @@ -26,6 +26,7 @@ module ocn_tracer_advection_std use mpas_pool_routines use mpas_io_units use mpas_threading + use ocn_diagnostics_variables use mpas_tracer_advection_helpers @@ -66,11 +67,13 @@ module ocn_tracer_advection_std subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdvCellsForEdge, advCellsForEdge, &!{{{ normalThicknessFlux, w, layerThickness, verticalCellSize, dt, meshPool, & tend, minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop, & - highOrderAdvectionMask, edgeSignOnCell) + highOrderAdvectionMask, edgeSignOnCell, tracerGroupName) real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: current tracer values real (kind=RKIND), dimension(:,:), intent(in) :: adv_coefs !< Input: Advection coefficients for 2nd order advection real (kind=RKIND), dimension(:,:), intent(in) :: adv_coefs_3rd !< Input: Advection coeffs for blending in 3rd/4th order + character (len=*), intent(in) :: & + tracerGroupName !< [in] Tracer name to check if active tracers integer, dimension(:), intent(in) :: nAdvCellsForEdge !< Input: Number of advection cells for each edge integer, dimension(:,:), intent(in) :: advCellsForEdge !< Input: List of advection cells for each edge real (kind=RKIND), dimension(:,:), intent(in) :: normalThicknessFlux !< Input: Thichness weighted velocitiy @@ -225,9 +228,41 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv end do ! k loop end do ! iCell loop !$omp end do + + if (tracerGroupName == 'IFHTracers') then + !$omp do schedule(runtime) private(invAreaCell1, i, iEdge, k) + do iCell = 1, nCells + invAreaCell1 = 1.0_RKIND / areaCell(iCell) + IFHTracersHorAdvTend(iTracer,:,iCell) = 0.0_RKIND + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = 1, maxLevelEdgeTop(iEdge) + IFHTracersHorAdvTend(iTracer,k,iCell) = IFHTracersHorAdvTend(iTracer,k,iCell) & + + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1 + end do + end do + end do + !$omp end do + end if + + if (tracerGroupName == 'IFVTracers') then + ! Accumulate the scaled high order vertical tendencies. + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsSolve + IFVTracersVertAdvTend(iTracer,:,iCell) = 0.0_RKIND + do k = 1,maxLevelCell(iCell) + IFVTracersVertAdvTend(iTracer,k,iCell) = IFVTracersVertAdvTend(iTracer,k,iCell) & + + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) & + - high_order_vert_flux(k, iCell)) + end do ! k loop + end do ! iCell loop + !$omp end do + end if + !$omp end parallel end do ! iTracer loop + deallocate(tracer_cur, high_order_horiz_flux, high_order_vert_flux) deallocate(verticalDivergenceFactor) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F index 54b3556cf5eb..7bdd14515005 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F @@ -1185,6 +1185,14 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, !$omp end do !$omp end parallel endif + elseif (trim(groupItr % memberName) == 'IFVTracers') then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + IFVTracersVertMixTend(:,:,iCell)=tracersGroup(:,:,iCell) + end do + !$omp end do + !$omp end parallel endif if ( associated(tracersGroup) ) then @@ -1214,6 +1222,15 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, !$omp end do !$omp end parallel endif + elseif (trim(groupItr % memberName) == 'IFVTracers') then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + IFVTracersVertMixTend(:,:,iCell) = & + (tracersGroup(:,:,iCell) - IFVTracersVertMixTend(:,:,iCell)) / dt + end do + !$omp end do + !$omp end parallel endif end if diff --git a/components/mpas-ocean/src/tracer_groups/Registry_IFH.xml b/components/mpas-ocean/src/tracer_groups/Registry_IFH.xml new file mode 100644 index 000000000000..704dc98b4e54 --- /dev/null +++ b/components/mpas-ocean/src/tracer_groups/Registry_IFH.xml @@ -0,0 +1,660 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/components/mpas-ocean/src/tracer_groups/Registry_IFV.xml b/components/mpas-ocean/src/tracer_groups/Registry_IFV.xml new file mode 100644 index 000000000000..13778d696088 --- /dev/null +++ b/components/mpas-ocean/src/tracer_groups/Registry_IFV.xml @@ -0,0 +1,183 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/components/mpas-ocean/src/tracer_groups/Registry_debugTracers.xml b/components/mpas-ocean/src/tracer_groups/Registry_debugTracers.xml index b2cab1dac8e6..253725b10773 100644 --- a/components/mpas-ocean/src/tracer_groups/Registry_debugTracers.xml +++ b/components/mpas-ocean/src/tracer_groups/Registry_debugTracers.xml @@ -11,6 +11,10 @@ description="Integer specifying number of layers at top to reset 2 at end of each timestep. Default is 20, chosen to be near a typical mixed layer depth of 200m." possible_values="Any positive integer value greater than or equal to 0." /> + 1 4 4 - - 96 - 72 + + + 72 run_coupling env_run.xml diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index f9d0547e5112..9d83ded23f77 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -137,36 +137,38 @@ module seq_diag_mct integer(in),parameter :: f_hlatf = 8 ! heat : latent, fusion, snow integer(in),parameter :: f_hioff = 9 ! heat : latent, fusion, frozen runoff integer(in),parameter :: f_hsen =10 ! heat : sensible - integer(in),parameter :: f_hh2ot =11 ! heat : water temperature - integer(in),parameter :: f_wfrz =12 ! water: freezing - integer(in),parameter :: f_wmelt =13 ! water: melting - integer(in),parameter :: f_wrain =14 ! water: precip, liquid - integer(in),parameter :: f_wsnow =15 ! water: precip, frozen - integer(in),parameter :: f_wevap =16 ! water: evaporation - integer(in),parameter :: f_wroff =17 ! water: runoff/flood - integer(in),parameter :: f_wioff =18 ! water: frozen runoff - integer(in),parameter :: f_wirrig =19 ! water: irrigation - integer(in),parameter :: f_wfrz_16O =20 ! water: freezing - integer(in),parameter :: f_wmelt_16O =21 ! water: melting - integer(in),parameter :: f_wrain_16O =22 ! water: precip, liquid - integer(in),parameter :: f_wsnow_16O =23 ! water: precip, frozen - integer(in),parameter :: f_wevap_16O =24 ! water: evaporation - integer(in),parameter :: f_wroff_16O =25 ! water: runoff/flood - integer(in),parameter :: f_wioff_16O =26 ! water: frozen runoff - integer(in),parameter :: f_wfrz_18O =27 ! water: freezing - integer(in),parameter :: f_wmelt_18O =28 ! water: melting - integer(in),parameter :: f_wrain_18O =29 ! water: precip, liquid - integer(in),parameter :: f_wsnow_18O =30 ! water: precip, frozen - integer(in),parameter :: f_wevap_18O =31 ! water: evaporation - integer(in),parameter :: f_wroff_18O =32 ! water: runoff/flood - integer(in),parameter :: f_wioff_18O =33 ! water: frozen runoff - integer(in),parameter :: f_wfrz_HDO =34 ! water: freezing - integer(in),parameter :: f_wmelt_HDO =35 ! water: melting - integer(in),parameter :: f_wrain_HDO =36 ! water: precip, liquid - integer(in),parameter :: f_wsnow_HDO =37 ! water: precip, frozen - integer(in),parameter :: f_wevap_HDO =38 ! water: evaporation - integer(in),parameter :: f_wroff_HDO =39 ! water: runoff/flood - integer(in),parameter :: f_wioff_HDO =40 ! water: frozen runoff + integer(in),parameter :: f_hberg =11 ! heat : data icebergs + integer(in),parameter :: f_hh2ot =12 ! heat : water temperature + integer(in),parameter :: f_wfrz =13 ! water: freezing + integer(in),parameter :: f_wmelt =14 ! water: melting + integer(in),parameter :: f_wrain =15 ! water: precip, liquid + integer(in),parameter :: f_wsnow =16 ! water: precip, frozen + integer(in),parameter :: f_wberg =17 ! water: data icebergs + integer(in),parameter :: f_wevap =18 ! water: evaporation + integer(in),parameter :: f_wroff =19 ! water: runoff/flood + integer(in),parameter :: f_wioff =20 ! water: frozen runoff + integer(in),parameter :: f_wirrig =21 ! water: irrigation + integer(in),parameter :: f_wfrz_16O =22 ! water: freezing + integer(in),parameter :: f_wmelt_16O =23 ! water: melting + integer(in),parameter :: f_wrain_16O =24 ! water: precip, liquid + integer(in),parameter :: f_wsnow_16O =25 ! water: precip, frozen + integer(in),parameter :: f_wevap_16O =26 ! water: evaporation + integer(in),parameter :: f_wroff_16O =27 ! water: runoff/flood + integer(in),parameter :: f_wioff_16O =28 ! water: frozen runoff + integer(in),parameter :: f_wfrz_18O =29 ! water: freezing + integer(in),parameter :: f_wmelt_18O =30 ! water: melting + integer(in),parameter :: f_wrain_18O =31 ! water: precip, liquid + integer(in),parameter :: f_wsnow_18O =32 ! water: precip, frozen + integer(in),parameter :: f_wevap_18O =33 ! water: evaporation + integer(in),parameter :: f_wroff_18O =34 ! water: runoff/flood + integer(in),parameter :: f_wioff_18O =35 ! water: frozen runoff + integer(in),parameter :: f_wfrz_HDO =36 ! water: freezing + integer(in),parameter :: f_wmelt_HDO =37 ! water: melting + integer(in),parameter :: f_wrain_HDO =38 ! water: precip, liquid + integer(in),parameter :: f_wsnow_HDO =39 ! water: precip, frozen + integer(in),parameter :: f_wevap_HDO =40 ! water: evaporation + integer(in),parameter :: f_wroff_HDO =41 ! water: runoff/flood + integer(in),parameter :: f_wioff_HDO =42 ! water: frozen runoff integer(in),parameter :: f_size = f_wioff_HDO ! Total array size of all elements integer(in),parameter :: f_a = f_area ! 1st index for area @@ -186,8 +188,9 @@ module seq_diag_mct (/' area',' hfreeze',' hmelt',' hnetsw',' hlwdn', & ' hlwup',' hlatvap',' hlatfus',' hiroff',' hsen', & - ' hh2otemp',' wfreeze',' wmelt',' wrain',' wsnow', & - ' wevap',' wrunoff',' wfrzrof',' wirrig', & + ' hberg',' hh2otemp',' wfreeze',' wmelt',' wrain', & + ' wsnow',' wberg',' wevap',' wrunoff',' wfrzrof', & + ' wirrig', & ' wfreeze_16O',' wmelt_16O',' wrain_16O',' wsnow_16O', & ' wevap_16O',' wrunoff_16O',' wfrzrof_16O', & ' wfreeze_18O',' wmelt_18O',' wrain_18O',' wsnow_18O', & @@ -305,6 +308,8 @@ module seq_diag_mct integer :: index_i2x_Fioi_melth integer :: index_i2x_Fioi_meltw + integer :: index_i2x_Fioi_bergh + integer :: index_i2x_Fioi_bergw integer :: index_i2x_Fioi_salt integer :: index_i2x_Faii_swnet integer :: index_i2x_Fioi_swpen @@ -1405,8 +1410,8 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa if (first_time) then index_x2o_Fioi_melth = mct_aVect_indexRA(x2o_o,'Fioi_melth') index_x2o_Fioi_meltw = mct_aVect_indexRA(x2o_o,'Fioi_meltw') - index_x2o_Fioi_bergh = mct_aVect_indexRA(x2o_o,'PFioi_bergh', perrWith='quiet') - index_x2o_Fioi_bergw = mct_aVect_indexRA(x2o_o,'PFioi_bergw', perrWith='quiet') + index_x2o_Fioi_bergh = mct_aVect_indexRA(x2o_o,'PFioi_bergh') + index_x2o_Fioi_bergw = mct_aVect_indexRA(x2o_o,'PFioi_bergw') index_x2o_Fioi_salt = mct_aVect_indexRA(x2o_o,'Fioi_salt') index_x2o_Foxx_swnet = mct_aVect_indexRA(x2o_o,'Foxx_swnet') index_x2o_Faxa_lwdn = mct_aVect_indexRA(x2o_o,'Faxa_lwdn') @@ -1461,27 +1466,18 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa ca_i = dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n) nf = f_area ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_o - if (index_x2o_Fioi_bergw == 0) then - nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_meltw,n) - else - nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + & - (ca_o+ca_i)*(x2o_o%rAttr(index_x2o_Fioi_meltw,n)+x2o_o%rAttr(index_x2o_Fioi_bergw,n)) - endif - - if (index_x2o_Fioi_bergh == 0) then - nf = f_hmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_melth,n) - else - nf = f_hmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + & - (ca_o+ca_i)*(x2o_o%rAttr(index_x2o_Fioi_melth,n)+x2o_o%rAttr(index_x2o_Fioi_bergh,n)) - endif - + nf = f_hmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_melth,n) nf = f_hswnet; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Foxx_swnet,n) nf = f_hlwdn ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Faxa_lwdn,n) + nf = f_hberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_bergh,n) + nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_meltw,n) nf = f_wrain ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Faxa_rain,n) nf = f_wsnow ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Faxa_snow,n) + nf = f_wberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_bergw,n) nf = f_wroff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Foxx_rofl,n) nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Foxx_rofi,n) + if ( flds_wiso_ocn )then nf = f_wmelt_16O; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + & @@ -1596,6 +1592,8 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i) if (present(do_i2x)) then index_i2x_Fioi_melth = mct_aVect_indexRA(i2x_i,'Fioi_melth') index_i2x_Fioi_meltw = mct_aVect_indexRA(i2x_i,'Fioi_meltw') + index_i2x_Fioi_bergh = mct_aVect_indexRA(i2x_i,'PFioi_bergh') + index_i2x_Fioi_bergw = mct_aVect_indexRA(i2x_i,'PFioi_bergw') index_i2x_Fioi_swpen = mct_aVect_indexRA(i2x_i,'Fioi_swpen') index_i2x_Faii_swnet = mct_aVect_indexRA(i2x_i,'Faii_swnet') index_i2x_Faii_lwup = mct_aVect_indexRA(i2x_i,'Faii_lwup') @@ -1626,12 +1624,14 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i) ca_i = dom_i%data%rAttr(kArea,n) * frac_i%rAttr(ki,n) nf = f_area ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i nf = f_hmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_melth,n) - nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_meltw,n) nf = f_hswnet; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_swnet,n) & - ca_i*i2x_i%rAttr(index_i2x_Fioi_swpen,n) nf = f_hlwup ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_lwup,n) nf = f_hlatv ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_lat,n) nf = f_hsen ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_sen,n) + nf = f_hberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*i2x_i%rAttr(index_i2x_Fioi_bergh,n) + nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_meltw,n) + nf = f_wberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*i2x_i%rAttr(index_i2x_Fioi_bergw,n) nf = f_wevap ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_evap,n) if ( flds_wiso_ice )then @@ -1698,6 +1698,7 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i) (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_frazil,n)) nf = f_hfrz ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - & (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_q,n)) + if ( flds_wiso_ice_x2i )then nf = f_wrain_16O; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + &