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
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 589 additions and 534 deletions
mk/bldxag.cfg
View file @ 66fb9523
......@@ -14,6 +14,7 @@ search_src 1
src::ioipsl $MAIN_DIR/ext/IOIPSL/src
src::nemo $CONFIG_DIR/$NEW_CONF/WORK
src::ppr_1d $MAIN_DIR/ext/PPR/src
bld::target nemo.exe
bld::exe_dep
......@@ -35,8 +36,10 @@ bld::tool::make %MK
# Pre-process code before analysing dependencies
bld::pp::ioipsl 1
bld::pp::nemo 1
bld::pp::ppr_1d 1
bld::tool::fppflags::nemo %FPPFLAGS
bld::tool::fppflags::ioipsl %FPPFLAGS
bld::tool::fppflags::ppr_1d %FPPFLAGS
# Ignore the following dependencies
bld::excl_dep inc::netcdf.inc
......
mk/bldxagxcdf.cfg
View file @ 66fb9523
......@@ -15,6 +15,7 @@ search_src 1
src::nocdf $MAIN_DIR/ext/DUMMY_NETCDF
src::ioipsl $MAIN_DIR/ext/IOIPSL/src
src::nemo $CONFIG_DIR/$NEW_CONF/WORK
src::ppr_1d $MAIN_DIR/ext/PPR/src
bld::target nemo.exe
bld::exe_dep
......@@ -37,9 +38,11 @@ bld::tool::make %MK
bld::pp::nocdf 1
bld::pp::ioipsl 1
bld::pp::nemo 1
bld::pp::ppr_1d 1
bld::tool::fppflags::nemo %FPPFLAGS
bld::tool::fppflags::nocdf %FPPFLAGS
bld::tool::fppflags::ioipsl %FPPFLAGS
bld::tool::fppflags::ppr_1d %FPPFLAGS
# Ignore the following dependencies
bld::excl_dep inc::VT.inc
......
sette/BATCH_TEMPLATE/batch-X64_AA_INTEL_OMPI
View file @ 66fb9523
......@@ -98,7 +98,7 @@ EOF
#
# Run SPMD case
#
time ./nemo
$MPIRUN --ntasks=TOTAL_NPROCS ./nemo
fi
#
......
sette/BATCH_TEMPLATE/batch-X64_JEANZAY
View file @ 66fb9523
......@@ -2,6 +2,7 @@
#SBATCH -A GROUP_IDRIS@cpu
#SBATCH --job-name=SETTE_JOB # nom du job
#SBATCH --partition=cpu_p1 # Nom de la partition d'exécution
#SBATCH --qos=qos_cpu-dev # quality of service
#SBATCH --ntasks=NPROCS # Nombre total de processus MPI
#SBATCH --ntasks-per-node=40 # Nombre de processus MPI par noeud
# /!\ Attention, la ligne suivante est trompeuse mais dans le vocabulaire
......
sette/all_functions.sh
View file @ 66fb9523
......@@ -172,7 +172,8 @@ clean_config() {
# define validation dir
set_valid_dir () {
if [ ${DETACHED_HEAD} == "no" ] ; then
REVISION_NB=`git -C ${MAIN_DIR} rev-list --abbrev-commit HEAD | head -1l`
branchname=$( git -C ${MAIN_DIR} branch --show-current )
REVISION_NB=$( git -C ${MAIN_DIR} rev-list --abbrev-commit origin/$branchname | head -1l )
else
REVISION_NB=${DETACHED_CMIT}
fi
......
sette/sette.sh
View file @ 66fb9523
......@@ -40,10 +40,10 @@ export USER_INPUT='yes' # Default: yes => request user input on decisions
#
# Check that git branch is usable
export DETACHED_HEAD="no"
git branch --show-current >& /dev/null
git -C ${MAIN_DIR} 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_SUB_VAL="$(git -C ${MAIN_DIR} branch --show-current)"
if [ -z $SETTE_SUB_VAL ] ; then
# Probabably on a detached HEAD (possibly testing an old commit).
# Verify this and try to recover original commit
......@@ -76,6 +76,7 @@ if [ $# -gt 0 ]; then
y) dry_run=1
echo "";;
b) NEMO_DEBUG="-b"
echo "-b: Nemo will be compiled with DEBUG options"
echo "";;
r) NO_REPORT=1
echo "";;
......@@ -166,7 +167,7 @@ if [ $# -gt 0 ]; then
echo ' for SETTE-built configurations (needed if sette.sh invocations may overlap)'
echo '-r to execute without waiting to run sette_rpt.sh at the end (useful for chaining sette.sh invocations)'
echo '-y to perform a dryrun to simply report what settings will be used'
echo '-d to compile Nemo with debug options (only if %DEBUG_FCFLAGS if defined in your arch file)'
echo '-b to compile Nemo with debug options (only if %DEBUG_FCFLAGS if defined in your arch file)'
echo '-c to clean each configuration'
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'
......
sette/sette_eval.sh
View file @ 66fb9523
......@@ -272,7 +272,7 @@ if [[ $? == 0 ]] ; then
branchname="Unknown"
fi
else
revision=`git rev-list --abbrev-commit origin | head -1l`
revision=$( git rev-list --abbrev-commit origin/$branchname | head -1l )
fi
else
branchname="Unknown"
......
sette/sette_rpt.sh
View file @ 66fb9523
......@@ -586,7 +586,7 @@ if [[ $? == 0 ]] ; then
branchname="Unknown"
fi
else
revision=`git rev-list --abbrev-commit HEAD | head -1l`
revision=$( git rev-list --abbrev-commit origin/$branchname | head -1l )
fi
else
branchname="Unknown"
......
sette/sette_test-cases.sh
View file @ 66fb9523
......@@ -108,7 +108,7 @@ if [ ${USING_MPMD} == "no" ]
fi
#
# Directory to run the tests
CONFIG_DIR0=${MAIN_DIR}/cfgs
CONFIG_DIR0=${MAIN_DIR}/tests
TOOLS_DIR=${MAIN_DIR}/tools
if [ -n "${CUSTOM_DIR}" ]; then
......
src/ABL/abl.F90
View file @ 66fb9523
......@@ -41,6 +41,8 @@ MODULE abl
!
INTEGER , PUBLIC :: nt_n, nt_a !: now / after indices (equal 1 or 2)
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 4.0 , NEMO Consortium (2011)
!! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $
......@@ -55,18 +57,18 @@ CONTAINS
INTEGER :: ierr
!!----------------------------------------------------------------------
!
ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime ), &
& v_abl (1:jpi,1:jpj,1:jpka,jptime ), &
& tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), &
& tke_abl (1:jpi,1:jpj,1:jpka,jptime ), &
& avm_abl (1:jpi,1:jpj,1:jpka ), &
& avt_abl (1:jpi,1:jpj,1:jpka ), &
& mxld_abl(1:jpi,1:jpj,1:jpka ), &
& mxlm_abl(1:jpi,1:jpj,1:jpka ), &
& fft_abl (1:jpi,1:jpj ), &
& pblh (1:jpi,1:jpj ), &
& msk_abl (1:jpi,1:jpj ), &
& rest_eq (1:jpi,1:jpj ), &
ALLOCATE( u_abl (A2D(0),1:jpka,jptime ), &
& v_abl (A2D(0),1:jpka,jptime ), &
& tq_abl (A2D(0),1:jpka,jptime,jptq), &
& tke_abl (A2D(0),1:jpka,jptime ), &
& avm_abl (A2D(0),1:jpka ), &
& avt_abl (A2D(0),1:jpka ), &
& mxld_abl(A2D(0),1:jpka ), &
& mxlm_abl(A2D(0),1:jpka ), &
& fft_abl (A2D(0) ), &
& pblh (A2D(1) ), & ! needed for optional smoothing
& msk_abl (A2D(0) ), &
& rest_eq (A2D(0) ), &
& e3t_abl (1:jpka), e3w_abl(1:jpka) , &
& ght_abl (1:jpka), ghw_abl(1:jpka) , STAT=ierr )
!
......
src/ABL/ablmod.F90
View file @ 66fb9523
......@@ -67,54 +67,50 @@ CONTAINS
!! - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step index
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pv_dta ! large-scale windj
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgu_dta ! large-scale hpgi
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pgv_dta ! large-scale hpgj
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pt_dta ! large-scale pot. temp.
REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pq_dta ! large-scale humidity
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp_dta ! sea-level pressure
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 - uoce||
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: plat ! latent heat flux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! ||tau||
!
INTEGER , INTENT(in ) :: kt ! time step index
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: psst ! sea-surface temperature [Celsius]
REAL(wp) , INTENT(in ), DIMENSION(A2D(nn_hls)) :: pssu ! sea-surface u (U-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(nn_hls)) :: pssv ! sea-surface v (V-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pssq ! sea-surface humidity
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pu_dta ! large-scale windi
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pv_dta ! large-scale windj
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pgu_dta ! large-scale hpgi
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pgv_dta ! large-scale hpgj
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pt_dta ! large-scale pot. temp.
REAL(wp) , INTENT(in ), DIMENSION(A2D(0),1:jpka) :: pq_dta ! large-scale humidity
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pslp_dta ! sea-level pressure
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pcd_du ! Cd x Du (T-point)
REAL(wp) , INTENT(inout), DIMENSION(A2D(0)) :: psen ! Ch x Du
REAL(wp) , INTENT(inout), DIMENSION(A2D(0)) :: pevp ! Ce x Du
REAL(wp) , INTENT(inout), DIMENSION(A2D(0)) :: pwndm ! ||uwnd - uoce||
REAL(wp) , INTENT( out), DIMENSION(A2D(0)) :: plat ! latent heat flux
REAL(wp) , INTENT( out), DIMENSION(A2D(nn_hls)) :: ptaui ! taux
REAL(wp) , INTENT( out), DIMENSION(A2D(nn_hls)) :: ptauj ! tauy
REAL(wp) , INTENT( out), DIMENSION(A2D(0)) :: ptaum ! ||tau||
#if defined key_si3
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptm_su ! ice-surface temperature [K]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu_ice ! ice-surface u (U-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du_ice ! Cd x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psen_ice ! Ch x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pevp_ice ! Ce x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! ||uwnd - uice||
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! ocean fraction
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui_ice ! ice-surface taux stress (U-point)
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: ptm_su ! ice-surface temperature [K]
REAL(wp) , INTENT(in ), DIMENSION(A2D(nn_hls)) :: pssu_ice ! ice-surface u (U-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(nn_hls)) :: pssv_ice ! ice-surface v (V-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pssq_ice ! ice-surface humidity
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pcd_du_ice ! Cd x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: psen_ice ! Ch x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pevp_ice ! Ce x Du over ice (T-point)
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pwndm_ice ! ||uwnd - uice||
REAL(wp) , INTENT(in ), DIMENSION(A2D(0)) :: pfrac_oce ! ocean fraction
REAL(wp) , INTENT( out), DIMENSION(A2D(nn_hls)) :: ptaui_ice ! ice-surface taux stress (T-point)
REAL(wp) , INTENT( out), DIMENSION(A2D(nn_hls)) :: ptauj_ice ! ice-surface tauy stress (T-point)
#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(A2D(0)) :: zwnd_i, zwnd_j
REAL(wp), DIMENSION(A2D(0)) :: zsspt, ztabs, zpre
REAL(wp), DIMENSION(A1Di(0),2:jpka) :: zCF
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: z_elem_a, z_elem_b, z_elem_c
!
REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_a
REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_b
REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_c
!
INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
LOGICAL :: SemiImp_Cor = .TRUE.
INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
LOGICAL :: SemiImp_Cor = .TRUE.
!
!!---------------------------------------------------------------------
!
......@@ -124,12 +120,12 @@ CONTAINS
WRITE(numout,*) '~~~~~~'
ENDIF
!
IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) )
IF( kt == nit000 ) ALLOCATE ( ustar2( A2D(0) ) )
IF( kt == nit000 ) ALLOCATE ( zrough( A2D(0) ) )
!! Compute ustar squared as Cd || Uatm-Uoce ||^2
!! needed for surface boundary condition of TKE
!! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zzoce = pCd_du (ji,jj) * pwndm (ji,jj)
#if defined key_si3
zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
......@@ -137,39 +133,40 @@ CONTAINS
#else
ustar2(ji,jj) = zzoce
#endif
!#LB: sorry Cdn_oce is gone:
!zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
zsspt(ji,jj) = theta_exner( psst(ji,jj)+rt0, pslp_dta(ji,jj) ) ! potential SST [K]
END_2D
!#LB: sorry Cdn_oce is gone:
!zrough(:,:) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(:,:), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
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(:,:) )
!zsspt(:,:) = theta_exner( psst(A2D(0))+rt0, pslp_dta(:,:) )
!
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (A1Di(0)) x (A1Dj(0))
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 2 *** Advance tracers to time n+1
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!-------------
DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj)
DO_1Dj(0,0) ! outer loop !--> tq_abl computed on (A1Di(0)) x (A1Dj(0))
!-------------
! Compute matrix elements for interior points
DO jk = 3, jpkam1
DO ji = 1, jpi ! vector opt.
DO_1Di(0,0)
z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
z_elem_c( ji, jk ) = - rDt_abl * Avt_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_1D
END DO
! Boundary conditions
DO ji = 1, jpi ! vector opt.
DO_1Di(0,0)
! Neumann at the bottom
z_elem_a( ji, 2 ) = 0._wp
z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )
......@@ -177,19 +174,18 @@ CONTAINS
z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
z_elem_c( ji, jpka ) = 0._wp
z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
END DO
END_1D
DO jtra = 1,jptq ! loop on active tracers
DO jk = 3, jpkam1
!DO ji = 2, jpim1
DO ji = 1,jpi !!GS: to be checked if needed
DO_1Di(0,0)
tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side
END DO
END_1D
END DO
IF(jtra == jp_ta) THEN
DO ji = 1,jpi ! surface boundary condition for temperature
DO_1Di(0,0) ! surface boundary condition for temperature
zztmp1 = psen(ji, jj)
zztmp2 = psen(ji, jj) * zsspt(ji, jj)
#if defined key_si3
......@@ -199,9 +195,9 @@ CONTAINS
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2
tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
END DO
END_1D
ELSE ! jp_qa
DO ji = 1,jpi ! surface boundary condition for humidity
DO_1Di(0,0) ! surface boundary condition for humidity
zztmp1 = pevp(ji, jj)
zztmp2 = pevp(ji, jj) * pssq(ji, jj)
#if defined key_si3
......@@ -211,31 +207,31 @@ CONTAINS
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2
tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra )
END DO
END_1D
END IF
!!
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1,jpi
DO_1Di(0,0)
zcff = 1._wp / z_elem_b( ji, 2 )
zCF ( ji, 2 ) = - zcff * z_elem_c( ji, 2 )
tq_abl( ji, jj, 2, nt_a, jtra ) = zcff * tq_abl( ji, jj, 2, nt_a, jtra )
END DO
END_1D
DO jk = 3, jpka
DO ji = 1,jpi
DO_1Di(0,0)
zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) )
zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) &
& - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
END DO
END_1D
END DO
!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
DO jk = jpkam1,2,-1
DO ji = 1,jpi
DO_1Di(0,0)
tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + &
& zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
END DO
END_1D
END DO
END DO !<-- loop on tracers
......@@ -254,7 +250,7 @@ CONTAINS
!-------------
!
! Advance u_abl & v_abl to time n+1
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2
u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( &
......@@ -269,13 +265,13 @@ CONTAINS
END_2D
!
!-------------
END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (A1Di(0)) x (A1Dj(0))
!-------------
!
IF( ln_geos_winds ) THEN
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) &
& - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk)
v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) &
......@@ -286,14 +282,14 @@ CONTAINS
END IF
!
IF( ln_hpgls_frc ) THEN
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
ENDDO
END_1D
ENDDO
ENDDO
END_1D
END IF
ELSE ! SemiImp_Cor = .FALSE.
......@@ -306,29 +302,25 @@ CONTAINS
!
IF( MOD( kt, 2 ) == 0 ) then
! Advance u_abl & v_abl to time n+1
DO jj = 1, jpj
DO ji = 1, jpi
zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) )
u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) )
v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a )
END DO
END DO
DO_2D( 0, 0, 0, 0 )
zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) )
u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) )
v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a )
END_2D
ELSE
DO jj = 1, jpj
DO ji = 1, jpi
zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) )
v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) )
u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a )
END DO
END DO
DO_2D( 0, 0, 0, 0 )
zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) )
v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) )
u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a )
END_2D
END IF
!
!-------------
END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (A1Di(0)) x (A1Dj(0))
!-------------
ENDIF ! ln_geos_winds
......@@ -340,18 +332,19 @@ CONTAINS
!
! Vertical diffusion for u_abl
!-------------
DO jj = 1, jpj ! outer loop
DO_1Dj(0,0) ! outer loop
!-------------
DO jk = 3, jpkam1
DO ji = 1, jpi
DO_1Di(0,0)
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_1D
END DO
DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
DO_1Di(0,0) ! boundary conditions
!++ 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 )
......@@ -381,29 +374,29 @@ CONTAINS
u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
!ENDIF
END DO
END_1D
!!
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1, jpi !DO ji = 2, jpi !!GS: TBI
DO_1Di(0,0)
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 )
END DO
END_1D
DO jk = 3, jpka
DO ji = 2, jpi
DO_1Di(0,0)
zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk ,nt_a) &
& - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
END DO
END_1D
END DO
DO jk = jpkam1,2,-1
DO ji = 2, jpi
DO_1Di(0,0)
u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
END DO
END_1D
END DO
!-------------
......@@ -413,18 +406,19 @@ CONTAINS
!
! Vertical diffusion for v_abl
!-------------
DO jj = 2, jpj ! outer loop
DO_1Dj(0,0) ! outer loop
!-------------
!
DO jk = 3, jpkam1
DO ji = 1, jpi
DO_1Di(0,0)
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_1D
END DO
DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
DO_1Di(0,0) ! boundary conditions
!++ 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 )
......@@ -454,29 +448,29 @@ CONTAINS
v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
!ENDIF
END DO
END_1D
!!
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1, jpi
DO_1Di(0,0)
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
END_1D
DO jk = 3, jpka
DO ji = 1, jpi
DO_1Di(0,0)
zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )
zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk ,nt_a) &
& - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
END DO
END_1D
END DO
DO jk = jpkam1,2,-1
DO ji = 1, jpi
DO_1Di(0,0)
v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
END DO
END_1D
END DO
!
!-------------
......@@ -491,7 +485,7 @@ CONTAINS
!-------------
DO jk = 2, jpka ! outer loop
!-------------
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
DO_2D( 0, 0, 0, 0)
zcff1 = pblh( ji, jj )
zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
......@@ -514,7 +508,7 @@ CONTAINS
!-------------
DO jk = 2, jpka ! outer loop
!-------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0)
zcff1 = pblh( ji, jj )
zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) )
zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )
......@@ -528,21 +522,17 @@ CONTAINS
tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) &
& + zcff * pq_dta( ji, jj, jk )
END_2D
!-------------
END DO ! end outer loop
!-------------
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 6 *** MPI exchanges & IOM outputs
! ! 6 *** IOM outputs
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
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...
!
#if defined key_xios
! 2D & first ABL level
IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(:,: ) )
IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(A2D(0) ) )
IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) )
IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) )
IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl", tq_abl(:,:,2,nt_a,jp_ta) )
......@@ -588,7 +578,7 @@ CONTAINS
! ! 7 *** Finalize flux computation
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
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 ) )
......@@ -599,15 +589,13 @@ CONTAINS
rhoa( ji, jj ) = zcff
END_2D
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
DO_2D( 0, 0, 0, 0 )
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 )
!
! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) &
& + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj)
zztmp = rhoa(ji,jj) * pcd_du(ji,jj)
......@@ -617,25 +605,18 @@ CONTAINS
zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
END_2D
! ... utau, vtau at U- and V_points, resp.
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
CALL iom_put( "taum_oce", ptaum )
! ... utau, vtau at T-points
DO_2D( 0, 0, 0, 0 )
zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
ptaui(ji,jj) = zwnd_i(ji,jj) * msk_abl(ji,jj)
ptauj(ji,jj) = zwnd_j(ji,jj) * msk_abl(ji,jj)
END_2D
!
CALL lbc_lnk( 'ablmod', ptaui(:,:), 'U', -1.0_wp, ptauj(:,:), 'V', -1.0_wp )
CALL iom_put( "taum_oce", ptaum )
CALL lbc_lnk( 'ablmod', ptaui(:,:) , 'T', -1.0_wp, ptauj(:,:) , 'T', -1.0_wp )
IF(sn_cfctl%l_prtctl) THEN
CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=umask, &
& tab2d_2=ptauj , clinfo2=' vtau : ', mask2=vmask )
CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=tmask, &
& tab2d_2=ptauj , clinfo2=' vtau : ', mask2=tmask )
CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' )
ENDIF
......@@ -643,37 +624,16 @@ 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 : ' )
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
DO_2D( 0, 0, 0, 0 )
zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
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) * 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) * rn_vfac )
DO_2D(0,0,0,0)
zztmp = rhoa(ji,jj) * pCd_du_ice(ji,jj)
ptaui_ice(ji,jj) = zztmp * ( u_abl(ji,jj,2,nt_a) - 0.5_wp * ( pssu_ice(ji,jj) + pssu_ice(ji-1,jj ) ) * rn_vfac )
ptauj_ice(ji,jj) = zztmp * ( v_abl(ji,jj,2,nt_a) - 0.5_wp * ( pssv_ice(ji,jj) + pssv_ice(ji ,jj-1) ) * rn_vfac )
END_2D
CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
CALL lbc_lnk( 'ablmod', ptaui_ice(:,:) , 'T', -1.0_wp, ptauj_ice(:,:) , 'T', -1.0_wp )
!
IF(sn_cfctl%l_prtctl) THEN
CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask, &
& tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=vmask )
CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=tmask, &
& tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=tmask )
END IF
#endif
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
......@@ -707,21 +667,22 @@ CONTAINS
!! - avmu, avmv : production of TKE by shear at u and v-points
!! (= Kz dz[Ub] * dz[Un] )
!! ---------------------------------------------------------------------
INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
INTEGER, DIMENSION(1:jpi ) :: ikbl
REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89
REAL(wp) :: zwndi, zwndj
REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2
REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2
REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF
REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a
REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b
REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c
LOGICAL :: ln_Patankar = .FALSE.
LOGICAL :: ln_dumpvar = .FALSE.
LOGICAL , DIMENSION(1:jpi ) :: ln_foundl
INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn
INTEGER, DIMENSION(A1Di(0)) :: ikbl
REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89
REAL(wp) :: zwndi, zwndj
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: zsh2
REAL(wp), DIMENSION(A2D(0) ,1:jpka) :: zbn2
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: zFC, zRH, zCF
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: z_elem_a
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: z_elem_b
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: z_elem_c
LOGICAL :: ln_Patankar = .FALSE.
LOGICAL :: ln_dumpvar = .FALSE.
LOGICAL , DIMENSION(A1Di(0)) :: ln_foundl
!
tind = nt_n
ziRic = 1._wp / rn_Ric
......@@ -730,33 +691,33 @@ CONTAINS
! ! Advance TKE equation to time n+1
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!-------------
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
!-------------
!
! Compute vertical shear
DO jk = 2, jpkam1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zcff = 1.0_wp / e3w_abl( jk )**2
zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
END DO
END_1D
END DO
!
! Compute brunt-vaisala frequency
DO jk = 2, jpkam1
DO ji = 1,jpi
DO_1Di( 0, 0 )
zcff = grav * itvref / e3w_abl( jk )
zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta)
zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) &
& - tq_abl( ji, jj, jk , tind, jp_ta) * tq_abl( ji, jj, jk , tind, jp_qa)
zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi)
END DO
END_1D
END DO
!
! Terms for the tridiagonal problem
DO jk = 2, jpkam1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point
zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
......@@ -773,10 +734,10 @@ CONTAINS
& - e3w_abl(jk) * rDt_abl * zbuoy
tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side
END IF
END DO
END_1D
END DO
DO ji = 1,jpi ! vector opt.
DO_1Di( 0, 0 )
zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min )
zetop = tke_min
......@@ -805,44 +766,41 @@ CONTAINS
!!
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1,jpi
DO_1Di( 0, 0 )
zcff = 1._wp / z_elem_b( ji, 1 )
zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 )
tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a )
END DO
END_1D
DO jk = 2, jpka
DO ji = 1,jpi
DO_1Di( 0, 0 )
zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
tke_abl(ji,jj,jk,nt_a) = zcff * ( tke_abl(ji,jj,jk ,nt_a) &
& - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
END DO
END_1D
END DO
DO jk = jpkam1,1,-1
DO ji = 1,jpi
DO_1Di( 0, 0 )
tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
END DO
END_1D
END DO
!!FL should not be needed because of Patankar procedure
tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
!!FL should not be needed because of Patankar procedure
tke_abl(A2D(0),1:jpka,nt_a) = MAX( tke_abl(A2D(0),1:jpka,nt_a), tke_min )
!!
!! Diagnose PBL height
!! ----------------------------------------------------------
!
! arrays zRH, zFC and zCF are available at this point
! and zFC(:, 1 ) = 0.
! diagnose PBL height based on zsh2 and zbn2
zFC ( : ,1) = 0._wp
ikbl( 1:jpi ) = 0
ikbl( A1Di(0) ) = 0
DO jk = 2,jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
zcff = ghw_abl( jk-1 )
zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
zcff = ghw_abl( jk )
......@@ -854,11 +812,11 @@ CONTAINS
- rn_Cek * ( fft_abl( ji, jj ) * fft_abl( ji, jj ) ) ) &
& )
IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
END DO
END_1D
END DO
!
! finalize the computation of the PBL height
DO ji = 1, jpi
DO_1Di( 0, 0 )
jk = ikbl(ji)
IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
pblh( ji, jj ) = ( ghw_abl( jk-1 ) * zFC( ji, jk ) &
......@@ -869,16 +827,16 @@ CONTAINS
ELSE
pblh( ji, jj ) = ghw_abl(jpka)
END IF
END DO
END_1D
!-------------
END DO
!-------------
!
! Optional : could add pblh smoothing if pblh is noisy horizontally ...
IF(ln_smth_pblh) THEN
CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
CALL smooth_pblh( pblh, msk_abl )
CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp )
CALL smooth_pblh( pblh, tmask(A2D(0),1) )
CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp )
ENDIF
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! Diagnostic mixing length computation
......@@ -889,55 +847,55 @@ CONTAINS
CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom
# define zlup zRH
# define zldw zFC
DO jj = 1, jpj ! outer loop
DO_1Dj(0,0) ! outer loop
!
DO ji = 1, jpi
DO_1Di(0,0)
mxld_abl( ji, jj, 1 ) = mxl_min
mxld_abl( ji, jj, jpka ) = mxl_min
mxlm_abl( ji, jj, 1 ) = mxl_min
mxlm_abl( ji, jj, jpka ) = mxl_min
zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc
zlup ( ji, jpka ) = mxl_min
END DO
END_1D
!
DO jk = 2, jpkam1
DO ji = 1, jpi
DO_1Di(0,0)
zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
mxlm_abl( ji, jj, jk ) = MAX( mxl_min, &
& SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
END DO
END_1D
END DO
!
! Limit mxl
DO jk = jpkam1,1,-1
DO ji = 1, jpi
DO_1Di(0,0)
zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
END DO
END_1D
END DO
!
DO jk = 2, jpka
DO ji = 1, jpi
DO_1Di(0,0)
zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
END DO
END_1D
END DO
!
! DO jk = 1, jpka
! DO ji = 1, jpi
! DO_1Di(0,0)
! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) )
! END DO
! END_1D
! END DO
!
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di(0,0)
! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp)
zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
END DO
END_1D
END DO
!
END DO
END_1D ! outer loop
# undef zlup
# undef zldw
!
......@@ -945,18 +903,18 @@ CONTAINS
CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
# define zlup zRH
# define zldw zFC
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
!
DO jk = 2, jpkam1
DO ji = 1,jpi
zcff = 1.0_wp / e3w_abl( jk )**2
DO_1Di( 0, 0 )
zcff = 1.0_wp / e3w_abl( jk )**2
zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
ENDDO
ENDDO
!
DO ji = 1, jpi
END_1D
END DO
!
DO_1Di( 0, 0 )
zcff = zrough(ji,jj) * rn_Lsfc
mxld_abl ( ji, jj, 1 ) = zcff
mxld_abl ( ji, jj, jpka ) = mxl_min
......@@ -964,42 +922,42 @@ CONTAINS
mxlm_abl ( ji, jj, jpka ) = mxl_min
zldw ( ji, 1 ) = zcff
zlup ( ji, jpka ) = mxl_min
END DO
END_1D
!
DO jk = 2, jpkam1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
zcff = 2.0_wp*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) &
& + SQRT(rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.0_wp*zbuoy ) )
mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
END DO
END_1D
END DO
!
! Limit mxl
DO jk = jpkam1,1,-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
END DO
END_1D
END DO
!
DO jk = 2, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
END DO
END_1D
END DO
!
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
!mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
!zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp)
zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
!mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
END DO
END_1D
END DO
!
END DO
!
END_1D ! outer loop
# undef zlup
# undef zldw
!
......@@ -1009,35 +967,35 @@ CONTAINS
# define zldw zFC
! zCF is used for matrix inversion
!
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
DO ji = 1, jpi
DO_1Di( 0, 0 )
zcff = zrough(ji,jj) * rn_Lsfc
zlup( ji, 1 ) = zcff
zldw( ji, 1 ) = zcff
zlup( ji, jpka ) = mxl_min
zldw( ji, jpka ) = mxl_min
END DO
END_1D
DO jk = 2,jpka-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
END DO
END_1D
END DO
!!
!! BL89 search for lup
!! ----------------------------------------------------------
DO jk=2,jpka-1
!
DO ji = 1, jpi
DO_1Di( 0, 0 )
zCF(ji,1:jpka) = 0._wp
zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
ln_foundl(ji ) = .false.
END DO
END_1D
!
DO jkup=jk+1,jpka-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
......@@ -1052,7 +1010,7 @@ CONTAINS
zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
ln_foundl(ji) = .true.
END IF
END DO
END_1D
END DO
!
END DO
......@@ -1061,14 +1019,14 @@ CONTAINS
!! ----------------------------------------------------------
DO jk=2,jpka-1
!
DO ji = 1, jpi
DO_1Di( 0, 0 )
zCF(ji,1:jpka) = 0._wp
zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
ln_foundl(ji ) = .false.
END DO
END_1D
!
DO jkdwn=jk-1,1,-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
......@@ -1083,21 +1041,21 @@ CONTAINS
zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
ln_foundl(ji) = .true.
END IF
END DO
END_1D
END DO
!
END DO
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
!zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp)
zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
END DO
END_1D
END DO
END DO
END_1D ! outer loop
# undef zlup
# undef zldw
!
......@@ -1107,46 +1065,46 @@ CONTAINS
# define zldw zFC
! zCF is used for matrix inversion
!
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
!
DO jk = 2, jpkam1
DO ji = 1,jpi
DO_1Di( 0, 0 )
zcff = 1.0_wp / e3w_abl( jk )**2
zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2
zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2
zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
ENDDO
ENDDO
END_1D
ENDDO
zsh2(:, 1) = zsh2( :, 2)
zsh2(:, jpka) = zsh2( :, jpkam1)
DO ji = 1, jpi
zcff = zrough(ji,jj) * rn_Lsfc
zlup( ji, 1 ) = zcff
zldw( ji, 1 ) = zcff
zlup( ji, jpka ) = mxl_min
zldw( ji, jpka ) = mxl_min
END DO
DO_1Di( 0, 0 )
zcff = zrough(ji,jj) * rn_Lsfc
zlup( ji, 1 ) = zcff
zldw( ji, 1 ) = zcff
zlup( ji, jpka ) = mxl_min
zldw( ji, jpka ) = mxl_min
END_1D
DO jk = 2,jpka-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1)
END DO
END_1D
END DO
!!
!! BL89 search for lup
!! ----------------------------------------------------------
DO jk=2,jpka-1
!
DO ji = 1, jpi
DO_1Di( 0, 0 )
zCF(ji,1:jpka) = 0._wp
zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
ln_foundl(ji ) = .false.
END DO
END_1D
!
DO jkup=jk+1,jpka-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall )
zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
......@@ -1165,7 +1123,7 @@ CONTAINS
zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk)
ln_foundl(ji) = .true.
END IF
END DO
END_1D
END DO
!
END DO
......@@ -1174,14 +1132,14 @@ CONTAINS
!! ----------------------------------------------------------
DO jk=2,jpka-1
!
DO ji = 1, jpi
DO_1Di( 0, 0 )
zCF(ji,1:jpka) = 0._wp
zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a )
ln_foundl(ji ) = .false.
END DO
END_1D
!
DO jkdwn=jk-1,1,-1
DO ji = 1, jpi
DO_1Di( 0, 0 )
zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall )
zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) &
......@@ -1200,18 +1158,18 @@ CONTAINS
zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn )
ln_foundl(ji) = .true.
END IF
END DO
END_1D
END DO
!
END DO
DO jk = 1, jpka
DO ji = 1, jpi
DO_1Di( 0, 0 )
!zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp)
zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min )
END DO
END_1D
END DO
END DO
......@@ -1225,10 +1183,10 @@ CONTAINS
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!-------------
DO jj = 1, jpj ! outer loop
DO_1Dj( 0, 0 ) ! outer loop
!-------------
DO jk = 1, jpka
DO ji = 1, jpi ! vector opt.
DO_1Di( 0, 0 )
zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) &
& * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z)
......@@ -1236,10 +1194,10 @@ CONTAINS
!!FL: MAX function probably useless because of the definition of mxl_min
Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak )
Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak )
END DO
END_1D
END DO
!-------------
END DO
END_1D ! outer loop
!-------------
!---------------------------------------------------------------------------------------------------
......@@ -1257,46 +1215,47 @@ CONTAINS
!! ** Purpose : 2D Hanning filter on atmospheric PBL height
!!
!! ---------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
INTEGER :: ji,jj
REAL(wp) :: smth_a, smth_b
REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY
REAL(wp) :: zumsk,zvmsk
!!
REAL(wp), DIMENSION(A2D(1)), INTENT(in ) :: msk
REAL(wp), DIMENSION(A2D(1)), INTENT(inout) :: pvar2d
INTEGER :: ji,jj
REAL(wp) :: smth_a, smth_b
REAL(wp), DIMENSION(A2D(1)) :: zdX,zdY,zFX,zFY
REAL(wp) :: zumsk,zvmsk
!!=========================================================
!!
!! Hanning filter
smth_a = 1._wp / 8._wp
smth_b = 1._wp / 4._wp
!
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
DO_2D( 1, 0, 1, 1 )
zumsk = msk(ji,jj) * msk(ji+1,jj)
zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk
END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
DO_2D( 1, 1, 1, 0 )
zvmsk = msk(ji,jj) * msk(ji,jj+1)
zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk
END_2D
DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
DO_2D( 0, 0, 1, 0 )
zFY ( ji, jj ) = zdY ( ji, jj ) &
& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
& + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & ! add () for NP repro
& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
END_2D
DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
DO_2D( 1, 0, 0, 0 )
zFX( ji, jj ) = zdX( ji, jj ) &
& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
& + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) & ! add () for NP repro
& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
END_2D
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
DO_2D( 0, 0, 0, 0 )
pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
& + msk(ji,jj) * smth_b * ( &
& zFX( ji, jj ) - zFX( ji-1, jj ) &
& +zFY( ji, jj ) - zFY( ji, jj-1 ) )
& ( zFX( ji, jj ) - zFX( ji-1, jj ) ) & ! add () for NP repro
& + ( zFY( ji, jj ) - zFY( ji, jj-1 ) ) )
END_2D
!---------------------------------------------------------------------------------------------------
......
src/ABL/ablrst.F90
View file @ 66fb9523
......@@ -202,12 +202,12 @@ CONTAINS
! --- mandatory fields --- !
CALL iom_get( numrar, jpdom_auto, 'u_abl', u_abl(:,:,:,nt_n ), cd_type = 'T', psgn = -1._wp )
CALL iom_get( numrar, jpdom_auto, 'v_abl', v_abl(:,:,:,nt_n ), cd_type = 'T', psgn = -1._wp )
CALL iom_get( numrar, jpdom_auto, 't_abl', tq_abl(:,:,:,nt_n,jp_ta), kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto, 't_abl', tq_abl(:,:,:,nt_n,jp_ta) ) !, kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto, 'q_abl', tq_abl(:,:,:,nt_n,jp_qa) )
CALL iom_get( numrar, jpdom_auto, 'tke_abl', tke_abl(:,:,:,nt_n ) )
CALL iom_get( numrar, jpdom_auto, 'avm_abl', avm_abl(:,:,: ), kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto, 'avm_abl', avm_abl(:,:,: ) ) !, kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto, 'avt_abl', avt_abl(:,:,: ) )
CALL iom_get( numrar, jpdom_auto,'mxld_abl',mxld_abl(:,:,: ), kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto,'mxld_abl',mxld_abl(:,:,: ) ) !, kfill = jpfillcopy )
CALL iom_get( numrar, jpdom_auto, 'pblh', pblh(:,: ), kfill = jpfillcopy )
IF(.NOT.lrxios) CALL iom_delay_rst( 'READ', 'ABL', numrar ) ! read only abl delayed global communication variables
......
src/ABL/sbcabl.F90
View file @ 66fb9523
......@@ -44,6 +44,8 @@ MODULE sbcabl
PUBLIC sbc_abl_init ! routine called in sbcmod module
PUBLIC sbc_abl ! routine called in sbcmod module
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 3.7 , NEMO-consortium (2014)
!! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
......@@ -259,7 +261,7 @@ CONTAINS
rest_eq(:,:) = 1._wp
END IF
! T-mask
msk_abl(:,:) = tmask(:,:,1)
msk_abl(:,:) = tmask(A2D(0),1)
!!-------------------------------------------------------------------------------------------
......@@ -307,8 +309,8 @@ CONTAINS
!! - Perform 1 time-step of the ABL model
!! - Finalize flux computation in blk_oce_2
!!
!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
!! - vtau : j-component of the stress at V-point (N/m2)
!! ** Outputs : - utau : i-component of the stress at T-point (N/m2)
!! - vtau : j-component of the stress at T-point (N/m2)
!! - taum : Wind stress module at T-point (N/m2)
!! - wndm : Wind speed module at T-point (m/s)
!! - qsr : Solar heat flux over the ocean (W/m2)
......@@ -318,9 +320,11 @@ 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
REAL(wp), DIMENSION(A2D(0)) :: zssq, zcd_du, zsen, zlat, zevp
#if defined key_si3
REAL(wp), DIMENSION(jpi,jpj) :: zssqi, zcd_dui, zseni, zevpi
!REAL(wp), DIMENSION(jpi,jpj) :: zssqi, zcd_dui, zseni, zevpi
REAL(wp), DIMENSION(A2D(0)) :: zssqi, zcd_dui, zseni, zevpi
#endif
INTEGER :: jbak, jbak_dta, ji, jj
!!---------------------------------------------------------------------
......@@ -339,7 +343,7 @@ CONTAINS
CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in
& tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in
& sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in
& sf(jp_slp )%fnow(:,:,1) , sst_m(A2D(0)), ssu_m(A2D(0)), ssv_m(A2D(0)), & ! <<= in
& sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in
& sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in
& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
......@@ -347,7 +351,7 @@ CONTAINS
#if defined key_si3
CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in
& tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in
& sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in
& sf(jp_slp)%fnow(:,:,1) , u_ice(A2D(0)), v_ice(A2D(0)), tm_su , & ! <<= in
& pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out
#endif
......@@ -355,7 +359,7 @@ CONTAINS
!! 3 - Advance ABL variables from now (n) to after (n+1)
!!-------------------------------------------------------------------------------------------
CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq, & ! <<= in
CALL abl_stp( kt, tsk_m(A2D(0)), ssu_m, ssv_m, zssq, & ! <<= in
& sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:), & ! <<= in
& sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), & ! <<= in
& sf(jp_slp )%fnow(:,:,1), & ! <<= in
......@@ -364,7 +368,7 @@ CONTAINS
& zlat, wndm, utau, vtau, taum & ! =>> out
#if defined key_si3
& , tm_su, u_ice, v_ice, zssqi, zcd_dui & ! <<= in
& , zseni, zevpi, wndm_ice, ato_i & ! <<= in
& , zseni, zevpi, wndm_ice, ato_i(A2D(0)) & ! <<= in
& , utau_ice, vtau_ice & ! =>> out
#endif
& )
......
src/ICE/ice.F90
View file @ 66fb9523
......@@ -451,6 +451,8 @@ MODULE ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_bot !: Bottom conduction flux (W/m2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface conduction flux (W/m2)
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
......@@ -464,71 +466,103 @@ CONTAINS
!!-----------------------------------------------------------------
INTEGER :: ice_alloc
!
INTEGER :: ierr(16), ii
INTEGER :: ierr(21), ii
!!-----------------------------------------------------------------
ierr(:) = 0
ii = 0
! ----------------- !
! == FULL ARRAYS == !
! ----------------- !
! * Ice global state variables
ii = ii + 1
ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , STAT=ierr(ii) )
ii = 1
ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , fraz_frac (jpi,jpj) , &
& strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , &
& delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , &
& aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) )
ii = ii + 1
ALLOCATE( h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , &
& v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , &
& s_i (jpi,jpj,jpl) , sv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl) , &
& a_ip (jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), &
& v_il (jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , &
& t_su (jpi,jpj,jpl) , t_s (jpi,jpj,nlay_s,jpl) , t_i(jpi,jpj,nlay_i,jpl) , sz_i(jpi,jpj,nlay_i,jpl) , &
& ato_i(jpi,jpj) , STAT = ierr(ii) )
ii = ii + 1
ALLOCATE( t_bo (jpi,jpj) , wfx_snw_sni(jpi,jpj) , &
& wfx_snw (jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) , &
& wfx_ice (jpi,jpj) , wfx_sub (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam (jpi,jpj) , &
& wfx_pnd (jpi,jpj) , &
& wfx_bog (jpi,jpj) , wfx_dyn (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , &
& wfx_res (jpi,jpj) , wfx_sni (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , &
& rn_amax_2d (jpi,jpj) , &
& qsb_ice_bot(jpi,jpj) , qlead (jpi,jpj) , &
& sfx_res (jpi,jpj) , sfx_bri (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) , &
& sfx_bog (jpi,jpj) , sfx_bom (jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , &
& hfx_res (jpi,jpj) , hfx_snw (jpi,jpj) , hfx_sub(jpi,jpj) , &
& qt_atm_oi (jpi,jpj) , qt_oce_ai (jpi,jpj) , fhld (jpi,jpj) , &
& hfx_sum (jpi,jpj) , hfx_bom (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , &
& hfx_opw (jpi,jpj) , hfx_thd (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , &
& hfx_err_dif(jpi,jpj) , wfx_err_sub(jpi,jpj) , STAT=ierr(ii) )
ALLOCATE( e_s(jpi,jpj,nlay_s,jpl) , e_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
! * Ice global state variables
! * Before values of global variables
ii = ii + 1
ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , &
& h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , &
& v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , &
& s_i (jpi,jpj,jpl) , sv_i (jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , &
& oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) )
ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) )
! * fluxes
ii = ii + 1
ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , &
& vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) , &
& et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) , &
& sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) , &
& om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), icb_mask(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( wfx_res(jpi,jpj) , sfx_res(jpi,jpj) , hfx_res(jpi,jpj) , STAT=ierr(ii) ) ! full arrays since it is used in conjonction with global variables
! * ice rheology
ii = ii+1
ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , &
& strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , &
& aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , &
& icb_mask (jpi,jpj) , STAT=ierr(ii) )
! * mean and total
ii = ii + 1
ALLOCATE( vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , & ! full arrays since they are used in rheology
& vt_ip(jpi,jpj) , vt_il(jpi,jpj) , at_ip(jpi,jpj) , STAT=ierr(ii) )
! * others
ii = ii + 1
ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
ALLOCATE( t_bo(jpi,jpj) , rn_amax_2d(jpi,jpj) , STAT=ierr(ii) )
! -------------------- !
! == REDUCED ARRAYS == !
! -------------------- !
! * Ice global state variables
ii = ii + 1
ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , sz_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
ALLOCATE( bv_i(A2D(0),jpl) , a_ip_frac(A2D(0),jpl) , a_ip_eff(A2D(0),jpl) , STAT=ierr(ii) )
! * Before values of global variables
ii = ii + 1
ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), &
& v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , &
& dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) )
ALLOCATE( at_i_b(A2D(0)) , h_i_b (A2D(0),jpl) , a_i_b(A2D(0),jpl) , v_i_b(A2D(0),jpl) , &
& v_s_b (A2D(0),jpl) , h_s_b (A2D(0),jpl) , &
& v_ip_b(A2D(0),jpl) , v_il_b(A2D(0),jpl) , &
& sv_i_b(A2D(0),jpl) , e_i_b (A2D(0),nlay_i,jpl) , e_s_b(A2D(0),nlay_s,jpl) , STAT=ierr(ii) )
! * fluxes
ii = ii + 1
ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) )
ALLOCATE( qsb_ice_bot(A2D(0)) , qlead (A2D(0)) , qt_atm_oi (A2D(0)) , qt_oce_ai (A2D(0)) , fhld (A2D(0)) , &
& wfx_snw_sni(A2D(0)) , wfx_snw (A2D(0)) , wfx_snw_dyn(A2D(0)) , wfx_snw_sum(A2D(0)) , wfx_snw_sub(A2D(0)) , &
& wfx_ice (A2D(0)) , wfx_sub (A2D(0)) , wfx_ice_sub(A2D(0)) , wfx_lam (A2D(0)) , &
& wfx_pnd (A2D(0)) , &
& wfx_bog (A2D(0)) , wfx_dyn (A2D(0)) , wfx_bom(A2D(0)) , wfx_sum(A2D(0)) , &
& wfx_sni (A2D(0)) , wfx_opw (A2D(0)) , wfx_spr(A2D(0)) , &
& &
& sfx_bri (A2D(0)) , sfx_dyn (A2D(0)) , sfx_sub(A2D(0)) , sfx_lam(A2D(0)) , &
& sfx_bog (A2D(0)) , sfx_bom (A2D(0)) , sfx_sum(A2D(0)) , sfx_sni(A2D(0)) , sfx_opw(A2D(0)) , &
& hfx_snw (A2D(0)) , hfx_sub (A2D(0)) , &
& hfx_sum (A2D(0)) , hfx_bom (A2D(0)) , hfx_bog(A2D(0)) , hfx_dif(A2D(0)) , &
& hfx_opw (A2D(0)) , hfx_thd (A2D(0)) , hfx_dyn(A2D(0)) , hfx_spr(A2D(0)) , &
& hfx_err_dif(A2D(0)) , wfx_err_sub(A2D(0)) , STAT=ierr(ii) )
ii = ii + 1
ALLOCATE( qtr_ice_bot(A2D(0),jpl) , cnd_ice(A2D(0),jpl) , t1_ice(A2D(0),jpl) , STAT=ierr(ii) )
! * ice rheology
ii = ii+1
ALLOCATE( delta_i(A2D(0)) , divu_i(A2D(0)) , shear_i(A2D(0)) , STAT=ierr(ii) )
! * Old values of global variables
! * mean and total
ii = ii + 1
ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), &
& v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , &
& a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , &
& STAT=ierr(ii) )
ALLOCATE( st_i (A2D(0)) , et_i (A2D(0)) , et_s(A2D(0)) , hm_i (A2D(0)) , &
& hm_ip(A2D(0)) , hm_il(A2D(0)) , tm_i(A2D(0)) , tm_s (A2D(0)) , &
& sm_i (A2D(0)) , hm_s (A2D(0)) , om_i(A2D(0)) , bvm_i(A2D(0)) , &
& tm_su(A2D(0)) , STAT=ierr(ii) )
! * others
ii = ii + 1
ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
ALLOCATE( tau_icebfr(A2D(0)) , dh_i_sum_2d(A2D(0),jpl) , dh_s_mlt_2d(A2D(0),jpl) , STAT=ierr(ii) )
ii = 1
ALLOCATE( ht_i_new (A2D(0)) , fraz_frac (A2D(0)) , STAT=ierr(ii) )
! * Ice thickness distribution variables
ii = ii + 1
......@@ -536,19 +570,19 @@ CONTAINS
! * Ice diagnostics
ii = ii + 1
ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), &
& diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), &
& diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), &
& diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( diag_trp_vi (A2D(0)) , diag_trp_vs (A2D(0)) , diag_trp_ei (A2D(0)) , &
& diag_trp_es (A2D(0)) , diag_trp_sv (A2D(0)) , diag_heat (A2D(0)) , &
& diag_sice (A2D(0)) , diag_vice (A2D(0)) , diag_vsnw (A2D(0)) , diag_aice(A2D(0)) , diag_vpnd(A2D(0)), &
& diag_adv_mass(A2D(0)) , diag_adv_salt(A2D(0)) , diag_adv_heat(A2D(0)) , STAT=ierr(ii) )
! * Ice conservation
ii = ii + 1
ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj), &
& diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) )
ALLOCATE( diag_v (A2D(0)) , diag_s (A2D(0)) , diag_t (A2D(0)), &
& diag_fv(A2D(0)) , diag_fs(A2D(0)) , diag_ft(A2D(0)), STAT=ierr(ii) )
! * SIMIP diagnostics
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) )
ALLOCATE( t_si(A2D(0),jpl) , tm_si(A2D(0)) , qcn_ice_bot(A2D(0),jpl) , qcn_ice_top(A2D(0),jpl) , STAT = ierr(ii) )
ice_alloc = MAXVAL( ierr(:) )
IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
......
src/ICE/ice1d.F90
View file @ 66fb9523
......@@ -134,9 +134,6 @@ MODULE ice1D
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_i_1d !: Ice enthalpy per unit volume
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_s_1d !: Snow enthalpy per unit volume
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: eh_i_old !: ice heat content (q*h, J.m-2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m)
! Conduction flux diagnostics (SIMIP)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qcn_ice_bot_1d
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qcn_ice_top_1d
......@@ -166,6 +163,9 @@ MODULE ice1D
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ib_2d
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_ib_2d
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_i_2d
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_s_2d
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
......@@ -218,8 +218,7 @@ CONTAINS
!
ii = ii + 1
ALLOCATE( t_s_1d (jpij,nlay_s) , t_i_1d (jpij,nlay_i) , sz_i_1d(jpij,nlay_i) , &
& e_i_1d (jpij,nlay_i) , e_s_1d (jpij,nlay_s) , &
& eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )
& e_i_1d (jpij,nlay_i) , e_s_1d (jpij,nlay_s) , STAT=ierr(ii) )
!
ii = ii + 1
ALLOCATE( qcn_ice_bot_1d(jpij) , qcn_ice_top_1d(jpij) , STAT=ierr(ii) )
......@@ -235,6 +234,9 @@ CONTAINS
& v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) , &
& a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) , &
& STAT=ierr(ii) )
!
ii = ii + 1
ALLOCATE( e_i_2d(jpij,nlay_i,jpl) , e_s_2d(jpij,nlay_s,jpl) , STAT=ierr(ii) )
ice1D_alloc = MAXVAL( ierr(:) )
IF( ice1D_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice1D_alloc: failed to allocate arrays.' )
......
src/ICE/icealb.F90
View file @ 66fb9523
......@@ -48,7 +48,7 @@ MODULE icealb
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
SUBROUTINE ice_alb( ld_pnd_alb, pt_su, ph_ice, ph_snw, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
!!----------------------------------------------------------------------
!! *** ROUTINE ice_alb ***
!!
......@@ -94,16 +94,16 @@ CONTAINS
!! Brandt et al. 2005, J. Climate, vol 18
!! Grenfell & Perovich 2004, JGR, vol 109
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pt_su ! ice surface temperature (Kelvin)
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_ice ! sea-ice thickness
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_snw ! snow depth
LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area)
REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_pnd ! melt pond depth
REAL(wp), INTENT(in ), DIMENSION(:,:) :: pcloud_fra ! cloud fraction
REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_ice ! albedo of ice
LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: pt_su ! ice surface temperature (Kelvin)
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_ice ! sea-ice thickness
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_snw ! snow depth
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: pafrac_pnd ! melt pond relative fraction (per unit ice area)
REAL(wp), INTENT(in ), DIMENSION(A2D(0),jpl) :: ph_pnd ! melt pond depth
REAL(wp), INTENT(in ), DIMENSION(A2D(0)) :: pcloud_fra ! cloud fraction
REAL(wp), INTENT( out), DIMENSION(A2D(0),jpl) :: palb_ice ! albedo of ice
!
REAL(wp), DIMENSION(jpi,jpj,jpl) :: za_s_fra ! ice fraction covered by snow
REAL(wp), DIMENSION(A2D(0),jpl) :: za_s_fra ! ice fraction covered by snow
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar
REAL(wp) :: z1_href_pnd ! inverse of the characteristic length scale (Lecomte et al. 2015)
......@@ -121,10 +121,10 @@ CONTAINS
z1_c3 = 1._wp / 0.02_wp
z1_c4 = 1._wp / 0.03_wp
!
CALL ice_var_snwfra( ph_snw, za_s_fra ) ! calculate ice fraction covered by snow
CALL ice_var_snwfra( ph_snw(:,:,:), za_s_fra(:,:,:) ) ! calculate ice fraction covered by snow
!
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! palb_ice used over the full domain in icesbc
DO_2D( 0, 0, 0, 0 ) ! palb_ice used over the full domain in icesbc
!
!---------------------------------------------!
!--- Specific snow, ice and pond fractions ---!
......@@ -164,10 +164,10 @@ CONTAINS
zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )
!
! !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * smask0(ji,jj)
!
zalb_cs = zalb_os - ( - 0.1010_wp * zalb_os * zalb_os &
& + 0.1933_wp * zalb_os - 0.0148_wp ) * tmask(ji,jj,1)
& + 0.1933_wp * zalb_os - 0.0148_wp ) * smask0(ji,jj)
!
! albedo depends on cloud fraction because of non-linear spectral effects
palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os
......
src/ICE/icecor.F90
View file @ 66fb9523
......@@ -71,7 +71,10 @@ CONTAINS
WHERE( a_i(:,:,:) >= epsi20 ) ; h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
ELSEWHERE ; h_i(:,:,:) = 0._wp
END WHERE
WHERE( h_i(:,:,:) < rn_himin ) a_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:) / rn_himin
IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
WHERE( h_i(:,:,:) < rn_himin ) a_ip(:,:,:) = a_ip(:,:,:) * h_i(:,:,:) / rn_himin
ENDIF
WHERE( h_i(:,:,:) < rn_himin ) a_i(:,:,:) = a_i (:,:,:) * h_i(:,:,:) / rn_himin
!
! !-----------------------------------------------------
! ! ice concentration should not exceed amax !
......@@ -83,7 +86,7 @@ CONTAINS
! !-----------------------------------------------------
! ! Rebin categories with thickness out of bounds !
! !-----------------------------------------------------
IF ( jpl > 1 ) CALL ice_itd_reb( kt )
IF( jpl > 1 ) CALL ice_itd_reb( kt )
!
! !-----------------------------------------------------
IF ( nn_icesal == 2 ) THEN ! salinity must stay in bounds [Simin,Simax] !
......
src/ICE/icectl.F90
View file @ 66fb9523
......@@ -83,25 +83,33 @@ CONTAINS
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
!!
INTEGER :: ji, jj, jl ! dummy loop index
REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat
REAL(wp), DIMENSION(jpi,jpj,10) :: ztmp3
REAL(wp), DIMENSION(jpi,jpj,jpl,8) :: ztmp4
REAL(wp), DIMENSION(10) :: zchk3
REAL(wp), DIMENSION(8) :: zchk4
REAL(wp), DIMENSION(A2D(0),10) :: ztmp3
REAL(wp), DIMENSION(A2D(0),jpl,8) :: ztmp4
REAL(wp), DIMENSION(10) :: zchk3
REAL(wp), DIMENSION(8) :: zchk4
!!-------------------------------------------------------------------
!
! -- quantities -- !
ztmp3(:,:,1) = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ! volume
ztmp3(:,:,2) = SUM( sv_i * rhoi, dim=3 ) * e1e2t ! salt
ztmp3(:,:,3) = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ! heat
!
! -- fluxes -- !
ztmp3(:,:,4) = ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd & ! mass
& + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t
ztmp3(:,:,5) = ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw & ! salt
& + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t
ztmp3(:,:,6) = ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & ! heat
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t
DO_2D( 0, 0, 0, 0 )
! -- quantities -- !
ztmp3(ji,jj,1) = SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos + &
& ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow ) * e1e2t(ji,jj) ! volume
ztmp3(ji,jj,2) = SUM( sv_i(ji,jj,:) * rhoi ) * e1e2t(ji,jj) ! salt
ztmp3(ji,jj,3) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) + & ! heat
& SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) ) * e1e2t(ji,jj)
!
! -- fluxes -- !
ztmp3(ji,jj,4) = ( wfx_bog (ji,jj) + wfx_bom (ji,jj) + wfx_sum (ji,jj) + wfx_sni (ji,jj) & ! mass
& + wfx_opw (ji,jj) + wfx_res (ji,jj) + wfx_dyn (ji,jj) + wfx_lam (ji,jj) + wfx_pnd(ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sub(ji,jj) &
& + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj) ) * e1e2t(ji,jj)
ztmp3(ji,jj,5) = ( sfx_bri(ji,jj) + sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj) & ! salt
& + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) ) * e1e2t(ji,jj)
ztmp3(ji,jj,6) = ( hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) & ! heat
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj) ) * e1e2t(ji,jj)
!
END_2D
!
! -- global sum -- !
zchk3(1:6) = glob_sum_vec( 'icectl', ztmp3(:,:,1:6) )
......@@ -123,25 +131,33 @@ CONTAINS
zdiag_heat = ( zchk3(3) - pdiag_t ) * r1_Dt_ice + ( zchk3(6) - pdiag_ft )
! -- max concentration diag -- !
ztmp3(:,:,7) = SUM( a_i, dim=3 )
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) )
DO_2D( 0, 0, 0, 0 )
ztmp3(ji,jj,7) = SUM( a_i(ji,jj,:) )
END_2D
zchk3(7) = glob_max( 'icectl', ztmp3(:,:,7) )
! -- advection scheme is conservative? -- !
ztmp3(:,:,8 ) = diag_adv_mass * e1e2t
ztmp3(:,:,9 ) = diag_adv_heat * e1e2t
ztmp3(:,:,10) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
zchk3(8:10) = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) )
DO_2D( 0, 0, 0, 0 )
ztmp3(ji,jj,8 ) = diag_adv_mass(ji,jj) * e1e2t(ji,jj)
ztmp3(ji,jj,9 ) = diag_adv_heat(ji,jj) * e1e2t(ji,jj)
ztmp3(ji,jj,10) = SUM( a_i(ji,jj,:) + epsi10 ) * e1e2t(ji,jj) ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
END_2D
zchk3(8:10) = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) )
! -- min diags -- !
ztmp4(:,:,:,1) = v_i
ztmp4(:,:,:,2) = v_s
ztmp4(:,:,:,3) = v_ip
ztmp4(:,:,:,4) = v_il
ztmp4(:,:,:,5) = a_i
ztmp4(:,:,:,6) = sv_i
ztmp4(:,:,:,7) = SUM( e_i, dim=3 )
ztmp4(:,:,:,8) = SUM( e_s, dim=3 )
zchk4(1:8) = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) )
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
ztmp4(ji,jj,jl,1) = v_i(ji,jj,jl)
ztmp4(ji,jj,jl,2) = v_s(ji,jj,jl)
ztmp4(ji,jj,jl,3) = v_ip(ji,jj,jl)
ztmp4(ji,jj,jl,4) = v_il(ji,jj,jl)
ztmp4(ji,jj,jl,5) = a_i(ji,jj,jl)
ztmp4(ji,jj,jl,6) = sv_i(ji,jj,jl)
ztmp4(ji,jj,jl,7) = SUM( e_i(ji,jj,:,jl) )
ztmp4(ji,jj,jl,8) = SUM( e_s(ji,jj,:,jl) )
END_2D
ENDDO
zchk4(1:8) = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) )
IF( lwp ) THEN
! check conservation issues
......@@ -188,17 +204,21 @@ CONTAINS
!!-------------------------------------------------------------------
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
!!
REAL(wp), DIMENSION(jpi,jpj,4) :: ztmp
REAL(wp), DIMENSION(4) :: zchk
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),4) :: ztmp
REAL(wp), DIMENSION(4) :: zchk
!!-------------------------------------------------------------------
ztmp(:,:,1) = ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ! mass diag
ztmp(:,:,2) = ( sfx + diag_sice - diag_adv_salt ) * e1e2t ! salt
ztmp(:,:,3) = ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ! heat
! equivalent to this:
!! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
ztmp(:,:,4) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
DO_2D( 0, 0, 0, 0 )
!
ztmp(ji,jj,1) = ( wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_pnd (ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) &
& + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj) ) * e1e2t(ji,jj) ! mass diag
ztmp(ji,jj,2) = ( sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) ) * e1e2t(ji,jj) ! salt
ztmp(ji,jj,3) = ( qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) ) * e1e2t(ji,jj) ! heat
! equivalent to this:
!! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
ztmp(ji,jj,4) = SUM( a_i(ji,jj,:) + epsi10 ) * e1e2t(ji,jj) ! ice area (+epsi10 to set a threshold > 0 when there is no ice)
END_2D
! global sums
zchk(1:4) = glob_sum_vec( 'icectl', ztmp(:,:,1:4) )
......@@ -226,11 +246,11 @@ CONTAINS
!!-------------------------------------------------------------------
INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end
CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
REAL(wp) , DIMENSION(jpi,jpj), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
REAL(wp) , DIMENSION(A2D(0)), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
!!
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass, zdiag_salt, zdiag_heat, &
& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
INTEGER :: jl, jk
REAL(wp), DIMENSION(A2D(0)) :: zdiag_mass, zdiag_salt, zdiag_heat, &
& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
INTEGER :: ji, jj, jl, jk
LOGICAL :: ll_stop_m = .FALSE.
LOGICAL :: ll_stop_s = .FALSE.
LOGICAL :: ll_stop_t = .FALSE.
......@@ -239,62 +259,79 @@ CONTAINS
!
IF( icount == 0 ) THEN
pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 )
pdiag_s = SUM( sv_i * rhoi , dim=3 )
pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 )
DO_2D( 0, 0, 0, 0 )
pdiag_v(ji,jj) = SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos + ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow )
pdiag_s(ji,jj) = SUM( sv_i(ji,jj,:) * rhoi )
pdiag_t(ji,jj) = SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) + SUM( SUM( e_s(ji,jj,:,:), dim=2 ) )
! mass flux
pdiag_fv = 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(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) &
& + wfx_opw(ji,jj) + wfx_res(ji,jj) + wfx_dyn(ji,jj) + wfx_lam(ji,jj) + wfx_pnd (ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) &
& + wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj)
! salt flux
pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam
pdiag_fs(ji,jj) = sfx_bri(ji,jj) + sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) &
& + sfx_opw(ji,jj) + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
! heat flux
pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr
pdiag_ft(ji,jj) = hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) &
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj)
END_2D
ELSEIF( icount == 1 ) THEN
! -- mass diag -- !
zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice &
& + ( 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
DO_2D( 0, 0, 0, 0 )
zdiag_mass(ji,jj) = ( SUM( v_i(ji,jj,:) * rhoi + v_s(ji,jj,:) * rhos &
& + ( v_ip(ji,jj,:) + v_il(ji,jj,:) ) * rhow ) - pdiag_v(ji,jj) ) * r1_Dt_ice &
& + ( wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) + wfx_opw(ji,jj) &
& + wfx_res(ji,jj) + wfx_dyn(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) &
& + wfx_snw_sni(ji,jj) + wfx_snw_sum(ji,jj) + wfx_snw_dyn(ji,jj) &
& + wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj) + wfx_spr(ji,jj) ) &
& - pdiag_fv(ji,jj)
END_2D
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
DO_2D( 0, 0, 0, 0 )
zdiag_salt(ji,jj) = ( SUM( sv_i(ji,jj,:) * rhoi ) - pdiag_s(ji,jj) ) * r1_Dt_ice &
& + ( sfx_bri(ji,jj) + sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) &
& + sfx_opw(ji,jj) + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) ) &
& - pdiag_fs(ji,jj)
END_2D
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
DO_2D( 0, 0, 0, 0 )
zdiag_heat(ji,jj) = ( SUM( SUM( e_i(ji,jj,:,:), dim=2 ) ) &
& + SUM( SUM( e_s(ji,jj,:,:), dim=2 ) ) - pdiag_t(ji,jj) ) * r1_Dt_ice &
& + ( hfx_sum(ji,jj) + hfx_bom(ji,jj) + hfx_bog(ji,jj) &
& + hfx_dif(ji,jj) + hfx_opw(ji,jj) + hfx_snw(ji,jj) &
& - hfx_thd(ji,jj) - hfx_dyn(ji,jj) - hfx_res(ji,jj) - hfx_sub(ji,jj) - hfx_spr(ji,jj) ) &
& - pdiag_ft(ji,jj)
END_2D
IF( MAXVAL( ABS(zdiag_heat) ) > rchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
!
! -- other diags -- !
! a_i < 0
zdiag_amin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( a_i(:,:,jl) < 0._wp ) zdiag_amin(:,:) = 1._wp
WHERE( a_i(A2D(0),jl) < 0._wp ) zdiag_amin(:,:) = 1._wp
ENDDO
! v_i < 0
zdiag_vmin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( v_i(:,:,jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp
WHERE( v_i(A2D(0),jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp
ENDDO
! s_i < 0
zdiag_smin(:,:) = 0._wp
DO jl = 1, jpl
WHERE( s_i(:,:,jl) < 0._wp ) zdiag_smin(:,:) = 1._wp
WHERE( s_i(A2D(0),jl) < 0._wp ) zdiag_smin(:,:) = 1._wp
ENDDO
! e_i < 0
zdiag_emin(:,:) = 0._wp
DO jl = 1, jpl
DO jk = 1, nlay_i
WHERE( e_i(:,:,jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp
WHERE( e_i(A2D(0),jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp
ENDDO
ENDDO
! a_i > amax
......@@ -528,8 +565,8 @@ CONTAINS
INTEGER :: jl, ji, jj
!!-------------------------------------------------------------------
DO ji = mi0(ki), mi1(ki)
DO jj = mj0(kj), mj1(kj)
DO ji = mi0(ki,nn_hls), mi1(ki,nn_hls)
DO jj = mj0(kj,nn_hls), mj1(kj,nn_hls)
WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title
......@@ -733,10 +770,10 @@ CONTAINS
CALL prt_ctl_info(' ')
CALL prt_ctl_info(' - Stresses : ')
CALL prt_ctl_info(' ~~~~~~~~~~ ')
CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = umask, &
& tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = vmask)
CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = umask, &
& tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = vmask)
CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', mask1 = tmask, &
& tab2d_2=vtau , clinfo2= ' vtau : ', mask2 = tmask)
CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', mask1 = tmask, &
& tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ', mask2 = tmask)
END SUBROUTINE ice_prt3D
......@@ -751,8 +788,9 @@ CONTAINS
!!-------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ice time-step index
!
REAL(wp), DIMENSION(jpi,jpj,6) :: ztmp
REAL(wp), DIMENSION(6) :: zchk
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),6) :: ztmp
REAL(wp), DIMENSION(6) :: zchk
!!-------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) THEN
......@@ -762,25 +800,27 @@ CONTAINS
ENDIF
!
! -- 2D budgets (must be close to 0) -- !
ztmp(:,:,1) = wfx_ice (:,:) + wfx_snw (:,:) + wfx_spr (:,:) + wfx_sub(:,:) + wfx_pnd(:,:) &
& + diag_vice(:,:) + diag_vsnw(:,:) + diag_vpnd(:,:) - diag_adv_mass(:,:)
ztmp(:,:,2) = sfx(:,:) + diag_sice(:,:) - diag_adv_salt(:,:)
ztmp(:,:,3) = qt_oce_ai(:,:) - qt_atm_oi(:,:) + diag_heat(:,:) - diag_adv_heat(:,:)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = wfx_ice (ji,jj) + wfx_snw (ji,jj) + wfx_spr (ji,jj) + wfx_sub(ji,jj) + wfx_pnd(ji,jj) &
& + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj)
ztmp(ji,jj,2) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj)
ztmp(ji,jj,3) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj)
END_2D
! write outputs
CALL iom_put( 'icedrift_mass', ztmp(:,:,1) )
CALL iom_put( 'icedrift_salt', ztmp(:,:,2) )
CALL iom_put( 'icedrift_heat', ztmp(:,:,3) )
! -- 1D budgets -- !
ztmp(:,:,1) = ztmp(:,:,1) * e1e2t * rDt_ice ! mass
ztmp(:,:,2) = ztmp(:,:,2) * e1e2t * rDt_ice * 1.e-3 ! salt
ztmp(:,:,3) = ztmp(:,:,3) * e1e2t ! heat
ztmp(:,:,4) = diag_adv_mass * e1e2t * rDt_ice
ztmp(:,:,5) = diag_adv_salt * e1e2t * rDt_ice * 1.e-3
ztmp(:,:,6) = diag_adv_heat * e1e2t
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = ztmp(ji,jj,1) * e1e2t(ji,jj) * rDt_ice ! mass
ztmp(ji,jj,2) = ztmp(ji,jj,2) * e1e2t(ji,jj) * rDt_ice * 1.e-3 ! salt
ztmp(ji,jj,3) = ztmp(ji,jj,3) * e1e2t(ji,jj) ! heat
ztmp(ji,jj,4) = diag_adv_mass(ji,jj) * e1e2t(ji,jj) * rDt_ice
ztmp(ji,jj,5) = diag_adv_salt(ji,jj) * e1e2t(ji,jj) * rDt_ice * 1.e-3
ztmp(ji,jj,6) = diag_adv_heat(ji,jj) * e1e2t(ji,jj)
END_2D
! global sums
zchk(1:6) = glob_sum_vec( 'icectl', ztmp(:,:,1:6) )
......
src/ICE/icedia.F90
View file @ 66fb9523
......@@ -33,6 +33,9 @@ MODULE icedia
PUBLIC ice_dia ! called by icestp.F90
PUBLIC ice_dia_init ! called in icestp.F90
!! * Substitutions
# include "do_loop_substitute.h90"
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
......@@ -48,7 +51,7 @@ CONTAINS
!!---------------------------------------------------------------------!
!! *** ROUTINE ice_dia_alloc ***
!!---------------------------------------------------------------------!
ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
ALLOCATE( vol_loc_ini(A2D(0)), sal_loc_ini(A2D(0)), tem_loc_ini(A2D(0)), STAT=ice_dia_alloc )
CALL mpp_sum ( 'icedia', ice_dia_alloc )
IF( ice_dia_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dia_alloc: failed to allocate arrays' )
......@@ -64,8 +67,9 @@ CONTAINS
!!---------------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
!!
REAL(wp), DIMENSION(jpi,jpj,16) :: ztmp
REAL(wp), DIMENSION(16) :: zbg
INTEGER :: ji, jj ! dummy loop index
REAL(wp), DIMENSION(A2D(0),16) :: ztmp
REAL(wp), DIMENSION(16) :: zbg
!!---------------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('ice_dia')
......@@ -85,31 +89,32 @@ CONTAINS
! 1 - Trends due to forcing !
! ---------------------------!
! they must be kept outside an IF(iom_use) because of the call to dia_rst below
ztmp(:,:,1) = - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ! freshwater flux ice/snow-ocean
ztmp(:,:,2) = - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ! freshwater flux ice/snow-atm
ztmp(:,:,3) = - sfx (:,:) * e1e2t(:,:) ! salt fluxes ice/snow-ocean
ztmp(:,:,4) = qt_atm_oi(:,:) * e1e2t(:,:) ! heat on top of ice-ocean
ztmp(:,:,5) = qt_oce_ai(:,:) * e1e2t(:,:) ! heat on top of ocean (and below ice)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) * e1e2t(ji,jj) ! freshwater flux ice/snow-ocean
ztmp(ji,jj,2) = - ( wfx_sub(ji,jj) + wfx_spr(ji,jj) ) * e1e2t(ji,jj) ! freshwater flux ice/snow-atm
ztmp(ji,jj,3) = - sfx (ji,jj) * e1e2t(ji,jj) ! salt fluxes ice/snow-ocean
ztmp(ji,jj,4) = qt_atm_oi(ji,jj) * e1e2t(ji,jj) ! heat on top of ice-ocean
ztmp(ji,jj,5) = qt_oce_ai(ji,jj) * e1e2t(ji,jj) ! heat on top of ocean (and below ice)
END_2D
! ----------------------- !
! 2 - Contents !
! ----------------------- !
IF( iom_use('ibgvol_tot' ) ) ztmp(:,:,6 ) = vt_i (:,:) * e1e2t(:,:) ! ice volume
IF( iom_use('sbgvol_tot' ) ) ztmp(:,:,7 ) = vt_s (:,:) * e1e2t(:,:) ! snow volume
IF( iom_use('ibgarea_tot') ) ztmp(:,:,8 ) = at_i (:,:) * e1e2t(:,:) ! area
IF( iom_use('ibgsalt_tot') ) ztmp(:,:,9 ) = st_i (:,:) * e1e2t(:,:) ! salt content
IF( iom_use('ibgheat_tot') ) ztmp(:,:,10) = et_i (:,:) * e1e2t(:,:) ! heat content
IF( iom_use('sbgheat_tot') ) ztmp(:,:,11) = et_s (:,:) * e1e2t(:,:) ! heat content
IF( iom_use('ipbgvol_tot') ) ztmp(:,:,12) = vt_ip(:,:) * e1e2t(:,:) ! ice pond volume
IF( iom_use('ilbgvol_tot') ) ztmp(:,:,13) = vt_il(:,:) * e1e2t(:,:) ! ice pond lid volume
IF( iom_use('ibgvol_tot' ) ) ztmp(:,:,6 ) = vt_i (A2D(0)) * e1e2t(A2D(0)) ! ice volume
IF( iom_use('sbgvol_tot' ) ) ztmp(:,:,7 ) = vt_s (A2D(0)) * e1e2t(A2D(0)) ! snow volume
IF( iom_use('ibgarea_tot') ) ztmp(:,:,8 ) = at_i (A2D(0)) * e1e2t(A2D(0)) ! area
IF( iom_use('ibgsalt_tot') ) ztmp(:,:,9 ) = st_i (:,:) * e1e2t(A2D(0)) ! salt content
IF( iom_use('ibgheat_tot') ) ztmp(:,:,10) = et_i (:,:) * e1e2t(A2D(0)) ! heat content
IF( iom_use('sbgheat_tot') ) ztmp(:,:,11) = et_s (:,:) * e1e2t(A2D(0)) ! heat content
IF( iom_use('ipbgvol_tot') ) ztmp(:,:,12) = vt_ip(A2D(0)) * e1e2t(A2D(0)) ! ice pond volume
IF( iom_use('ilbgvol_tot') ) ztmp(:,:,13) = vt_il(A2D(0)) * e1e2t(A2D(0)) ! ice pond lid volume
! ---------------------------------- !
! 3 - Content variations and drifts !
! ---------------------------------- !
IF( iom_use('ibgvolume') ) ztmp(:,:,14) = ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ! freshwater trend
IF( iom_use('ibgsaltco') ) ztmp(:,:,15) = ( rhoi*st_i(:,:) - sal_loc_ini(:,:) ) * e1e2t(:,:) ! salt content trend
IF( iom_use('ibgvolume') ) ztmp(:,:,14) = ( rhoi*vt_i(A2D(0)) + rhos*vt_s(A2D(0)) - vol_loc_ini(:,:) ) * e1e2t(A2D(0)) ! freshwater trend
IF( iom_use('ibgsaltco') ) ztmp(:,:,15) = ( rhoi*st_i(:,:) - sal_loc_ini(:,:) ) * e1e2t(A2D(0)) ! salt content trend
IF( iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) &
& ztmp(:,:,16) = ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(:,:) ! heat content trend
& ztmp(:,:,16) = ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(A2D(0)) ! heat content trend
! global sum
zbg(1:16) = glob_sum_vec( 'icedia', ztmp(:,:,1:16) )
......@@ -261,9 +266,9 @@ CONTAINS
frc_tembot = 0._wp
frc_sal = 0._wp
! record initial ice volume, salt and temp
vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:) ! ice/snow volume (kg/m2)
tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J)
sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2)
vol_loc_ini(:,:) = rhoi * vt_i(A2D(0)) + rhos * vt_s(A2D(0)) ! ice/snow volume (kg/m2)
tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J)
sal_loc_ini(:,:) = rhoi * st_i(:,:) ! ice salt content (pss*kg/m2)
ENDIF
!
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
......
src/ICE/icedyn.F90
View file @ 66fb9523
......@@ -93,8 +93,8 @@ CONTAINS
!
! retrieve thickness from volume for landfast param. and UMx advection scheme
WHERE( a_i(:,:,:) >= epsi20 )
h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:)
ELSEWHERE
h_i(:,:,:) = 0._wp
h_s(:,:,:) = 0._wp
......@@ -162,15 +162,12 @@ CONTAINS
CASE ( np_dynADV1D , np_dynADV2D )
ALLOCATE( zdivu_i(jpi,jpj) )
ALLOCATE( zdivu_i(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
zdivu_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)
zdivu_i(ji,jj) = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) & ! add () for NP repro
& + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) ) * r1_e1e2t(ji,jj)
END_2D
CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.0_wp )
! output
CALL iom_put( 'icediv' , zdivu_i )
DEALLOCATE( zdivu_i )
END SELECT
......