Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 441 additions and 373 deletions
sette/param.default
View file @ 4a9f4b18
......@@ -65,8 +65,9 @@ LC_MESSAGES=en_US
#
# LIST OF CONFIGURATIONS
# to be updated if you added an new configuration in sette_test-cases.sh or sette_reference-configuration.sh
TEST_CONFIG_AVAILABLE=(ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 AGRIF_DEMO WED025 GYRE_PISCES ORCA2_SAS_ICE ORCA2_ICE_OBS SWG ICE_AGRIF OVERFLOW LOCK_EXCHANGE VORTEX ISOMIP+)
if [ -z "${TEST_CONFIGS}" ]; then
export TEST_CONFIGS=(${SETTE_TEST_CONFIGS[@]:-"ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 AGRIF WED025 GYRE_PISCES SAS ORCA2_ICE_OBS SWG ICE_AGRIF OVERFLOW LOCK_EXCHANGE VORTEX ISOMIP+"})
export TEST_CONFIGS=(${SETTE_TEST_CONFIGS[@]:-${TEST_CONFIG_AVAILABLE[@]}})
fi
#
# TYPES OF TESTS TO PERFORM
......
sette/sette.sh
View file @ 4a9f4b18
......@@ -12,17 +12,17 @@ NO_REPORT=0
#
export USING_TIMING='yes' # Default: yes => set ln_timing=.true. ; use -T to disable
export USING_ICEBERGS='yes' # Default: yes => set ln_icebergs=.true. ; use -i to disable
export USING_EXTRA_HALO='yes' # Default: yes => set nn_hls=2 ; use -e to set nn_hls=1
export USING_EXTRA_HALO='no' # Default: no => set nn_hls=1 ; use "--halo 2" to set nn_hls=2
export USING_COLLECTIVES='yes' # Default: yes => set nn_comm=2 ; use -C to set nn_comm=1
export USING_NOGATHER='yes' # Default: yes => set ln_nnogather=.true.; use -N to set ln_nnogather=.false.
export USING_TILING='yes' # Default: yes => set ln_tile=.true. ; use -t to disable
export USING_TILING='no' # Default: no => set ln_tile=.false. ; use "--tiling yes" to set ln_tile=.true.
# Note: yes also ensures nn_hls=2 but -t will not alter nn_hls
#
# controls for some common compile-time keys:
#
export USING_QCO='yes' # Default: yes => add key_qco ; use -q to delete key_qco
export USING_RK3='no' # Default: yes => add key_RK3 & key_qco ; use -Q to delete key_RK3
export USING_LOOP_FUSION='yes' # Default: yes => add key_loop_fusion ; use -F to delete key_loop_fusion
export USING_LOOP_FUSION='no' # Default: no => del key_loop_fusion ; use "--loop_fusion yes" to add key_loop_fusion
export USING_XIOS='yes' # Default: yes => add key_xios ; use -X to delete key_xios
# Note: changing USING_XIOS may require a change in arch file
#
......@@ -62,22 +62,40 @@ else
export SETTE_THIS_BRANCH="Unknown"
fi
# small functions
getyn() {
case $( echo $1 | tr '[:upper:]' '[:lower:]' ) in
0|1) [ $1 -eq 0 ] && echo 'no' || echo 'yes' ;;
'no'|'yes') echo $1 ;;
*) echo 'err' ;;
esac
}
prtopt() {
[ "$3" == "err" ] && echo -e "\033[0;31m\n \"$2\" BAD OPTION for ${1}\033[0m\n" && exit 2
echo "$1 $2: $4"
echo ""
}
yn2tf (){
[ "$1" == "no" ] && echo false || echo true
}
yn2not (){
[ "$1" == "no" ] && echo "not " || echo ""
}
# Parse command-line arguments
if [ $# -gt 0 ]; then
while getopts n:x:v:g:cdrshTqQteiACFNXu option; do
case $option in
c) export SETTE_CLEAN_CONFIGS='yes'
while [ ${#1} -gt 0 ]; do
case $1 in
-c) export SETTE_CLEAN_CONFIGS='yes'
export SETTE_SYNC_CONFIGS='yes'
echo "-c: Configuration ${SETTE_TEST_CONFIGS[@]} will be cleaned; this option enforces also synchronisation"
echo "";;
d) dry_run=1
-d) dry_run=1
echo "";;
r) NO_REPORT=1
-r) NO_REPORT=1
echo "";;
s) export SETTE_SYNC_CONFIGS='yes'
-s) export SETTE_SYNC_CONFIGS='yes'
echo "-s: MY_SRC and EXP00 in ${SETTE_TEST_CONFIGS[@]} will be synchronised with the MY_SRC and EXPREF from the reference configuration"
echo "";;
n) OPTSTR="$OPTARG"
-n) OPTSTR=$2
OPTSTR="${OPTSTR/ORCA2_SAS_ICE/SAS}" # Permit either shortened (expected) or full name for SAS
OPTSTR="${OPTSTR/AGRIF_DEMO/AGRIF}" # Permit either shortened (expected) or full name for AGRIF
export SETTE_TEST_CONFIGS=(${OPTSTR})
......@@ -87,66 +105,89 @@ if [ $# -gt 0 ]; then
else
echo "-n: Configuration ${SETTE_TEST_CONFIGS[@]} will be tested if it is available"
fi
echo "";;
g) case $OPTARG in
[0-9,a-z,A-Z] ) echo "-g: Using ${SETTE_STG}${OPTARG} as the configuration suffix";;
echo ""
shift ;;
-g) case $2 in
[0-9,a-z,A-Z] ) echo "-g: Using ${SETTE_STG}${2} as the configuration suffix";;
* ) echo "-g only accepts a single, alphanumeric character. Processing halted"; exit 42;;
esac
export SETTE_STG=${SETTE_STG}${OPTARG}
echo "";;
x) export SETTE_TEST_TYPES=(${OPTARG})
export SETTE_STG=${SETTE_STG}${2}
echo ""
shift ;;
-x) export SETTE_TEST_TYPES=(${2}) # SETTE_TEST_TYPES is an array
echo "-x: ${SETTE_TEST_TYPES[@]} tests requested"
echo "";;
v) export SETTE_SUB_VAL=($OPTARG)
echo ""
shift ;;
-v) export SETTE_SUB_VAL=(${2}) # SETTE_SUB_VAL is an array
echo "-v: $SETTE_SUB_VAL validation sub-directory requested"
echo "";;
T) export USING_TIMING='no'
echo ""
shift ;;
-T) export USING_TIMING='no'
echo "-T: ln_timing will be set to false"
echo "";;
t) export USING_TILING='no'
--tiling)
export USING_TILING=$( getyn $2 )
prtopt $1 $2 $USING_TILING "ln_tile will be set to "$( yn2tf $USING_TILING )
shift ;;
-t) export USING_TILING='no'
echo "-t: ln_tile will be set to false"
echo "";;
e) export USING_EXTRA_HALO='no'
--halo)
case $2 in
1|2) USING_HALO_SIZE=$2 ;;
* ) USING_HALO_SIZE='err' ;;
esac
prtopt "$1" "$2" $USING_HALO_SIZE "nn_hls will be set to $2"
[ "${USING_HALO_SIZE}" == "2" ] && export USING_EXTRA_HALO='yes'
shift ;;
-e) export USING_EXTRA_HALO='no'
echo "-e: nn_hls will be set to 1"
echo "";;
i) export USING_ICEBERGS='no'
-i) export USING_ICEBERGS='no'
echo "-i: ln_icebergs will be set to false"
echo "";;
C) export USING_COLLECTIVES='no'
-C) export USING_COLLECTIVES='no'
echo "-C: nn_comm will be set to 1"
echo "";;
N) export USING_NOGATHER='no'
-N) export USING_NOGATHER='no'
echo "-N: ln_nnogather will be set to false"
echo "";;
q) export USING_QCO='no'
-q) export USING_QCO='no'
echo "-q: key_qco and key_linssh will NOT be activated"
echo "";;
Q) export USING_RK3='no'
-Q) export USING_RK3='no'
echo "-Q: key_qco and key_RK3 will not be activated"
echo " This is the curent default for now since RK3 is not ready"
echo "";;
F) export USING_LOOP_FUSION='no'
--loop_fusion)
export USING_LOOP_FUSION=$( getyn $2 )
prtopt $1 $2 $USING_LOOP_FUSION "key_loop_fusion will $( yn2not $USING_LOOP_FUSION )be activated"
shift ;;
-F) export USING_LOOP_FUSION='no'
echo "-F: key_loop_fusion will not be activated"
echo "";;
X) export USING_XIOS='no'
-X) export USING_XIOS='no'
echo "-X: key_xios will not be activated"
echo "";;
A) export USING_MPMD='no'
-A) export USING_MPMD='no'
echo "-A: Tasks will be run in attached (SPMD) mode"
echo "";;
u) export USER_INPUT='no'
-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)'
-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)'
echo '-e set nn_hls=1 (default: nn_hls=2)'
echo '-t set ln_tile false in all tests that support it (default: false)'
echo '"--tiling $opt" with opt=yes|1|no|0, set ln_tile true/false in all tests that support it (default: false)'
echo '-e set nn_hls=1 (default: nn_hls=1)'
echo '"--halo $opt", with opt=1|2, set nn_hls=1 or 2 (default: nn_hls=1)'
echo '-i set ln_icebergs false (default: true)'
echo '-C set nn_comm=1 (default: nn_comm=2 ==> use MPI3 collective comms)'
echo '-N set ln_nnogather false for ORCA2 configurations (default: true)'
echo '-q to remove the key_qco key (default: added)'
echo '-X to remove the key_xios key (default: added)'
echo '-F to remove the key_loop_fusion key (default: added)'
echo '-F to remove the key_loop_fusion key (default: removed)'
echo '"--loop_fusion $opt" with opt=yes|1|no|0, to add/remove the key_loop_fusion key (default: removed)'
echo '-Q to remove the key_RK3 key (currently a null-op since key_RK3 is not used)'
echo '-A to run tests in attached (SPMD) mode (default: MPMD with key_xios)'
echo '-n "CFG1_to_test CFG2_to_test ..." to test some specific configurations'
......@@ -161,10 +202,9 @@ if [ $# -gt 0 ]; then
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))
fi
esac
shift
done
#
# Option dependency tests
#
......@@ -172,16 +212,16 @@ if [ ${USING_TILING} == "yes" ] ; then
if [ ${USING_EXTRA_HALO} == "no" ] ; then
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
read -p 'Tiling requires the extra halo but you did not used "--halo 2" to select it. Would you like to select it? (y/n)?: ' yn
case $yn in
[Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
[Yy]* ) echo 'Ok, selecting the "--halo 2" 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."
echo 'Tiling requires the extra halo but you did not used "--halo 2" to select it. Tiling will not be used.'
USING_TILING="no"
fi
fi
......@@ -190,16 +230,16 @@ if [ ${USING_LOOP_FUSION} == "yes" ] ; then
if [ ${USING_EXTRA_HALO} == "no" ] ; then
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
read -p 'Loop fusion requires the extra halo but you did not used "--halo 2" to select it. Would you like to select it? (y/n)?: ' yn
case $yn in
[Yy]* ) echo "Ok, ignoring the -e option"; USING_EXTRA_HALO="yes"; break;;
[Yy]* ) echo 'Ok, selecting the "--halo 2" 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."
echo 'Loop fusion requires the extra halo but you did not used "--halo 2" to select it. Loop fusion will not be used.'
USING_LOOP_FUSION="no"
fi
fi
......@@ -253,9 +293,11 @@ if [ ! -d $NEMO_VALIDATION_DIR/$SETTE_SUB_VAL ] && [ ${dry_run} -eq 0 ] ; then
fi
export NEMO_VALIDATION_DIR=$NEMO_VALIDATION_DIR/$SETTE_SUB_VAL
TEST_CONFIGS=(${TEST_CONFIGS[@]/ORCA2_SAS_ICE/SAS}) # Shortening of 'ORCA2_SAS_ICE' to 'SAS'
TEST_CONFIGS=(${TEST_CONFIGS[@]/AGRIF_DEMO/AGRIF}) # Shortening of 'AGRIF_DEMO' to 'AGRIF'
if [ ${#SETTE_TEST_CONFIGS[@]} -eq 0 ]; then
echo "=================================="
echo "Configurations $TEST_CONFIGS will be tested if they are available"
echo "Configurations ${TEST_CONFIGS[@]} will be tested if they are available"
fi
echo "Carrying out the following tests : ${TEST_TYPES[@]}"
echo "requested by the command : "$cmd $cmdargs
......
sette/sette_eval.sh
View file @ 4a9f4b18
......@@ -155,6 +155,9 @@ function runcmpres(){
MAIN_DIR=$(dirname $SETTE_DIR)
quiet=0
. ./param.cfg
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]:-${TEST_CONFIGS[@]}} # workaround for some dated param.cfgs files
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]/ SAS / ORCA2_SAS_ICE } # workaround for some dated param.cfgs files
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]/ AGRIF / AGRIF_DEMO } # workaround for some dated param.cfgs files
USER_INPUT='yes' # Default: yes => request user input on decisions.
mach=${COMPILER}
......@@ -327,8 +330,7 @@ fi
echo "REFERENCE directory : $NEMO_VALID_REF at rev $NEMO_REV_REF"
echo ''
fi
checklist=(GYRE_PISCES ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 ORCA2_SAS_ICE ORCA2_ICE_OBS AGRIF_DEMO WED025 ISOMIP+ VORTEX ICE_AGRIF OVERFLOW LOCK_EXCHANGE SWG)
for repro_test in ${checklist[@]}
for repro_test in ${TEST_CONFIGS_AVAILABLE[@]}
do
runcmpres $NEMO_VALID $repro_test $NEMO_VALID_REF $NEMO_REV_REF $quiet
done
......
sette/sette_reference-configurations.sh
View file @ 4a9f4b18
......@@ -983,7 +983,7 @@ if [ ${config} == "ORCA2_ICE_OBS" ] ; then
#
. ./makenemo -m ${CMP_NAM} -n ${SETTE_CONFIG} -r ORCA2_ICE_PISCES -d "OCE ICE" -j ${CMPL_CORES} add_key "key_asminc ${ADD_KEYS}" del_key "key_top ${DEL_KEYS}"
fi
if [ ${config} == "ORCA2_ICE_OBS" ] && [ ${DO_RESTART} == "1" ] ; then
if [ ${config} == "ORCA2_ICE_OBS" ] && [ ${DO_REPRO} == "1" ] ; then
## Reproducibility tests
export TEST_NAME="REPRO_4_8"
cd ${SETTE_DIR}
......@@ -1422,8 +1422,7 @@ if [ ${config} == "AGRIF" ] && [ ${DO_CORRUPT} == "1" ] ; then
sync_config AGRIF_DEMO ${SETTE_CONFIG} 'cfgs'
clean_config AGRIF_DEMO ${SETTE_CONFIG} 'cfgs'
#
# AGRIF_DEMO does not yet support nn_hls=2 => key_loop_fusion can not be used
. ./makenemo -m ${CMP_NAM} -n ${SETTE_CONFIG} -r AGRIF_DEMO -j ${CMPL_CORES} add_key "${ADD_KEYS/key_loop_fusion}" del_key "key_agrif ${DEL_KEYS}"
. ./makenemo -m ${CMP_NAM} -n ${SETTE_CONFIG} -r AGRIF_DEMO -j ${CMPL_CORES} add_key "${ADD_KEYS}" del_key "key_agrif ${DEL_KEYS}"
cd ${SETTE_DIR}
. ./prepare_exe_dir.sh
set_valid_dir
......
sette/sette_rpt.sh
View file @ 4a9f4b18
......@@ -267,10 +267,12 @@ function runcmpres(){
#
if [ -d $vdir/$mach/$dorv/$nam ]; then
f1s=$vdir/$mach/$dorv/$nam/LONG/run.stat
f1t=$vdir/$mach/$dorv/$nam/LONG/tracer.stat
f2s=$vdirref/$mach/$dorvref/$nam/LONG/run.stat
f2t=$vdirref/$mach/$dorvref/$nam/LONG/tracer.stat
# Selection of the test run used for the comparison (LONG or one of the reproducibility-test runs)
TESTD=$(ls -1 ${vdir}/${mach}/${dorv}/${nam}/ | grep -m 1 -e '^LONG$' -e '^REPRO_'); TESTD=${TESTD:-LONG}
f1s=$vdir/$mach/$dorv/${nam}/${TESTD}/run.stat
f1t=$vdir/$mach/$dorv/${nam}/${TESTD}/tracer.stat
f2s=$vdirref/$mach/$dorvref/${nam}/${TESTD}/run.stat
f2t=$vdirref/$mach/$dorvref/${nam}/${TESTD}/tracer.stat
if [ ! -f $f1s ] && [ ! -f $f1t ] ; then
printf "%-20s %s\n" $nam " incomplete test";
return;
......@@ -286,11 +288,11 @@ function runcmpres(){
cmp -s $f1s $f2s
if [ $? == 0 ]; then
if [ $pass == 0 ]; then
printf "%-20s %s %s\n" $nam " run.stat files are identical "
printf "%-20s %s (%s)\n" $nam " run.stat files are identical " ${TESTD}
fi
else
get_ktdiff $f1s $f2s
printf "%-20s %s %s %-5s %s\n" $nam " run.stat files are DIFFERENT (results are different after " $ktdiff " time steps)"
printf "%-20s %s %s %-5s (%s)\n" $nam " run.stat files are DIFFERENT (results are different after " $ktdiff " time steps) " ${TESTD}
#
# Offer view of differences on the second pass
#
......@@ -310,11 +312,11 @@ function runcmpres(){
cmp -s $f1t $f2t
if [ $? == 0 ]; then
if [ $pass == 0 ]; then
printf "%-20s %s %s\n" $nam " tracer.stat files are identical "
printf "%-20s %s (%s)\n" $nam " tracer.stat files are identical " ${TESTD}
fi
else
get_ktdiff2 $f1t $f2t
printf "%-20s %s %s %-5s %s\n" $nam " tracer.stat files are DIFFERENT (results are different after " $ktdiff " time steps) "
printf "%-20s %s %s %-5s (%s)\n" $nam " tracer.stat files are DIFFERENT (results are different after " $ktdiff " time steps) " ${TESTD}
#
# Offer view of differences on the second pass
#
......@@ -351,8 +353,10 @@ function runcmptim(){
#
if [ -d $vdir/$mach/$dorv/$nam ]; then
f1a=$vdir/$mach/$dorv/$nam/LONG/timing.output
f2a=$vdirref/$mach/$dorvref/$nam/LONG/timing.output
# Selection of the test run used for the comparison (LONG or one of the reproducibility-test runs)
TESTD=$(ls -1 ${vdir}/${mach}/${dorv}/${nam}/ | grep -m 1 -e '^LONG$' -e '^REPRO_'); TESTD=${TESTD:-LONG}
f1a=$vdir/$mach/$dorv/${nam}/${TESTD}/timing.output
f2a=$vdirref/$mach/$dorvref/${nam}/${TESTD}/timing.output
#
# Report average CPU time differences (if available)
#
......@@ -471,6 +475,9 @@ function identictest(){
SETTE_DIR=$(cd $(dirname "$0"); pwd)
MAIN_DIR=$(dirname $SETTE_DIR)
. ./param.cfg
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]:-${TEST_CONFIGS[@]}} # Workaround for some dated param.cfgs files
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]/ SAS / ORCA2_SAS_ICE } # Workaround for some dated param.cfgs files
TEST_CONFIGS_AVAILABLE=${TEST_CONFIGS_AVAILABLE[@]/ AGRIF / AGRIF_DEMO } # Workaround for some dated param.cfgs files
if [ -z $USER_INPUT ] ; then USER_INPUT='yes' ; fi # Default: yes => request user input on decisions.
# (but may br inherited/imported from sette.sh)
......@@ -584,7 +591,7 @@ if [[ $? == 0 ]] ; then
branchname="Unknown"
fi
else
revision=`git rev-list --abbrev-commit origin | head -1l`
revision=`git rev-parse --short HEAD 2> /dev/null`
fi
else
branchname="Unknown"
......@@ -634,17 +641,17 @@ do
# Restartability test
echo ""
echo " !----restart----! "
for restart_test in GYRE_PISCES ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 ORCA2_SAS_ICE AGRIF_DEMO WED025 ISOMIP+ OVERFLOW LOCK_EXCHANGE VORTEX ICE_AGRIF SWG
for restart_test in ${TEST_CONFIGS_AVAILABLE[@]}
do
resttest $NEMO_VALID $restart_test $pass
[ "${restart_test}" != "ORCA2_ICE_OBS" ] && resttest $NEMO_VALID $restart_test $pass
done
#
# Reproducibility tests
echo ""
echo " !----repro----! "
for repro_test in GYRE_PISCES ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 ORCA2_SAS_ICE ORCA2_ICE_OBS AGRIF_DEMO WED025 ISOMIP+ VORTEX ICE_AGRIF SWG
for repro_test in ${TEST_CONFIGS_AVAILABLE[@]}
do
reprotest $NEMO_VALID $repro_test $pass
[ "${repro_test}" != "OVERFLOW" -a "${repro_test}" != "LOCK_EXCHANGE" ] && reprotest $NEMO_VALID $repro_test $pass
done
# AGRIF special check to ensure results are unchanged with and without key_agrif
......@@ -669,14 +676,13 @@ do
echo 'and'
echo "REFERENCE directory : $NEMO_VALID_REF at rev $NEMO_REV_REF"
echo ''
checklist=(GYRE_PISCES ORCA2_ICE_PISCES ORCA2_OFF_PISCES AMM12 ORCA2_SAS_ICE AGRIF_DEMO WED025 ISOMIP+ VORTEX ICE_AGRIF OVERFLOW LOCK_EXCHANGE SWG)
for repro_test in ${checklist[@]}
for runcmp_test in ${TEST_CONFIGS_AVAILABLE[@]}
do
runcmpres $NEMO_VALID $repro_test $NEMO_VALID_REF $NEMO_REV_REF $pass
runcmpres $NEMO_VALID $runcmp_test $NEMO_VALID_REF $NEMO_REV_REF $pass
done
echo ''
echo 'Report timing differences between REFERENCE and VALID (if available) :'
for repro_test in ${checklist[@]}
for repro_test in ${TEST_CONFIGS_AVAILABLE[@]}
do
runcmptim $NEMO_VALID $repro_test $NEMO_VALID_REF $NEMO_REV_REF $pass
done
......
src/ICE/ice.F90
View file @ 4a9f4b18
......@@ -241,6 +241,7 @@ MODULE ice
REAL(wp), PUBLIC :: rn_hpnd !: prescribed pond depth (0<rn_hpnd<1)
LOGICAL, PUBLIC :: ln_pnd_lids !: Allow ponds to have frozen lids
LOGICAL , PUBLIC :: ln_pnd_alb !: melt ponds affect albedo
INTEGER , PUBLIC :: nn_pnd_brsal !: brine salinity formulation 0 = Consistent expression with SI3 (linear liquidus) ; 1 = used in GOSI9
! !!** ice-diagnostics namelist (namdia) **
LOGICAL , PUBLIC :: ln_icediachk !: flag for ice diag (T) or not (F)
......@@ -452,6 +453,11 @@ MODULE ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface conduction flux (W/m2)
!
!!----------------------------------------------------------------------
!! * Only for atmospheric coupling
!!----------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
!
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
......@@ -550,6 +556,10 @@ CONTAINS
ii = ii + 1
ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
! * For atmospheric coupling
ii = ii + 1
ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) )
ice_alloc = MAXVAL( ierr(:) )
IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
!
......
src/ICE/icedyn_adv_umx.F90
View file @ 4a9f4b18
......@@ -32,7 +32,7 @@ MODULE icedyn_adv_umx
PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90
!
INTEGER, PARAMETER :: np_advS = 1 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1
INTEGER, PARAMETER :: np_advS = 2 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1
! or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
! or dVS/dt = -div( uV * uS / u ) => np_advS = 3
INTEGER, PARAMETER :: np_limiter = 1 ! limiter: 1 = nonosc
......@@ -1157,7 +1157,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
!!----------------------------------------------------------------------
zbig = 1.e+40_wp
zbig = 1.e+20_wp ! works ok with simple/double precison
! antidiffusive flux : high order minus low order
! --------------------------------------------------
......
src/ICE/icedyn_rdgrft.F90
View file @ 4a9f4b18
......@@ -51,6 +51,7 @@ MODULE icedyn_rdgrft
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrexp ! e-folding ridge thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging
......@@ -59,7 +60,7 @@ MODULE icedyn_rdgrft
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d
!
REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi)
REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipliyer: (hi/hraft)
REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multiplier: (hi/hraft)
!
! ** namelist (namdyn_rdgrft) **
LOGICAL :: ln_str_R75 ! ice strength parameterization: Rothrock 75
......@@ -67,13 +68,13 @@ MODULE icedyn_rdgrft
LOGICAL :: ln_str_CST ! ice strength parameterization: Constant
REAL(wp) :: rn_str ! constant value of ice strength
LOGICAL :: ln_str_smooth ! ice strength spatial smoothing
LOGICAL :: ln_distf_lin ! redistribution of ridged ice: linear (Hibler 1980)
LOGICAL :: ln_distf_exp ! redistribution of ridged ice: exponential
LOGICAL :: ln_distf_lin ! redistribution of ridged ice: linear (Hibler, 1980)
LOGICAL :: ln_distf_exp ! redistribution of ridged ice: exponential (Lipscomb et al., 2007)
REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5)
REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging
LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975))
LOGICAL :: ln_partf_lin ! participation function: linear (Thorndike et al., 1975)
REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging
LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007))
LOGICAL :: ln_partf_exp ! participation function: exponential (Lipscomb et al., 2007)
REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function
LOGICAL :: ln_ridging ! ridging of ice or not
REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice
......@@ -99,9 +100,9 @@ CONTAINS
!!-------------------------------------------------------------------
!! *** ROUTINE ice_dyn_rdgrft_alloc ***
!!-------------------------------------------------------------------
ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , &
& apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
& hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , &
ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij), &
& apartf(jpij,0:jpl), hrmin (jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl), &
& hrmax (jpij,jpl) , hrexp (jpij,jpl), hi_hrdg(jpij,jpl) , araft(jpij,jpl) , &
& ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
......@@ -301,8 +302,9 @@ CONTAINS
!! ** Method : Compute the thickness distribution of the ice and open water
!! participating in ridging and of the resulting ridges.
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net
REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i
REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i
REAL(wp), DIMENSION(:) , INTENT(in), OPTIONAL :: pclosing_net
!!
INTEGER :: ji, jl ! dummy loop indices
REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar
......@@ -353,7 +355,7 @@ CONTAINS
zGsum(1:npti,-1) = 0._wp
zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
DO jl = 1, jpl
zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok (and not jpl)
zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is correct (and not jpl)
END DO
!
IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975)
......@@ -371,7 +373,7 @@ CONTAINS
END DO
END DO
!
ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007)
ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al., 2007)
!
zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
DO jl = -1, jpl
......@@ -420,8 +422,10 @@ CONTAINS
! 2) Transfer function
!-----------------------------------------------------------------
! If assuming ridged ice is uniformly distributed between hrmin and
! hrmax (ln_distf_lin):
!
! Compute max and min ridged ice thickness for each ridging category.
! Assume ridged ice is uniformly distributed between hrmin and hrmax.
!
! This parameterization is a modified version of Hibler (1980).
! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
......@@ -436,9 +440,37 @@ CONTAINS
! These modifications have the effect of reducing the ice strength
! (relative to the Hibler formulation) when very thick ice is ridging.
!
! zaksum = net area removed/ total area removed
! where total area removed = area of ice that ridges
! net area removed = total area removed - area of new ridges
!-----------------------------------------------------------------
! If assuming ridged ice ITD is a negative exponential
! (ln_distf_exp) and following CICE implementation:
!
! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin
!
! where hrmin is the minimum thickness of ridging ice and
! hrexp is the e-folding thickness.
!
! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
! That is, the minimum ridge thickness results from rafting,
! unless the ice is thicker than maxraft.
!
! Also, assume that hrexp = mu_rdg*sqrt(hi).
! The parameter mu_rdg is tuned to give e-folding scales mostly
! in the range 2-4 m as observed by upward-looking sonar.
!
! Values of mu_rdg in the right column give ice strengths
! roughly equal to values of Hstar in the left column
! (within ~10 kN/m for typical ITDs):
!
! Hstar mu_rdg
!
! 25 3.0
! 50 4.0
! 75 5.0
! 100 6.0
!
! zaksum = net area removed/ total area participating
! where total area participating = area of ice that ridges
! net area removed = total area participating - area of new ridges
!-----------------------------------------------------------------
zfac = 1._wp / hi_hrft
zaksum(1:npti) = apartf(1:npti,0)
......@@ -448,9 +480,16 @@ CONTAINS
IF ( apartf(ji,jl) > 0._wp ) THEN
zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
hraft (ji,jl) = zhi(ji,jl) * zfac
hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
!
IF( ln_distf_lin ) THEN
hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
ELSEIF( ln_distf_exp ) THEN
hrexp (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) )
hi_hrdg(ji,jl) = zhi(ji,jl) / ( hrmin(ji,jl) + hrexp(ji,jl) )
!!clem: set a mini for zhi??
ENDIF
!
! Normalization factor : zaksum, ensures mass conservation
zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) &
......@@ -458,45 +497,51 @@ CONTAINS
ELSE
hrmin (ji,jl) = 0._wp
hrmax (ji,jl) = 0._wp
hrexp (ji,jl) = 0._wp
hraft (ji,jl) = 0._wp
hi_hrdg(ji,jl) = 1._wp
!!clem zaksum(ji,jl) = 0._wp
ENDIF
END DO
END DO
!
! 3) closing_gross
!-----------------
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
ELSEWHERE ; closing_gross(1:npti) = 0._wp
END WHERE
IF( PRESENT( pclosing_net ) ) THEN
!
! 3) closing_gross
!-----------------
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
ELSEWHERE ; closing_gross(1:npti) = 0._wp
END WHERE
! correction to closing rate if excessive ice removal
!----------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category would be removed
! Reduce the opening rate in proportion
DO jl = 1, jpl
DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
ENDIF
END DO
END DO
! correction to closing rate if excessive ice removal
!----------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category would be removed
! Reduce the opening rate in proportion
DO jl = 1, jpl
! 4) correction to opening if excessive open water removal
!---------------------------------------------------------
! Reduce the closing rate if more than 100% of the open water would be removed
! Reduce the opening rate in proportion
DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
ENDIF
END DO
END DO
! 4) correction to opening if excessive open water removal
!---------------------------------------------------------
! Reduce the closing rate if more than 100% of the open water would be removed
! Reduce the opening rate in proportion
DO ji = 1, npti
zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
ENDIF
END DO
!
ENDIF
!
END SUBROUTINE rdgrft_prep
......@@ -513,6 +558,7 @@ CONTAINS
!
INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2
REAL(wp) :: expL, expR ! exponentials involving hL, hR
REAL(wp) :: vsw ! vol of water trapped into ridges
REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted
REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1
......@@ -694,18 +740,50 @@ CONTAINS
IF( ll_shift(ji) ) THEN
! Compute the fraction of ridged ice area and volume going to thickness category jl2
IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
IF( ln_distf_lin ) THEN ! Hibler (1980) linear formulation
!
itest_rdg(ji) = 1 ! test for conservation
ELSE
farea = 0._wp
fvol(ji) = 0._wp
ENDIF
IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
!
itest_rdg(ji) = 1 ! test for conservation
ELSE
farea = 0._wp
fvol(ji) = 0._wp
ENDIF
!
ELSEIF( ln_distf_exp ) THEN ! Lipscomb et al. (2007) exponential formulation
!
IF( jl2 < jpl ) THEN
!
IF( hrmin(ji,jl1) <= hi_max(jl2) ) THEN
hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
hR = hi_max(jl2)
expL = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
expR = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
farea = expL - expR
fvol(ji) = ( ( hL + hrexp(ji,jl1) ) * expL &
- ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
ELSE
farea = 0._wp
fvol(ji) = 0._wp
END IF
!
ELSE ! jl2 = jpl
!
hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
expL = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
farea = expL
fvol(ji) = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
!
END IF ! jl2 < jpl
!
itest_rdg(ji) = 1 ! test for conservation => clem: I am not sure about that
!
END IF ! ridge redistribution
! Compute the fraction of rafted ice area and volume going to thickness category jl2
IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN
zswitch(ji) = 1._wp
......@@ -780,18 +858,20 @@ CONTAINS
!!
!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
!!
!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
!! dissipated per unit area removed from the ice pack under compression,
!! and assumed proportional to the change in potential energy caused
!! by ridging. Note that only Hibler's formulation is stable and that
!! ice strength has to be smoothed
!! by ridging. Note that ice strength using Hibler's formulation must be
!! smoothed.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_3 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zworka ! temporary array used here
!!clem
LOGICAL :: ln_str_R75
REAL(wp) :: zhi, zcp, rn_pe_rdg
!!
LOGICAL :: ln_str_R75
REAL(wp) :: zhi, zcp
REAL(wp) :: h2rdg ! mean value of h^2 for new ridge
REAL(wp), PARAMETER :: zmax_strength = 200.e3_wp ! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m
REAL(wp), DIMENSION(jpij) :: zstrength, zaksum ! strength in 1D
!!----------------------------------------------------------------------
! prepare the mask
......@@ -804,7 +884,7 @@ CONTAINS
CASE ( np_strr75 ) !== Rothrock(1975)'s method ==!
! these 2 param should be defined once for all at the 1st time step
! this should be defined once for all at the 1st time step
zcp = 0.5_wp * grav * (rho0-rhoi) * rhoi * r1_rho0 ! proport const for PE
!
strength(:,:) = 0._wp
......@@ -824,7 +904,7 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti) , strength )
CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
!
zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module
DO jl = 1, jpl
......@@ -841,16 +921,31 @@ CONTAINS
DO ji = 1, npti
!
IF( apartf(ji,jl) > 0._wp ) THEN
!
IF( ln_distf_lin ) THEN ! Uniform redistribution of ridged ice
h2rdg = z1_3 * ( hrmax(ji,jl) * hrmax(ji,jl) + & ! (a**3-b**3)/(a-b) = a*a+ab+b*b
& hrmin(ji,jl) * hrmin(ji,jl) + &
& hrmax(ji,jl) * hrmin(ji,jl) )
!
ELSEIF( ln_distf_exp ) THEN ! Exponential redistribution of ridged ice
h2rdg = hrmin(ji,jl) * hrmin(ji,jl) &
& + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) &
& + 2._wp * hrexp(ji,jl) * hrexp(ji,jl)
END IF
!
IF( a_i_2d(ji,jl) > epsi10 ) THEN ; zhi = v_i_2d(ji,jl) / a_i_2d(ji,jl)
ELSE ; zhi = 0._wp
ENDIF
zstrength(ji) = zstrength(ji) - apartf(ji,jl) * zhi * zhi ! PE loss from deforming ice
zstrength(ji) = zstrength(ji) + 2._wp * araft (ji,jl) * zhi * zhi ! PE gain from rafting ice
zstrength(ji) = zstrength(ji) + aridge(ji,jl) * hi_hrdg(ji,jl) * z1_3 * & ! PE gain from ridging ice
& ( hrmax(ji,jl) * hrmax(ji,jl) + & ! (a**3-b**3)/(a-b) = a*a+ab+b*b
& hrmin(ji,jl) * hrmin(ji,jl) + &
& hrmax(ji,jl) * hrmin(ji,jl) )
!!$ zstrength(ji) = zstrength(ji) - apartf(ji,jl) * zhi * zhi ! PE loss from deforming ice
!!$ zstrength(ji) = zstrength(ji) + 2._wp * araft (ji,jl) * zhi * zhi ! PE gain from rafting ice
!!$ zstrength(ji) = zstrength(ji) + aridge(ji,jl) * hi_hrdg(ji,jl) * z1_3 * & ! PE gain from ridging ice
!!$ & ( hrmax(ji,jl) * hrmax(ji,jl) + & ! (a**3-b**3)/(a-b) = a*a+ab+b*b
!!$ & hrmin(ji,jl) * hrmin(ji,jl) + &
!!$ & hrmax(ji,jl) * hrmin(ji,jl) )
zstrength(ji) = zstrength(ji) - apartf(ji,jl) * zhi * zhi ! PE loss
zstrength(ji) = zstrength(ji) + 2._wp * araft(ji,jl) * zhi * zhi ! PE gain (rafting)
zstrength(ji) = zstrength(ji) + aridge(ji,jl) * h2rdg * hi_hrdg(ji,jl) ! PE gain (ridging)
ENDIF
!
END DO
......@@ -858,6 +953,11 @@ CONTAINS
!
zstrength(1:npti) = rn_pe_rdg * zcp * zstrength(1:npti) / zaksum(1:npti)
!
! Enforce a maximum for strength
! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m
WHERE( zstrength(1:npti) > zmax_strength ) ; zstrength(1:npti) = zmax_strength
ENDWHERE
!
CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength )
!
ENDIF
......@@ -896,24 +996,22 @@ CONTAINS
!!-----------------------------------------------------------------------
!! *** ROUTINE ice_dyn_1d2d ***
!!
!! ** Purpose : move arrays from 1d to 2d and the reverse
!! ** Purpose : move arrays between 1d <=> 2d forms and
!! between 2d <=> 3d forms
!!-----------------------------------------------------------------------
INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
INTEGER, INTENT(in) :: kn ! 1= from 2D/3D to 1D/2D
! 2= from 1D/2D to 2D/3D
!
INTEGER :: jl, jk ! dummy loop indices
!!-----------------------------------------------------------------------
!
SELECT CASE( kn )
! !---------------------!
CASE( 1 ) !== from 2D to 1D ==!
! !---------------------!
! !---------------------------!
CASE( 1 ) !== from 2D/3D to 1D/2D ==!
! !---------------------------!
! fields used but not modified
CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
! the following fields are modified in this routine
!!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
!!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
!!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
......@@ -935,9 +1033,9 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
!
! !---------------------!
CASE( 2 ) !== from 1D to 2D ==!
! !---------------------!
! !---------------------------!
CASE( 2 ) !== from 1D/2D to 2D/3D ==!
! !---------------------------!
CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
......@@ -1009,7 +1107,7 @@ CONTAINS
WRITE(numout,*) ' ice strength value rn_str = ', rn_str
WRITE(numout,*) ' spatial smoothing of the strength ln_str_smooth= ', ln_str_smooth
WRITE(numout,*) ' redistribution of ridged ice: linear (Hibler 1980) ln_distf_lin = ', ln_distf_lin
WRITE(numout,*) ' redistribution of ridged ice: exponential ln_distf_exp = ', ln_distf_exp
WRITE(numout,*) ' redistribution of ridged ice: exponential(Lipscomb 2007) ln_distf_exp = ', ln_distf_exp
WRITE(numout,*) ' e-folding scale of ridged ice rn_murdg = ', rn_murdg
WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg
WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
......@@ -1034,11 +1132,13 @@ CONTAINS
IF( ln_str_CST ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strcst ; ENDIF
IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_rdgrft_init: one and only one ice strength option has to be defined ' )
!
IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN
CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' )
ENDIF
!
IF ( ( ln_distf_lin .AND. ln_distf_exp ) .OR. ( .NOT.ln_distf_lin .AND. .NOT.ln_distf_exp ) ) THEN
CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one redistribution function (ln_distf_lin or ln_distf_exp)' )
ENDIF
!!clem
IF( ln_distf_exp ) CALL ctl_stop( 'ice_dyn_rdgrft_init: exponential redistribution function not coded yet (ln_distf_exp)' )
!
IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
......
src/ICE/icedyn_rhg_vp.F90
View file @ 4a9f4b18
......@@ -1115,7 +1115,7 @@ CONTAINS
END_2D
!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1._wp, zs2, 'T', 1._wp, zs12, 'T', 1._wp )
ENDIF
......@@ -1132,7 +1132,7 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1._wp )
ENDIF
......@@ -1169,8 +1169,8 @@ CONTAINS
END_2D
!
CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1._wp, ztauy_oi, 'V', -1._wp, ztaux_ai, 'U', -1._wp, ztauy_ai, 'V', -1._wp ) !, &
! & ztaux_bi, 'U', -1._wp, ztauy_bi, 'V', -1._wp )
!
CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
......@@ -1200,7 +1200,7 @@ CONTAINS
zsig_II(ji,jj) = 0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )
END_2D
!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1._wp, zsig_II, 'T', 1._wp)
IF( iom_use('normstr') ) CALL iom_put( 'normstr' , zsig_I(:,:) * zmsk00(:,:) ) ! Normal stress
IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
......@@ -1241,7 +1241,7 @@ CONTAINS
END_2D
!
! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1._wp, zsig2_p, 'T', 1._wp)
!
CALL iom_put( 'sig1_pnorm' , zsig1_p )
CALL iom_put( 'sig2_pnorm' , zsig2_p )
......@@ -1265,8 +1265,8 @@ CONTAINS
& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
END_2D
!
CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
& zCorU, 'U', -1., zCorV, 'V', -1. )
CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1._wp, zspgV, 'V', -1._wp, &
& zCorU, 'U', -1._wp, zCorV, 'V', -1._wp )
!
CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
......@@ -1296,7 +1296,7 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1._wp, zfV, 'V', -1._wp )
CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
......@@ -1326,9 +1326,9 @@ CONTAINS
END_2D
CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
& zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
& zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. )
CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1._wp, zdiag_ymtrp_ice, 'V', -1._wp, &
& zdiag_xmtrp_snw, 'U', -1._wp, zdiag_ymtrp_snw, 'V', -1._wp, &
& zdiag_xatrp , 'U', -1._wp, zdiag_yatrp , 'V', -1._wp )
CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
......@@ -1612,14 +1612,14 @@ CONTAINS
END_2D
IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. )
IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1._wp, zv_res , 'V', 1._wp )
DO_2D( 0, 0, 0, 0 ) !clem check bounds
pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) )
END_2D
CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1._wp )
ELSE
......
src/ICE/icerst.F90
View file @ 4a9f4b18
......@@ -26,6 +26,7 @@ MODULE icerst
!
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
USE ioipsl , ONLY : ju2ymds ! for calendar
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
......@@ -51,6 +52,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! number of iteration
!
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(len=20) :: clkt ! ocean time-step define as a character
CHARACTER(len=50) :: clname ! ice output restart file name
CHARACTER(len=256) :: clpath ! full path to ice output restart file
......@@ -67,8 +71,15 @@ CONTAINS
& .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN
IF( nitrst <= nitend .AND. nitrst > 0 ) THEN
! beware of the format used to write kt (default is i8.8, that should be large enough...)
IF( nitrst > 99999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
IF ( ln_rstdate ) THEN
zfjulday = fjulday + (2*nn_fsbc+1)*rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( nitrst > 99999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
ENDIF
ENDIF
! create the file
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
......@@ -313,6 +324,11 @@ CONTAINS
ENDIF
ENDIF
! If this is a coupled model we need to pick up a_i for use as a_i_last_couple
IF (ln_cpl) then
a_i_last_couple = a_i
ENDIF
IF(.NOT.lrxios) CALL iom_delay_rst( 'READ', 'ICE', numrir ) ! read only ice delayed global communication variables
! ! ---------------------------------- !
ELSE ! == case of a simplified restart == !
......
src/ICE/icesbc.F90
View file @ 4a9f4b18
......@@ -322,23 +322,21 @@ CONTAINS
! qlead is the energy received from the atm. in the leads.
! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2)
! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2)
IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN
! upper bound for fhld: fhld should be equal to zqld
! but we have to make sure that this heat will not make the sst drop below the freezing point
! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos
! The following formulation is ok for both normal conditions and supercooling
fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
& - qsb_ice_bot(ji,jj) )
qlead(ji,jj) = 0._wp
ELSE
IF( ( zqld - zqfr ) < 0._wp .OR. at_i(ji,jj) < epsi10 ) THEN
fhld (ji,jj) = 0._wp
! upper bound for qlead: qlead should be equal to zqld
! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point.
! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb)
! and freezing point is reached if zqfr = zqld - qsb*a/dt
! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr
! The energy for this cooling down is zqfr and freezing point is reached if zqfr = zqld
! so the max heat that can be pulled out of the ocean is zqld - zqfr
! The following formulation is ok for both normal conditions and supercooling
qlead(ji,jj) = MIN( 0._wp , zqld - zqfr )
ELSE
! upper bound for fhld: fhld should be equal to zqld
! but we have to make sure that this heat will not make the sst drop below the freezing point
! so the max heat that can be pulled out of the ocean is zqld - zqfr_pos
! The following formulation is ok for both normal conditions and supercooling
qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr )
fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
qlead(ji,jj) = 0._wp
ENDIF
!
! If ice is landfast and ice concentration reaches its max
......
src/ICE/icestp.F90
View file @ 4a9f4b18
......@@ -275,7 +275,7 @@ CONTAINS
CALL ice_istate( nit000, Kbb, Kmm, Kaa ) ! start from rest or read a file
ENDIF
CALL ice_var_glo2eqv
CALL ice_var_agg(1)
CALL ice_var_agg(2)
!
CALL ice_dyn_init ! set ice dynamics parameters
!
......
src/ICE/icetab.F90
View file @ 4a9f4b18
......@@ -32,14 +32,15 @@ MODULE icetab
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d )
SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab2d, tab3d )
!!----------------------------------------------------------------------
!! *** ROUTINE tab_2d_1d ***
!! *** ROUTINE tab_3d_2d ***
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab1d ! output 1D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab3d ! input 3D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab2d ! output 2D field
!
INTEGER :: jl, jn, jid, jjd
!!----------------------------------------------------------------------
......@@ -47,7 +48,7 @@ CONTAINS
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
tab1d(jn,jl) = tab2d(jid,jjd,jl)
tab2d(jn,jl) = tab3d(jid,jjd,jl)
END DO
END DO
END SUBROUTINE tab_3d_2d
......@@ -72,14 +73,14 @@ CONTAINS
END SUBROUTINE tab_2d_1d
SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d )
SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab2d, tab3d )
!!----------------------------------------------------------------------
!! *** ROUTINE tab_2d_1d ***
!! *** ROUTINE tab_2d_3d ***
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab2d ! output 2D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab3d ! output 3D field
!
INTEGER :: jl, jn, jid, jjd
!!----------------------------------------------------------------------
......@@ -87,7 +88,7 @@ CONTAINS
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
tab2d(jid,jjd,jl) = tab1d(jn,jl)
tab3d(jid,jjd,jl) = tab2d(jn,jl)
END DO
END DO
END SUBROUTINE tab_2d_3d
......
src/ICE/icethd.F90
View file @ 4a9f4b18
......@@ -342,11 +342,22 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
!
! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
! Here we make sure that we don't divide by very small, but physically
! meaningless, products of sea ice thicknesses/snow depths and sea ice
! concentration
DO jk = 1, nlay_i
WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
ELSEWHERE
e_i_1d(1:npti,jk) = 0._wp
ENDWHERE
END DO
DO jk = 1, nlay_s
WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
ELSEWHERE
e_s_1d(1:npti,jk) = 0._wp
ENDWHERE
END DO
!
! !---------------------!
......
src/ICE/icethd_dh.F90
View file @ 4a9f4b18
......@@ -192,8 +192,8 @@ CONTAINS
! Snow melting
! ------------
! If heat still available (zq_top > 0)
! then all snw precip has been melted and we need to melt more snow
! Melt snow layers, starting with newly fallen snow layer 0
! and moving downward, until zq_top=0
DO jk = 0, nlay_s
DO ji = 1, npti
IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
......@@ -216,10 +216,10 @@ CONTAINS
END DO
END DO
! Snow sublimation
!-----------------
! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
! Snow sublimation and deposition
!--------------------------------
! when evap_ice_1d > 0 (upwards) snow sublimates and snow thickness decreases
! when evap_ice_1d < 0 (downwards) deposition occurs and snow thickness increases
zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0
zevap_rema(1:npti) = 0._wp
DO ji = 1, npti
......@@ -486,29 +486,6 @@ CONTAINS
END DO
END DO
! Snow load on ice
! -----------------
! When snow load exceeds Archimede's limit and sst is positive,
! snow-ice formation (next bloc) can lead to negative ice enthalpy.
! Therefore we consider here that this excess of snow falls into the ocean
zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos
DO jk = 0, nlay_s
DO ji = 1, npti
IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN
! snow layer thickness that falls into the ocean
zdum = MIN( zdeltah(ji) , zh_s(ji,jk) )
! mass & energy loss to the ocean
hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! mass flux
! update thickness and energy
h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum )
zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum )
! update snow thickness that still has to fall
zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum )
ENDIF
END DO
END DO
! Snow-Ice formation
! ------------------
! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,
......
src/ICE/icethd_pnd.F90
View file @ 4a9f4b18
......@@ -140,10 +140,10 @@ CONTAINS
!------------------------------------
! Diagnostics
!------------------------------------
CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting
CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid
CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage
CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow
IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting
IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid
IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage
IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow
!
DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf )
DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d )
......@@ -544,7 +544,7 @@ CONTAINS
! a_ip -> apond
! a_ip_frac -> apnd
CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' )
!CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' )
!---------------------------------------------------------------
! Initialise
......@@ -644,12 +644,6 @@ CONTAINS
!--------------------------
! Pond lid growth and melt
!--------------------------
! Mean surface temperature
zTavg = 0._wp
DO jl = 1, jpl
zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
END DO
zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
DO jl = 1, jpl-1
......@@ -692,8 +686,8 @@ CONTAINS
! differential growth of base of surface floating ice layer
zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0
zomega = rcnd_i * zdTice / zrhoi_L
zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
- v_il(ji,jj,jl) / a_i(ji,jj,jl)
zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_ip(ji,jj,jl) )**2 ) &
- v_il(ji,jj,jl) / a_ip(ji,jj,jl)
zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
IF ( zdvice > epsi10 ) THEN
......@@ -1319,7 +1313,9 @@ CONTAINS
!-----------------------------------------------------------------
! brine salinity and liquid fraction
!-----------------------------------------------------------------
SELECT CASE( nn_pnd_brsal )
CASE( 0 )
DO k = 1, nlay_i
Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
......@@ -1328,6 +1324,16 @@ CONTAINS
phi(k) = salin(k) / Sbr
END DO
CASE( 1 )
DO k = 1, nlay_i
Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3
phi(k) = salin(k) / Sbr
END DO
END SELECT
!-----------------------------------------------------------------
! permeability
......@@ -1354,7 +1360,7 @@ CONTAINS
NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, &
& ln_pnd_CST , rn_apnd, rn_hpnd, &
& ln_pnd_TOPO, &
& ln_pnd_lids, ln_pnd_alb
& ln_pnd_lids, ln_pnd_alb, nn_pnd_brsal
!!-------------------------------------------------------------------
!
READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
......@@ -1379,6 +1385,7 @@ CONTAINS
WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids
WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb
WRITE(numout,*) ' Brine salinity formulation nn_pnd_brsal = ', nn_pnd_brsal
ENDIF
!
! !== set the choice of ice pond scheme ==!
......
src/ICE/icethd_zdf_bl99.F90
View file @ 4a9f4b18
......@@ -939,7 +939,7 @@ CONTAINS
ELSE
cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl
ENDIF
t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
t1_ice_1d(ji) = isnow_comb(ji) * t_s_1d(ji,1) + ( 1._wp - isnow_comb(ji) ) * t_i_1d(ji,1)
END DO
!
IF( k_cnd == np_cnd_EMU ) THEN
......
src/ICE/iceupdate.F90
View file @ 4a9f4b18
......@@ -31,7 +31,7 @@ MODULE iceupdate
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
USE timing ! Timing
USE oce
IMPLICIT NONE
PRIVATE
......@@ -383,8 +383,8 @@ CONTAINS
zat_v = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji ,jj+1 ) * tmask(ji ,jj+1,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji ,jj+1,1) )
! ! linearized quadratic drag formulation
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - zflagi * pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - zflagi * pv_oce(ji,jj) )
! ! stresses at the ocean surface
utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
......@@ -446,12 +446,15 @@ CONTAINS
CALL iom_get( numrir, jpdom_auto, 'snwice_mass_b', snwice_mass_b )
ELSE ! start from rest
IF(lwp) WRITE(numout,*) ' ==>> previous run without snow-ice mass output then set it'
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) &
& + rhow * (vt_ip(:,:) + vt_il(:,:)) )
snwice_mass_b(:,:) = snwice_mass(:,:)
ENDIF
ELSE !* Start from rest
!JC: I think this is useless with what is now done in ice_istate
IF(lwp) WRITE(numout,*) ' ==>> start from rest: set the snow-ice mass'
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) &
& + rhow * (vt_ip(:,:) + vt_il(:,:)) )
snwice_mass_b(:,:) = snwice_mass(:,:)
ENDIF
!
......
src/ICE/icewri.F90
View file @ 4a9f4b18
......@@ -103,7 +103,7 @@ CONTAINS
IF( iom_use('icethic' ) ) CALL iom_put( 'icethic', hm_i * zmsk00 ) ! ice thickness
IF( iom_use('snwthic' ) ) CALL iom_put( 'snwthic', hm_s * zmsk00 ) ! snw thickness
IF( iom_use('icebrv' ) ) CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 ) ! brine volume
IF( iom_use('iceage' ) ) CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) ) ! ice age
IF( iom_use('iceage' ) ) CALL iom_put( 'iceage' , om_i / rday * zmsk15 ) ! ice age
IF( iom_use('icehnew' ) ) CALL iom_put( 'icehnew', ht_i_new ) ! new ice thickness formed in the leads
IF( iom_use('snwvolu' ) ) CALL iom_put( 'snwvolu', vt_s * zmsksn ) ! snow volume
IF( iom_use('icefrb' ) ) THEN ! Ice freeboard
......@@ -118,14 +118,15 @@ CONTAINS
IF( iom_use('icehlid' ) ) CALL iom_put( 'icehlid', hm_il * zmsk00 ) ! melt pond lid depth
IF( iom_use('icevlid' ) ) CALL iom_put( 'icevlid', vt_il * zmsk00 ) ! melt pond lid total volume per unit area
! salt
IF( iom_use('icesalt' ) ) CALL iom_put( 'icesalt', sm_i * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity
IF( iom_use('icesalt' ) ) CALL iom_put( 'icesalt', sm_i * zmsk00 ) ! mean ice salinity
IF( iom_use('icesalm' ) ) CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 ) ! Mass of salt in sea ice per cell area
IF( iom_use('iceepnd' ) ) CALL iom_put( 'iceepnd', SUM( a_ip_eff * a_i, dim=3 ) * zmsk00 ) ! melt pond total effective fraction per cell area
! heat
IF( iom_use('icetemp' ) ) CALL iom_put( 'icetemp', ( tm_i - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! ice mean temperature
IF( iom_use('snwtemp' ) ) CALL iom_put( 'snwtemp', ( tm_s - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) ) ! snw mean temperature
IF( iom_use('icettop' ) ) CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the ice surface
IF( iom_use('icetbot' ) ) CALL iom_put( 'icetbot', ( t_bo - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the ice bottom
IF( iom_use('icetsni' ) ) CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the snow-ice interface
IF( iom_use('icetemp' ) ) CALL iom_put( 'icetemp', ( tm_i - rt0 ) * zmsk00 ) ! ice mean temperature
IF( iom_use('snwtemp' ) ) CALL iom_put( 'snwtemp', ( tm_s - rt0 ) * zmsksn ) ! snw mean temperature
IF( iom_use('icettop' ) ) CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 ) ! temperature at the ice surface
IF( iom_use('icetbot' ) ) CALL iom_put( 'icetbot', ( t_bo - rt0 ) * zmsk00 ) ! temperature at the ice bottom
IF( iom_use('icetsni' ) ) CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 ) ! temperature at the snow-ice interface
IF( iom_use('icehc' ) ) CALL iom_put( 'icehc' , -et_i * zmsk00 ) ! ice heat content
IF( iom_use('snwhc' ) ) CALL iom_put( 'snwhc' , -et_s * zmsksn ) ! snow heat content
! momentum
......@@ -169,16 +170,16 @@ CONTAINS
! --- category-dependent fields --- !
IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0%
IF( iom_use('iceconc_cat' ) ) CALL iom_put( 'iceconc_cat' , a_i * zmsk00l ) ! area for categories
IF( iom_use('icethic_cat' ) ) CALL iom_put( 'icethic_cat' , h_i * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories
IF( iom_use('snwthic_cat' ) ) CALL iom_put( 'snwthic_cat' , h_s * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories
IF( iom_use('icesalt_cat' ) ) CALL iom_put( 'icesalt_cat' , s_i * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories
IF( iom_use('iceage_cat' ) ) CALL iom_put( 'iceage_cat' , o_i / rday * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age
IF( iom_use('icethic_cat' ) ) CALL iom_put( 'icethic_cat' , h_i * zmsk00l ) ! thickness for categories
IF( iom_use('snwthic_cat' ) ) CALL iom_put( 'snwthic_cat' , h_s * zmsksnl ) ! snow depth for categories
IF( iom_use('icesalt_cat' ) ) CALL iom_put( 'icesalt_cat' , s_i * zmsk00l ) ! salinity for categories
IF( iom_use('iceage_cat' ) ) CALL iom_put( 'iceage_cat' , o_i / rday * zmsk00l ) ! ice age
IF( iom_use('icetemp_cat' ) ) CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) &
& * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature
& * zmsk00l ) ! ice temperature
IF( iom_use('snwtemp_cat' ) ) CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) &
& * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature
IF( iom_use('icettop_cat' ) ) CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature
IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
& * zmsksnl ) ! snow temperature
IF( iom_use('icettop_cat' ) ) CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l ) ! surface temperature
IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l ) ! brine volume
IF( iom_use('iceapnd_cat' ) ) CALL iom_put( 'iceapnd_cat' , a_ip * zmsk00l ) ! melt pond frac for categories
IF( iom_use('icevpnd_cat' ) ) CALL iom_put( 'icevpnd_cat' , v_ip * zmsk00l ) ! melt pond volume for categories
IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories
......
src/NST/README.rst deleted 100644 → 0
View file @ 266206a5
**************
Embedded zooms
**************
.. todo::
.. contents::
:local:
Overview
========
AGRIF (Adaptive Grid Refinement In Fortran) is a library that
allows the seamless space and time refinement over rectangular regions in NEMO.
Refinement factors can be odd or even (usually lower than 5 to maintain stability).
Interaction between grid is "two-ways" in the sense that
the parent grid feeds the child grid open boundaries and
the child grid provides volume averages of prognostic variables once
a given number of time step is completed.
These pages provide guidelines how to use AGRIF in NEMO.
For a more technical description of the library itself, please refer to AGRIF_.
Compilation
===========
Activating AGRIF requires to append the cpp key ``key_agrif`` at compilation time:
.. code-block:: sh
./makenemo [...] add_key 'key_agrif'
Although this is transparent to users,
the way the code is processed during compilation is different from the standard case:
a preprocessing stage (the so called ``conv`` program) translates the actual code so that
saved arrays may be switched in memory space from one domain to an other.
Definition of the grid hierarchy
================================
An additional text file :file:`AGRIF_FixedGrids.in` is required at run time.
This is where the grid hierarchy is defined.
An example of such a file, here taken from the ``ICEDYN`` test case, is given below
.. literalinclude:: ../../../tests/ICE_AGRIF/EXPREF/AGRIF_FixedGrids.in
The first line indicates the number of zooms (1).
The second line contains the starting and ending indices in both directions on the root grid
(``imin=34 imax=63 jmin=34 jmax=63``) followed by the space and time refinement factors (3 3 3).
The last line is the number of child grid nested in the refined region (0).
A more complex example with telescoping grids can be found below and
in the :file:`AGRIF_DEMO` reference configuration directory.
.. todo::
Add some plots here with grid staggering and positioning?
When creating the nested domain, one must keep in mind that
the child domain is shifted toward north-east and
depends on the number of ghost cells as illustrated by
the *attempted* drawing below for ``nbghostcells=1`` and ``nbghostcells=3``.
The grid refinement is 3 and ``nxfin`` is the number of child grid points in i-direction.
.. image:: _static/agrif_grid_position.jpg
Note that rectangular regions must be defined so that they are connected to a single parent grid.
Hence, defining overlapping grids with the same refinement ratio will not work properly,
boundary data exchange and update being only performed between root and child grids.
Use of east-west periodic or north-fold boundary conditions is not allowed in child grids either.
Defining for instance a circumpolar zoom in a global model is therefore not possible.
Preprocessing
=============
Knowing the refinement factors and area,
a ``NESTING`` pre-processing tool may help to create needed input files
(mesh file, restart, climatological and forcing files).
The key is to ensure volume matching near the child grid interface,
a step done by invoking the :file:`Agrif_create_bathy.exe` program.
You may use the namelists provided in the :file:`NESTING` directory as a guide.
These correspond to the namelists used to create ``AGRIF_DEMO`` inputs.
Namelist options
================
Each child grid expects to read its own namelist so that different numerical choices can be made
(these should be stored in the form :file:`1_namelist_cfg`, :file:`2_namelist_cfg`, etc...
according to their rank in the grid hierarchy).
Consistent time steps and number of steps with the chosen time refinement have to be provided.
Specific to AGRIF is the following block:
.. literalinclude:: ../../namelists/namagrif
:language: fortran
where sponge layer coefficients have to be chosen according to the child grid mesh size.
The sponge area is hard coded in NEMO and applies on the following grid points:
2 x refinement factor (from ``i=1+nbghostcells+1`` to ``i=1+nbghostcells+sponge_area``)
.. rubric:: References
.. bibliography:: zooms.bib
:all:
:style: unsrt
:labelprefix: A
:keyprefix: a-