diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..403ab42569278283f6742ea0944cfbac508d4a8a
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,16 @@
+# configurations & testcases
+
+*/work_cfgs.txt
+cfgs/*/BLD
+cfgs/*/EXP00
+cfgs/*/WORK
+cfgs/*_ST
+tests/*/BLD
+tests/*/EXP00
+tests/*/WORK
+tests/*_ST
+
+# mk
+mk/arch*
+mk/cpp.*
+mk/full_key_list.txt
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
deleted file mode 100644
index 8ea34149a7fcc014b553eb18148b401a84d79d71..0000000000000000000000000000000000000000
--- a/.gitlab-ci.yml
+++ /dev/null
@@ -1,61 +0,0 @@
-# This file is a template, and might need editing before it works on your project.
-# To contribute improvements to CI/CD templates, please follow the Development guide at:
-# https://docs.gitlab.com/ee/development/cicd/templates.html
-# This specific template is located at:
-# https://gitlab.com/gitlab-org/gitlab/-/blob/master/lib/gitlab/ci/templates/Getting-Started.gitlab-ci.yml
-
-# This is a sample GitLab CI/CD configuration file that should run without any modifications.
-# It demonstrates a basic 3 stage CI/CD pipeline. Instead of real tests or scripts,
-# it uses echo commands to simulate the pipeline execution.
-#
-# A pipeline is composed of independent jobs that run scripts, grouped into stages.
-# Stages run in sequential order, but jobs within stages run in parallel.
-#
-# For more information, see: https://docs.gitlab.com/ee/ci/yaml/index.html#stages
-variables:
- DATADIR: "NOT_USED"
-
-stages: # List of stages for jobs, and their order of execution
- - test
-
-sette_test1:
- stage: test
- script:
-# test will failed if any of the line below finished on an exit status not 0
- - echo "start sette1"
- - cd sette
- - echo "./sette_ci_dummy.sh -A (to run it activate it in .gitlab-ci.yml file)" >> log_sette1
- tags:
- - test-ci
- artifacts:
- paths:
- - sette/log_sette1
- when: always
-
-sette_test2:
- stage: test
- script:
-# test will failed if any of the line below finished on an exit status not 0
- - echo "start sette2"
- - cd sette
- - echo "./sette.sh -A -n \"GYRE_PISCES\""
- tags:
- - test-ci
- artifacts:
- paths:
- - sette/log_sette2
- when: always
-
-sette_test3:
- stage: test
- script:
-# test will failed if any of the line below finished on an exit status not 0
- - echo "start sette3"
- - cd sette
- - echo "./sette_ci_dummy.sh" >> log_sette3
- tags:
- - test-ci
- artifacts:
- paths:
- - sette/log_sette3
- when: always
diff --git a/arch/NOC/arch-X86_ARCHER2-Gnu.fcm b/arch/NOC/arch-X86_ARCHER2-Gnu.fcm
new file mode 100644
index 0000000000000000000000000000000000000000..eaf9fbc0c6a8900550ef7ed29956d33662f314a2
--- /dev/null
+++ b/arch/NOC/arch-X86_ARCHER2-Gnu.fcm
@@ -0,0 +1,62 @@
+# compiler options for Archer CRAY XC-30 (using crayftn compiler)
+#
+# NCDF_HOME root directory containing lib and include subdirectories for netcdf4
+# HDF5_HOME root directory containing lib and include subdirectories for HDF5
+# XIOS_HOME root directory containing lib for XIOS
+# OASIS_HOME root directory containing lib for OASIS
+#
+# NCDF_INC netcdf4 include file
+# NCDF_LIB netcdf4 library
+# XIOS_INC xios include file (taken into accound only if key_xios is activated)
+# XIOS_LIB xios library (taken into accound only if key_xios is activated)
+# OASIS_INC oasis include file (taken into accound only if key_oasis3 is activated)
+# OASIS_LIB oasis library (taken into accound only if key_oasis3 is activated)
+#
+# FC Fortran compiler command
+# FCFLAGS Fortran compiler flags
+# FFLAGS Fortran 77 compiler flags
+# LD linker
+# LDFLAGS linker flags, e.g. -L<lib dir> if you have libraries
+# FPPFLAGS pre-processing flags
+# AR assembler
+# ARFLAGS assembler flags
+# MK make
+# USER_INC complete list of include files
+# USER_LIB complete list of libraries to pass to the linker
+# CC C compiler used to compile conv for AGRIF
+# CFLAGS compiler flags used with CC
+#
+# Note that:
+# - unix variables "$..." are accpeted and will be evaluated before calling fcm.
+# - fcm variables are starting with a % (and not a $)
+#
+%NCDF_HOME $NETCDF_DIR
+%HDF5_HOME $HDF5_DIR
+%XIOS_HOME /work/n01/shared/nemo/xios-trunk-gnu
+#OASIS_HOME
+
+%NCDF_INC -I%NCDF_HOME/include -I%HDF5_HOME/include
+%NCDF_LIB -L%HDF5_HOME/lib -L%NCDF_HOME/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 -lz
+%XIOS_INC -I%XIOS_HOME/inc
+%XIOS_LIB -L%XIOS_HOME/lib -lxios
+#OASIS_INC -I%OASIS_HOME/build/lib/mct -I%OASIS_HOME/build/lib/psmile.MPI1
+#OASIS_LIB -L%OASIS_HOME/lib -lpsmile.MPI1 -lmct -lmpeu -lscrip
+
+%CPP cpp -Dkey_nosignedzero
+%FC ftn
+%FCFLAGS -O2 -cpp -fallow-argument-mismatch -fdefault-real-8 -fcray-pointer -ffree-line-length-none
+%FFLAGS -O2 -cpp -fallow-argument-mismatch -fdefault-real-8 -fcray-pointer -ffree-line-length-none
+%LD CC
+%FPPFLAGS -P -traditional
+%LDFLAGS -lmpichf90
+%AR ar
+%ARFLAGS -r
+%MK gmake
+%USER_INC %XIOS_INC %NCDF_INC
+%USER_LIB %XIOS_LIB %NCDF_LIB
+#USER_INC %XIOS_INC %OASIS_INC %NCDF_INC
+#USER_LIB %XIOS_LIB %OASIS_LIB %NCDF_LIB
+
+%CC cc -Wl,"--allow-multiple-definition"
+%CFLAGS -O2 -Wl,"--allow-multiple-definition"
+bld::tool::fc_modsearch -J
diff --git a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
index a417d038a6b8e8d61bf6c5d5b02a1f59af86eaea..24a0667fdb27b1e5c4e673aae7216b6da817e35a 100644
--- a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
+++ b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
@@ -399,19 +399,19 @@
!-----------------------------------------------------------------------
&namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T)
!-----------------------------------------------------------------------
- nn_zpyc = 2 ! pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
- ln_mevar = .true. ! variable (T) or constant (F) mixing efficiency
+ ln_mevar = .false. ! variable (T) or constant (F) mixing efficiency
ln_tsdiff = .true. ! account for differential T/S mixing (T) or not (F)
- cn_dir = './' ! root directory for the iwm data location
+ cn_dir = './' ! root directory for the iwm data location
!___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename !
- sn_mpb = 'int_wave_mix' , -12. , 'mixing_power_bot' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_mpp = 'int_wave_mix' , -12. , 'mixing_power_pyc' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_mpc = 'int_wave_mix' , -12. , 'mixing_power_cri' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_dsb = 'int_wave_mix' , -12. , 'decay_scale_bot' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_dsc = 'int_wave_mix' , -12. , 'decay_scale_cri' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpb = 'zdfiwm_forcing' , -12. , 'power_bot' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpc = 'zdfiwm_forcing' , -12. , 'power_cri' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpn = 'zdfiwm_forcing' , -12. , 'power_nsq' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mps = 'zdfiwm_forcing' , -12. , 'power_sho' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_dsb = 'zdfiwm_forcing' , -12. , 'scale_bot' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_dsc = 'zdfiwm_forcing' , -12. , 'scale_cri' , .false. , .true. , 'yearly' , '' , '' , ''
/
!!======================================================================
!! *** Diagnostics namelists *** !!
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index a397a232e5846115fa1d867bf29b0c81c6b57ca1..5a2400ff42a09b9cc5cbc31b098730e0de359497 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -1307,19 +1307,19 @@
!-----------------------------------------------------------------------
&namzdf_iwm ! internal wave-driven mixing parameterization (ln_zdfiwm =T)
!-----------------------------------------------------------------------
- nn_zpyc = 1 ! pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
- ln_mevar = .true. ! variable (T) or constant (F) mixing efficiency
+ ln_mevar = .false. ! variable (T) or constant (F) mixing efficiency
ln_tsdiff = .true. ! account for differential T/S mixing (T) or not (F)
cn_dir = './' ! root directory for the iwm data location
!___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename !
- sn_mpb = 'NOT USED' , -12. , 'mixing_power_bot' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_mpp = 'NOT USED' , -12. , 'mixing_power_pyc' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_mpc = 'NOT USED' , -12. , 'mixing_power_cri' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_dsb = 'NOT USED' , -12. , 'decay_scale_bot' , .false. , .true. , 'yearly' , '' , '' , ''
- sn_dsc = 'NOT USED' , -12. , 'decay_scale_cri' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpb = 'NOT USED' , -12. , 'power_bot' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpc = 'NOT USED' , -12. , 'power_cri' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mpn = 'NOT USED' , -12. , 'power_nsq' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_mps = 'NOT USED' , -12. , 'power_sho' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_dsb = 'NOT USED' , -12. , 'scale_bot' , .false. , .true. , 'yearly' , '' , '' , ''
+ sn_dsc = 'NOT USED' , -12. , 'scale_cri' , .false. , .true. , 'yearly' , '' , '' , ''
/
!!======================================================================
!! *** Diagnostics namelists *** !!
diff --git a/doc/latex/.svnignore b/doc/latex/.gitignore
similarity index 100%
rename from doc/latex/.svnignore
rename to doc/latex/.gitignore
diff --git a/sette/.gitignore b/sette/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..2bbe23174acbad247dbfbd4524b134af0f5d3896
--- /dev/null
+++ b/sette/.gitignore
@@ -0,0 +1,3 @@
+param.cfg
+job_batch_template
+output.sette
diff --git a/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Cray b/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Cray
index 33dfbf16fa861c4d4a53eef51439de3b38d1046b..a3223fdf3bfba6fb0d5279e7274399fb1c492007 100644
--- a/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Cray
+++ b/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Cray
@@ -1,6 +1,6 @@
#!/bin/bash
#
# A batch script will be generated using:
-# /work/n01/shared/acc/mkslurm_settejob -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template
+# /work/n01/shared/acc/mkslurm_settejob_4.2 -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template
# by prepare_job.sh
#
diff --git a/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Gnu b/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Gnu
new file mode 100644
index 0000000000000000000000000000000000000000..88ef23b40562e07b5e2cc4c664c57bb03589c41d
--- /dev/null
+++ b/sette/BATCH_TEMPLATE/batch-X86_ARCHER2-Gnu
@@ -0,0 +1,6 @@
+#!/bin/bash
+#
+# A batch script will be generated using:
+# /work/n01/shared/nemo/mkslurm_settejob_4.2_Gnu -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template
+# by prepare_job.sh
+#
diff --git a/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI b/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI
new file mode 120000
index 0000000000000000000000000000000000000000..a2ec08987081cb561dc7c6f12d0482f00ed755a0
--- /dev/null
+++ b/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI
@@ -0,0 +1 @@
+batch-X64_KARA_GCC_OMPI
\ No newline at end of file
diff --git a/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI_DEBUG b/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI_DEBUG
new file mode 120000
index 0000000000000000000000000000000000000000..a2ec08987081cb561dc7c6f12d0482f00ed755a0
--- /dev/null
+++ b/sette/BATCH_TEMPLATE/batch-mpmd-X64_KARA_GCC_OMPI_DEBUG
@@ -0,0 +1 @@
+batch-X64_KARA_GCC_OMPI
\ No newline at end of file
diff --git a/sette/all_functions.sh b/sette/all_functions.sh
index 8e817db7cb2fb3ebfed00aa7a674db142e8dd2d9..e8bc906cdcb328724bf82b9cad6a246d57fe495f 100755
--- a/sette/all_functions.sh
+++ b/sette/all_functions.sh
@@ -171,8 +171,8 @@ clean_config() {
# define validation dir
set_valid_dir () {
- REVISION_NB=`git rev-list --abbrev-commit origin | head -1l`
- REV_DATE0="`git log -1 | grep Date | sed -e 's/.*Date: *//' -e's/ +.*$//'`"
+ REVISION_NB=`git -C ${MAIN_DIR} rev-list --abbrev-commit origin | head -1l`
+ REV_DATE0="`git -C ${MAIN_DIR} log -1 | grep Date | sed -e 's/.*Date: *//' -e's/ +.*$//'`"
REV_DATE=`${DATE_CONV}"${REV_DATE0}" +"%y%j"`
REVISION_NB=${REV_DATE}_${REVISION_NB}
if [ ${#REVISION_NB} -eq 0 ]
@@ -185,8 +185,8 @@ set_valid_dir () {
else
echo "value of revision number of NEMOGCM: ${REVISION_NB}"
fi
- localchanges=`git status --short -uno | wc -l`
- if [ $localchanges > 0 ] ; then
+ localchanges=`git -C ${MAIN_DIR} status --short -uno | wc -l`
+ if [[ $localchanges > 0 ]] ; then
REVISION_NB=${REVISION_NB}+
fi
# remove last _ST followed by zero or more alphanumeric characters
diff --git a/sette/param.cfg b/sette/param.default
similarity index 97%
rename from sette/param.cfg
rename to sette/param.default
index 4d1651b432b4fafb0dad41bae4b2062d5e6509d2..adac7711765dfcc7dea674260de3457acf438abd 100644
--- a/sette/param.cfg
+++ b/sette/param.default
@@ -37,6 +37,8 @@ FORCING_DIR=${SETTE_FORCING_DIR:-$WORKDIR/FORCING}
NEMO_VALIDATION_DIR=${SETTE_NEMO_VALIDATION_DIR:-$MAIN_DIR}/NEMO_VALIDATION
# input files storing (namelist, iodef ...) (DO NOT CHANGE)
INPUT_DIR=${CONFIG_DIR}/${NEW_CONF}/EXP00
+# optional custom SETTE tests directory
+#export CUSTOM_DIR=/path/to/custom/sette/tests
# ------------------------------------------------------------------------------------------
#
# RUN setup
diff --git a/sette/prepare_exe_dir.sh b/sette/prepare_exe_dir.sh
index 62856fac01496e40e4502104facd77d5c337b289..1be9693fb337178140688848929f353fe14e1b6c 100755
--- a/sette/prepare_exe_dir.sh
+++ b/sette/prepare_exe_dir.sh
@@ -64,12 +64,15 @@ set -o posix
#-
-cd ${CONFIG_DIR}
-mkdir -p ${NEW_CONF}/${TEST_NAME}
-
# PREPARE EXEC_DIR
#==================
-export EXE_DIR=${CONFIG_DIR}/${NEW_CONF}/${TEST_NAME}
+if [ -z "${CUSTOM_DIR}" ]; then
+ export EXE_DIR=${CONFIG_DIR}/${NEW_CONF}/${TEST_NAME}
+else
+ NEMO_REV=$( git rev-parse --short HEAD 2> /dev/null )
+ export EXE_DIR=${CUSTOM_DIR}/${SETTE_SUB_VAL}_${NEMO_REV}/${NEW_CONF}/${TEST_NAME}
+fi
+mkdir -p ${EXE_DIR}
cp -RL ${CONFIG_DIR}/${NEW_CONF}/EXP00/* ${EXE_DIR}/.
#cat ${SETTE_DIR}/iodef_sette.xml | sed -e"s;DEF_SHARED;${CONFIG_DIR0}/SHARED;" > ${EXE_DIR}/iodef.xml
@@ -80,6 +83,7 @@ COMP_KEYS="`cat ${CONFIG_DIR}/${NEW_CONF}/cpp_${NEW_CONF}.fcm | sed -e 's/.*fppk
echo "Summary of sette environment" > ./sette_config
echo "----------------------------" >> ./sette_config
echo "requested by the command : "$cmd $cmdargs >> ./sette_config
+echo "on branch : "$SETTE_THIS_BRANCH >> ./sette_config
printf "%-33s : %s\n" USING_TIMING $USING_TIMING >> ./sette_config
printf "%-33s : %s\n" USING_ICEBERGS $USING_ICEBERGS >> ./sette_config
printf "%-33s : %s\n" USING_EXTRA_HALO $USING_EXTRA_HALO >> ./sette_config
@@ -90,6 +94,7 @@ printf "%-33s : %s\n" USING_LOOP_FUSION $USING_LOOP_FUSION >> ./sette_con
printf "%-33s : %s\n" USING_XIOS $USING_XIOS >> ./sette_config
printf "%-33s : %s\n" USING_MPMD $USING_MPMD >> ./sette_config
printf "%-33s : %s\n" USING_RK3 $USING_RK3 >> ./sette_config
+printf "%-33s : %s\n" USER_INPUT $USER_INPUT >> ./sette_config
printf "%-33s : %s\n" "Common compile keys added" "$ADD_KEYS" >> ./sette_config
printf "%-33s : %s\n" "Common compile keys deleted" "$DEL_KEYS" >> ./sette_config
printf "%-33s : %s\n" "Compile keys actually used" "${COMP_KEYS}" >> ./sette_config
diff --git a/sette/prepare_job.sh b/sette/prepare_job.sh
index d4339d1e027853ec065e0f9d984afb267a74b643..2cc3471461edb949448d94171373de24e76d55a6 100755
--- a/sette/prepare_job.sh
+++ b/sette/prepare_job.sh
@@ -198,9 +198,12 @@ fi
NB_NODES=$( echo $NB_PROC $NXIO_PROC | awk '{printf("%d",($1 + $2 ) / 16 + 1 )}')
fi
;;
- X86_ARCHER2*)
+ X86_ARCHER2-Cray)
MK_TEMPLATE=$( /work/n01/shared/nemo/mkslurm_settejob_4.2 -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template )
;;
+ X86_ARCHER2-Gnu)
+ MK_TEMPLATE=$( /work/n01/shared/nemo/mkslurm_settejob_4.2_Gnu -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template )
+ ;;
XC40_METO*) #Setup for Met Office XC40 with any compiler
# ocean cores are packed 32 to a node
# If we need more than one node then have to use parallel queue and XIOS must have a node to itself
diff --git a/sette/sette.sh b/sette/sette.sh
index 266a0b7020039499a6760cc37014e6606f3d873e..8724abd3184d66ccb858476b6d8c3c481eb90bcd 100755
--- a/sette/sette.sh
+++ b/sette/sette.sh
@@ -30,11 +30,26 @@ export USING_XIOS='yes' # Default: yes => add key_xios ; use -X
#
export USING_MPMD='yes' # Default: yes => run with detached XIOS servers ; use -A to run in attached (SPMD) mode
# Note: yes also ensures key_xios but -A will not remove it
-export SETTE_SUB_VAL="MAIN" # Default subdirectory below NEMO_VALIDATION_DIR
+export USER_INPUT='yes' # Default: yes => request user input on decisions. For example:
+ # 1. regarding mismatched options
+ # 2. regardin incompatible options
+ # 3. regarding creation of directories
+#
+# Check that git branch is usable
+git branch --show-current >& /dev/null
+if [[ $? == 0 ]] ; then
+ # subdirectory below NEMO_VALIDATION_DIR defaults to branchname
+ export SETTE_SUB_VAL="$(git branch --show-current)"
+ export SETTE_THIS_BRANCH=${SETTE_SUB_VAL}
+else
+ # subdirectory below NEMO_VALIDATION_DIR defaults to "MAIN"
+ export SETTE_SUB_VAL="MAIN"
+ export SETTE_THIS_BRANCH="Unknown"
+fi
# Parse command-line arguments
if [ $# -gt 0 ]; then
- while getopts n:x:v:g:cdrshTqQteiACFNX option; do
+ while getopts n:x:v:g:cdrshTqQteiACFNXu option; do
case $option in
c) export SETTE_CLEAN_CONFIGS='yes'
export SETTE_SYNC_CONFIGS='yes'
@@ -104,6 +119,9 @@ if [ $# -gt 0 ]; then
A) export USING_MPMD='no'
echo "-A: Tasks will be run in attached (SPMD) mode"
echo "";;
+ u) export USER_INPUT='no'
+ echo "-u: sette.sh will not expect any user interaction == no safety net!"
+ echo "";;
h | *) echo 'sette.sh with no arguments (in this case all configuration will be tested with default options)'
echo '-T to set ln_timing false for all non-AGRIF configurations (default: true)'
echo '-t set ln_tile false in all tests that support it (default: true)'
@@ -125,7 +143,9 @@ if [ $# -gt 0 ]; then
echo '-r to execute without waiting to run sette_rpt.sh at the end (useful for chaining sette.sh invocations)'
echo '-d to perform a dryrun to simply report what settings will be used'
echo '-c to clean each configuration'
- echo '-s to synchronise the sette MY_SRC and EXP00 with the reference MY_SRC and EXPREF'; exit 42 ;;
+ echo '-s to synchronise the sette MY_SRC and EXP00 with the reference MY_SRC and EXPREF'
+ echo '-u to run sette.sh without any user interaction. This means no checks on creating'
+ echo ' directories etc. i.e. no safety net!' ; exit 42 ;;
esac
done
shift $((OPTIND - 1))
@@ -135,26 +155,38 @@ fi
#
if [ ${USING_TILING} == "yes" ] ; then
if [ ${USING_EXTRA_HALO} == "no" ] ; then
- while true; do
- read -p "Tiling requires the extra halo but you have used -e to deselect it. Would you like to reselect it? (y/n)?: " yn
- case $yn in
- [Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
- [Nn]* ) echo "Ok, exiting instead"; exit 42;;
- * ) echo "Please answer yes or no.";;
- esac
- done
+ if [ ${USER_INPUT} == "yes" ] ; then
+ while true; do
+ read -p "Tiling requires the extra halo but you have used -e to deselect it. Would you like to reselect it? (y/n)?: " yn
+ case $yn in
+ [Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
+ [Nn]* ) echo "Ok, exiting instead"; exit 42;;
+ * ) echo "Please answer yes or no.";;
+ esac
+ done
+ else
+ # Without user input, the best option is to disable tiling
+ echo "Tiling requires the extra halo but you have used -e to deselect it. Tiling will not be used."
+ USING_TILING="no"
+ fi
fi
fi
if [ ${USING_LOOP_FUSION} == "yes" ] ; then
if [ ${USING_EXTRA_HALO} == "no" ] ; then
- while true; do
- read -p "Loop fusion requires the extra halo but you have used -e to deselect it. Would you like to reselect it? (y/n)?: " yn
- case $yn in
- [Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
- [Nn]* ) echo "Ok, exiting instead"; exit 42;;
- * ) echo "Please answer yes or no.";;
- esac
- done
+ if [ ${USER_INPUT} == "yes" ] ; then
+ while true; do
+ read -p "Loop fusion requires the extra halo but you have used -e to deselect it. Would you like to reselect it? (y/n)?: " yn
+ case $yn in
+ [Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
+ [Nn]* ) echo "Ok, exiting instead"; exit 42;;
+ * ) echo "Please answer yes or no.";;
+ esac
+ done
+ else
+ # Without user input, the best option is to disable loop fusion
+ echo "Loop fusion requires the extra halo but you have used -e to deselect it. Loop fusion will not be used."
+ USING_LOOP_FUSION="no"
+ fi
fi
fi
#
@@ -182,14 +214,20 @@ if [ ${USING_RK3} == "no" ] ; then export DEL_KEYS="${DEL_KEYS}key_RK3 " ; fi
#
if [ ! -d $NEMO_VALIDATION_DIR ] ; then
if [ ${dry_run} -eq 0 ] ; then
- while true; do
- read -p "$NEMO_VALIDATION_DIR does not exist. Do you wish to create it? " yn
- case $yn in
- [Yy]* ) echo "Ok, creating $NEMO_VALIDATION_DIR"; mkdir $NEMO_VALIDATION_DIR; break;;
- [Nn]* ) echo "Ok, exiting instead"; exit 42;;
- * ) echo "Please answer yes or no.";;
- esac
- done
+ if [ ${USER_INPUT} == "yes" ] ; then
+ while true; do
+ read -p "$NEMO_VALIDATION_DIR does not exist. Do you wish to create it? " yn
+ case $yn in
+ [Yy]* ) echo "Ok, creating $NEMO_VALIDATION_DIR"; mkdir $NEMO_VALIDATION_DIR; break;;
+ [Nn]* ) echo "Ok, exiting instead"; exit 42;;
+ * ) echo "Please answer yes or no.";;
+ esac
+ done
+ else
+ # Without user input, carry on regardless
+ echo "$NEMO_VALIDATION_DIR does not exist. It will be created"
+ mkdir $NEMO_VALIDATION_DIR
+ fi
else
echo "$NEMO_VALIDATION_DIR does not exist"
echo "but this is a dry run so it will not be created"
@@ -206,6 +244,7 @@ if [ ${#SETTE_TEST_CONFIGS[@]} -eq 0 ]; then
fi
echo "Carrying out the following tests : ${TEST_TYPES[@]}"
echo "requested by the command : "$cmd $cmdargs
+echo "on branch : "$SETTE_THIS_BRANCH
printf "%-33s : %s\n" USING_TIMING $USING_TIMING
printf "%-33s : %s\n" USING_ICEBERGS $USING_ICEBERGS
printf "%-33s : %s\n" USING_EXTRA_HALO $USING_EXTRA_HALO
@@ -217,6 +256,7 @@ printf "%-33s : %s\n" USING_LOOP_FUSION $USING_LOOP_FUSION
printf "%-33s : %s\n" USING_XIOS $USING_XIOS
printf "%-33s : %s\n" USING_MPMD $USING_MPMD
printf "%-33s : %s\n" USING_RK3 $USING_RK3
+printf "%-33s : %s\n" USER_INPUT $USER_INPUT
printf "%-33s : %s\n" "Common compile keys to be added" "$ADD_KEYS"
printf "%-33s : %s\n" "Common compile keys to be deleted" "$DEL_KEYS"
echo "Validation records to appear under: "$NEMO_VALIDATION_DIR
diff --git a/sette/sette_eval.sh b/sette/sette_eval.sh
index f942f177fa5a75170a13a8a49feacf9b9c39e81e..93d70715a693565cdc967d4a057c7f6bb7b7e81f 100755
--- a/sette/sette_eval.sh
+++ b/sette/sette_eval.sh
@@ -155,11 +155,12 @@ function runcmpres(){
MAIN_DIR=$(dirname $SETTE_DIR)
quiet=0
. ./param.cfg
+ USER_INPUT='yes' # Default: yes => request user input on decisions.
mach=${COMPILER}
# overwrite revision (later) or compiler
if [ $# -gt 0 ]; then
- while getopts r:R:c:v:V:T:qh option; do
+ while getopts r:R:c:v:V:T:quh option; do
case $option in
c) mach=$OPTARG;;
r) rev=$OPTARG;;
@@ -174,6 +175,7 @@ function runcmpres(){
fi
;;
T) TESTD_ROOT=$OPTARG;;
+ u) USER_INPUT='no';;
h | *) echo ''
echo 'sette_eval.sh : '
echo ' display result for the latest change'
@@ -190,6 +192,7 @@ function runcmpres(){
echo ' -V sub_dir2 :'
echo ' 2nd validation sub-directory below NEMO_VALIDATION_DIR'
echo ' if set the comparison is between two subdirectory trees beneath NEMO_VALIDATION_DIR'
+ echo ' -u to run sette_eval.sh without any user interaction'
echo ' -q : Activate quiet mode - only the number of differing results is returned'
echo ''
exit 42;;
@@ -200,24 +203,41 @@ function runcmpres(){
# if $1 (remaining arguments)
if [[ ! -z $1 ]] ; then rev=$1 ; fi
+ # Check that git branch is usable
+ git branch --show-current >&/dev/null
+ if [[ $? == 0 ]] ; then
+ # subdirectory below NEMO_VALIDATION_DIR defaults to branchname
+ NAM_MAIN="$(git branch --show-current)"
+ else
+ # subdirectory below NEMO_VALIDATION_DIR defaults to "MAIN"
+ NAM_MAIN="MAIN"
+ fi
if [ ! -z $SETTE_SUB_VAL ] ; then
export NEMO_VALIDATION_DIR=$NEMO_VALIDATION_DIR/$SETTE_SUB_VAL
- if [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] && [ ${quiet} -eq 0 ] ; then
+ if [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] && [ ${USER_INPUT} == "yes" ] ; then
while true; do
read -p "$NEMO_VALIDATION_REF/$SETTE_SUB_VAL exists. Do you wish to use it as a reference? " yn
case $yn in
- [Yy]* ) export $NEMO_VALIDATION_REF/$SETTE_SUB_VAL; break;;
- [Nn]* ) echo "Ok, continuing with ${NEMO_VALIDATION_REF}/MAIN as the reference directory"
- export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/MAIN
+ [Yy]* ) export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/$SETTE_SUB_VAL; break;;
+ [Nn]* ) echo "Ok, continuing with ${NEMO_VALIDATION_REF}/${NAM_MAIN} as the reference directory"
+ export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/${NAM_MAIN}
break
;;
* ) echo "Please answer yes or no.";;
esac
done
+ elif [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] ; then
+ # No user input: make a best guess as to intent
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/$SETTE_SUB_VAL
+ elif [ -z $SETTE_SUB_VAL2 ] ; then
+ # No user input: default to branchname or MAIN
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/${NAM_MAIN}
fi
else
- export NEMO_VALIDATION_DIR=${NEMO_VALIDATION_DIR}/MAIN
- export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/MAIN
+ export NEMO_VALIDATION_DIR=${NEMO_VALIDATION_DIR}/${NAM_MAIN}
+ if [ -z $SETTE_SUB_VAL2 ] ; then
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/${NAM_MAIN}
+ fi
fi
NEMO_VALID=${NEMO_VALIDATION_DIR}
NEMO_VALID_REF=${NEMO_VALIDATION_REF}
@@ -240,14 +260,12 @@ rev_date0=`git log -1 | grep Date | sed -e 's/.*Date: *//' -e's/ +.*$//'`
rev_date=`${DATE_CONV}"${rev_date0}" +"%y%j"`
revision=${rev_date}_${revision}
branchname=`git branch --show-current`
-if [ ${quiet} -eq 0 ] ; then
- if [ $localchanges > 0 ] ; then
- echo "Current code is : $branchname @ $revision ( with local changes )"
- lastchange=${revision}+
- else
- echo "Current code is : $branchname @ $revision"
- lastchange=$revision
- fi
+if [[ $localchanges > 0 ]] ; then
+ if [ ${quiet} -eq 0 ] ; then echo "Current code is : $branchname @ $revision ( with local changes )" ; fi
+ lastchange=${revision}+
+else
+ if [ ${quiet} -eq 0 ] ; then echo "Current code is : $branchname @ $revision" ; fi
+ lastchange=$revision
fi
# by default use the current lastchanged revision
@@ -257,7 +275,7 @@ if [ ${quiet} -eq 0 ] ; then
echo ""
echo "SETTE evaluation for : "
echo ""
- if [ $localchanges > 0 ] ; then
+ if [[ $localchanges > 0 ]] ; then
echo " $branchname @ $revision (with local changes)"
else
echo " $branchname @ $revision"
diff --git a/sette/sette_list_avail_rev.sh b/sette/sette_list_avail_rev.sh
index 90684e707bf47b0224668dce374c67f8cbb43fc1..594f4538309789dbaed1ea3d8fb6caf369ce867d 100755
--- a/sette/sette_list_avail_rev.sh
+++ b/sette/sette_list_avail_rev.sh
@@ -48,9 +48,9 @@ lst_rev () {
for rev in $ALLLST
do
if [ -d ${VALSUB}/$rev/${CONFIG} ] ; then
- printf "%-6s " $rev
+ printf "%-14s " $rev
else
- printf "%-5s " "----- "
+ printf "%-14s " "------------ "
fi
done
}
@@ -71,7 +71,7 @@ lst_rev () {
echo ""
printf " List of all avail. rev. in :"${NEMO_VALID}"\n"
printf " is : "
- for dir in `echo $DIRLIST`; do printf "%-6s " $dir ; done
+ for dir in `echo $DIRLIST`; do printf "%-14s " $dir ; done
printf "\n"
# start checking configuration revision
diff --git a/sette/sette_rpt.sh b/sette/sette_rpt.sh
index 5d922bf8fcf7eb1becb345ef5e292f55dcb8cb27..16998096f82cb555f23696ec72672277500a8656 100755
--- a/sette/sette_rpt.sh
+++ b/sette/sette_rpt.sh
@@ -429,11 +429,12 @@ function identictest(){
#
get_dorv
#
- rep=`ls -1rt $vdir/$mach/$dorv/$nam/ | tail -1l`
- f1s=${vdir}/${mach}/${dorv}/${nam}/${rep}/run.stat
- f2s=${vdir}/${mach}/${dorv2}/${nam2}/${rep}/run.stat
+ if [ -d $vdir/$mach/$dorv/$nam ] ; then
+ rep=`ls -1rt $vdir/$mach/$dorv/$nam/ | tail -1l`
+ f1s=${vdir}/${mach}/${dorv}/${nam}/${rep}/run.stat
+ f2s=${vdir}/${mach}/${dorv2}/${nam2}/${rep}/run.stat
#
- if [ -f $f1s ] && [ -f $f2s ] ; then
+ if [ -f $f1s ] && [ -f $f2s ] ; then
cmp -s $f1s $f2s
if [ $? == 0 ]; then
if [ $pass == 0 ]; then
@@ -453,8 +454,11 @@ function identictest(){
read y
fi
fi
- else
+ else
printf "%-27s %-27s %s\n" $nam $nam2 " incomplete test"
+ fi
+ else
+ printf "%-27s %-27s %s\n" " " " " " non-existent test directory"
fi
}
########################### END of function definitions #################################
@@ -467,11 +471,13 @@ function identictest(){
SETTE_DIR=$(cd $(dirname "$0"); pwd)
MAIN_DIR=$(dirname $SETTE_DIR)
. ./param.cfg
+ if [ -z $USER_INPUT ] ; then USER_INPUT='yes' ; fi # Default: yes => request user input on decisions.
+ # (but may br inherited/imported from sette.sh)
mach=${COMPILER}
# overwrite revision (later) or compiler
if [ $# -gt 0 ]; then
- while getopts r:R:c:v:V:h option; do
+ while getopts r:R:c:v:V:uh option; do
case $option in
c) mach=$OPTARG;;
r) rev=$OPTARG;;
@@ -484,6 +490,7 @@ function identictest(){
echo "Requested comparison subdirectory: ${NEMO_VALIDATION_DIR}/${SETTE_SUB_VAL2} does not exist"
fi
;;
+ u) USER_INPUT='no';;
h | *) echo ''
echo 'sette_rpt.sh : '
echo ' display result for the latest change'
@@ -498,6 +505,7 @@ function identictest(){
echo ' -V sub_dir2 :'
echo ' 2nd validation sub-directory below NEMO_VALIDATION_DIR'
echo ' if set the comparison is between two subdirectory trees beneath NEMO_VALIDATION_DIR'
+ echo ' -u to run sette_rpt.sh without any user interaction'
echo ''
exit 42;;
esac
@@ -507,24 +515,41 @@ function identictest(){
# if $1 (remaining arguments)
if [[ ! -z $1 ]] ; then rev=$1 ; fi
+ # Check that git branch is usable
+ git branch --show-current >&/dev/null
+ if [[ $? == 0 ]] ; then
+ # subdirectory below NEMO_VALIDATION_DIR defaults to branchname
+ NAM_MAIN="$(git branch --show-current)"
+ else
+ # subdirectory below NEMO_VALIDATION_DIR defaults to "MAIN"
+ NAM_MAIN="MAIN"
+ fi
if [ ! -z $SETTE_SUB_VAL ] ; then
export NEMO_VALIDATION_DIR=$NEMO_VALIDATION_DIR/$SETTE_SUB_VAL
- if [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] ; then
+ if [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] && [ ${USER_INPUT} == "yes" ] ; then
while true; do
read -p "$NEMO_VALIDATION_REF/$SETTE_SUB_VAL exists. Do you wish to use it as a reference? " yn
case $yn in
- [Yy]* ) export $NEMO_VALIDATION_REF/$SETTE_SUB_VAL; break;;
- [Nn]* ) echo "Ok, continuing with ${NEMO_VALIDATION_REF}/MAIN as the reference directory"
- export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/MAIN
+ [Yy]* ) export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/$SETTE_SUB_VAL; break;;
+ [Nn]* ) echo "Ok, continuing with ${NEMO_VALIDATION_REF}/${NAM_MAIN} as the reference directory"
+ export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/${NAM_MAIN}
break
;;
* ) echo "Please answer yes or no.";;
esac
done
+ elif [ -d $NEMO_VALIDATION_REF/$SETTE_SUB_VAL ] && [ -z $SETTE_SUB_VAL2 ] ; then
+ # No user input: make a best guess as to intent
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/$SETTE_SUB_VAL
+ elif [ -z $SETTE_SUB_VAL2 ] ; then
+ # No user input: default to branchname or MAIN
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/${NAM_MAIN}
fi
else
- export NEMO_VALIDATION_DIR=${NEMO_VALIDATION_DIR}/MAIN
- export NEMO_VALIDATION_REF=${NEMO_VALIDATION_REF}/MAIN
+ export NEMO_VALIDATION_DIR=${NEMO_VALIDATION_DIR}/${NAM_MAIN}
+ if [ -z $SETTE_SUB_VAL2 ] ; then
+ export NEMO_VALIDATION_REF=$NEMO_VALIDATION_REF/${NAM_MAIN}
+ fi
fi
NEMO_VALID=${NEMO_VALIDATION_DIR}
NEMO_VALID_REF=${NEMO_VALIDATION_REF}
@@ -547,7 +572,7 @@ rev_date0=`git log -1 | grep Date | sed -e 's/.*Date: *//' -e's/ +.*$//'`
rev_date=`${DATE_CONV}"${rev_date0}" +"%y%j"`
revision=${rev_date}_${revision}
branchname=`git branch --show-current`
-if [ $localchanges > 0 ] ; then
+if [[ $localchanges > 0 ]] ; then
echo "Current code is : $branchname @ $revision ( with local changes )"
lastchange=${revision}+
else
@@ -561,7 +586,7 @@ lastchange=${rev:-$lastchange}
echo ""
echo "SETTE validation report generated for : "
echo ""
-if [ $localchanges > 0 ] ; then
+if [[ $localchanges > 0 ]] ; then
echo " $branchname @ $revision (with local changes)"
else
echo " $branchname @ $revision"
diff --git a/sette/super_sette_waitq.sh b/sette/super_sette_waitq.sh
new file mode 100755
index 0000000000000000000000000000000000000000..78a7efd93f314e7ad8f558fd399d993af4fb6944
--- /dev/null
+++ b/sette/super_sette_waitq.sh
@@ -0,0 +1,65 @@
+#!/bin/bash
+# set -vx
+# Simple script to robustly run a full suite of SETTE tests
+#
+########################################
+function wait_on_q()
+{
+# SETTE testing on ARCHER2 uses the test queue in which users are limited to
+# 16 queued jobs (including a maximum of 4 running). To prevent sette testing
+# breaching this limit each configuration is processed separately and
+# processing only begins when the user has no more than 12 jobs already queued.
+# The supposition here is that each config forks no more than 4 tests - may need
+# re-revaluating if used for PHYSICS tests.
+#
+# This function checks the queue usage and waits if necessary until the queue
+# has drained sufficiently for the next test.
+NRUN=999
+NIT=0
+BATCH_STAT="squeue -u $USER"
+BATCH_NAME=sette
+echo "Checking queues"
+while [[ $NRUN -gt 12 && $nit -le 1080 ]]; do
+ nit=$((nit+1))
+ NRUN=$( ${BATCH_STAT} | grep ${BATCH_NAME} | wc -l )
+ if [[ $NRUN -gt 12 ]]; then
+ printf "%-3d %s\r" $NRUN 'nemo_sette runs still in queue or running ...';
+ else
+ printf "%-50s\n" "Queues sufficiently drained"
+ return 99
+ fi
+ sleep 10
+done
+echo "Something has gone wrong. Excessive wait time has been exceeded"
+exit
+}
+#
+########################################
+# Start of main script
+########################################
+FULLSET=( ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 AGRIF WED025 GYRE_PISCES SAS ORCA2_ICE_OBS SWG ICE_AGRIF OVERFLOW LOCK_EXCHANGE VORTEX ISOMIP+ )
+#
+GROUP_SETS=( "-r" )
+#
+# These groups sets correspond to the following test regimes:
+#
+# A. complete sets with various combinations of options:
+#
+ printf "%-93s %s\n" "Full tests - <branch_name> (using *_ST config dirs) : " "${GROUP_SETS[0]}"
+#
+# A. Full tests
+for gs in 0
+do
+ for n in `seq 0 1 $(( ${#FULLSET[@]} - 1 ))`
+ do
+ confstr="${FULLSET[$n]}"
+ # compile seperately since the final link sometimes has a bus error on ARCHER2 (which never happens on the 2nd attempt)
+ echo ./sette.sh ${GROUP_SETS[$gs]} -x "COMPILE" -n "$confstr"
+ ./sette.sh ${GROUP_SETS[$gs]} -x "COMPILE" -n "$confstr"
+ # Now run the test (and finish linking if necessary)
+ echo ./sette.sh ${GROUP_SETS[$gs]} -x "RESTART REPRO CORRUPT" -n "$confstr"
+ ./sette.sh ${GROUP_SETS[$gs]} -x "RESTART REPRO CORRUPT" -n "$confstr"
+ wait_on_q
+ done
+done
+exit
diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index 5a4dd65ae4bdf2a92cb3409cccca846cb994f908..b5b66b8a6d129be51693c33859efd507f543476d 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -82,7 +82,7 @@ CONTAINS
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du ! Cd x Du (T-point)
REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du
REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du
- REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd||
+ REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd - uoce||
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: plat ! latent heat flux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
@@ -103,6 +103,8 @@ CONTAINS
#endif
!
REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j
+ REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zsspt
+ REAL(wp), DIMENSION(1:jpi,1:jpj ) :: ztabs, zpre
REAL(wp), DIMENSION(1:jpi ,2:jpka) :: zCF
!
REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_a
@@ -141,6 +143,9 @@ CONTAINS
zrough(:,:) = z0_from_Cd( ght_abl(2), pCd_du(:,:) / MAX( pwndm(:,:), 0.5_wp ) ) ! #LB: z0_from_Cd is define in sbc_phy.F90...
+ ! sea surface potential temperature [K]
+ zsspt(:,:) = theta_exner( psst(:,:)+rt0, pslp_dta(:,:) )
+
!
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
@@ -186,7 +191,7 @@ CONTAINS
IF(jtra == jp_ta) THEN
DO ji = 1,jpi ! surface boundary condition for temperature
zztmp1 = psen(ji, jj)
- zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
+ zztmp2 = psen(ji, jj) * zsspt(ji, jj)
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
@@ -340,26 +345,26 @@ CONTAINS
DO jk = 3, jpkam1
DO ji = 1, jpi
- z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
- z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
- z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
+ z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
+ z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
+ z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
END DO
END DO
- DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
+ DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
!++ Surface boundary condition
z_elem_a( ji, 2 ) = 0._wp
z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
!
- zztmp1 = pcd_du(ji, jj)
- zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )
+ zztmp1 = pcd_du(ji, jj)
+ zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji, jj ) ) * rn_vfac
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
- zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) )
+ zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj ) ) * rn_vfac
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
#endif
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
- u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
+ u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
! idealized test cases only
!IF( ln_topbc_neumann ) THEN
@@ -380,11 +385,10 @@ CONTAINS
!!
!! Matrix inversion
!! ----------------------------------------------------------
- !DO ji = 2, jpi
- DO ji = 1, jpi !!GS: TBI
+ DO ji = 1, jpi !DO ji = 2, jpi !!GS: TBI
zcff = 1._wp / z_elem_b( ji, 2 )
- zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
- u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a)
+ zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
+ u_abl (ji,jj,2,nt_a) = zcff * u_abl ( ji, jj, 2, nt_a )
END DO
DO jk = 3, jpka
@@ -414,9 +418,9 @@ CONTAINS
!
DO jk = 3, jpkam1
DO ji = 1, jpi
- z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
- z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
- z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
+ z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
+ z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
+ z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
END DO
END DO
@@ -426,10 +430,10 @@ CONTAINS
z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
!
zztmp1 = pcd_du(ji, jj)
- zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )
+ zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
- zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) )
+ zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
#endif
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
@@ -455,9 +459,9 @@ CONTAINS
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1, jpi
- zcff = 1._wp / z_elem_b( ji, 2 )
- zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
- v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
+ zcff = 1._wp / z_elem_b( ji, 2 )
+ zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
+ v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
END DO
DO jk = 3, jpka
@@ -494,8 +498,8 @@ CONTAINS
zmsk = msk_abl(ji,jj)
zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 &
& + jp_alp1_dyn * zsig + jp_alp0_dyn
- zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
- ! rn_Dt = rDt_abl / nn_fsbc
+ zcff = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl ! zcff = 1 for masked points
+ ! rn_Dt = rDt_abl / nn_fsbc
zcff = zcff * rest_eq(ji,jj)
u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) &
& + zcff * pu_dta( ji, jj, jk )
@@ -517,8 +521,8 @@ CONTAINS
zmsk = msk_abl(ji,jj)
zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 &
& + jp_alp1_tra * zsig + jp_alp0_tra
- zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
- ! rn_Dt = rDt_abl / nn_fsbc
+ zcff = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl ! zcff = 1 for masked points
+ ! rn_Dt = rDt_abl / nn_fsbc
tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
& + zcff * pt_dta( ji, jj, jk )
@@ -534,7 +538,7 @@ CONTAINS
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
CALL lbc_lnk( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp )
- CALL lbc_lnk( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp , tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp , kfillmode = jpfillnothing ) ! ++++ this should not be needed...
+ !CALL lbc_lnk( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp, tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp, kfillmode = jpfillnothing ) ! ++++ this should not be needed...
!
#if defined key_xios
! 2D & first ABL level
@@ -553,8 +557,8 @@ CONTAINS
IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )
END IF
IF( ln_hpgls_frc ) THEN
- IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
- IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
+ IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", -pgv_dta(:,:,2)/MAX( ABS( fft_abl(:,:)), 2.5e-5_wp)*fft_abl(:,:)/ABS(fft_abl(:,:)) )
+ IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgu_dta(:,:,2)/MAX( ABS( fft_abl(:,:)), 2.5e-5_wp)*fft_abl(:,:)/ABS(fft_abl(:,:)) )
END IF
! 3D (all ABL levels)
IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) )
@@ -576,8 +580,8 @@ CONTAINS
IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) )
END IF
IF( ln_hpgls_frc ) THEN
- IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
- IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
+ IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", -pgv_dta(:,:,2:jpka)/MAX( ABS( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:))), 2.5e-5_wp) * RESHAPE( fft_abl(:,:)/ABS(fft_abl(:,:)), (/jpi,jpj,jpka-1/), fft_abl(:,:)/ABS(fft_abl(:,:))) )
+ IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgu_dta(:,:,2:jpka)/MAX( ABS( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:))), 2.5e-5_wp) * RESHAPE( fft_abl(:,:)/ABS(fft_abl(:,:)), (/jpi,jpj,jpka-1/), fft_abl(:,:)/ABS(fft_abl(:,:))) )
END IF
#endif
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -585,18 +589,19 @@ CONTAINS
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )
- zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )
- zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
- psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention
- pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
+ ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )
+ zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )
+ zpre( ji, jj ) = pres_temp( zhumi, pslp_dta(ji,jj), ght_abl(2), ptpot=ztemp, pta=ztabs( ji, jj ) )
+ zcff = rho_air( ztabs( ji, jj ), zhumi, zpre( ji, jj ) )
+ psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( zsspt(ji,jj) - ztemp ) !GS: negative sign to respect aerobulk convention
+ pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
plat( ji, jj ) = - L_vap( psst(ji,jj) ) * pevp( ji, jj )
rhoa( ji, jj ) = zcff
END_2D
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
- zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) )
- zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) )
+ zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) ) * rn_vfac
+ zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) ) * rn_vfac
END_2D
!
CALL lbc_lnk( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp )
@@ -638,16 +643,17 @@ CONTAINS
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
- DO_2D( 0, 0, 0, 0 )
- ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
- & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
- ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
- & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
- END_2D
- CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
- !
- IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
- & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
+ !DO_2D( 0, 0, 0, 0 )
+ ! ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
+ ! & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
+ ! ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
+ ! & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
+ !END_2D
+ !CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
+ !!
+ !IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
+ ! & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
+
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
@@ -658,10 +664,10 @@ CONTAINS
ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
& + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
- & * ( zztmp1 - pssu_ice(ji,jj) )
+ & * ( zztmp1 - pssu_ice(ji,jj) * rn_vfac )
ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
& + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
- & * ( zztmp2 - pssv_ice(ji,jj) )
+ & * ( zztmp2 - pssv_ice(ji,jj) * rn_vfac )
END_2D
CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
!
diff --git a/src/ABL/par_abl.F90 b/src/ABL/par_abl.F90
index 432dda3b7991e40ef9da15dce6d43ed88cabeca4..35dcda6f381faf91b4e6bd1ba9a7db8020dde3e5 100644
--- a/src/ABL/par_abl.F90
+++ b/src/ABL/par_abl.F90
@@ -60,6 +60,7 @@ MODULE par_abl
REAL(wp), PUBLIC :: rn_ltra_min !: namelist parameter
REAL(wp), PUBLIC :: rn_ltra_max !: namelist parameter
REAL(wp), PUBLIC :: rn_Ric !: critical Richardson number
+ REAL(wp), PUBLIC :: rn_vfac !: multiplicative factor for ocean/ice velocity
!!---------------------------------------------------------------------
!! ABL parameters for the vertical profile of the restoring term
diff --git a/src/ABL/sbcabl.F90 b/src/ABL/sbcabl.F90
index 9c56d719e8c7cccfe96884a072a12c1a4afabf7e..f087614e39166824aa5084b19c44237ccc957963 100644
--- a/src/ABL/sbcabl.F90
+++ b/src/ABL/sbcabl.F90
@@ -71,7 +71,7 @@ CONTAINS
& ln_hpgls_frc, ln_geos_winds, nn_dyn_restore, &
& rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max, &
& nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, &
- & ln_smth_pblh
+ & rn_vfac, ln_smth_pblh
!!---------------------------------------------------------------------
! Namelist namsbc_abl in reference namelist : ABL parameters
@@ -217,8 +217,8 @@ CONTAINS
& + jp_alp1_dyn * jp_bmax + jp_alp0_dyn ) * rDt_abl
IF(lwp) THEN
IF(nn_dyn_restore > 0) THEN
- WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff
- WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1
+ WRITE(numout,*) ' Minimum value for ABL dynamics restoring = ',zcff
+ WRITE(numout,*) ' Maximum value for ABL dynamics restoring = ',zcff1
! Check that restoring coefficients are between 0 and 1
IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) &
& CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
@@ -235,8 +235,8 @@ CONTAINS
zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2 &
& + jp_alp1_tra * jp_bmax + jp_alp0_tra ) * rDt_abl
IF(lwp) THEN
- WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff
- WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1
+ WRITE(numout,*) ' Minimum value for ABL tracers restoring = ',zcff
+ WRITE(numout,*) ' Maximum value for ABL tracers restoring = ',zcff1
! Check that restoring coefficients are between 0 and 1
IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) &
& CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
@@ -369,6 +369,7 @@ CONTAINS
& , utau_ice, vtau_ice & ! =>> out
#endif
& )
+
!!-------------------------------------------------------------------------------------------
!! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
!! time swap is done in abl_stp
diff --git a/src/ICE/icectl.F90 b/src/ICE/icectl.F90
index 568b86fbcdd31f58e0eb4ee91a80d8019532cbcd..e66fe77985a8370c07564c29a8c63fb9b154e4c8 100644
--- a/src/ICE/icectl.F90
+++ b/src/ICE/icectl.F90
@@ -47,9 +47,9 @@ MODULE icectl
! thresold rates for conservation
! these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk
- REAL(wp), PARAMETER :: zchk_m = 2.5e-7 ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
- REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
- REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
+ REAL(wp), PARAMETER :: rchk_m = 2.5e-7 ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
+ REAL(wp), PARAMETER :: rchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
+ REAL(wp), PARAMETER :: rchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
! for drift outputs
CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename
@@ -75,7 +75,7 @@ CONTAINS
!!
!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
!! It prints in ocean.output if there is a violation of conservation at each time-step
- !! The thresholds (zchk_m, zchk_s, zchk_t) determine violations
+ !! The thresholds (rchk_m, rchk_s, rchk_t) determine violations
!! For salt and heat thresholds, ice is considered to have a salinity of 10
!! and a heat content of 3e5 J/kg (=latent heat of fusion)
!!-------------------------------------------------------------------
@@ -145,11 +145,11 @@ CONTAINS
IF( lwp ) THEN
! check conservation issues
- IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zchk3(10) ) &
+ IF( ABS(zdiag_mass) > rchk_m * rn_icechk_glo * zchk3(10) ) &
& WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rDt_ice
- IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zchk3(10) ) &
+ IF( ABS(zdiag_salt) > rchk_s * rn_icechk_glo * zchk3(10) ) &
& WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rDt_ice
- IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zchk3(10) ) &
+ IF( ABS(zdiag_heat) > rchk_t * rn_icechk_glo * zchk3(10) ) &
& WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice
! check negative values
IF( zchk4(1) < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zchk4(1)
@@ -164,9 +164,9 @@ CONTAINS
IF( zchk3(7)>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
& WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zchk3(7)
! check if advection scheme is conservative
- IF( ABS(zchk3(8)) > zchk_m * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
+ IF( ABS(zchk3(8)) > rchk_m * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
& WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zchk3(8) * rDt_ice
- IF( ABS(zchk3(9)) > zchk_t * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
+ IF( ABS(zchk3(9)) > rchk_t * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
& WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zchk3(9) * rDt_ice
ENDIF
!
@@ -182,7 +182,7 @@ CONTAINS
!!
!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
!! It prints in ocean.output if there is a violation of conservation at each time-step
- !! The thresholds (zchk_m, zchk_s, zchk_t) determine the violations
+ !! The thresholds (rchk_m, rchk_s, rchk_t) determine the violations
!! For salt and heat thresholds, ice is considered to have a salinity of 10
!! and a heat content of 3e5 J/kg (=latent heat of fusion)
!!-------------------------------------------------------------------
@@ -204,11 +204,11 @@ CONTAINS
zchk(1:4) = glob_sum_vec( 'icectl', ztmp(:,:,1:4) )
IF( lwp ) THEN
- IF( ABS(zchk(1)) > zchk_m * rn_icechk_glo * zchk(4) ) &
+ IF( ABS(zchk(1)) > rchk_m * rn_icechk_glo * zchk(4) ) &
& WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zchk(1) * rDt_ice
- IF( ABS(zchk(2)) > zchk_s * rn_icechk_glo * zchk(4) ) &
+ IF( ABS(zchk(2)) > rchk_s * rn_icechk_glo * zchk(4) ) &
& WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zchk(2) * rDt_ice
- IF( ABS(zchk(3)) > zchk_t * rn_icechk_glo * zchk(4) ) &
+ IF( ABS(zchk(3)) > rchk_t * rn_icechk_glo * zchk(4) ) &
& WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zchk(3) * rDt_ice
ENDIF
!
@@ -259,20 +259,20 @@ CONTAINS
& + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) &
& - pdiag_fv
- IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel ) ll_stop_m = .TRUE.
+ IF( MAXVAL( ABS(zdiag_mass) ) > rchk_m * rn_icechk_cel ) ll_stop_m = .TRUE.
!
! -- salt diag -- !
zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_Dt_ice &
& + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &
& - pdiag_fs
- IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel ) ll_stop_s = .TRUE.
+ IF( MAXVAL( ABS(zdiag_salt) ) > rchk_s * rn_icechk_cel ) ll_stop_s = .TRUE.
!
! -- heat diag -- !
zdiag_heat = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_Dt_ice &
& + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) &
& - pdiag_ft
- IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
+ IF( MAXVAL( ABS(zdiag_heat) ) > rchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
!
! -- other diags -- !
! a_i < 0
diff --git a/src/ICE/icedia.F90 b/src/ICE/icedia.F90
index 290d77f6935484a257de6336c057159e301709db..154d89de1e4311b194a31b1e645aa65cf523a6d9 100644
--- a/src/ICE/icedia.F90
+++ b/src/ICE/icedia.F90
@@ -33,7 +33,7 @@ MODULE icedia
PUBLIC ice_dia ! called by icestp.F90
PUBLIC ice_dia_init ! called in icestp.F90
- REAL(wp), SAVE :: z1_e1e2 ! inverse of the ocean area
+ REAL(wp), SAVE :: r1_area ! inverse of the ocean area
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
REAL(wp) :: frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot ! global forcing trends
@@ -76,7 +76,7 @@ CONTAINS
ENDIF
IF( kt == nit000 ) THEN
- z1_e1e2 = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
+ r1_area = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
ENDIF
ztmp(:,:,:) = 0._wp ! should be better coded
@@ -152,8 +152,8 @@ CONTAINS
CALL iom_put( 'ibgfrcsal' , frc_sal ) ! sal forcing (psu*km3 equivalent ocean water)
CALL iom_put( 'ibgfrctemtop' , frc_temtop ) ! heat on top of ice/snw/ocean (1.e20 J)
CALL iom_put( 'ibgfrctembot' , frc_tembot ) ! heat on top of ocean(below ice) (1.e20 J)
- CALL iom_put( 'ibgfrchfxtop' , frc_temtop * z1_e1e2 * 1.e-20 * kt*rn_Dt ) ! heat on top of ice/snw/ocean (W/m2)
- CALL iom_put( 'ibgfrchfxbot' , frc_tembot * z1_e1e2 * 1.e-20 * kt*rn_Dt ) ! heat on top of ocean(below ice) (W/m2)
+ CALL iom_put( 'ibgfrchfxtop' , frc_temtop * r1_area * 1.e-20 * kt*rn_Dt ) ! heat on top of ice/snw/ocean (W/m2)
+ CALL iom_put( 'ibgfrchfxbot' , frc_tembot * r1_area * 1.e-20 * kt*rn_Dt ) ! heat on top of ocean(below ice) (W/m2)
CALL iom_put( 'ibgvol_tot' , zbg(6) )
CALL iom_put( 'sbgvol_tot' , zbg(7) )
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index f0a4aab4f92b0173bc2e211e3371dc1586f66036..e5113ac90de67aa1d107f9c1aea3275f2c71f9a2 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -44,8 +44,8 @@ MODULE icedyn_adv_umx
! then interpolate T at u/v points using the upstream scheme
LOGICAL :: ll_prelim = .FALSE. ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
!
- REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
- REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
+ REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
+ REAL(wp) :: r1_120 = 1._wp / 120._wp ! =1/120
!
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: imsk_small, jmsk_small
!
@@ -936,7 +936,7 @@ CONTAINS
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -950,7 +950,7 @@ CONTAINS
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -967,9 +967,9 @@ CONTAINS
zdx4 = zdx2 * zdx2
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
- & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
+ & + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1072,7 +1072,7 @@ CONTAINS
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1085,7 +1085,7 @@ CONTAINS
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1102,9 +1102,9 @@ CONTAINS
zdy4 = zdy2 * zdy2
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
- & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
+ & + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
END_2D
END DO
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index 2a8437e06da2a9a3b3c09806d64f8c3087bbd8b8..306fc5e5c31cc80c0dc63b6d7068eeea184ad1c1 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -788,12 +788,17 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_3 ! local scalars
- REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
+ REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zworka ! temporary array used here
!!clem
LOGICAL :: ln_str_R75
REAL(wp) :: zhi, zcp, rn_pe_rdg
REAL(wp), DIMENSION(jpij) :: zstrength, zaksum ! strength in 1D
!!----------------------------------------------------------------------
+ ! prepare the mask
+ at_i(:,:) = SUM( a_i, dim=3 )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ zmsk(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice
+ END_2D
!
SELECT CASE( nice_str ) !--- Set which ice strength is chosen
@@ -805,7 +810,6 @@ CONTAINS
strength(:,:) = 0._wp
!
! Identify grid cells with ice
- at_i(:,:) = SUM( a_i, dim=3 )
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( at_i(ji,jj) > epsi10 ) THEN
@@ -859,10 +863,10 @@ CONTAINS
ENDIF
!
CASE ( np_strh79 ) !== Hibler(1979)'s method ==!
- strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
+ strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - at_i(:,:) ) ) * zmsk(:,:)
!
CASE ( np_strcst ) !== Constant strength ==!
- strength(:,:) = rn_str
+ strength(:,:) = rn_str * zmsk(:,:)
!
END SELECT
!
@@ -879,7 +883,7 @@ CONTAINS
END_2D
DO_2D( 0, 0, 0, 0 )
- strength(ji,jj) = zworka(ji,jj)
+ strength(ji,jj) = zworka(ji,jj) * zmsk(ji,jj)
END_2D
CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
!
diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90
index 21c9f53a9cb5809ee3c9581b7d4519b7cbf76c53..2efe027c4f8e0cc917a521ec1f35a20e6f226c70 100644
--- a/src/ICE/icedyn_rhg_evp.F90
+++ b/src/ICE/icedyn_rhg_evp.F90
@@ -158,7 +158,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast)
REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast)
!
- REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15
+ REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zmsk00, zmsk15
REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence
@@ -189,6 +189,7 @@ CONTAINS
! for diagnostics and convergence tests
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice
+ zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice
END_2D
IF( nn_rhg_chkcvg > 0 ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -383,10 +384,10 @@ CONTAINS
zdt2 = zdt * zdt
! delta at T points
- zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
+ zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! zmsk is for reducing cpu
! P/delta at T points
- zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+ zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) * zmsk(ji,jj) ! zmsk is for reducing cpu
END_2D
CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
@@ -421,8 +422,10 @@ CONTAINS
ENDIF
! stress at T points (zkt/=0 if landfast)
- zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
- zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+ zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) &
+ & * z1_alph1 * zmsk(ji,jj) ! zmsk is for reducing cpu
+ zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) &
+ & * z1_alph2 * zmsk(ji,jj) ! zmsk is for reducing cpu
END_2D
@@ -449,7 +452,8 @@ CONTAINS
zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) )
! stress at F points (zkt/=0 if landfast)
- zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+ zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) &
+ & * z1_alph2 * zmsk(ji,jj) ! zmsk is for reducing cpu
END_2D
@@ -746,15 +750,15 @@ CONTAINS
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! shear at T points
- pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+ pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj)
! divergence at T points
pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
- & ) * r1_e1e2t(ji,jj)
+ & ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
! delta at T points
- zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
+ zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
@@ -851,8 +855,8 @@ CONTAINS
zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
END_2D
!
- CALL iom_put( 'sig1_pnorm' , zsig1_p )
- CALL iom_put( 'sig2_pnorm' , zsig2_p )
+ CALL iom_put( 'sig1_pnorm' , zsig1_p * zmsk00 )
+ CALL iom_put( 'sig2_pnorm' , zsig2_p * zmsk00 )
DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90
index 9484a039044977c5590f17cc351db80aa003dd58..a973312f04364781465403618f175cf52c6170ed 100644
--- a/src/ICE/iceupdate.F90
+++ b/src/ICE/iceupdate.F90
@@ -339,7 +339,7 @@ CONTAINS
REAL(wp) :: zat_v, zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
!!---------------------------------------------------------------------
- IF( ln_timing ) CALL timing_start('ice_update')
+ IF( ln_timing ) CALL timing_start('iceupdate')
IF( kt == nit000 .AND. lwp ) THEN
WRITE(numout,*)
@@ -391,7 +391,7 @@ CONTAINS
END_2D
CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp ) ! lateral boundary condition
!
- IF( ln_timing ) CALL timing_stop('ice_update')
+ IF( ln_timing ) CALL timing_stop('iceupdate')
!
END SUBROUTINE ice_update_tau
diff --git a/src/OCE/DOM/dommsk.F90 b/src/OCE/DOM/dommsk.F90
index 5c7b801af14a99f641668266aa3920383cb586cc..f153b531a8c835830cde1cad784e3a34609bee49 100644
--- a/src/OCE/DOM/dommsk.F90
+++ b/src/OCE/DOM/dommsk.F90
@@ -160,9 +160,6 @@ CONTAINS
! In case of a coarsened grid, account her for possibly aditionnal
! masked points; these have been read in the mesh file and stored in mbku, mbkv, mbkf
DO_2D( 0, 0, 0, 0 )
- IF (mbku(ji,jj)<=1 ) umask(ji,jj,:) = 0._wp
- IF (mbkv(ji,jj)<=1 ) vmask(ji,jj,:) = 0._wp
- IF (mbkf(ji,jj)<=1 ) fmask(ji,jj,:) = 0._wp
IF ( MAXVAL(umask(ji,jj,:))/=0._wp ) umask(ji,jj,mbku(ji,jj)+1:jpk) = 0._wp
IF ( MAXVAL(vmask(ji,jj,:))/=0._wp ) vmask(ji,jj,mbkv(ji,jj)+1:jpk) = 0._wp
IF ( MAXVAL(fmask(ji,jj,:))/=0._wp ) fmask(ji,jj,mbkf(ji,jj)+1:jpk) = 0._wp
diff --git a/src/OCE/DOM/istate.F90 b/src/OCE/DOM/istate.F90
index ad97be1bbdb732ed2494d4d9c89f8e4d1f46f1f9..a8ea487a05a3e0f8ddae46ca9a39ed1d9c9c270c 100644
--- a/src/OCE/DOM/istate.F90
+++ b/src/OCE/DOM/istate.F90
@@ -79,10 +79,12 @@ CONTAINS
CALL dta_tsd_init ! Initialisation of T & S input data
IF( ln_c1d) CALL dta_uvd_init ! Initialisation of U & V input data (c1d only)
- rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
- rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
- ts (:,:,:,:,Kaa) = 0._wp ! set one for all to 0 at level jpk
- rab_b(:,:,:,: ) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
+ ts (:,:,:,:,Kaa) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
+ IF ( ALLOCATED( rhd ) ) THEN ! SWE, for example, will not have allocated these
+ rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
+ rn2b (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
+ rab_b(:,:,:,: ) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
+ ENDIF
#if defined key_agrif
uu (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization
vv (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index 67cfd73dc8747d71516eccb6c33f4dd32292ec70..87e1d746557589acb22e01749b03156bbc5172cb 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -518,7 +518,7 @@ CONTAINS
WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv
WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia
WRITE(numout,*) ' coefficients :'
- WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aht_ijk_t
+ WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t
WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s'
WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m'
WRITE(numout,*)
diff --git a/src/OCE/SBC/sbc_phy.F90 b/src/OCE/SBC/sbc_phy.F90
index fd527e3885953f11ba3f474d8dff42f97b0f4d2b..cabec70fcfda4a096d045c62979b86cb029f02a8 100644
--- a/src/OCE/SBC/sbc_phy.F90
+++ b/src/OCE/SBC/sbc_phy.F90
@@ -22,20 +22,25 @@ MODULE sbc_phy
USE phycst ! physical constants
IMPLICIT NONE
- !PRIVATE
- PUBLIC !! Haleluja that was the solution
+ PUBLIC !! Haleluja that was the solution for AGRIF
INTEGER , PARAMETER, PUBLIC :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument)
!! (mainly removed from sbcblk.F90)
- REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
- REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
- REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...)
- REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant
+ REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
+ REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+ REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...)
+ REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant
+ REAL(wp), PARAMETER, PUBLIC :: R_gas = 8.314510_wp !: Universal molar gas constant [J/mol/K]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_dryair = 28.9647e-3_wp !: dry air molar mass / molecular weight [kg/mol]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_water = 18.0153e-3_wp !: water molar mass / molecular weight [kg/mol]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_ratio = rmm_water / rmm_dryair
+ REAL(wp), PARAMETER, PUBLIC :: rgamma_dry = R_gas / ( rmm_dryair * rCp_dry ) !: Poisson constant for dry air
+ REAL(wp), PARAMETER, PUBLIC :: rpref = 1.e5_wp !: reference air pressure for exner function [Pa]
!
REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3]
REAL(wp), PARAMETER, PUBLIC :: rho0_w = 1025._wp !: Density of sea-water (ECMWF->1025) [kg/m^3]
@@ -82,6 +87,14 @@ MODULE sbc_phy
MODULE PROCEDURE virt_temp_vctr, virt_temp_sclr
END INTERFACE virt_temp
+ INTERFACE pres_temp
+ MODULE PROCEDURE pres_temp_vctr, pres_temp_sclr
+ END INTERFACE pres_temp
+
+ INTERFACE theta_exner
+ MODULE PROCEDURE theta_exner_vctr, theta_exner_sclr
+ END INTERFACE theta_exner
+
INTERFACE visc_air
MODULE PROCEDURE visc_air_vctr, visc_air_sclr
END INTERFACE visc_air
@@ -97,6 +110,7 @@ MODULE sbc_phy
INTERFACE e_sat_ice
MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr
END INTERFACE e_sat_ice
+
INTERFACE de_sat_dt_ice
MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr
END INTERFACE de_sat_dt_ice
@@ -147,6 +161,8 @@ MODULE sbc_phy
PUBLIC virt_temp
+ PUBLIC pres_temp
+ PUBLIC theta_exner
PUBLIC rho_air
PUBLIC visc_air
PUBLIC L_vap
@@ -157,7 +173,7 @@ MODULE sbc_phy
PUBLIC q_sat
PUBLIC q_air_rh
PUBLIC dq_sat_dt_ice
- !:
+ !
PUBLIC update_qnsol_tau
PUBLIC alpha_sw
PUBLIC bulk_formula
@@ -204,14 +220,140 @@ CONTAINS
!! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
!
END FUNCTION virt_temp_sclr
- !!
+
FUNCTION virt_temp_vctr( pta, pqa )
+
REAL(wp), DIMENSION(jpi,jpj) :: virt_temp_vctr !: virtual temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air [kg/kg]
+
virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
+
END FUNCTION virt_temp_vctr
- !===============================================================================================
+
+
+ FUNCTION pres_temp_sclr( pqspe, pslp, pz, ptpot, pta, l_ice )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION pres_temp ***
+ !!
+ !! ** Purpose : compute air pressure using barometric equation
+ !! from either potential or absolute air temperature
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp) :: pres_temp_sclr ! air pressure [Pa]
+ REAL(wp), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
+ REAL(wp), INTENT(in ) :: pslp ! sea-level pressure [Pa]
+ REAL(wp), INTENT(in ) :: pz ! height above surface [m]
+ REAL(wp), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
+ REAL(wp), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
+ REAL(wp) :: ztpot, zta, zpa, zxm, zmask, zqsat
+ INTEGER :: it, niter = 3 ! iteration indice and number
+ LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
+ LOGICAL :: lice ! sea-ice presence
+
+ IF( PRESENT(ptpot) ) THEN
+ zmask = 1._wp
+ ztpot = ptpot
+ zta = 0._wp
+ ELSE
+ zmask = 0._wp
+ ztpot = 0._wp
+ zta = pta
+ ENDIF
+
+ lice = .FALSE.
+ IF( PRESENT(l_ice) ) lice = l_ice
+
+ zpa = pslp ! air pressure first guess [Pa]
+ DO it = 1, niter
+ zta = ztpot * ( zpa / rpref )**rgamma_dry * zmask + (1._wp - zmask) * zta
+ zqsat = q_sat( zta, zpa, l_ice=lice ) ! saturation specific humidity [kg/kg]
+ zxm = (1._wp - pqspe/zqsat) * rmm_dryair + pqspe/zqsat * rmm_water ! moist air molar mass [kg/mol]
+ zpa = pslp * EXP( -grav * zxm * pz / ( R_gas * zta ) )
+ END DO
+
+ pres_temp_sclr = zpa
+ IF(( PRESENT(pta) ).AND.( PRESENT(ptpot) )) pta = zta
+
+ END FUNCTION pres_temp_sclr
+
+
+ FUNCTION pres_temp_vctr( pqspe, pslp, pz, ptpot, pta, l_ice )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION pres_temp ***
+ !!
+ !! ** Purpose : compute air pressure using barometric equation
+ !! from either potential or absolute air temperature
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp), DIMENSION(jpi,jpj) :: pres_temp_vctr ! air pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pslp ! sea-level pressure [Pa]
+ REAL(wp), INTENT(in ) :: pz ! height above surface [m]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
+ INTEGER :: ji, jj ! loop indices
+ LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
+ LOGICAL :: lice ! sea-ice presence
+
+ lice = .FALSE.
+ IF( PRESENT(l_ice) ) lice = l_ice
+
+ IF( PRESENT(ptpot) ) THEN
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, ptpot=ptpot(ji,jj), pta=pta(ji,jj), l_ice=lice )
+ END_2D
+ ELSE
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, pta=pta(ji,jj), l_ice=lice )
+ END_2D
+ ENDIF
+
+ END FUNCTION pres_temp_vctr
+
+
+ FUNCTION theta_exner_sclr( pta, ppa )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION theta_exner ***
+ !!
+ !! ** Purpose : compute air/surface potential temperature from absolute temperature
+ !! and pressure using Exner function
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp) :: theta_exner_sclr ! air/surface potential temperature [K]
+ REAL(wp), INTENT(in) :: pta ! air/surface absolute temperature [K]
+ REAL(wp), INTENT(in) :: ppa ! air/surface pressure [Pa]
+
+ theta_exner_sclr = pta * ( rpref / ppa ) ** rgamma_dry
+
+ END FUNCTION theta_exner_sclr
+
+ FUNCTION theta_exner_vctr( pta, ppa )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION theta_exner ***
+ !!
+ !! ** Purpose : compute air/surface potential temperature from absolute temperature
+ !! and pressure using Exner function
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp), DIMENSION(jpi,jpj) :: theta_exner_vctr ! air/surface potential temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta ! air/surface absolute temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! air/surface pressure [Pa]
+ INTEGER :: ji, jj ! loop indices
+
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ theta_exner_vctr(ji,jj) = theta_exner_sclr( pta(ji,jj), ppa(ji,jj) )
+ END_2D
+
+ END FUNCTION theta_exner_vctr
FUNCTION rho_air_vctr( ptak, pqa, ppa )
@@ -227,7 +369,9 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! pressure in [Pa]
REAL(wp), DIMENSION(jpi,jpj) :: rho_air_vctr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
+
rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
+
END FUNCTION rho_air_vctr
FUNCTION rho_air_sclr( ptak, pqa, ppa )
@@ -244,9 +388,8 @@ CONTAINS
REAL(wp) :: rho_air_sclr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
- END FUNCTION rho_air_sclr
-
+ END FUNCTION rho_air_sclr
FUNCTION visc_air_sclr(ptak)
@@ -268,12 +411,15 @@ CONTAINS
END FUNCTION visc_air_sclr
FUNCTION visc_air_vctr(ptak)
+
REAL(wp), DIMENSION(jpi,jpj) :: visc_air_vctr ! kinetic viscosity (m^2/s)
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)
INTEGER :: ji, jj ! dummy loop indices
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
+ visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
END_2D
+
END FUNCTION visc_air_vctr
@@ -318,10 +464,12 @@ CONTAINS
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp), DIMENSION(jpi,jpj) :: cp_air_vctr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
+
cp_air_vctr = rCp_dry + rCp_vap * pqa
+
END FUNCTION cp_air_vctr
FUNCTION cp_air_sclr( pqa )
@@ -335,12 +483,12 @@ CONTAINS
REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp) :: cp_air_sclr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
+
cp_air_sclr = rCp_dry + rCp_vap * pqa
- END FUNCTION cp_air_sclr
+ END FUNCTION cp_air_sclr
- !===============================================================================================
FUNCTION gamma_moist_sclr( ptak, pqa )
!!----------------------------------------------------------------------------------
!! ** Purpose : Compute the moist adiabatic lapse-rate.
@@ -365,17 +513,19 @@ CONTAINS
gamma_moist_sclr = grav * ( 1._wp + zLvap*zwa*ziRT ) / ( rCp_dry + zLvap*zLvap*zwa*reps0*ziRT/zta )
!!
END FUNCTION gamma_moist_sclr
- !!
+
FUNCTION gamma_moist_vctr( ptak, pqa )
+
REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
+ gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
END_2D
+
END FUNCTION gamma_moist_vctr
- !===============================================================================================
FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
@@ -398,17 +548,15 @@ CONTAINS
!!-------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- !
- zqa = (1._wp + rctv0*pqa(ji,jj))
- !
- ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
- ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
- ! or
- ! b/ -u* [ theta* + 0.61 theta q* ]
- !
- One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
- & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
- !
+ zqa = (1._wp + rctv0*pqa(ji,jj))
+ !
+ ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
+ ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
+ ! or
+ ! b/ -u* [ theta* + 0.61 theta q* ]
+ !
+ One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
+ & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
END_2D
!
One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
@@ -416,7 +564,6 @@ CONTAINS
END FUNCTION One_on_L
- !===============================================================================================
FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
!!----------------------------------------------------------------------------------
!! Bulk Richardson number according to "wide-spread equation"...
@@ -427,7 +574,7 @@ CONTAINS
!!----------------------------------------------------------------------------------
REAL(wp) :: Ri_bulk_sclr
REAL(wp), INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
- REAL(wp), INTENT(in) :: psst ! SST [K]
+ REAL(wp), INTENT(in) :: psst ! potential SST [K]
REAL(wp), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K]
REAL(wp), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg]
REAL(wp), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg]
@@ -437,29 +584,20 @@ CONTAINS
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv ! local scalars
+ REAL(wp) :: ztptv
!!-------------------------------------------------------------------
- IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
+ IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
!
- zsstv = virt_temp_sclr( psst, pssq ) ! virtual SST (absolute==potential because z=0!)
+ zsstv = virt_temp_sclr( psst, pssq ) ! virtual potential SST
+ ztptv = virt_temp_sclr( ptha, pqa ) ! virtual potential air temperature
+ zdthv = ztptv - zsstv ! air-sea delta of "virtual potential temperature"
!
- zdthv = virt_temp_sclr( ptha, pqa ) - zsstv ! air-sea delta of "virtual potential temperature"
- !
- !! ztv: estimate of the ABSOLUTE virtual temp. within the layer
- IF( l_ptqa_l_prvd ) THEN
- ztv = virt_temp_sclr( pta_layer, pqa_layer )
- ELSE
- zqa = 0.5_wp*( pqa + pssq ) ! ~ mean q within the layer...
- zta = 0.5_wp*( psst + ptha - gamma_moist(ptha, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
- zta = 0.5_wp*( psst + ptha - gamma_moist( zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
- zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer
- ztv = 0.5_wp*( zsstv + virt_temp_sclr( ptha-zgamma*pz, pqa ) )
- END IF
- !
- Ri_bulk_sclr = grav*zdthv*pz / ( ztv*pub*pub ) ! the usual definition of Ri_bulk_sclr
+ Ri_bulk_sclr = grav * zdthv * pz / ( ztptv * pub * pub ) ! the usual definition of Ri_bulk_sclr
!
END FUNCTION Ri_bulk_sclr
- !!
+
FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
+
REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk_vctr
REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K]
@@ -472,19 +610,22 @@ CONTAINS
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
INTEGER :: ji, jj
- IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+
+ IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
IF( l_ptqa_l_prvd ) THEN
- Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
- & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
+ & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) )
+ END_2D
ELSE
- Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
+ END_2D
END IF
- END_2D
+
END FUNCTION Ri_bulk_vctr
- !===============================================================================================
- !===============================================================================================
+
FUNCTION e_sat_sclr( ptak )
!!----------------------------------------------------------------------------------
!! *** FUNCTION e_sat_sclr ***
@@ -509,7 +650,7 @@ CONTAINS
& + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
!
END FUNCTION e_sat_sclr
- !!
+
FUNCTION e_sat_vctr(ptak)
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_vctr !: vapour pressure at saturation [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: temperature (K)
@@ -518,9 +659,8 @@ CONTAINS
e_sat_vctr(ji,jj) = e_sat_sclr(ptak(ji,jj))
END_2D
END FUNCTION e_sat_vctr
- !===============================================================================================
- !===============================================================================================
+
FUNCTION e_sat_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! Same as "e_sat" but over ice rather than water!
@@ -536,8 +676,9 @@ CONTAINS
zle = rAg_i*(ztmp - 1._wp) + rBg_i*LOG10(ztmp) + rCg_i*(1._wp - zta/rtt0) + rDg_i
!!
e_sat_ice_sclr = 100._wp * 10._wp**zle
+
END FUNCTION e_sat_ice_sclr
- !!
+
FUNCTION e_sat_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
@@ -545,10 +686,12 @@ CONTAINS
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
+ e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
END_2D
+
END FUNCTION e_sat_ice_vctr
- !!
+
+
FUNCTION de_sat_dt_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! d [ e_sat_ice ] / dT (derivative / temperature)
@@ -566,7 +709,7 @@ CONTAINS
!!
de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta)
END FUNCTION de_sat_dt_ice_sclr
- !!
+
FUNCTION de_sat_dt_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: de_sat_dt_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
@@ -574,15 +717,12 @@ CONTAINS
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
+ de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
END_2D
- END FUNCTION de_sat_dt_ice_vctr
-
+ END FUNCTION de_sat_dt_ice_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION q_sat_sclr( pta, ppa, l_ice )
!!---------------------------------------------------------------------------------
!! *** FUNCTION q_sat_sclr ***
@@ -606,9 +746,11 @@ CONTAINS
ze_s = e_sat( pta ) ! Vapour pressure at saturation (Goff) :
END IF
q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s)
+
END FUNCTION q_sat_sclr
- !!
+
FUNCTION q_sat_vctr( pta, ppa, l_ice )
+
REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
@@ -619,13 +761,12 @@ CONTAINS
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
+ q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
END_2D
+
END FUNCTION q_sat_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION dq_sat_dt_ice_sclr( pta, ppa )
!!---------------------------------------------------------------------------------
!! *** FUNCTION dq_sat_dt_ice_sclr ***
@@ -646,21 +787,21 @@ CONTAINS
dq_sat_dt_ice_sclr = reps0*ppa*zde_s_dt / ( ztmp*ztmp )
!
END FUNCTION dq_sat_dt_ice_sclr
- !!
+
FUNCTION dq_sat_dt_ice_vctr( pta, ppa )
+
REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
+ dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
END_2D
+
END FUNCTION dq_sat_dt_ice_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION q_air_rh(prha, ptak, ppa)
!!----------------------------------------------------------------------------------
!! Specific humidity of air out of Relative Humidity
@@ -677,14 +818,14 @@ CONTAINS
!!----------------------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
- q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
+ ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
+ q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
END_2D
!
END FUNCTION q_air_rh
- SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, &
+ SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, prhoa, &
& pQns, pTau, &
& Qlat)
!!----------------------------------------------------------------------------------
@@ -704,6 +845,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prlw ! downwelling longwave radiative flux [W/m^2]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! air density [kg/m3]
!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
@@ -714,50 +856,53 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
- zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
- zz0 = pust(ji,jj)/pUb(ji,jj)
- zCd = zz0*zz0
- zCh = zz0*ptst(ji,jj)/zdt
- zCe = zz0*pqst(ji,jj)/zdq
- CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
- & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), &
- & pTau(ji,jj), zQsen, zQlat )
+ zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
+ zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
+ zz0 = pust(ji,jj)/pUb(ji,jj)
+ zCd = zz0*zz0
+ zCh = zz0*ptst(ji,jj)/zdt
+ zCe = zz0*pqst(ji,jj)/zdq
- zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
+ CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
+ & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
+ & pTau(ji,jj), zQsen, zQlat )
- pQns(ji,jj) = zQlat + zQsen + zQlw
+ zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
+
+ pQns(ji,jj) = zQlat + zQsen + zQlw
+
+ IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
- IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
END_2D
+
END SUBROUTINE UPDATE_QNSOL_TAU
SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
- & pwnd, pUb, ppa, &
+ & pwnd, pUb, ppa, prhoa, &
& pTau, pQsen, pQlat, &
- & pEvap, prhoa, pfact_evap )
+ & pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
- REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
- REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
- REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
- REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
- REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
+ REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
+ REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
+ REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
+ REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
+ REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), INTENT(in) :: pCd
REAL(wp), INTENT(in) :: pCh
REAL(wp), INTENT(in) :: pCe
- REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
- REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
- REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
+ REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+ REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
- REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
@@ -766,16 +911,7 @@ CONTAINS
zfact_evap = 1._wp
IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
- !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
- ztaa = pTa ! first guess...
- DO jq = 1, 4
- zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) !#LB: why not "0.5*(pqs+pqa)" rather then "pqa" ???
- ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder...
- END DO
- zrho = rho_air(ztaa, pqa, ppa)
- zrho = rho_air(ztaa, pqa, ppa-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-
- zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10
+ zUrho = pUb*MAX(prhoa, 1._wp) ! rho*U10
pTau = zUrho * pCd * pwnd ! Wind stress module
@@ -784,34 +920,33 @@ CONTAINS
pQlat = L_vap(pTs) * zevap
IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
- IF( PRESENT(prhoa) ) prhoa = zrho
END SUBROUTINE BULK_FORMULA_SCLR
- !!
+
SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
- & pwnd, pUb, ppa, &
+ & pwnd, pUb, ppa, prhoa, &
& pTau, pQsen, pQlat, &
- & pEvap, prhoa, pfact_evap )
+ & pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
- REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
+ REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCh
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCe
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
@@ -822,15 +957,15 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
- & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &
- & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), &
- & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &
- & pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap )
+ CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+ & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &
+ & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
+ & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &
+ & pEvap=zevap, pfact_evap=zfact_evap )
- IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
- IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
+ IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
END_2D
+
END SUBROUTINE BULK_FORMULA_VCTR
@@ -846,6 +981,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
+
END FUNCTION alpha_sw_vctr
FUNCTION alpha_sw_sclr( psst )
@@ -860,10 +996,10 @@ CONTAINS
REAL(wp), INTENT(in) :: psst ! sea-water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
+
END FUNCTION alpha_sw_sclr
- !===============================================================================================
FUNCTION qlw_net_sclr( pdwlw, pts, l_ice )
!!---------------------------------------------------------------------------------
!! *** FUNCTION qlw_net_sclr ***
@@ -886,9 +1022,11 @@ CONTAINS
END IF
zt2 = pts*pts
qlw_net_sclr = zemiss*( pdwlw - stefan*zt2*zt2) ! zemiss used both as the IR albedo and IR emissivity...
+
END FUNCTION qlw_net_sclr
- !!
+
FUNCTION qlw_net_vctr( pdwlw, pts, l_ice )
+
REAL(wp), DIMENSION(jpi,jpj) :: qlw_net_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts !: surface temperature [K]
@@ -899,12 +1037,14 @@ CONTAINS
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice )
+ qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice )
END_2D
+
END FUNCTION qlw_net_vctr
- !===============================================================================================
+
FUNCTION z0_from_Cd( pzu, pCd, ppsi )
+
REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd !: roughness length [m]
REAL(wp) , INTENT(in) :: pzu !: reference height zu [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd !: (neutral or non-neutral) drag coefficient []
@@ -920,9 +1060,12 @@ CONTAINS
!! Cd provided is the neutral-stability Cd, aka CdN :
z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) ) !LB: ok, double-checked!
END IF
+
END FUNCTION z0_from_Cd
+
FUNCTION Cd_from_z0( pzu, pz0, ppsi )
+
REAL(wp), DIMENSION(jpi,jpj) :: Cd_from_z0 !: (neutral or non-neutral) drag coefficient []
REAL(wp) , INTENT(in) :: pzu !: reference height zu [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 !: roughness length [m]
@@ -939,6 +1082,7 @@ CONTAINS
Cd_from_z0 = 1._wp / LOG( pzu / pz0(:,:) )
END IF
Cd_from_z0 = vkarmn2 * Cd_from_z0 * Cd_from_z0
+
END FUNCTION Cd_from_z0
@@ -964,17 +1108,20 @@ CONTAINS
& + zstab * 1._wp / ( 1._wp + ram_louis * zts ) ! Stable Eq.(A7)
!
END FUNCTION f_m_louis_sclr
- !!
+
FUNCTION f_m_louis_vctr( pzu, pRib, pCdn, pz0 )
+
REAL(wp), DIMENSION(jpi,jpj) :: f_m_louis_vctr
REAL(wp), INTENT(in) :: pzu
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCdn
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) )
+ f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) )
END_2D
+
END FUNCTION f_m_louis_vctr
@@ -1000,19 +1147,23 @@ CONTAINS
& + zstab * 1._wp / ( 1._wp + rah_louis * zts ) ! Stable Eq.(A7) !#LB: in paper it's "ram_louis" and not "rah_louis" typo or what????
!
END FUNCTION f_h_louis_sclr
- !!
+
FUNCTION f_h_louis_vctr( pzu, pRib, pChn, pz0 )
+
REAL(wp), DIMENSION(jpi,jpj) :: f_h_louis_vctr
REAL(wp), INTENT(in) :: pzu
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pChn
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) )
+ f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) )
END_2D
+
END FUNCTION f_h_louis_vctr
+
FUNCTION UN10_from_ustar( pzu, pUzu, pus, ppsi )
!!----------------------------------------------------------------------------------
!! Provides the neutral-stability wind speed at 10 m
@@ -1122,5 +1273,5 @@ CONTAINS
END FUNCTION z0tq_LKB
- !!======================================================================
+
END MODULE sbc_phy
diff --git a/src/OCE/SBC/sbcblk.F90 b/src/OCE/SBC/sbcblk.F90
index 677005d82cf129028624c4d6b12c3a3a39b211f8..a88f532dc52a56a7ac794b9dcb4433e7a34701e8 100644
--- a/src/OCE/SBC/sbcblk.F90
+++ b/src/OCE/SBC/sbcblk.F90
@@ -79,7 +79,7 @@ MODULE sbcblk
INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
- INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
+ INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity (kg/kg)
INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
@@ -170,7 +170,7 @@ MODULE sbcblk
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: sbcblk.F90 14718 2021-04-16 09:43:50Z clem $
+ !! $Id: sbcblk.F90 15551 2021-11-28 20:19:36Z gsamson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -274,8 +274,8 @@ CONTAINS
& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )
IF( ln_ANDREAS ) &
& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' )
- IF( nn_fsbc /= 1 ) &
- & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
+ !IF( nn_fsbc /= 1 ) &
+ ! & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
END IF
IF( ln_skin_wl ) THEN
@@ -502,7 +502,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
+ REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp, zpre, ztheta
REAL(wp) :: ztst
LOGICAL :: llerr
!!----------------------------------------------------------------------
@@ -539,7 +539,7 @@ CONTAINS
! ! compute the surface ocean fluxes using bulk formulea
IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
- ! Specific humidity of air at z=rn_zqt !
+ ! Specific humidity of air at z=rn_zqt
SELECT CASE( nhumi )
CASE( np_humi_sph )
q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity!
@@ -552,16 +552,15 @@ CONTAINS
& sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file
END SELECT
- ! POTENTIAL temperature of air at z=rn_zqt
- ! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
- ! (most reanalysis products provide absolute temp., not potential temp.)
+ ! Potential temperature of air at z=rn_zqt (most reanalysis products provide absolute temp., not potential temp.)
IF( ln_tair_pot ) THEN
! temperature read into file is already potential temperature, do nothing...
theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
ELSE
! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
- theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
+ zpre(:,:) = pres_temp( q_air_zt(:,:), sf(jp_slp)%fnow(:,:,1), rn_zu, pta=sf(jp_tair)%fnow(:,:,1) )
+ theta_air_zt(:,:) = theta_exner( sf(jp_tair)%fnow(:,:,1), zpre(:,:) )
ENDIF
!
CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in
@@ -585,15 +584,12 @@ CONTAINS
ELSE
qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
ENDIF
- tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ????
- !tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ?
-
- qatm_ice(:,:) = q_air_zt(:,:) !#LB:
-
- tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
- sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
- wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
- wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
+ tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature (theta_air_zt) instead ????
+ qatm_ice(:,:) = q_air_zt(:,:)
+ tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
+ sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
+ wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
+ wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
ENDIF
#endif
!
@@ -651,6 +647,8 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean
+ REAL(wp), DIMENSION(jpi,jpj) :: zsspt ! potential sea-surface temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, ztabs ! air pressure [Pa] & absolute temperature [K]
REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2
!!---------------------------------------------------------------------
!
@@ -658,6 +656,9 @@ CONTAINS
! ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K)
ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
+ ! sea surface potential temperature [K]
+ zsspt(:,:) = theta_exner( ptsk(:,:), pslp(:,:) )
+
! --- cloud cover --- !
cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
@@ -704,7 +705,7 @@ CONTAINS
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! Backup "bulk SST" and associated spec. hum.
- zztmp1(:,:) = ptsk(:,:)
+ zztmp1(:,:) = zsspt(:,:)
zztmp2(:,:) = pssq(:,:)
ENDIF
@@ -713,33 +714,33 @@ CONTAINS
SELECT CASE( nblk )
CASE( np_NCAR )
- CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_ncar ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
CASE( np_COARE_3p0 )
- CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_COARE_3p6 )
- CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ECMWF )
- CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ANDREAS )
- CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_andreas ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
@@ -760,34 +761,43 @@ CONTAINS
IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
- !! ptsk and pssq have been updated!!!
- !!
- !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
+ !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
WHERE ( fr_i(:,:) > 0.001_wp )
! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
- ptsk(:,:) = zztmp1(:,:)
- pssq(:,:) = zztmp2(:,:)
+ zsspt(:,:) = zztmp1(:,:)
+ pssq(:,:) = zztmp2(:,:)
END WHERE
+ ! apply potential temperature increment to abolute SST
+ ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
END IF
! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90
! -------------------------------------------------------------
IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zztmp = zU_zu(ji,jj)
wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
psen(ji,jj) = zztmp * zch_oce(ji,jj)
pevp(ji,jj) = zztmp * zce_oce(ji,jj)
- rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
+ zpre(ji,jj) = pres_temp( pqair(ji,jj), pslp(ji,jj), rn_zu, ptpot=ptair(ji,jj), pta=ztabs(ji,jj) )
+ rhoa(ji,jj) = rho_air( ztabs(ji,jj), pqair(ji,jj), zpre(ji,jj) )
END_2D
+
ELSE !== BLK formulation ==! turbulent fluxes computation
- CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
+
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ zpre(ji,jj) = pres_temp( q_zu(ji,jj), pslp(ji,jj), rn_zu, ptpot=theta_zu(ji,jj), pta=ztabs(ji,jj) )
+ rhoa(ji,jj) = rho_air( ztabs(ji,jj), q_zu(ji,jj), zpre(ji,jj) )
+ END_2D
+
+ CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
- & wndm(:,:), zU_zu(:,:), pslp(:,:), &
+ & wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
- & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
+ & pEvap=pevp(:,:), pfact_evap=rn_efac )
psen(:,:) = psen(:,:) * tmask(:,:,1)
plat(:,:) = plat(:,:) * tmask(:,:,1)
@@ -795,19 +805,19 @@ CONTAINS
pevp(:,:) = pevp(:,:) * tmask(:,:,1)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- IF( wndm(ji,jj) > 0._wp ) THEN
- zztmp = taum(ji,jj) / wndm(ji,jj)
+ IF( wndm(ji,jj) > 0._wp ) THEN
+ zztmp = taum(ji,jj) / wndm(ji,jj)
#if defined key_cyclone
- ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
- ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+ ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+ ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
#else
- ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
- ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
+ ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
+ ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
#endif
- ELSE
- ztau_i(ji,jj) = 0._wp
- ztau_j(ji,jj) = 0._wp
- ENDIF
+ ELSE
+ ztau_i(ji,jj) = 0._wp
+ ztau_j(ji,jj) = 0._wp
+ ENDIF
END_2D
IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
@@ -850,13 +860,13 @@ CONTAINS
CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ')
ENDIF
!
- ENDIF !IF( ln_abl )
+ ENDIF ! ln_blk / ln_abl
ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
- CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
+ CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference
ENDIF
!
END SUBROUTINE blk_oce_1
@@ -975,8 +985,8 @@ CONTAINS
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
- REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
- REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s]
+ REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric potential temperature at T-point [K]
+ REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric specific humidity at T-point [kg/kg]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
@@ -990,11 +1000,9 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
- REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array
+ REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
!!---------------------------------------------------------------------
!
- ! LB: ptsui is in K !!!
- !
! ------------------------------------------------------------ !
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
@@ -1003,9 +1011,10 @@ CONTAINS
wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
END_2D
!
- ! Make ice-atm. drag dependent on ice concentration
-
+ ! potential sea-ice surface temperature [K]
+ zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
+ ! sea-ice <-> atmosphere bulk transfer coefficients
SELECT CASE( nblk_ice )
CASE( np_ice_cst )
@@ -1019,17 +1028,17 @@ CONTAINS
CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, &
+ CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lu12 )
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
+ CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
+ CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
END SELECT
@@ -1049,8 +1058,8 @@ CONTAINS
! supress moving ice in wind stress computation as we don't know how to do it properly...
DO_2D( 0, 1, 0, 1 ) ! at T point
zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
- putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
- pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
+ putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
+ pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
END_2D
!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
@@ -1073,15 +1082,15 @@ CONTAINS
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
ELSE ! ln_abl
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
- pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
- pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
+ pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
+ pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
+ pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
END_2D
- !#LB:
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
- !#LB.
- ENDIF !IF( ln_blk )
+
+ ENDIF ! ln_blk / ln_abl
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
!
@@ -1112,7 +1121,7 @@ CONTAINS
REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
!!
INTEGER :: ji, jj, jl ! dummy loop indices
- REAL(wp) :: zst, zst3, zsq ! local variable
+ REAL(wp) :: zst, zst3, zsq, zsipt ! local variable
REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - -
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
@@ -1120,13 +1129,12 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
+ REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
REAL(wp), DIMENSION(jpi,jpj) :: ztri
REAL(wp), DIMENSION(jpi,jpj) :: zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
!!---------------------------------------------------------------------
!
zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars
- !
-
zztmp = 1. / ( 1. - albo )
dqla_ice(:,:,:) = 0._wp
@@ -1140,52 +1148,53 @@ CONTAINS
! ! ========================== !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
- zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
-
- ! ----------------------------!
- ! I Radiative FLUXES !
- ! ----------------------------!
- ! Short Wave (sw)
- qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
-
- ! Long Wave (lw)
- zst3 = zst * zst * zst
- z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
- ! lw sensitivity
- z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
-
- ! ----------------------------!
- ! II Turbulent FLUXES !
- ! ----------------------------!
-
- ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
-
- ! Common term in bulk F. equations...
- zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
-
- ! Sensible Heat
- zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
- z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj))
- z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
-
- ! Latent Heat
- zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
- qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
- IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
- ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
- !#LB: without this unjustified "condensation sensure":
- !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
- !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
-
-
- ! ----------------------------!
- ! III Total FLUXES !
- ! ----------------------------!
- ! Downward Non Solar flux
- qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
- ! Total non solar heat flux sensitivity for ice
- dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
+ zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
+ zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
+ zsipt = theta_exner( zst, pslp(ji,jj) ) ! potential sea-ice surface temperature [K]
+
+ ! ----------------------------!
+ ! I Radiative FLUXES !
+ ! ----------------------------!
+ ! Short Wave (sw)
+ qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
+
+ ! Long Wave (lw)
+ zst3 = zst * zst * zst
+ z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
+ ! lw sensitivity
+ z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
+
+ ! ----------------------------!
+ ! II Turbulent FLUXES !
+ ! ----------------------------!
+
+ ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
+
+ ! Common term in bulk F. equations...
+ zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
+
+ ! Sensible Heat
+ zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
+ z_qsb (ji,jj,jl) = zztmp1 * (zsipt - theta_zu_i(ji,jj))
+ z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
+
+ ! Latent Heat
+ zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
+ qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
+ IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
+ ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
+ !#LB: without this unjustified "condensation sensure":
+ !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
+ !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
+
+
+ ! ----------------------------!
+ ! III Total FLUXES !
+ ! ----------------------------!
+ ! Downward Non Solar flux
+ qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
+ ! Total non solar heat flux sensitivity for ice
+ dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
END_2D
!
diff --git a/src/OCE/SBC/sbcblk_algo_coare3p0.F90 b/src/OCE/SBC/sbcblk_algo_coare3p0.F90
index 1dfd38e8d8a57fddf314d1769c6f6ba20bb46d2e..2ad244fcb4d788a921e25c502cfd966d479db70f 100644
--- a/src/OCE/SBC/sbcblk_algo_coare3p0.F90
+++ b/src/OCE/SBC/sbcblk_algo_coare3p0.F90
@@ -188,6 +188,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
@@ -320,11 +321,16 @@ CONTAINS
q_zu = q_zt - q_star/vkarmn*ztmp1
ENDIF
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
+
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
@@ -335,7 +341,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u
!! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
diff --git a/src/OCE/SBC/sbcblk_algo_coare3p6.F90 b/src/OCE/SBC/sbcblk_algo_coare3p6.F90
index 8255e81c09360dd54e2cac843a89d58ad6126787..cb7fff12b99999e0c740d5d8a23393f8dcfec3e7 100644
--- a/src/OCE/SBC/sbcblk_algo_coare3p6.F90
+++ b/src/OCE/SBC/sbcblk_algo_coare3p6.F90
@@ -178,6 +178,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
@@ -310,11 +311,16 @@ CONTAINS
q_zu = q_zt - q_star/vkarmn*ztmp1
ENDIF
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
+
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
@@ -325,7 +331,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u
!! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
diff --git a/src/OCE/SBC/sbcblk_algo_ecmwf.F90 b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
index 7d1c886312ffe33a5401ce714359469434a0ca89..d86c997565a355c17a2eee34b0f6b1ab72c06b92 100644
--- a/src/OCE/SBC/sbcblk_algo_ecmwf.F90
+++ b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
@@ -183,6 +183,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air
REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length...
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
+ REAL(wp), DIMENSION(jpi,jpj) :: zrhoa, zpre, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
!
@@ -346,12 +347,16 @@ CONTAINS
func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0
CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1
@@ -363,7 +368,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2
CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
!! Updating T_s and q_s !!!
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index d6e020564ec0dd19b742e10ec737caefd4b324a5..05070b410a823177a2029e61964eb95a29268c14 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -53,7 +53,7 @@ MODULE sbccpl
USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut
#endif
- USE sbc_phy, ONLY : pp_cldf
+ USE sbc_phy, ONLY : pp_cldf, rpref
IMPLICIT NONE
PRIVATE
@@ -219,9 +219,6 @@ MODULE sbccpl
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
#endif
- REAL(wp) :: rpref = 101000._wp ! reference atmospheric pressure[N/m2]
- REAL(wp) :: r1_grau ! = 1.e0 / (grav * rho0)
-
INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) :: nrcvinfo ! OASIS info argument
!! Substitution
@@ -229,7 +226,7 @@ MODULE sbccpl
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: sbccpl.F90 15004 2021-06-16 10:33:18Z mathiot $
+ !! $Id: sbccpl.F90 15551 2021-11-28 20:19:36Z gsamson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -1188,6 +1185,7 @@ CONTAINS
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
REAL(wp) :: zzx, zzy ! temporary variables
+ REAL(wp) :: r1_grau ! = 1.e0 / (grav * rho0)
REAL(wp), DIMENSION(jpi,jpj) :: ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
!!----------------------------------------------------------------------
!
diff --git a/src/OCE/ZDF/zdfiwm.F90 b/src/OCE/ZDF/zdfiwm.F90
index a6b5aaeec1135d8edca9a775c4b5e98947f41631..30a2dfd71a5ece2e41e15b323ac6f09b165ee8ec 100644
--- a/src/OCE/ZDF/zdfiwm.F90
+++ b/src/OCE/ZDF/zdfiwm.F90
@@ -8,6 +8,8 @@ MODULE zdfiwm
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 3.6 ! 2016-03 (C. de Lavergne) New param: internal wave-driven mixing
!! 4.0 ! 2017-04 (G. Madec) renamed module, remove the old param. and the CPP keys
+ !! 4.0 ! 2020-12 (C. de Lavergne) Update param to match published one
+ !! 4.0 ! 2021-09 (C. de Lavergne) Add energy from trapped and shallow internal tides
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -36,24 +38,25 @@ MODULE zdfiwm
PUBLIC zdf_iwm_init ! called in nemogcm module
! !!* Namelist namzdf_iwm : internal wave-driven mixing *
- INTEGER :: nn_zpyc ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
LOGICAL :: ln_mevar ! variable (=T) or constant (=F) mixing efficiency
LOGICAL :: ln_tsdiff ! account for differential T/S wave-driven mixing (=T) or not (=F)
REAL(wp):: r1_6 = 1._wp / 6._wp
+ REAL(wp):: rnu = 1.4e-6_wp ! molecular kinematic viscosity
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_iwm ! power available from high-mode wave breaking (W/m2)
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: epyc_iwm ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_iwm ! power available from low-mode, critical slope wave breaking (W/m2)
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_iwm ! WKB decay scale for high-mode energy dissipation (m)
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_iwm ! decay scale for low-mode critical slope dissipation (m)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_iwm ! bottom-intensified dissipation above abyssal hills (W/m2)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_iwm ! bottom-intensified dissipation at topographic slopes (W/m2)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ensq_iwm ! dissipation scaling with squared buoyancy frequency (W/m2)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esho_iwm ! dissipation due to shoaling internal tides (W/m2)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_iwm ! decay scale for abyssal hill dissipation (m)
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_iwm ! inverse decay scale for topographic slope dissipation (m-1)
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: zdfiwm.F90 14882 2021-05-18 16:32:47Z gsamson $
+ !! $Id: zdfiwm.F90 15533 2021-11-24 12:07:20Z cdllod $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -62,8 +65,8 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** FUNCTION zdf_iwm_alloc ***
!!----------------------------------------------------------------------
- ALLOCATE( ebot_iwm(jpi,jpj), epyc_iwm(jpi,jpj), ecri_iwm(jpi,jpj) , &
- & hbot_iwm(jpi,jpj), hcri_iwm(jpi,jpj) , STAT=zdf_iwm_alloc )
+ ALLOCATE( ebot_iwm(jpi,jpj), ecri_iwm(jpi,jpj), ensq_iwm(jpi,jpj) , &
+ & esho_iwm(jpi,jpj), hbot_iwm(jpi,jpj), hcri_iwm(jpi,jpj) , STAT=zdf_iwm_alloc )
!
CALL mpp_sum ( 'zdfiwm', zdf_iwm_alloc )
IF( zdf_iwm_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_iwm_alloc: failed to allocate arrays' )
@@ -78,7 +81,7 @@ CONTAINS
!! breaking internal waves.
!!
!! ** Method : - internal wave-driven vertical mixing is given by:
- !! Kz_wave = min( 100 cm2/s, f( Reb = zemx_iwm /( Nu * N^2 ) )
+ !! Kz_wave = min( f( Reb = zemx_iwm / (Nu * N^2) ), 100 cm2/s )
!! where zemx_iwm is the 3D space distribution of the wave-breaking
!! energy and Nu the molecular kinematic viscosity.
!! The function f(Reb) is linear (constant mixing efficiency)
@@ -86,52 +89,53 @@ CONTAINS
!!
!! - Compute zemx_iwm, the 3D power density that allows to compute
!! Reb and therefrom the wave-induced vertical diffusivity.
- !! This is divided into three components:
- !! 1. Bottom-intensified low-mode dissipation at critical slopes
+ !! This is divided into four components:
+ !! 1. Bottom-intensified dissipation at topographic slopes, expressed
+ !! as an exponential decay above the bottom.
!! zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm )
!! / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm
!! where hcri_iwm is the characteristic length scale of the bottom
- !! intensification, ecri_iwm a map of available power, and H the ocean depth.
- !! 2. Pycnocline-intensified low-mode dissipation
- !! zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc )
- !! / SUM( sqrt(rn2(z))^nn_zpyc * e3w[z) )
- !! where epyc_iwm is a map of available power, and nn_zpyc
- !! is the chosen stratification-dependence of the internal wave
- !! energy dissipation.
- !! 3. WKB-height dependent high mode dissipation
- !! zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
- !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w[z) )
- !! where hbot_iwm is the characteristic length scale of the WKB bottom
- !! intensification, ebot_iwm is a map of available power, and z_wkb is the
- !! WKB-stretched height above bottom defined as
- !! z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w[z'>=z) )
- !! / SUM( sqrt(rn2(z')) * e3w[z') )
+ !! intensification, ecri_iwm a static 2D map of available power, and
+ !! H the ocean depth.
+ !! 2. Bottom-intensified dissipation above abyssal hills, expressed
+ !! as an algebraic decay above bottom.
+ !! zemx_iwm(z) = ( ebot_iwm / rho0 ) * ( 1 + hbot_iwm/H )
+ !! / ( 1 + (H-z)/hbot_iwm )^2
+ !! where hbot_iwm is the characteristic length scale of the bottom
+ !! intensification and ebot_iwm is a static 2D map of available power.
+ !! 3. Dissipation scaling in the vertical with the squared buoyancy
+ !! frequency (N^2).
+ !! zemx_iwm(z) = ( ensq_iwm / rho0 ) * rn2(z)
+ !! / ZSUM( rn2 * e3w )
+ !! where ensq_iwm is a static 2D map of available power.
+ !! 4. Dissipation due to shoaling internal tides, scaling in the
+ !! vertical with the buoyancy frequency (N).
+ !! zemx_iwm(z) = ( esho_iwm / rho0 ) * sqrt(rn2(z))
+ !! / ZSUM( sqrt(rn2) * e3w )
+ !! where esho_iwm is a static 2D map of available power.
!!
- !! - update the model vertical eddy viscosity and diffusivity:
- !! avt = avt + av_wave
+ !! - update the model vertical eddy viscosity and diffusivity:
+ !! avt = avt + av_wave
+ !! avs = avs + av_wave
!! avm = avm + av_wave
!!
!! - if namelist parameter ln_tsdiff = T, account for differential mixing:
- !! avs = avt + av_wave * diffusivity_ratio(Reb)
+ !! avs = avs + av_wave * diffusivity_ratio(Reb)
!!
- !! ** Action : - avt, avs, avm, increased by tide internal wave-driven mixing
+ !! ** Action : - avt, avs, avm, increased by internal wave-driven mixing
!!
- !! References : de Lavergne et al. 2015, JPO; 2016, in prep.
+ !! References : de Lavergne et al. JAMES 2020, https://doi.org/10.1029/2020MS002065
+ !! de Lavergne et al. JPO 2016, https://doi.org/10.1175/JPO-D-14-0259.1
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time step
- INTEGER , INTENT(in ) :: Kmm ! time level index
+ INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points)
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avt, p_avs ! tracer Kz (w-points)
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), SAVE :: zztmp
- REAL(wp) :: ztmp1, ztmp2 ! scalar workspace
+ !
REAL(wp), DIMENSION(A2D(nn_hls)) :: zfact ! Used for vertical structure
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zhdep ! Ocean depth
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwkb ! WKB-stretched height above bottom
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zweight ! Weight for high mode vertical distribution
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_t ! Molecular kinematic viscosity (T grid)
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: znu_w ! Molecular kinematic viscosity (W grid)
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zReb ! Turbulence intensity parameter
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zemx_iwm ! local energy density available for mixing (W/kg)
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T)
@@ -140,170 +144,140 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D - - - -
!!----------------------------------------------------------------------
!
- !
- ! Set to zero the 1st and last vertical levels of appropriate variables
- IF( iom_use("emix_iwm") ) THEN
- zemx_iwm(:,:,:) = 0._wp
- ENDIF
- IF( iom_use("av_ratio") ) THEN
- zav_ratio(:,:,:) = 0._wp
- ENDIF
- IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN
- zav_wave(:,:,:) = 0._wp
- ENDIF
+ ! !* Initialize appropriately certain variables
+ zav_ratio(:,:,:) = 1._wp * wmask(:,:,:) ! important to set it to 1 here
+ IF( iom_use("emix_iwm") ) zemx_iwm (:,:,:) = 0._wp
+ IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) zav_wave (:,:,:) = 0._wp
!
! ! ----------------------------- !
! ! Internal wave-driven mixing ! (compute zav_wave)
! ! ----------------------------- !
!
- ! !* Critical slope mixing: distribute energy over the time-varying ocean depth,
- ! using an exponential decay from the seafloor.
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level
- zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
- zfact(ji,jj) = rho0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) ) )
- IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
- END_2D
-!!gm gde3w ==>>> check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
- IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization
- zemx_iwm(ji,jj,jk) = 0._wp
- ELSE
- zemx_iwm(ji,jj,jk) = zfact(ji,jj) * ( EXP( ( gde3w(ji,jj,jk ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) &
- & - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) ) &
- & / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
+ ! !* 'cri' component: distribute energy over the time-varying
+ ! !* ocean depth using an exponential decay from the seafloor.
+ DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level
+ IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ecri_iwm(ji,jj) * r1_rho0 / ( 1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) ) )
+ ELSE ; zfact(ji,jj) = 0._wp
ENDIF
+ END_2D
+
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
+ zemx_iwm(ji,jj,jk) = zfact(ji,jj) * ( EXP( ( gdept(ji,jj,jk ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) &
+ & - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) &
+ & ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm)
END_3D
-!!gm delta(gde3w) = e3t(:,:,:,Kmm) !! Please verify the grid-point position w versus t-point
-!!gm it seems to me that only 1/hcri_iwm is used ==> compute it one for all
+ !* 'bot' component: distribute energy over the time-varying
+ !* ocean depth using an algebraic decay above the seafloor.
+ DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! part independent of the level
+ IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ebot_iwm(ji,jj) * ( 1._wp + hbot_iwm(ji,jj) / ht(ji,jj) ) * r1_rho0
+ ELSE ; zfact(ji,jj) = 0._wp
+ ENDIF
+ END_2D
- ! !* Pycnocline-intensified mixing: distribute energy over the time-varying
- ! !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
- ! ! (NB: N2 is masked, so no use of wmask here)
- SELECT CASE ( nn_zpyc )
- !
- CASE ( 1 ) ! Dissipation scales as N (recommended)
- !
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zfact(ji,jj) = 0._wp
- END_2D
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level
- zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk)
- END_3D
- !
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
- END_2D
- !
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
- zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk)
- END_3D
- !
- CASE ( 2 ) ! Dissipation scales as N^2
- !
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zfact(ji,jj) = 0._wp
- END_2D
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level
- zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk)
- END_3D
- !
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
- END_2D
- !
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk)
- END_3D
- !
- END SELECT
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
+ zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + &
+ & zfact(ji,jj) * ( 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk ,Kmm) ) / hbot_iwm(ji,jj) ) &
+ & - 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) ) &
+ & ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm)
+ END_3D
- ! !* WKB-height dependent mixing: distribute energy over the time-varying
- ! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
- !
+ !* 'nsq' component: distribute energy over the time-varying
+ !* ocean depth as proportional to rn2
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zwkb(ji,jj,1) = 0._wp
+ zfact(ji,jj) = 0._wp
END_2D
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- zwkb(ji,jj,jk) = zwkb(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk)
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level
+ zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) )
END_3D
- DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zfact(ji,jj) = zwkb(ji,jj,jpkm1)
- END_2D
!
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) &
- & * wmask(ji,jj,jk) / zfact(ji,jj)
- END_3D
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1)
+ IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = ensq_iwm(ji,jj) * r1_rho0 / zfact(ji,jj)
END_2D
!
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN ! optimization: EXP coast a lot
- zweight(ji,jj,jk) = 0._wp
- ELSE
- zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj) &
- & * ( EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) ) )
- ENDIF
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
+ zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) )
END_3D
- !
+
+ !* 'sho' component: distribute energy over the time-varying
+ !* ocean depth as proportional to sqrt(rn2)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zfact(ji,jj) = 0._wp
END_2D
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level
- zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk)
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! part independent of the level
+ zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) )
END_3D
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+ IF( zfact(ji,jj) /= 0._wp ) zfact(ji,jj) = esho_iwm(ji,jj) * r1_rho0 / zfact(ji,jj)
END_2D
!
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
- zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk) &
- & / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
-!!gm use of e3t(ji,jj,:,Kmm) just above?
- END_3D
- !
-!!gm this is to be replaced by just a constant value znu=1.e-6 m2/s
- ! Calculate molecular kinematic viscosity
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
- znu_t(ji,jj,jk) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm) &
- & + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) &
- & + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm) ) * tmask(ji,jj,jk) * r1_rho0
- END_3D
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk)
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! complete with the level-dependent part
+ zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) )
END_3D
-!!gm end
- !
+
! Calculate turbulence intensity parameter Reb
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) )
+ zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, rnu * rn2(ji,jj,jk) )
END_3D
!
! Define internal wave-induced diffusivity
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- zav_wave(ji,jj,jk) = znu_w(ji,jj,jk) * zReb(ji,jj,jk) * r1_6 ! This corresponds to a constant mixing efficiency of 1/6
+ zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * r1_6 * rnu ! This corresponds to a constant mixing efficiency of 1/6
END_3D
!
- IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+ IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224) regimes
IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
- zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+ zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( zReb(ji,jj,jk) )
ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
- zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+ zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
ENDIF
END_3D
ENDIF
!
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Bound diffusivity by molecular value and 100 cm2/s
- zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk)
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Bound diffusivity by molecular value and 100 cm2/s
+ zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk)
+ END_3D
+ !
+ ! ! ----------------------- !
+ ! ! Update mixing coefs !
+ ! ! ----------------------- !
+ !
+ IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature
+ DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb (else it is set to 1)
+ zav_ratio(ji,jj,jk) = ( 0.505_wp + &
+ & 0.495_wp * TANH( 0.92_wp * ( LOG10( MAX( 1.e-20, zReb(ji,jj,jk) * 5._wp * r1_6 ) ) - 0.60_wp ) ) &
+ & ) * wmask(ji,jj,jk)
+ END_3D
+ ENDIF
+ CALL iom_put( "av_ratio", zav_ratio )
+ !
+ DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing
+ p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)
+ p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
+ p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
END_3D
+ ! !* output internal wave-driven mixing coefficient
+ CALL iom_put( "av_wave", zav_wave )
+ !* output useful diagnostics: Kz*N^2 ,
+ ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
+ IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
+ ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
+ z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp ! Initialisation for iom_put
+ DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+ z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
+ z2d(ji,jj) = z2d(ji,jj) + rho0 * e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)
+ END_3D
+ CALL iom_put( "bflx_iwm", z3d )
+ CALL iom_put( "pcmap_iwm", z2d )
+ DEALLOCATE( z2d , z3d )
+ ENDIF
+ CALL iom_put( "emix_iwm", zemx_iwm )
+
!
IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave
IF( .NOT. l_istiled .OR. ntile == 1 ) zztmp = 0._wp ! Do only on the first tile
-!!gm used of glosum 3D....
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj) &
& * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
@@ -311,7 +285,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
CALL mpp_sum( 'zdfiwm', zztmp )
- zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing
+ zztmp = rho0 * zztmp ! Global integral of rho0 * Kz * N^2 = power contributing to mixing
!
IF(lwp) THEN
WRITE(numout,*)
@@ -323,58 +297,6 @@ CONTAINS
ENDIF
ENDIF
- ! ! ----------------------- !
- ! ! Update mixing coefs !
- ! ! ----------------------- !
- !
- IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature
- ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
- DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Calculate S/T diffusivity ratio as a function of Reb
- ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
- IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
- zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) )
- ELSE
- zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk)
- ENDIF
- END_3D
- CALL iom_put( "av_ratio", zav_ratio )
- DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing
- p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)
- p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
- p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
- END_3D
- !
- ELSE !* update momentum & tracer diffusivity with wave-driven mixing
- DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk)
- p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
- p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
- END_3D
- ENDIF
-
- ! !* output internal wave-driven mixing coefficient
- CALL iom_put( "av_wave", zav_wave )
- !* output useful diagnostics: Kz*N^2 ,
-!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5)
- ! vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
- IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
- ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
- ! Initialisation for iom_put
- z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp
-
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
- z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)
- END_3D
- DO_2D( 0, 0, 0, 0 )
- z2d(ji,jj) = rho0 * z2d(ji,jj)
- END_2D
- CALL iom_put( "bflx_iwm", z3d )
- CALL iom_put( "pcmap_iwm", z2d )
- DEALLOCATE( z2d , z3d )
- ENDIF
- CALL iom_put( "emix_iwm", zemx_iwm )
-
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
!
END SUBROUTINE zdf_iwm
@@ -384,48 +306,50 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_iwm_init ***
!!
- !! ** Purpose : Initialization of the wave-driven vertical mixing, reading
- !! of input power maps and decay length scales in netcdf files.
+ !! ** Purpose : Initialization of the internal wave-driven vertical mixing, reading
+ !! of input power maps and decay length scales in a netcdf file.
!!
!! ** Method : - Read the namzdf_iwm namelist and check the parameters
!!
- !! - Read the input data in NetCDF files :
- !! power available from high-mode wave breaking (mixing_power_bot.nc)
- !! power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
- !! power available from critical slope wave-breaking (mixing_power_cri.nc)
- !! WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
- !! decay scale for critical slope wave-breaking (decay_scale_cri.nc)
+ !! - Read the input data in a NetCDF file (zdfiwm_forcing.nc) with variables:
+ !! 'power_bot' bottom-intensified dissipation above abyssal hills
+ !! 'power_cri' bottom-intensified dissipation at topographic slopes
+ !! 'power_nsq' dissipation scaling with squared buoyancy frequency
+ !! 'power_sho' dissipation due to shoaling internal tides
+ !! 'scale_bot' decay scale for abyssal hill dissipation
+ !! 'scale_cri' decay scale for topographic-slope dissipation
!!
!! ** input : - Namlist namzdf_iwm
- !! - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
- !! decay_scale_bot.nc decay_scale_cri.nc
+ !! - NetCDF file : zdfiwm_forcing.nc
!!
!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
- !! - Define ebot_iwm, epyc_iwm, ecri_iwm, hbot_iwm, hcri_iwm
+ !! - Define ebot_iwm, ecri_iwm, ensq_iwm, esho_iwm, hbot_iwm, hcri_iwm
!!
- !! References : de Lavergne et al. JPO, 2015 ; de Lavergne PhD 2016
- !! de Lavergne et al. in prep., 2017
+ !! References : de Lavergne et al. JAMES 2020, https://doi.org/10.1029/2020MS002065
!!----------------------------------------------------------------------
INTEGER :: ifpr ! dummy loop indices
INTEGER :: inum ! local integer
INTEGER :: ios
- REAL(wp) :: zbot, zpyc, zcri ! local scalars
!
CHARACTER(len=256) :: cn_dir ! Root directory for location of ssr files
- INTEGER, PARAMETER :: jpiwm = 5 ! maximum number of files to read
+ INTEGER, PARAMETER :: jpiwm = 6 ! maximum number of variables to read
INTEGER, PARAMETER :: jp_mpb = 1
- INTEGER, PARAMETER :: jp_mpp = 2
- INTEGER, PARAMETER :: jp_mpc = 3
- INTEGER, PARAMETER :: jp_dsb = 4
- INTEGER, PARAMETER :: jp_dsc = 5
+ INTEGER, PARAMETER :: jp_mpc = 2
+ INTEGER, PARAMETER :: jp_mpn = 3
+ INTEGER, PARAMETER :: jp_mps = 4
+ INTEGER, PARAMETER :: jp_dsb = 5
+ INTEGER, PARAMETER :: jp_dsc = 6
+ !
+ TYPE(FLD_N), DIMENSION(jpiwm) :: slf_iwm ! array of namelist informations
+ TYPE(FLD_N) :: sn_mpb, sn_mpc, sn_mpn, sn_mps ! information about Mixing Power field to be read
+ TYPE(FLD_N) :: sn_dsb, sn_dsc ! information about Decay Scale field to be read
+ TYPE(FLD ), DIMENSION(jpiwm) :: sf_iwm ! structure of input fields (file informations, fields read)
!
- TYPE(FLD_N), DIMENSION(jpiwm) :: slf_iwm ! array of namelist informations
- TYPE(FLD_N) :: sn_mpb, sn_mpp, sn_mpc ! informations about Mixing Power field to be read
- TYPE(FLD_N) :: sn_dsb, sn_dsc ! informations about Decay Scale field to be read
- TYPE(FLD ), DIMENSION(jpiwm) :: sf_iwm ! structure of input fields (file informations, fields read)
+ REAL(wp), DIMENSION(jpi,jpj,4) :: ztmp
+ REAL(wp), DIMENSION(4) :: zdia
!
- NAMELIST/namzdf_iwm/ nn_zpyc, ln_mevar, ln_tsdiff, &
- & cn_dir, sn_mpb, sn_mpp, sn_mpc, sn_dsb, sn_dsc
+ NAMELIST/namzdf_iwm/ ln_mevar, ln_tsdiff, &
+ & cn_dir, sn_mpb, sn_mpc, sn_mpn, sn_mps, sn_dsb, sn_dsc
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namzdf_iwm, IOSTAT = ios, ERR = 901)
@@ -440,17 +364,16 @@ CONTAINS
WRITE(numout,*) 'zdf_iwm_init : internal wave-driven mixing'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namzdf_iwm : set wave-driven mixing parameters'
- WRITE(numout,*) ' Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
WRITE(numout,*) ' Variable (T) or constant (F) mixing efficiency = ', ln_mevar
WRITE(numout,*) ' Differential internal wave-driven mixing (T) or not (F) = ', ln_tsdiff
ENDIF
- ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
+ ! This internal-wave-driven mixing parameterization elevates avt and avm in the interior, and
! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should
! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
- avmb(:) = 1.4e-6_wp ! viscous molecular value
+ avmb(:) = rnu ! molecular value
avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_iwm)
- avtb_2d(:,:) = 1.e0_wp ! uniform
+ avtb_2d(:,:) = 1._wp ! uniform
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) ' Force the background value applied to avm & avt in TKE to be everywhere ', &
@@ -461,41 +384,50 @@ CONTAINS
IF( zdf_iwm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' )
!
! store namelist information in an array
- slf_iwm(jp_mpb) = sn_mpb ; slf_iwm(jp_mpp) = sn_mpp ; slf_iwm(jp_mpc) = sn_mpc
+ slf_iwm(jp_mpb) = sn_mpb ; slf_iwm(jp_mpc) = sn_mpc ; slf_iwm(jp_mpn) = sn_mpn ; slf_iwm(jp_mps) = sn_mps
slf_iwm(jp_dsb) = sn_dsb ; slf_iwm(jp_dsc) = sn_dsc
!
DO ifpr= 1, jpiwm
ALLOCATE( sf_iwm(ifpr)%fnow(jpi,jpj,1) )
- IF( slf_iwm(ifpr)%ln_tint )ALLOCATE( sf_iwm(ifpr)%fdta(jpi,jpj,1,2) )
+ IF( slf_iwm(ifpr)%ln_tint ) ALLOCATE( sf_iwm(ifpr)%fdta(jpi,jpj,1,2) )
END DO
! fill sf_iwm with sf_iwm and control print
CALL fld_fill( sf_iwm, slf_iwm , cn_dir, 'zdfiwm_init', 'iwm input file', 'namiwm' )
- ! ! hard-coded default definition (to be defined in namelist ?)
- sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-6
- sf_iwm(jp_mpp)%fnow(:,:,1) = 1.e-6
- sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10
- sf_iwm(jp_dsb)%fnow(:,:,1) = 100.
- sf_iwm(jp_dsc)%fnow(:,:,1) = 100.
+ ! ! hard-coded default values
+ sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-10_wp
+ sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10_wp
+ sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-5_wp
+ sf_iwm(jp_mps)%fnow(:,:,1) = 1.e-10_wp
+ sf_iwm(jp_dsb)%fnow(:,:,1) = 100._wp
+ sf_iwm(jp_dsc)%fnow(:,:,1) = 100._wp
! ! read necessary fields
CALL fld_read( nit000, 1, sf_iwm )
- ebot_iwm(:,:) = sf_iwm(1)%fnow(:,:,1) * ssmask(:,:) ! energy flux for high-mode wave breaking [W/m2]
- epyc_iwm(:,:) = sf_iwm(2)%fnow(:,:,1) * ssmask(:,:) ! energy flux for pynocline-intensified wave breaking [W/m2]
- ecri_iwm(:,:) = sf_iwm(3)%fnow(:,:,1) * ssmask(:,:) ! energy flux for critical slope wave breaking [W/m2]
- hbot_iwm(:,:) = sf_iwm(4)%fnow(:,:,1) ! spatially variable decay scale for high-mode wave breaking [m]
- hcri_iwm(:,:) = sf_iwm(5)%fnow(:,:,1) ! spatially variable decay scale for critical slope wave breaking [m]
+ ebot_iwm(:,:) = sf_iwm(1)%fnow(:,:,1) * ssmask(:,:) ! energy flux for dissipation above abyssal hills [W/m2]
+ ecri_iwm(:,:) = sf_iwm(2)%fnow(:,:,1) * ssmask(:,:) ! energy flux for dissipation at topographic slopes [W/m2]
+ ensq_iwm(:,:) = sf_iwm(3)%fnow(:,:,1) * ssmask(:,:) ! energy flux for dissipation scaling with N^2 [W/m2]
+ esho_iwm(:,:) = sf_iwm(4)%fnow(:,:,1) * ssmask(:,:) ! energy flux for dissipation due to shoaling [W/m2]
+ hbot_iwm(:,:) = sf_iwm(5)%fnow(:,:,1) ! spatially variable decay scale for abyssal hill dissipation [m]
+ hcri_iwm(:,:) = sf_iwm(6)%fnow(:,:,1) ! spatially variable decay scale for topographic-slope [m]
+
+ hcri_iwm(:,:) = 1._wp / hcri_iwm(:,:) ! only the inverse height is used, hence calculated here once for all
+
+ ! diags
+ ztmp(:,:,1) = e1e2t(:,:) * ebot_iwm(:,:)
+ ztmp(:,:,2) = e1e2t(:,:) * ecri_iwm(:,:)
+ ztmp(:,:,3) = e1e2t(:,:) * ensq_iwm(:,:)
+ ztmp(:,:,4) = e1e2t(:,:) * esho_iwm(:,:)
- zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) )
- zpyc = glob_sum( 'zdfiwm', e1e2t(:,:) * epyc_iwm(:,:) )
- zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) )
+ zdia(1:4) = glob_sum_vec( 'zdfiwm', ztmp(:,:,1:4) )
IF(lwp) THEN
- WRITE(numout,*) ' High-mode wave-breaking energy: ', zbot * 1.e-12_wp, 'TW'
- WRITE(numout,*) ' Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
- WRITE(numout,*) ' Critical slope wave-breaking energy: ', zcri * 1.e-12_wp, 'TW'
+ WRITE(numout,*) ' Dissipation above abyssal hills: ', zdia(1) * 1.e-12_wp, 'TW'
+ WRITE(numout,*) ' Dissipation along topographic slopes: ', zdia(2) * 1.e-12_wp, 'TW'
+ WRITE(numout,*) ' Dissipation scaling with N^2: ', zdia(3) * 1.e-12_wp, 'TW'
+ WRITE(numout,*) ' Dissipation due to shoaling: ', zdia(4) * 1.e-12_wp, 'TW'
ENDIF
!
END SUBROUTINE zdf_iwm_init
diff --git a/tests/ICE_AGRIF/EXPREF/1_namelist_cfg b/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
index ffbec2ea8d93641b7948c63504218313540ba16c..d1028790e36deb9e89b4a3d8779b815cfddbe823 100644
--- a/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
+++ b/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
@@ -104,6 +104,7 @@
!-----------------------------------------------------------------------
&namagrif ! AGRIF zoom ("key_agrif")
!-----------------------------------------------------------------------
+ ln_init_chfrpar = .true. ! initialize child grids from parent
/
!!======================================================================
!! *** Top/Bottom boundary condition *** !!