From 40d24f1163ddd085c6e0cd61d202ca549413e761 Mon Sep 17 00:00:00 2001
From: Clement Rousset <clement.rousset@locean.ipsl.fr>
Date: Sat, 1 Oct 2022 19:19:26 +0000
Subject: [PATCH] Resolve "Summer 2022 work- ISF halo cleanup"
---
cfgs/SHARED/field_def_nemo-oce.xml | 58 ++--
src/OCE/ASM/asminc.F90 | 6 +-
src/OCE/DOM/closea.F90 | 16 +-
src/OCE/DYN/divhor.F90 | 17 +-
src/OCE/DYN/sshwzv.F90 | 2 +
src/OCE/ISF/isf_oce.F90 | 171 ++++--------
src/OCE/ISF/isfcav.F90 | 171 +++++++-----
src/OCE/ISF/isfcavgam.F90 | 188 +++++--------
src/OCE/ISF/isfcavmlt.F90 | 181 ++++++-------
src/OCE/ISF/isfcpl.F90 | 407 ++++++++++++++++-------------
src/OCE/ISF/isfdiags.F90 | 51 ++--
src/OCE/ISF/isfdynatf.F90 | 73 +++---
src/OCE/ISF/isfhdiv.F90 | 68 ++---
src/OCE/ISF/isfload.F90 | 4 +-
src/OCE/ISF/isfpar.F90 | 94 +++----
src/OCE/ISF/isfparmlt.F90 | 167 ++++++------
src/OCE/ISF/isfrst.F90 | 27 +-
src/OCE/ISF/isfstp.F90 | 136 ++++++----
src/OCE/ISF/isftbl.F90 | 215 ++++++---------
src/OCE/ISF/isfutils.F90 | 53 ++--
src/OCE/SBC/sbc_oce.F90 | 5 +-
src/OCE/SBC/sbccpl.F90 | 12 +-
src/OCE/SBC/sbcmod.F90 | 6 +-
src/OCE/SBC/sbcrnf.F90 | 52 ++--
src/OCE/TRA/eosbn2.F90 | 230 +++++++++++++---
src/OCE/TRA/traisf.F90 | 59 ++---
src/OCE/step.F90 | 4 +-
src/OCE/stpmlf.F90 | 2 +
src/OCE/stprk3.F90 | 2 +
src/OFF/dtadyn.F90 | 32 +--
tests/ISOMIP+/MY_SRC/isf_oce.F90 | 173 +++++-------
31 files changed, 1365 insertions(+), 1317 deletions(-)
diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index bc6b20728..6cd9bdb8e 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -241,39 +241,39 @@ that are available in the tidal-forcing implementation (see
<!-- * variable related to ice shelf forcing * -->
<!-- * fwf * -->
- <field id="fwfisf_cav" long_name="Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" />
- <field id="fwfisf_par" long_name="Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" />
- <field id="fwfisf3d_cav" long_name="3d Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_3D" />
- <field id="fwfisf3d_par" long_name="3d Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_3D" />
+ <field id="fwfisf_cav" long_name="Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_2D_inner" />
+ <field id="fwfisf_par" long_name="Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_2D_inner" />
+ <field id="fwfisf3d_cav" long_name="3d Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_3D_inner" />
+ <field id="fwfisf3d_par" long_name="3d Ice shelf fresh water flux ( from isf to oce )" unit="kg/m2/s" grid_ref="grid_T_3D_inner" />
<!-- * heat fluxes * -->
- <field id="qoceisf_cav" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" />
- <field id="qoceisf_par" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" />
- <field id="qlatisf_cav" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" />
- <field id="qlatisf_par" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" />
- <field id="qhcisf_cav" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" />
- <field id="qhcisf_par" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" />
- <field id="qoceisf3d_cav" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qoceisf3d_par" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qlatisf3d_cav" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qlatisf3d_par" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qhcisf3d_cav" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qhcisf3d_par" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D" />
- <field id="qconisf" long_name="Conductive heat flux through the ice shelf ( from isf to oce )" unit="W/m2" />
+ <field id="qoceisf_cav" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qoceisf_par" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qlatisf_cav" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qlatisf_par" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qhcisf_cav" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qhcisf_par" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
+ <field id="qoceisf3d_cav" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qoceisf3d_par" long_name="Ice shelf ocean heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qlatisf3d_cav" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qlatisf3d_par" long_name="Ice shelf latent heat flux ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qhcisf3d_cav" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qhcisf3d_par" long_name="Ice shelf heat content flux of injected water ( from isf to oce )" unit="W/m2" grid_ref="grid_T_3D_inner" />
+ <field id="qconisf" long_name="Conductive heat flux through the ice shelf ( from isf to oce )" unit="W/m2" grid_ref="grid_T_2D_inner" />
<!-- top boundary layer properties -->
- <field id="isftfrz_cav" long_name="freezing point temperature at ocean/isf interface" unit="degC" />
- <field id="isftfrz_par" long_name="freezing point temperature in the parametrization boundary layer" unit="degC" />
- <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting" unit="degC" />
- <field id="isfthermald_par" long_name="thermal driving of ice shelf melting" unit="degC" />
- <field id="isfgammat" long_name="Ice shelf heat-transfert velocity" unit="m/s" />
- <field id="isfgammas" long_name="Ice shelf salt-transfert velocity" unit="m/s" />
- <field id="ttbl_cav" long_name="temperature in Losch tbl" unit="degC" />
- <field id="ttbl_par" long_name="temperature in the parametrisation boundary layer" unit="degC" />
- <field id="stbl" long_name="salinity in the Losh tbl" unit="1e-3" />
- <field id="utbl" long_name="zonal current in the Losh tbl at T point" unit="m/s" />
- <field id="vtbl" long_name="merid current in the Losh tbl at T point" unit="m/s" />
- <field id="isfustar" long_name="ustar at T point used in ice shelf melting" unit="m/s" />
+ <field id="isftfrz_cav" long_name="freezing point temperature at ocean/isf interface" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="isftfrz_par" long_name="freezing point temperature in the parametrization boundary layer" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="isfthermald_par" long_name="thermal driving of ice shelf melting" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="isfgammat" long_name="Ice shelf heat-transfert velocity" unit="m/s" grid_ref="grid_T_2D_inner" />
+ <field id="isfgammas" long_name="Ice shelf salt-transfert velocity" unit="m/s" grid_ref="grid_T_2D_inner" />
+ <field id="ttbl_cav" long_name="temperature in Losch tbl" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="ttbl_par" long_name="temperature in the parametrisation boundary layer" unit="degC" grid_ref="grid_T_2D_inner" />
+ <field id="stbl" long_name="salinity in the Losh tbl" unit="1e-3" grid_ref="grid_T_2D_inner" />
+ <field id="utbl" long_name="zonal current in the Losh tbl at T point" unit="m/s" grid_ref="grid_T_2D_inner" />
+ <field id="vtbl" long_name="merid current in the Losh tbl at T point" unit="m/s" grid_ref="grid_T_2D_inner" />
+ <field id="isfustar" long_name="ustar at T point used in ice shelf melting" unit="m/s" grid_ref="grid_T_2D_inner" />
</field_group> <!-- grid_T -->
diff --git a/src/OCE/ASM/asminc.F90 b/src/OCE/ASM/asminc.F90
index 54371a62c..82f63a2f1 100644
--- a/src/OCE/ASM/asminc.F90
+++ b/src/OCE/ASM/asminc.F90
@@ -821,12 +821,12 @@ CONTAINS
CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments)
!
IF( ln_linssh ) THEN
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+ DO_2D_OVR( 0, 0, 0, 0 )
phdivn(ji,jj,1) = phdivn(ji,jj,1) - ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) * tmask(ji,jj,1)
END_2D
ELSE
- ALLOCATE( ztim(A2D(nn_hls)) )
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+ ALLOCATE( ztim(A2D(0)) )
+ DO_2D_OVR( 0, 0, 0, 0 )
ztim(ji,jj) = ssh_iau(ji,jj) / ( ht(ji,jj) + 1.0 - ssmask(ji,jj) )
DO jk = 1, jpkm1
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ztim(ji,jj) * tmask(ji,jj,jk)
diff --git a/src/OCE/DOM/closea.F90 b/src/OCE/DOM/closea.F90
index eb523c63c..3505bb0f7 100644
--- a/src/OCE/DOM/closea.F90
+++ b/src/OCE/DOM/closea.F90
@@ -51,6 +51,8 @@ MODULE closea
INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csrnf , mask_csgrprnf !: mask of integers defining closed seas rnf mappings
INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csemp , mask_csgrpemp !: mask of integers defining closed seas empmr mappings
+ !! * Substitutions
+# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: closea.F90 13558 2020-10-02 15:30:22Z smasson $
@@ -173,17 +175,17 @@ CONTAINS
!! ** Action : update (p_)mskrnf (set 1 at closed sea outflow)
!!----------------------------------------------------------------------
!! subroutine parameter
- REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: p_rnfmsk ! river runoff mask (rnfmsk array)
- !!
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: p_rnfmsk ! river runoff mask (rnfmsk array)
!! local variables
- REAL(wp), DIMENSION(jpi,jpj) :: zmsk
+ INTEGER :: ji, jj
+ INTEGER :: zmsk
!!----------------------------------------------------------------------
!
! zmsk > 0 where cs river mouth defined (case rnf and emp)
- zmsk(:,:) = ( mask_csgrprnf (:,:) + mask_csgrpemp(:,:) ) * mask_opnsea(:,:)
- WHERE( zmsk(:,:) > 0 )
- p_rnfmsk(:,:) = 1.0_wp
- END WHERE
+ DO_2D( 0, 0, 0, 0 )
+ zmsk = ( mask_csgrprnf(ji,jj) + mask_csgrpemp(ji,jj) ) * mask_opnsea(ji,jj)
+ IF( zmsk > 0 ) p_rnfmsk(ji,jj) = 1.0_wp
+ END_2D
!
END SUBROUTINE clo_rnf
diff --git a/src/OCE/DYN/divhor.F90 b/src/OCE/DYN/divhor.F90
index 2cdf53d18..795277715 100644
--- a/src/OCE/DYN/divhor.F90
+++ b/src/OCE/DYN/divhor.F90
@@ -78,12 +78,11 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence '
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
- hdiv (:,:,:) = 0._wp ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step
+ hdiv(:,:,:) = 0._wp ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step
+ pe3divUh(:,:,jpk) = 0._wp
ENDIF
!
- pe3divUh(:,:,:) = 0._wp !!gm to be applied to the halos only
- !
- DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
+ DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 )
hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * puu(ji ,jj,jk) &
& - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk) &
& + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * pvv(ji,jj ,jk) &
@@ -100,7 +99,7 @@ CONTAINS
!
IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== + ice-shelf mass exchange ==!
!
- IF( nn_hls==1 ) CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp ) ! (no sign change)
+ CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp ) ! (no sign change)
!
!!gm Patch before suppression of hdiv from all modules that use it
! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== e3t * Horizontal divergence ==!
@@ -149,7 +148,7 @@ CONTAINS
END_3D
ENDIF
!
- DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 ) !== Horizontal divergence ==!
+ DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 ) !== Horizontal divergence ==!
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
hdiv(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) &
@@ -167,10 +166,10 @@ CONTAINS
IF( ln_sshinc .AND. ln_asmiau ) CALL ssh_asm_div( kt, Kbb, Kmm, hdiv ) !== SSH assimilation ==! (update hdiv field)
!
#endif
- IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== ice shelf ==! (update hdiv field)
+ IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== ice shelf ==! (update hdiv field)
!
- IF( nn_hls==1 ) CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp ) ! (no sign change)
- ! ! needed for ww in sshwzv
+ CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp ) ! needed for ww in sshwzv (no sign change)
+ !
IF( ln_timing ) CALL timing_stop('div_hor')
!
END SUBROUTINE div_hor_old
diff --git a/src/OCE/DYN/sshwzv.F90 b/src/OCE/DYN/sshwzv.F90
index 71452dc72..52027ff1b 100644
--- a/src/OCE/DYN/sshwzv.F90
+++ b/src/OCE/DYN/sshwzv.F90
@@ -198,6 +198,7 @@ CONTAINS
! ! Is it problematic to have a wrong vertical velocity in boundary cells?
! ! Same question holds for hdiv. Perhaps just for security
! ! clem: yes it is a problem because ww is used in many other places where we need the halos
+ ! ! clem: I do not understand this lbc. Not needed to me
!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
! computation of w
@@ -339,6 +340,7 @@ CONTAINS
IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
! ! Is it problematic to have a wrong vertical velocity in boundary cells?
! ! Same question holds for hdiv. Perhaps just for security
+ ! ! clem: I do not understand this lbc. Not needed to me
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( ze3div(ji,jj,jk) + zhdiv(ji,jj,jk) &
& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
diff --git a/src/OCE/ISF/isf_oce.F90 b/src/OCE/ISF/isf_oce.F90
index eb62d83e9..0ac9e8a5e 100644
--- a/src/OCE/ISF/isf_oce.F90
+++ b/src/OCE/ISF/isf_oce.F90
@@ -13,8 +13,7 @@ MODULE isf_oce
!! isf : define and allocate ice shelf variables
!!----------------------------------------------------------------------
- USE par_oce , ONLY: jpi, jpj, jpk
- USE in_out_manager, ONLY: wp, jpts ! I/O manager
+ USE par_oce
USE lib_mpp , ONLY: ctl_stop, mpp_sum ! MPP library
USE fldread ! read input fields
@@ -22,7 +21,7 @@ MODULE isf_oce
PRIVATE
- PUBLIC isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl, isf_dealloc_cpl
+ PUBLIC isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl
!
!-------------------------------------------------------
! 0 : namelist parameter
@@ -79,33 +78,35 @@ MODULE isf_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_oasis
!
! 2.2 -------- ice shelf cavity melt namelist parameter -------------
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
- TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
+ TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
!
- REAL(wp) , PUBLIC :: risf_lamb1, risf_lamb2, risf_lamb3 ! freezing point linearization coeficient
+ REAL(wp) , PUBLIC :: risf_lamb1, risf_lamb2, risf_lamb3 !: freezing point linearization coeficient
!
! 2.3 -------- ice shelf param. melt namelist parameter -------------
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_par !:
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_par , misfkb_par !:
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_par, rfrac_tbl_par !:
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_par , fwfisf_par_b !: before and now net fwf from the ice shelf [kg/m2/s]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
- TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfpar_fwf !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_par !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_par , misfkb_par !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_par, rfrac_tbl_par !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_par , fwfisf_par_b !: before and now net fwf from the ice shelf [kg/m2/s]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
+ TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfpar_fwf !:
!
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf0_tbl_par !: thickness of tbl (initial value) [m]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf0_tbl_par !: thickness of tbl (initial value) [m]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !:
!
! 2.4 -------- coupling namelist parameter -------------
INTEGER , PUBLIC :: nstp_iscpl !:
REAL(wp), PUBLIC :: rdt_iscpl !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfcpl_ssh, risfcpl_cons_ssh, risfcpl_cons_ssh_b !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risfcpl_vol, risfcpl_cons_vol, risfcpl_cons_vol_b !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: risfcpl_tsc, risfcpl_cons_tsc, risfcpl_cons_tsc_b !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfcpl_ssh, risfcpl_cons_ssh !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risfcpl_vol, risfcpl_cons_vol !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: risfcpl_tsc, risfcpl_cons_tsc !:
!
+ !! * Substitutions
+# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
@@ -125,20 +126,12 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0 ! set to zero if no array to be allocated
- !
- ALLOCATE(risfLeff(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE(misfkt_par(jpi,jpj), misfkb_par(jpi,jpj), STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( rfrac_tbl_par(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( rhisf_tbl_par(jpi,jpj), rhisf0_tbl_par(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( mskisf_par(jpi,jpj), STAT=ialloc)
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( misfkt_par (A2D(0)) , misfkb_par (A2D(0)) , rfrac_tbl_par(A2D(0)) , &
+ & rhisf_tbl_par(A2D(0)) , rhisf0_tbl_par(A2D(0)) , &
+ & risfLeff (A2D(0)) , mskisf_par (A2D(0)) , STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
@@ -159,14 +152,10 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0 ! set to zero if no array to be allocated
- !
- ALLOCATE(misfkt_cav(jpi,jpj), misfkb_cav(jpi,jpj), STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( rfrac_tbl_cav(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( rhisf_tbl_cav(jpi,jpj), STAT=ialloc)
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( misfkt_cav(A2D(0)), misfkb_cav(A2D(0)), rfrac_tbl_cav(A2D(0)) , rhisf_tbl_cav(A2D(0)) , STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
@@ -185,46 +174,25 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0
- !
- ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts) , risfcpl_vol(jpi,jpj,jpk) , STAT=ialloc )
+ ! ----------------- !
+ ! == FULL ARRAYS == !
+ ! ----------------- !
+ ALLOCATE( risfcpl_ssh (jpi,jpj) , risfcpl_vol (jpi,jpj,jpk) , &
+ & risfcpl_cons_ssh(jpi,jpj) , risfcpl_cons_vol(jpi,jpj,jpk) , STAT=ialloc )
+ ierr = ierr + ialloc
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( risfcpl_tsc (A2D(0),jpk,jpts) , &
+ & risfcpl_cons_tsc(A2D(0),jpk,jpts) , STAT=ialloc )
ierr = ierr + ialloc
- !
- risfcpl_tsc(:,:,:,:) = 0._wp ; risfcpl_vol(:,:,:) = 0._wp ; risfcpl_ssh(:,:) = 0._wp
-
- IF ( ln_isfcpl_cons ) THEN
- ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- risfcpl_cons_tsc(:,:,:,:) = 0._wp ; risfcpl_cons_vol(:,:,:) = 0._wp ; risfcpl_cons_ssh(:,:) = 0._wp
- !
- END IF
!
CALL mpp_sum ( 'isf', ierr )
IF( ierr /= 0 ) CALL ctl_stop('STOP','isfcpl: failed to allocate arrays.')
!
END SUBROUTINE isf_alloc_cpl
-
- SUBROUTINE isf_dealloc_cpl()
- !!---------------------------------------------------------------------
- !! *** ROUTINE isf_dealloc_cpl ***
- !!
- !! ** Purpose : de-allocate useless public 3d array used for ice sheet coupling
- !!
- !!----------------------------------------------------------------------
- INTEGER :: ierr, ialloc
- !!----------------------------------------------------------------------
- ierr = 0
- !
- DEALLOCATE( risfcpl_ssh , risfcpl_tsc , risfcpl_vol , STAT=ialloc )
- ierr = ierr + ialloc
- !
- CALL mpp_sum ( 'isf', ierr )
- IF( ierr /= 0 ) CALL ctl_stop('STOP','isfcpl: failed to deallocate arrays.')
- !
- END SUBROUTINE isf_dealloc_cpl
-
-
+
SUBROUTINE isf_alloc()
!!---------------------------------------------------------------------
!! *** ROUTINE isf_alloc ***
@@ -237,52 +205,29 @@ CONTAINS
!
ierr = 0 ! set to zero if no array to be allocated
!
- ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , &
- & fwfisf_oasis(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
+ ! ----------------- !
+ ! == FULL ARRAYS == !
+ ! ----------------- !
+ ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , risfload(jpi,jpj) , &
#if ! defined key_RK3
- ! MLF : need to allocate before arrays
- ALLOCATE( fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
+ & fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , & ! MLF : need to allocate before arrays
#endif
- !
- ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
+ & STAT=ialloc )
ierr = ierr + ialloc
!
- ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc )
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( fwfisf_oasis(A2D(0)) , risf_par_tsc (A2D(0),jpts) , risf_cav_tsc (A2D(0),jpts) , mskisf_cav(A2D(0)) , &
+#if ! defined key_RK3
+ & risf_par_tsc_b(A2D(0),jpts) , risf_cav_tsc_b(A2D(0),jpts) , & ! MLF : need to allocate before arrays
+#endif
+ & STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
!
- ! initalisation of fwf and tsc array to 0
- risfload (:,:) = 0._wp
- fwfisf_oasis(:,:) = 0._wp
-#if defined key_RK3
- fwfisf_par (:,:) = 0._wp
- fwfisf_cav (:,:) = 0._wp
- risf_cav_tsc(:,:,:) = 0._wp
- risf_par_tsc(:,:,:) = 0._wp
-#else
- fwfisf_par (:,:) = 0._wp ; fwfisf_par_b (:,:) = 0._wp
- fwfisf_cav (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp
- risf_cav_tsc(:,:,:) = 0._wp ; risf_cav_tsc_b(:,:,:) = 0._wp
- risf_par_tsc(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp
-#endif
- !
END SUBROUTINE isf_alloc
!!======================================================================
diff --git a/src/OCE/ISF/isfcav.F90 b/src/OCE/ISF/isfcav.F90
index 27ca9912e..790324c4e 100644
--- a/src/OCE/ISF/isfcav.F90
+++ b/src/OCE/ISF/isfcav.F90
@@ -14,16 +14,17 @@ MODULE isfcav
!!----------------------------------------------------------------------
USE isf_oce ! ice shelf public variables
!
- USE isfrst , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
+ USE isfrst , ONLY: isfrst_read ! ice shelf restart read/write subroutine
USE isfutils , ONLY: debug ! ice shelf debug subroutine
USE isftbl , ONLY: isf_tbl ! ice shelf top boundary layer properties subroutine
USE isfcavmlt, ONLY: isfcav_mlt ! ice shelf melt formulation subroutine
USE isfcavgam, ONLY: isfcav_gammats ! ice shelf melt exchange coeficient subroutine
USE isfdiags , ONLY: isf_diags_flx ! ice shelf diags subroutine
+ USE zdfdrg , ONLY: rCd0_top, r_ke0_top ! vertical physics: top/bottom drag coef.
!
USE oce , ONLY: ts, uu, vv, rn2 ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
- USE par_oce , ONLY: jpi,jpj ! ocean space and time domain
+ USE par_oce ! ocean space and time domain
USE phycst , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp ! physical constants
USE eosbn2 , ONLY: ln_teos10 ! use ln_teos10 or not
!
@@ -49,7 +50,7 @@ MODULE isfcav
!!----------------------------------------------------------------------
CONTAINS
- SUBROUTINE isf_cav( kt, Kmm, ptsc, pqfwf )
+ SUBROUTINE isf_cav( kt, Kmm, ptsc, pfwf )
!!---------------------------------------------------------------------
!! *** ROUTINE isf_cav ***
!!
@@ -67,37 +68,76 @@ CONTAINS
!! ** Convention : all fluxes are from isf to oce
!!
!!---------------------------------------------------------------------
- !!-------------------------- OUT --------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pqfwf ! ice shelf fwf
- REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc ! T & S ice shelf cavity contents
- !!-------------------------- IN --------------------------------------
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
INTEGER, INTENT(in) :: kt ! ocean time step
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0),jpts), INTENT(inout) :: ptsc ! T & S ice shelf cavity contents
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(inout) :: pfwf ! ice shelf fwf
!!---------------------------------------------------------------------
- LOGICAL :: lit
- INTEGER :: nit, ji, jj, ikt
- REAL(wp) :: zerr
- REAL(wp) :: zcoef, zdku, zdkv
- REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh ! heat fluxes
- REAL(wp), DIMENSION(jpi,jpj) :: zqh_b, zRc !
- REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas ! exchange coeficient
- REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl ! temp. and sal. in top boundary layer
+ LOGICAL :: lit
+ INTEGER :: nit, ji, jj, ikt
+ REAL(wp) :: zerr
+ REAL(wp) :: zcoef, zdku, zdkv
+ REAL(wp), DIMENSION(A2D(0)) :: zqlat, zqoce, zqhc, zqh ! heat fluxes
+ REAL(wp), DIMENSION(A2D(0)) :: zqh_b, zRc !
+ REAL(wp), DIMENSION(A2D(0)) :: zgammat, zgammas ! exchange coeficient
+ REAL(wp), DIMENSION(A2D(0)) :: zttbl, zstbl ! temp. and sal. in top boundary layer
+ REAL(wp), DIMENSION(A2D(0)) :: zustar ! u*
+ REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity
!!---------------------------------------------------------------------
!
- ! compute T/S/U/V for the top boundary layer
- CALL isf_tbl(Kmm, ts(:,:,:,jp_tem,Kmm), zttbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
- CALL isf_tbl(Kmm, ts(:,:,:,jp_sal,Kmm), zstbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
+ !===============================
+ ! 1.: compute T and S in the tbl
+ !===============================
+ CALL isf_tbl( Kmm, ts(A2D(0),:,jp_tem,Kmm), 'T', misfkt_cav, rhisf_tbl_cav, & ! <<== in
+ & zttbl, & ! ==>> out
+ & misfkb_cav, rfrac_tbl_cav ) ! <<== in (optional)
+ !
+ CALL isf_tbl( Kmm, ts(A2D(0),:,jp_sal,Kmm), 'T', misfkt_cav, rhisf_tbl_cav, & ! <<== in
+ & zstbl, & ! ==>> out
+ & misfkb_cav, rfrac_tbl_cav ) ! <<== in (optional)
+ !
+ ! output T/S for the top boundary layer
+ CALL iom_put( 'ttbl_cav', zttbl(:,:) * mskisf_cav(:,:) )
+ CALL iom_put( 'stbl' , zstbl(:,:) * mskisf_cav(:,:) )
!
- ! output T/S/U/V for the top boundary layer
- CALL iom_put('ttbl_cav',zttbl(:,:) * mskisf_cav(:,:))
- CALL iom_put('stbl' ,zstbl(:,:) * mskisf_cav(:,:))
+ !==========================================
+ ! 2.: compute velocity in the tbl if needed
+ !==========================================
+ IF ( TRIM(cn_gammablk) == 'vel_stab' .OR. TRIM(cn_gammablk) == 'vel' ) THEN
+ ! compute velocity in tbl
+ CALL isf_tbl( Kmm, uu(A2D(0),:,Kmm), 'U', miku(A2D(0)), rhisf_tbl_cav, & ! <<== in
+ & zutbl(A2D(0)) ) ! ==>> out
+ !
+ CALL isf_tbl( Kmm, vv(A2D(0),:,Kmm), 'V', mikv(A2D(0)), rhisf_tbl_cav, & ! <<== in
+ & zvtbl(A2D(0)) ) ! ==>> out
+ !
+ ! compute ustar (AD15 eq. 27)
+ CALL lbc_lnk( 'isfcav', zutbl, 'U', -1._wp, zvtbl, 'V', -1._wp )
+ DO_2D( 0, 0, 0, 0 )
+ zustar(ji,jj) = SQRT( rCd0_top(ji,jj) * ( 0.25_wp * ( (zutbl(ji,jj)+zutbl(ji-1,jj))*(zutbl(ji,jj)+zutbl(ji-1,jj)) &
+ & + (zvtbl(ji,jj)+zvtbl(ji,jj-1))*(zvtbl(ji,jj)+zvtbl(ji,jj-1)) ) &
+ & + r_ke0_top ) ) * mskisf_cav(ji,jj)
+ END_2D
+ !
+ ! output
+ CALL iom_put( 'isfustar', zustar(:,:) )
+ CALL iom_put( 'utbl' , zutbl(A2D(0)) )
+ CALL iom_put( 'vtbl' , zvtbl(A2D(0)) )
+ !
+ ELSE
+ zustar(:,:) = 0._wp ! not used
+ ENDIF
!
- ! initialisation
+ !==============================
+ ! 3.: compute ice shelf melting
+ !==============================
+ ! --- initialisation --- !
IF ( TRIM(cn_gammablk) == 'vel_stab' ) THEN
- zqoce(:,:) = -pqfwf(:,:) * rLfusisf !
- zqh_b(:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
-
+ !
DO_2D( 0, 0, 0, 0 )
+ zqoce(ji,jj) = -pfwf(ji,jj) * rLfusisf
+ zqh_b(ji,jj) = ptsc(ji,jj,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
+
ikt = mikt(ji,jj)
! compute Rc number (as done in zdfric.F90)
!!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation
@@ -108,82 +148,83 @@ CONTAINS
zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) &
& -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) )
! ! richardson number (minimum value set to zero)
- zRc(ji,jj) = MAX(rn2(ji,jj,ikt+1), 1.e-20_wp) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
+ zRc(ji,jj) = MAX( rn2(ji,jj,ikt+1), 1.e-20_wp ) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
END_2D
- CALL lbc_lnk( 'isfmlt', zRc, 'T', 1._wp )
+ !
ENDIF
!
- ! compute ice shelf melting
+ ! --- iteration loop --- !
nit = 1 ; lit = .TRUE.
DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine
!
! compute gammat everywhere (2d)
! useless if melt specified
IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
- CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce , pqfwf, zRc, &
- & zgammat, zgammas )
+ CALL isfcav_gammats( Kmm, zustar, zttbl, zstbl, zqoce, pfwf, zRc, & ! <<== in
+ & zgammat, zgammas ) ! ==>> out
END IF
!
! compute tfrz, latent heat and melt (2d)
- CALL isfcav_mlt(kt, zgammat, zgammas, zttbl, zstbl, &
- & zqhc , zqoce, pqfwf )
+ CALL isfcav_mlt(kt, zgammat, zgammas, zttbl, zstbl, & ! <<== in
+ & zqhc , zqoce, pfwf ) ! ==>> out
!
! define if we need to iterate
SELECT CASE ( cn_gammablk )
- CASE ( 'spe','vel' )
+ CASE ( 'spe', 'vel' )
! no convergence needed
lit = .FALSE.
CASE ( 'vel_stab' )
! compute error between 2 iterations
zerr = 0._wp
DO_2D( 0, 0, 0, 0 )
- zerr = MAX( zerr, ABS(zqhc(ji,jj)+zqoce(ji,jj) - zqh_b(ji,jj)) )
+ zerr = MAX( zerr, ABS( zqhc(ji,jj) + zqoce(ji,jj) - zqh_b(ji,jj) ) )
END_2D
CALL mpp_max( 'isfcav', zerr ) ! max over the global domain
!
! define if iteration needed
- IF (nit >= 100) THEN ! too much iteration
+ IF( nit >= 100 ) THEN ! too much iteration
CALL ctl_stop( 'STOP', 'isf_cav: vel_stab gamma formulation had too many iterations ...' )
- ELSE IF ( zerr <= 0.01_wp ) THEN ! convergence is achieve
+ ELSEIF( zerr <= 0.01_wp ) THEN ! convergence is achieve
lit = .FALSE.
- ELSE ! converge is not yet achieve
+ ELSE ! converge is not yet achieve
nit = nit + 1
- zqh_b(:,:) = zqhc(:,:)+zqoce(:,:)
+ zqh_b(:,:) = zqhc(:,:) + zqoce(:,:)
END IF
END SELECT
!
END DO
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
! compute heat and water flux ( > 0 from isf to oce)
- pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_cav(ji,jj)
+ pfwf (ji,jj) = pfwf (ji,jj) * mskisf_cav(ji,jj)
zqoce(ji,jj) = zqoce(ji,jj) * mskisf_cav(ji,jj)
- zqhc (ji,jj) = zqhc(ji,jj) * mskisf_cav(ji,jj)
+ zqhc (ji,jj) = zqhc (ji,jj) * mskisf_cav(ji,jj)
!
! compute heat content flux ( > 0 from isf to oce)
- zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2)
+ zqlat(ji,jj) = - pfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2)
!
! total heat flux ( > 0 from isf to oce)
- zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
+ zqh(ji,jj) = ( zqhc(ji,jj) + zqoce(ji,jj) )
!
! set temperature content
ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
END_2D
- CALL lbc_lnk( 'isfmlt', pqfwf, 'T', 1.0_wp)
!
! output fluxes
- CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
+ CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pfwf, zqoce, zqlat, zqhc )
!
-#if ! defined key_RK3
- ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
- IF (lrst_oce) CALL isfrst_write(kt, 'cav', ptsc, pqfwf)
-#endif
+ ! --- output and debug --- !
+ !
+ CALL iom_put( 'isfgammat', zgammat(:,:) )
+ CALL iom_put( 'isfgammas', zgammas(:,:) )
!
IF ( ln_isfdebug ) THEN
IF(lwp) WRITE(numout,*) ''
+ CALL debug( 'isfcav gammaT:', zgammat(:,:) )
+ CALL debug( 'isfcav gammaS:', zgammas(:,:) )
CALL debug('isf_cav: ptsc T',ptsc(:,:,1))
CALL debug('isf_cav: ptsc S',ptsc(:,:,2))
- CALL debug('isf_cav: pqfwf fwf',pqfwf(:,:))
+ CALL debug('isf_cav: pfwf fwf',pfwf(:,:))
IF(lwp) WRITE(numout,*) ''
END IF
!
@@ -196,27 +237,31 @@ CONTAINS
!! ** Purpose : initialisation of variable needed to compute melt under an ice shelf
!!
!!----------------------------------------------------------------------
- INTEGER :: ierr
+ INTEGER :: ierr
+ INTEGER :: ji, jj ! dummy loop indices
!!---------------------------------------------------------------------
!
!==============
! 0: allocation
!==============
- !
CALL isf_alloc_cav()
!
!==================
! 1: initialisation
!==================
- !
- ! top and bottom level of the 'top boundary layer'
- misfkt_cav(:,:) = mikt(:,:) ; misfkb_cav(:,:) = 1
- !
- ! thickness of 'tbl' and fraction of bottom cell affected by 'tbl'
- rhisf_tbl_cav(:,:) = 0.0_wp ; rfrac_tbl_cav(:,:) = 0.0_wp
- !
- ! cavity mask
- mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
+ DO_2D( 0, 0, 0, 0 )
+ ! top and bottom level of the 'top boundary layer'
+ misfkt_cav(ji,jj) = mikt(ji,jj)
+ misfkb_cav(ji,jj) = 1
+ !
+ ! thickness of 'tbl' and fraction of bottom cell affected by 'tbl'
+ rhisf_tbl_cav(ji,jj) = 0.0_wp
+ rfrac_tbl_cav(ji,jj) = 0.0_wp
+ !
+ ! cavity mask
+ mskisf_cav(ji,jj) = ( 1._wp - tmask(ji,jj,1) ) * ssmask(ji,jj)
+ !
+ END_2D
!
#if ! defined key_RK3
!================
@@ -239,7 +284,7 @@ CONTAINS
CASE( 'spe' )
ALLOCATE( sf_isfcav_fwf(1), STAT=ierr )
- ALLOCATE( sf_isfcav_fwf(1)%fnow(jpi,jpj,1), sf_isfcav_fwf(1)%fdta(jpi,jpj,1,2) )
+ ALLOCATE( sf_isfcav_fwf(1)%fnow(A2D(0),1), sf_isfcav_fwf(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_isfdir, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
IF(lwp) WRITE(numout,*)
diff --git a/src/OCE/ISF/isfcavgam.F90 b/src/OCE/ISF/isfcavgam.F90
index 6c0ac2a4d..fbab4647e 100644
--- a/src/OCE/ISF/isfcavgam.F90
+++ b/src/OCE/ISF/isfcavgam.F90
@@ -16,10 +16,10 @@ MODULE isfcavgam
USE oce , ONLY: uu, vv ! ocean dynamics
USE phycst , ONLY: grav, vkarmn ! physical constant
USE eosbn2 , ONLY: eos_rab ! equation of state
- USE zdfdrg , ONLY: rCd0_top, r_ke0_top ! vertical physics: top/bottom drag coef.
USE iom , ONLY: iom_put !
USE lib_mpp , ONLY: ctl_stop
+ USE par_oce ! ocean space and time domain
USE dom_oce ! ocean space and time domain
USE in_out_manager ! I/O manager
!
@@ -39,88 +39,45 @@ MODULE isfcavgam
!!----------------------------------------------------------------------
CONTAINS
!
- !!-----------------------------------------------------------------------------------------------------
- !! PUBLIC SUBROUTINES
- !!-----------------------------------------------------------------------------------------------------
- !
- SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs )
+ SUBROUTINE isfcav_gammats( Kmm, pustar, pttbl, pstbl, pqoce, pfwf, pRc, pgt, pgs )
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange for heat and fwf flux
!!
!! ** Method : select the gamma formulation
!! 3 method available (cst, vel and vel_stab)
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt , pgs ! gamma t and gamma s
- !!-------------------------- IN -------------------------------------
- INTEGER :: Kmm ! ocean time level index
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! isf heat and fwf
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number
- !!---------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity
+ INTEGER :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pustar ! u*
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pqoce, pfwf ! isf heat and fwf
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pRc ! Richardson number
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pgt , pgs ! gamma t and gamma s
!!---------------------------------------------------------------------
!
- !==========================================
- ! 1.: compute velocity in the tbl if needed
- !==========================================
+ ! --- compute gamma --- !
!
- SELECT CASE ( cn_gammablk )
- CASE ( 'spe' )
- ! gamma is constant (specified in namelist)
- ! nothing to do
- CASE ('vel', 'vel_stab')
- ! compute velocity in tbl
- CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
- CALL isf_tbl(Kmm, vv(:,:,:,Kmm) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav)
- !
- ! mask velocity in tbl with ice shelf mask
- zutbl(:,:) = zutbl(:,:) * mskisf_cav(:,:)
- zvtbl(:,:) = zvtbl(:,:) * mskisf_cav(:,:)
+ SELECT CASE( cn_gammablk )
+ CASE( 'spe' ) ! gamma is constant (specified in namelist)
!
- ! output
- CALL iom_put('utbl',zutbl(:,:))
- CALL iom_put('vtbl',zvtbl(:,:))
- CASE DEFAULT
- CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
- END SELECT
- !
- !==========================================
- ! 2.: compute gamma
- !==========================================
- !
- SELECT CASE ( cn_gammablk )
- CASE ( 'spe' ) ! gamma is constant (specified in namelist)
pgt(:,:) = rn_gammat0
pgs(:,:) = rn_gammas0
- CASE ( 'vel' ) ! gamma is proportional to u*
- CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs )
- CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
- CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs )
+ !
+ CASE( 'vel' ) ! gamma is proportional to u*
+ !
+ CALL gammats_vel( pustar, pgt, pgs )
+ !
+ CASE( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
+ !
+ CALL gammats_vel_stab( Kmm, pustar, pttbl, pstbl, pqoce, pfwf, pRc, pgt, pgs )
+ !
CASE DEFAULT
CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
END SELECT
!
- !==========================================
- ! 3.: output and debug
- !==========================================
- !
- CALL iom_put('isfgammat', pgt(:,:))
- CALL iom_put('isfgammas', pgs(:,:))
- !
- IF (ln_isfdebug) THEN
- CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
- CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
- END IF
- !
END SUBROUTINE isfcav_gammats
!
- !!-----------------------------------------------------------------------------------------------------
- !! PRIVATE SUBROUTINES
- !!-----------------------------------------------------------------------------------------------------
- !
- SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, & ! <<== in
- & pgt, pgs ) ! ==>> out gammats [m/s]
+ SUBROUTINE gammats_vel( pustar, & ! <<== in
+ & pgt, pgs ) ! ==>> out gammats [m/s]
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange coefficient
!!
@@ -128,33 +85,22 @@ CONTAINS
!!
!! ** Reference : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt, pgs ! gammat and gammas [m/s]
- !!-------------------------- IN -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pCd ! drag coefficient
- REAL(wp), INTENT(in ) :: pke2 ! background velocity
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pustar ! u*
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pgt, pgs ! gammat and gammas [m/s]
!!---------------------------------------------------------------------
- INTEGER :: ji, jj ! loop index
- REAL(wp), DIMENSION(jpi,jpj) :: zustar
+ INTEGER :: ji, jj ! loop index
!!---------------------------------------------------------------------
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- ! compute ustar (AD15 eq. 27)
- zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj)
- !
- ! Compute gammats
- pgt(ji,jj) = zustar(ji,jj) * rn_gammat0
- pgs(ji,jj) = zustar(ji,jj) * rn_gammas0
+ DO_2D( 0, 0, 0, 0 )
+ pgt(ji,jj) = pustar(ji,jj) * rn_gammat0
+ pgs(ji,jj) = pustar(ji,jj) * rn_gammas0
END_2D
!
- ! output ustar
- CALL iom_put('isfustar',zustar(:,:))
!
END SUBROUTINE gammats_vel
- SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, pRc, & ! <<== in
- & pgt , pgs ) ! ==>> out gammats [m/s]
+ SUBROUTINE gammats_vel_stab( Kmm, pustar, pttbl, pstbl, pqoce, pfwf, pRc, & ! <<== in
+ & pgt , pgs ) ! ==>> out gammats [m/s]
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange coefficient
!!
@@ -162,52 +108,38 @@ CONTAINS
!!
!! ** Reference : Holland and Jenkins, 1999, JPO, p1787-1800
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt, pgs ! gammat and gammas
- !!-------------------------- IN -------------------------------------
- INTEGER :: Kmm ! ocean time level index
- REAL(wp), INTENT(in ) :: pke2 ! background velocity squared
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! surface heat flux and fwf flux
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pCd ! drag coeficient
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! tracer in the losch top boundary layer
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number
+ INTEGER :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pustar ! u*
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pttbl, pstbl ! tracer in the losch top boundary layer
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pqoce, pfwf ! surface heat flux and fwf flux
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pRc ! Richardson number
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pgt, pgs ! gammat and gammas
!!---------------------------------------------------------------------
- INTEGER :: ji, jj ! loop index
- INTEGER :: ikt ! local integer
- REAL(wp) :: zdku, zdkv ! U, V shear
- REAL(wp) :: zPr, zSc ! Prandtl and Scmidth number
- REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point
- REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness
- REAL(wp) :: zhmax ! limitation of mol
- REAL(wp) :: zetastar ! stability parameter
- REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence
- REAL(wp) :: zcoef ! temporary coef
- REAL(wp) :: zdep
- REAL(wp) :: zeps = 1.0e-20_wp
- REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant
- REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
- REAL(wp), DIMENSION(2) :: zts, zab
- REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
+ INTEGER :: ji, jj ! loop index
+ INTEGER :: ikt ! local integer
+ REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point
+ REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness
+ REAL(wp) :: zhmax ! limitation of mol
+ REAL(wp) :: zetastar ! stability parameter
+ REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence
+ REAL(wp) :: zdep
+ REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant
+ REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
+ REAL(wp), PARAMETER :: zPr = 13.8_wp
+ REAL(wp), PARAMETER :: zSc = 2432.0_wp
+ REAL(wp), DIMENSION(jpts) :: zts, zab
!!---------------------------------------------------------------------
!
- ! compute Pr and Sc number (eq ??)
- zPr = 13.8_wp
- zSc = 2432.0_wp
- !
! compute gamma mole (eq ??)
zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
!
! compute gamma
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
ikt = mikt(ji,jj)
- ! compute ustar
- zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )
-
- IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think
+ IF( pustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think
pgt(ji,jj) = rn_gammat0
pgs(ji,jj) = rn_gammas0
ELSE
@@ -219,34 +151,32 @@ CONTAINS
CALL eos_rab( zts, zdep, zab, Kmm )
!
! compute length scale (Eq ??)
- zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
+ zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pfwf(ji,jj) )
!
! compute Monin Obukov Length
! Maximum boundary layer depth (Eq ??)
zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
!
! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
- zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+ zmob = pustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + 1.e-20_wp ))
zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
!
! compute eta* (stability parameter) (Eq ??)
- zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * zustar(ji,jj) &
+ zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * pustar(ji,jj) &
& / MAX( 1.e-20, ABS(ff_t(ji,jj)) * zmols * pRc(ji,jj) ) )))
!
! compute the sublayer thickness (Eq ??)
- zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) )
+ zhnu = 5 * znu / MAX( 1.e-20, pustar(ji,jj) )
!
! compute gamma turb (Eq ??)
- zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) &
- & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+ zgturb = 1._wp / vkarmn * LOG(pustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10_wp, ABS(ff_t(ji,jj)) * zhnu )) &
+ & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
!
! compute gammats
- pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
- pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
+ pgt(ji,jj) = pustar(ji,jj) / (zgturb + zgmolet)
+ pgs(ji,jj) = pustar(ji,jj) / (zgturb + zgmoles)
END IF
END_2D
- ! output ustar
- CALL iom_put('isfustar',zustar(:,:))
END SUBROUTINE gammats_vel_stab
diff --git a/src/OCE/ISF/isfcavmlt.F90 b/src/OCE/ISF/isfcavmlt.F90
index 665d75bbc..7bf1a4976 100644
--- a/src/OCE/ISF/isfcavmlt.F90
+++ b/src/OCE/ISF/isfcavmlt.F90
@@ -16,6 +16,7 @@ MODULE isfcavmlt
USE isftbl , ONLY: isf_tbl ! ice shelf depth average
USE isfutils,ONLY: debug ! debug subroutine
+ USE par_oce ! ocean space and time domain
USE dom_oce ! ocean space and time domain
USE phycst , ONLY: rcp, rho0, rho0_rcp ! physical constants
USE eosbn2 , ONLY: eos_fzp ! equation of state
@@ -40,40 +41,44 @@ MODULE isfcavmlt
!!----------------------------------------------------------------------
CONTAINS
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
- SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
- & pqhc, pqoce, pqfwf )
+ SUBROUTINE isfcav_mlt( kt, pgt, pgs , pttbl, pstbl, &
+ & pqhc, pqoce, pfwf )
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt ***
!!
!! ** Purpose : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER, INTENT(in) :: kt
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt , pgs ! gamma t and gamma s
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
+ !!---------------------------------------------------------------------
+ INTEGER, INTENT(in) :: kt
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pgt , pgs ! gamma t and gamma s
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! heat and fwf fluxes
!!---------------------------------------------------------------------
!
! compute latent heat and melt (2d)
SELECT CASE ( cn_isfcav_mlt )
CASE ( 'spe' ) ! ice shelf melt specified (read input file, and heat fluxes derived from
+ !
CALL isfcav_mlt_spe( kt, pstbl, &
- & pqhc, pqoce, pqfwf )
+ & pqhc, pqoce, pfwf )
+ !
CASE ( '2eq' ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
+ !
CALL isfcav_mlt_2eq( pgt, pttbl, pstbl, &
- & pqhc , pqoce, pqfwf )
+ & pqhc , pqoce, pfwf )
+ !
CASE ( '3eq' ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+ !
CALL isfcav_mlt_3eq( pgt, pgs , pttbl, pstbl, &
- & pqhc, pqoce, pqfwf )
+ & pqhc, pqoce, pfwf )
+ !
CASE ( 'oasis' ) ! fwf pass trough oasis
- CALL isfcav_mlt_oasis( kt, pstbl, &
- & pqhc, pqoce, pqfwf )
+ !
+ CALL isfcav_mlt_oasis( kt, pstbl, &
+ & pqhc, pqoce, pfwf )
+ !
CASE DEFAULT
CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
END SELECT
@@ -82,18 +87,15 @@ CONTAINS
IF(lwp) WRITE(numout,*) ''
CALL debug( 'isfcav_mlt qhc :', pqhc (:,:) )
CALL debug( 'isfcav_mlt qoce :', pqoce(:,:) )
- CALL debug( 'isfcav_mlt qfwf :', pqfwf(:,:) )
+ CALL debug( 'isfcav_mlt fwf :', pfwf(:,:) )
IF(lwp) WRITE(numout,*) ''
END IF
!
END SUBROUTINE isfcav_mlt
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
SUBROUTINE isfcav_mlt_spe(kt, pstbl, & ! <<== in
- & pqhc , pqoce, pqfwf ) ! ==>> out
+ & pqhc , pqoce, pfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_spe ***
@@ -102,13 +104,12 @@ CONTAINS
!! - compute ocea-ice heat flux (assuming it is equal to latent heat)
!! - compute heat content flux
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat content, latent heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER , INTENT(in ) :: kt ! current time step
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pstbl ! salinity in tbl
+ INTEGER , INTENT(in ) :: kt ! current time step
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pstbl ! salinity in tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! heat content, latent heat and fwf fluxes
!!--------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! tbl freezing temperature
+ INTEGER :: ji, jj ! dummy loop indices
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! tbl freezing temperature
!!--------------------------------------------------------------------
!
! Compute freezing temperature
@@ -119,9 +120,11 @@ CONTAINS
!
! define fwf and qoce
! ocean heat flux is assume to be equal to the latent heat
- pqfwf(:,:) = sf_isfcav_fwf(1)%fnow(:,:,1) ! fwf ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
+ DO_2D( 0, 0, 0, 0 )
+ pfwf(ji,jj) = sf_isfcav_fwf(1)%fnow(ji,jj,1) ! fwf ( > 0 from isf to oce)
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ END_2D
!
! output freezing point at the interface
CALL iom_put('isftfrz_cav', ztfrz(:,:) * mskisf_cav(:,:) )
@@ -129,7 +132,7 @@ CONTAINS
END SUBROUTINE isfcav_mlt_spe
SUBROUTINE isfcav_mlt_2eq(pgt , pttbl, pstbl, & ! <<== in
- & pqhc, pqoce, pqfwf ) ! ==>> out
+ & pqhc, pqoce, pfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_2eq ***
@@ -138,43 +141,45 @@ CONTAINS
!!
!! ** Method : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
!! From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
- !! qfwf = Gammat * rho0 * Cp * ( Tw - Tfrz ) / Lf
+ !! fwf = Gammat * rho0 * Cp * ( Tw - Tfrz ) / Lf
!! qhoce = qlat
- !! qhc = qfwf * Cp * Tfrz
+ !! qhc = fwf * Cp * Tfrz
!!
!! ** Reference : Hunter, J. R.: Specification for test models of ice shelf cavities,
!! Tech. Rep. June, Antarctic Climate & Ecosystems Cooperative Research Centre, available at:
!! http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! hean content, ocean-ice heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt ! temperature exchange coeficient
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! temperature and salinity in top boundary layer
!!--------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing temperature
- REAL(wp), DIMENSION(jpi,jpj) :: zthd ! thermal driving
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pgt ! temperature exchange coeficient
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pttbl, pstbl ! temperature and salinity in top boundary layer
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! hean content, ocean-ice heat and fwf fluxes
+ !!--------------------------------------------------------------------
+ INTEGER :: ji, jj ! dummy loop indices
+ REAL(wp) :: zthd ! thermal driving
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! freezing temperature
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
!
- ! thermal driving
- zthd (:,:) = ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:)
- !
- ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
- pqfwf(:,:) = pgt(:,:) * rho0_rcp * zthd(:,:) / rLfusisf ! fresh water flux ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocea-ice flux ( > 0 from isf to oce)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
- !
+ DO_2D( 0, 0, 0, 0 )
+ ! thermal driving
+ zthd = ( pttbl(ji,jj) - ztfrz(ji,jj) ) * mskisf_cav(ji,jj)
+ !
+ ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
+ pfwf (ji,jj) = pgt (ji,jj) * rho0_rcp * zthd / rLfusisf ! fresh water flux ( > 0 from isf to oce)
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocea-ice flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ !
+ END_2D
! output thermal driving and freezinpoint at the ice shelf interface
- CALL iom_put('isfthermald_cav', zthd )
- CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
+ CALL iom_put('isfthermald_cav', ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:) )
+ CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_2eq
- SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, & ! <<== in
- & pqhc, pqoce, pqfwf ) ! ==>> out
+ SUBROUTINE isfcav_mlt_3eq( pgt, pgs , pttbl, pstbl, & ! <<== in
+ & pqhc, pqoce, pfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_3eq ***
@@ -194,25 +199,24 @@ CONTAINS
!! MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1),
!! Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016.
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! latent heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt , pgs ! heat/salt exchange coeficient
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! mean temperature and salinity in top boundary layer
+ !!--------------------------------------------------------------------
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pgt , pgs ! heat/salt exchange coeficient
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pttbl, pstbl ! mean temperature and salinity in top boundary layer
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! latent heat and fwf fluxes
!!--------------------------------------------------------------------
REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 ! dummy local scalar for quadratic equation resolution
REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac ! dummy local scalar for quadratic equation resolution
REAL(wp) :: zeps = 1.e-20
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing point
- REAL(wp), DIMENSION(jpi,jpj) :: zqcon ! conductive flux through the ice shelf
- REAL(wp), DIMENSION(jpi,jpj) :: zthd ! thermal driving
+ REAL(wp) :: zthd ! thermal driving
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! freezing point
+ REAL(wp), DIMENSION(A2D(0)) :: zqcnd ! conductive flux through the ice shelf
!
INTEGER :: ji, jj ! dummy loop indices
!!--------------------------------------------------------------------
!
! compute upward heat flux zhtflx and upward water flux zwflx
! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
!
! compute coeficient to solve the 2nd order equation
zeps1 = rho0_rcp * pgt(ji,jj)
@@ -237,28 +241,28 @@ CONTAINS
ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
!
! thermal driving
- zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) )
+ zthd = ( pttbl(ji,jj) - ztfrz(ji,jj) )
!
! compute the upward water and heat flux (eq. 24 and eq. 26)
- pqfwf(ji,jj) = - rho0 * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux ( > 0 from isf to oce)
- pqoce(ji,jj) = - rho0_rcp * pgt(ji,jj) * zthd (ji,jj) ! ocean-ice heat flux ( > 0 from isf to oce)
- pqhc (ji,jj) = rcp * pqfwf(ji,jj) * ztfrz(ji,jj) ! heat content flux ( > 0 from isf to oce)
+ pfwf(ji,jj) = - rho0 * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux ( > 0 from isf to oce)
+ pqoce(ji,jj) = - rho0_rcp * pgt(ji,jj) * zthd ! ocean-ice heat flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = rcp * pfwf(ji,jj) * ztfrz(ji,jj) ! heat content flux ( > 0 from isf to oce)
!
- zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )
+ zqcnd(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )
!
END_2D
!
! output conductive heat flux through the ice
- CALL iom_put('qconisf', zqcon(:,:) * mskisf_cav(:,:) )
+ CALL iom_put('qconisf', zqcnd(:,:) * mskisf_cav(:,:) )
!
! output thermal driving and freezing point at the interface
- CALL iom_put('isfthermald_cav', zthd (:,:) * mskisf_cav(:,:) )
- CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
+ CALL iom_put('isfthermald_cav', ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:) )
+ CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_3eq
- SUBROUTINE isfcav_mlt_oasis(kt, pstbl, & ! <<== in
- & pqhc , pqoce, pqfwf ) ! ==>> out
+ SUBROUTINE isfcav_mlt_oasis(kt, pstbl, & ! <<== in
+ & pqhc , pqoce, pfwf ) ! ==>> out
!!----------------------------------------------------------------------
!! *** ROUTINE isfcav_mlt_oasis ***
!!
@@ -269,44 +273,43 @@ CONTAINS
!! - scale fwf and compute heat fluxes
!!
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat content, latent heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER , INTENT(in ) :: kt ! current time step
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pstbl ! salinity in tbl
+ INTEGER , INTENT(in ) :: kt ! current time step
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in ) :: pstbl ! salinity in tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! heat content, latent heat and fwf fluxes
!!--------------------------------------------------------------------
- REAL(wp) :: zfwf_fld, zfwf_oasis ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! tbl freezing temperature
- REAL(wp), DIMENSION(jpi,jpj) :: zfwf ! 2d fwf map after scaling
+ REAL(wp) :: zfwf_fld, zfwf_oasis ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! tbl freezing temperature
+ INTEGER :: ji, jj
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
- CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+ CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(A2D(0)) )
!
! read input file of fwf from isf to oce
CALL fld_read ( kt, 1, sf_isfcav_fwf )
!
! ice shelf 2d map
- zfwf(:,:) = sf_isfcav_fwf(1)%fnow(:,:,1)
+ pfwf(:,:) = sf_isfcav_fwf(1)%fnow(:,:,1)
!
! compute glob sum from input file
! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
- zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
+ zfwf_fld = glob_sum( 'isfcav_mlt', e1e2t(A2D(0)) * pfwf(:,:) )
!
! compute glob sum from atm->oce ice shelf fwf
! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
- zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
+ zfwf_oasis = glob_sum( 'isfcav_mlt', e1e2t(A2D(0)) * fwfisf_oasis(:,:) )
!
! scale fwf
- zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
+ pfwf(:,:) = pfwf(:,:) * zfwf_oasis / zfwf_fld
!
- ! define fwf and qoce
+ ! define qoce
! ocean heat flux is assume to be equal to the latent heat
- pqfwf(:,:) = zfwf(:,:) ! fwf ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
+ DO_2D( 0, 0, 0, 0 )
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ END_2D
!
- CALL iom_put('isftfrz_cav', ztfrz * mskisf_cav(:,:) )
+ CALL iom_put('isftfrz_cav', ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_oasis
diff --git a/src/OCE/ISF/isfcpl.F90 b/src/OCE/ISF/isfcpl.F90
index 4989dc850..574875e3d 100644
--- a/src/OCE/ISF/isfcpl.F90
+++ b/src/OCE/ISF/isfcpl.F90
@@ -11,6 +11,7 @@ MODULE isfcpl
!!----------------------------------------------------------------------
!! isfrst : read/write iceshelf variables in/from restart
!!----------------------------------------------------------------------
+ USE par_oce ! ocean space and time domain
USE oce ! ocean dynamics and tracers
#if defined key_qco
USE domqco , ONLY : dom_qco_zgr ! vertical scale factor interpolation
@@ -32,7 +33,6 @@ MODULE isfcpl
PRIVATE
PUBLIC isfcpl_rst_write, isfcpl_init ! iceshelf restart read and write
- PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons ! iceshelf correction for ssh, tra, dyn and conservation
TYPE isfcons
INTEGER :: ii ! i global
@@ -79,6 +79,13 @@ CONTAINS
! allocation and initialisation to 0
CALL isf_alloc_cpl()
!
+ risfcpl_tsc (:,:,:,:) = 0._wp
+ risfcpl_cons_tsc(:,:,:,:) = 0._wp
+ risfcpl_vol (:,:,:) = 0._wp
+ risfcpl_cons_vol(:,:,:) = 0._wp
+ risfcpl_ssh (:,:) = 0._wp
+ risfcpl_cons_ssh(:,:) = 0._wp
+ !
! check presence of variable needed for coupling
! iom_varid return 0 if not found
id = 1
@@ -115,11 +122,11 @@ CONTAINS
!
! all before fields set to now values
ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)
- uu (:,:,:,Kbb) = uu (:,:,:,Kmm)
- vv (:,:,:,Kbb) = vv (:,:,:,Kmm)
+ uu (:,:,:,Kbb) = uu (:,:,:,Kmm)
+ vv (:,:,:,Kbb) = vv (:,:,:,Kmm)
ssh (:,:,Kbb) = ssh (:,:,Kmm)
#if ! defined key_qco && ! defined key_linssh
- e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+ e3t (:,:,:,Kbb) = e3t (:,:,:,Kmm)
#endif
END SUBROUTINE isfcpl_init
@@ -135,23 +142,30 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! ocean time level index
!!----------------------------------------------------------------------
INTEGER :: jk ! loop index
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v, zgdepw ! for qco substitution
+ REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmp ! for qco substitution
!!----------------------------------------------------------------------
!
+ CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask )
+ CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask )
+ !
+ DO jk = 1, jpk
+ ztmp(:,:,jk) = e3t(:,:,jk,Kmm)
+ END DO
+ CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , ztmp )
DO jk = 1, jpk
- ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
- ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
- ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+ ztmp(:,:,jk) = e3u(:,:,jk,Kmm)
+ END DO
+ CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , ztmp )
+ DO jk = 1, jpk
+ ztmp(:,:,jk) = e3v(:,:,jk,Kmm)
+ END DO
+ CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , ztmp )
!
- zgdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
+ DO jk = 1, jpk
+ ztmp(:,:,jk) = gdepw(:,:,jk,Kmm)
END DO
+ CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', ztmp )
!
- CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask )
- CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask )
- CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , ze3t )
- CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , ze3u )
- CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , ze3v )
- CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', zgdepw )
!
END SUBROUTINE isfcpl_rst_write
@@ -166,13 +180,14 @@ CONTAINS
!!
!!----------------------------------------------------------------------
!!
- INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
+ INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
!!----------------------------------------------------------------------
- INTEGER :: ji, jj, jd, jk !! loop index
- INTEGER :: jip1, jim1, jjp1, jjm1
+ INTEGER :: ji, jj, jd, jk !! loop index
+ INTEGER :: jip1, jim1, jjp1, jjm1
+ INTEGER :: ibnd
!!
- REAL(wp):: zsummsk
- REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
+ REAL(wp):: zdssmask, zsummsk
+ REAL(wp), DIMENSION(jpi,jpj) :: zssmask0, zssmask_b, zssh
!!----------------------------------------------------------------------
!
CALL iom_get( numror, jpdom_auto, 'ssmask' , zssmask_b ) ! need to extrapolate T/S
@@ -185,14 +200,21 @@ CONTAINS
!
DO jd = 1, nn_drown
!
- zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
- DO_2D( 0, 0, 0, 0 )
+ IF ( MOD(jd,2) == 0 ) THEN ! even
+ ibnd = nn_hls-1
+ ELSEIF( MOD(jd,2) == 1 ) THEN ! odd
+ ibnd = 0
+ ENDIF
+
+ DO_2D( ibnd, ibnd, ibnd, ibnd )
+ zdssmask = ssmask(ji,jj) - zssmask0(ji,jj)
+
jip1=ji+1 ; jim1=ji-1
jjp1=jj+1 ; jjm1=jj-1
!
zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
!
- IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
+ IF( zdssmask == 1._wp .AND. zsummsk /= 0._wp ) THEN
ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj) &
& + zssh(jim1,jj)*zssmask0(jim1,jj) &
& + zssh(ji,jjp1)*zssmask0(ji,jjp1) &
@@ -200,11 +222,16 @@ CONTAINS
zssmask_b(ji,jj) = 1._wp
ENDIF
END_2D
- CALL lbc_lnk( 'isfcpl', ssh(:,:,Kmm), 'T', 1.0_wp, zssmask_b(:,:), 'T', 1.0_wp )
- !
- zssh(:,:) = ssh(:,:,Kmm)
- zssmask0(:,:) = zssmask_b(:,:)
!
+ ! update ssh and mask
+ DO_2D( ibnd, ibnd, ibnd, ibnd )
+ zssh (ji,jj) = ssh (ji,jj,Kmm)
+ zssmask0(ji,jj) = zssmask_b(ji,jj)
+ END_2D
+
+ IF( ibnd == 0 ) THEN
+ CALL lbc_lnk( 'isfcpl', zssh(:,:), 'T', 1.0_wp, zssmask0(:,:), 'T', 1.0_wp )
+ ENDIF
!
END DO
!
@@ -244,19 +271,16 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: Kmm ! ocean time level index
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
+ REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b, ztmask0
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
!REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b !! depth w before
!!
- INTEGER :: ji, jj, jk, jd !! loop index
- INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
- !!
- REAL(wp):: zsummsk
- REAL(wp):: zdz, zdzm1, zdzp1
+ INTEGER :: ji, jj, jk, jd !! loop index
+ INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
+ INTEGER :: ibnd
!!
- REAL(wp), DIMENSION(jpi,jpj) :: zdmask
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d
- REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
+ REAL(wp):: zdmask, zsummsk
+ REAL(wp):: zdz, zdzm1, zdzp1
!!----------------------------------------------------------------------
!
CALL iom_get( numror, jpdom_auto, 'tmask' , ztmask_b ) ! need to extrapolate T/S
@@ -301,69 +325,77 @@ CONTAINS
! END DO
! END IF
- zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
- ztmask0(:,:,:) = ztmask_b(:,:,:)
- ztmask1(:,:,:) = ztmask_b(:,:,:)
+ zts0 (:,:,:,:) = ts (:,:,:,:,Kmm)
+ ztmask0(:,:,:) = ztmask_b(:,:,:)
!
! iterate the extrapolation processes nn_drown times
DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
- DO jk = 1,jpk-1
- !
+
+ IF ( MOD(jd,2) == 0 ) THEN ! even
+ ibnd = nn_hls-1
+ ELSEIF( MOD(jd,2) == 1 ) THEN ! odd
+ ibnd = 0
+ ENDIF
+
+ DO_3D( ibnd, ibnd, ibnd, ibnd, 1, jpkm1 )
! define new wet cell
- zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
+ zdmask = tmask(ji,jj,jk) - ztmask0(ji,jj,jk)
+
+ jip1=ji+1; jim1=ji-1;
+ jjp1=jj+1; jjm1=jj-1;
!
- DO_2D( 0, 0, 0, 0 )
- jip1=ji+1; jim1=ji-1;
- jjp1=jj+1; jjm1=jj-1;
- !
- ! check if a wet neigbourg cell is present
- zsummsk = ztmask0(jip1,jj ,jk) + ztmask0(jim1,jj ,jk) &
- + ztmask0(ji ,jjp1,jk) + ztmask0(ji ,jjm1,jk)
+ ! check if a wet neigbourg cell is present
+ zsummsk = ztmask0(jip1,jj ,jk) + ztmask0(jim1,jj ,jk) &
+ & + ztmask0(ji ,jjp1,jk) + ztmask0(ji ,jjm1,jk)
+ !
+ ! if neigbourg wet cell available at the same level
+ IF( zdmask == 1._wp .AND. zsummsk /= 0._wp ) THEN
!
- ! if neigbourg wet cell available at the same level
- IF ( zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
- !
- ! horizontal basic extrapolation
- ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj ,jk,1) * ztmask0(jip1,jj ,jk) &
+ ! horizontal basic extrapolation
+ ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj ,jk,1) * ztmask0(jip1,jj ,jk) &
& + zts0(jim1,jj ,jk,1) * ztmask0(jim1,jj ,jk) &
& + zts0(ji ,jjp1,jk,1) * ztmask0(ji ,jjp1,jk) &
& + zts0(ji ,jjm1,jk,1) * ztmask0(ji ,jjm1,jk) ) / zsummsk
- ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj ,jk,2) * ztmask0(jip1,jj ,jk) &
+ ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj ,jk,2) * ztmask0(jip1,jj ,jk) &
& + zts0(jim1,jj ,jk,2) * ztmask0(jim1,jj ,jk) &
& + zts0(ji ,jjp1,jk,2) * ztmask0(ji ,jjp1,jk) &
& + zts0(ji ,jjm1,jk,2) * ztmask0(ji ,jjm1,jk) ) / zsummsk
- !
- ! update mask for next pass
- ztmask1(ji,jj,jk)=1
- !
+ !
+ ! update mask for next pass
+ ztmask_b(ji,jj,jk)=1
+ !
! in case no neigbourg wet cell available at the same level
! check if a wet cell is available below
- ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
- !
- ! vertical extrapolation if horizontal extrapolation failed
- jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
- !
- ! check if a wet neigbourg cell is present
- zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
- IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
- ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
+ ELSEIF( zdmask == 1._wp .AND. zsummsk == 0._wp ) THEN
+ !
+ ! vertical extrapolation if horizontal extrapolation failed
+ jkm1=MAX(1,jk-1) ; jkp1=MIN(jpk,jk+1)
+ !
+ ! check if a wet neigbourg cell is present
+ zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
+ IF( zdmask == 1._wp .AND. zsummsk /= 0._wp ) THEN
+ ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
& + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
- ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
+ ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
& + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
- !
- ! update mask for next pass
- ztmask1(ji,jj,jk)=1._wp
- END IF
+ !
+ ! update mask for next pass
+ ztmask_b(ji,jj,jk)=1._wp
END IF
- END_2D
- END DO
- !
- CALL lbc_lnk( 'isfcpl', ts(:,:,:,jp_tem,Kmm), 'T', 1.0_wp, ts(:,:,:,jp_sal,Kmm), 'T', 1.0_wp, ztmask1, 'T', 1.0_wp)
- !
+ END IF
+ !
+ END_3D
+
! update temperature and salinity and mask
- zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
- ztmask0(:,:,:) = ztmask1(:,:,:)
+ DO_3D( ibnd, ibnd, ibnd, ibnd, 1, jpk )
+ zts0 (ji,jj,jk,1) = ts (ji,jj,jk,1,Kmm)
+ zts0 (ji,jj,jk,2) = ts (ji,jj,jk,2,Kmm)
+ ztmask0(ji,jj,jk) = ztmask_b(ji,jj,jk)
+ END_3D
!
+ IF( ibnd == 0 ) THEN
+ CALL lbc_lnk( 'isfcpl', zts0(:,:,:,1), 'T', 1.0_wp, zts0(:,:,:,2), 'T', 1.0_wp, ztmask0(:,:,:), 'T', 1.0_wp )
+ ENDIF
!
END DO ! nn_drown
!
@@ -374,7 +406,7 @@ CONTAINS
! sanity check
! -----------------------------------------------------------------------------------------
! case we open a cell but no neigbour cells available to get an estimate of T and S
- DO_3D( 0, 0, 0, 0, 1,jpk-1 )
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp) &
& CALL ctl_stop('STOP', 'failing to fill all new weet cell, &
& try increase nn_drown or activate XXXX &
@@ -455,25 +487,29 @@ CONTAINS
!
END_2D
!
- CALL lbc_lnk( 'isfcpl', risfcpl_vol, 'T', 1.0_wp )
- !
! 3.0: set total correction (div, tr(:,:,:,:,Krhs), ssh)
!
! 3.1: mask volume flux divergence correction
- risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ risfcpl_vol(ji,jj,jk) = risfcpl_vol(ji,jj,jk) * tmask(ji,jj,jk)
+ END_3D
!
! 3.2: get 3d tr(:,:,:,:,Krhs) increment to apply at the first time step
! temperature and salt content flux computed using local ts(:,:,:,:,Kmm)
! (very simple advection scheme)
! (>0 out)
- risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_tem,Kmm)
- risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_sal,Kmm)
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ risfcpl_tsc(ji,jj,jk,jp_tem) = -risfcpl_vol(ji,jj,jk) * ts(ji,jj,jk,jp_tem,Kmm)
+ risfcpl_tsc(ji,jj,jk,jp_sal) = -risfcpl_vol(ji,jj,jk) * ts(ji,jj,jk,jp_sal,Kmm)
+ END_3D
!
! 3.3: ssh correction (for dynspg_ts)
- risfcpl_ssh(:,:) = 0.0
- DO jk = 1,jpk
- risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
- END DO
+ DO_2D( 0, 0, 0, 0 )
+ risfcpl_ssh(ji,jj) = 0.0
+ END_2D
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ risfcpl_ssh(ji,jj) = risfcpl_ssh(ji,jj) + risfcpl_vol(ji,jj,jk) * r1_e1e2t(ji,jj)
+ END_3D
!
END SUBROUTINE isfcpl_vol
@@ -507,21 +543,17 @@ CONTAINS
REAL(wp) :: z1_sum, z1_rdtiscpl
REAL(wp) :: zdtem, zdsal, zdvol, zratio ! tem, sal, vol increment
REAL(wp) :: zlon , zlat ! target location
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b ! scale factor before
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b ! scale factor before
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b ! scale factor before
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztmask_b ! mask before
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ze3t_b ! scale factor before
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztmp ! scale factor before
!!----------------------------------------------------------------------
!==============================================================================
! 1.0: initialisation
!==============================================================================
-
! get restart variable
CALL iom_get( numror, jpdom_auto, 'tmask' , ztmask_b(:,:,:) ) ! need to extrapolate T/S
- CALL iom_get( numror, jpdom_auto, 'e3t_n' , ze3t_b(:,:,:) )
- CALL iom_get( numror, jpdom_auto, 'tn' , zt_b(:,:,:) )
- CALL iom_get( numror, jpdom_auto, 'sn' , zs_b(:,:,:) )
+
! compute run length
nstp_iscpl = nitend - nit000 + 1
@@ -540,32 +572,36 @@ CONTAINS
! 2.0: diagnose the heat, salt and volume input and compute the correction variable
! for case where we wet a cell or cell still wet (no change in cell status)
!==============================================================================
+ CALL iom_get( numror, jpdom_auto, 'e3t_n' , ze3t_b(:,:,:) )
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ ! volume diff
+ zdvol = e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
+ & - ze3t_b(ji,jj,jk ) * ztmask_b(ji,jj,jk)
+ ! volume differences in each cell (>0 means correction is an outward flux)
+ ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
+ risfcpl_cons_vol(ji,jj,jk) = ( zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk) ) * z1_rdtiscpl
+ END_3D
- DO jk = 1,jpk-1
- DO jj = Njs0,Nje0
- DO ji = Nis0,Nie0
-
- ! volume diff
- zdvol = e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
- & - ze3t_b(ji,jj,jk ) * ztmask_b(ji,jj,jk)
-
- ! heat diff
- zdtem = ts (ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
- - zt_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
-
- ! salt diff
- zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
- - zs_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+ CALL iom_get( numror, jpdom_auto, 'tn' , ztmp(:,:,:) )
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ ! heat diff
+ zdtem = ts (ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
+ & - ztmp(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+ ! heat differences in each cell (>0 means correction is an outward flux)
+ ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
+ risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
+ END_3D
- ! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
- ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
- risfcpl_cons_vol(ji,jj,jk) = ( zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk) ) * z1_rdtiscpl
- risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
- risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
+ CALL iom_get( numror, jpdom_auto, 'sn' , ztmp(:,:,:) )
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ ! salt diff
+ zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
+ & - ztmp(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
- END DO
- END DO
- END DO
+ ! salt differences in each cell (>0 means correction is an outward flux)
+ ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
+ risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
+ END_3D
!
!==============================================================================
! 3.0: diagnose the heat, salt and volume input and compute the correction variable
@@ -575,16 +611,15 @@ CONTAINS
! compute the total number of point receiving a correction increment for each processor
! local
nisfl(:)=0
- DO jk = 1,jpk-1
- DO jj = Njs0,Nje0
- DO ji = Nis0,Nie0
- jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
- IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
- nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
- ENDIF
- ENDDO
- ENDDO
- ENDDO
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ jip1=MIN(ji+1,jpi)
+ jim1=MAX(ji-1,1 )
+ jjp1=MIN(jj+1,jpj)
+ jjm1=MAX(jj-1,1 )
+ IF( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
+ nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
+ ENDIF
+ END_3D
!
! global
CALL mpp_sum('isfcpl',nisfl )
@@ -597,46 +632,42 @@ CONTAINS
! start computing the correction and fill zisfpts
! local
jisf = 0
- DO jk = 1,jpk-1
- DO jj = Njs0,Nje0
- DO ji = Nis0,Nie0
- IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
-
- jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
-
- zdvol = risfcpl_cons_vol(ji,jj,jk )
- zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
- zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
-
- IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
- ! spread correction amoung neigbourg wet cells (horizontal direction first)
- ! as it is a rude correction corner and lateral cell have the same weight
- !
- z1_sum = 1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
- !
- ! lateral cells
- IF (tmask(jip1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(jim1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(ji ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(ji ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
- !
- ! corner cells
- IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
- IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
- !
- ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
- ! spread correction amoung neigbourg wet cells (vertical direction)
- CALL update_isfpts(zisfpts, jisf, ji , jj , jk+1, zdvol, zdsal, zdtem, 1.0_wp, 0)
- ELSE
- ! need to find where to put correction in later on
- CALL update_isfpts(zisfpts, jisf, ji , jj , jk , zdvol, zdsal, zdtem, 1.0_wp, 1)
- END IF
- END IF
- END DO
- END DO
- END DO
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
+
+ jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
+
+ zdvol = risfcpl_cons_vol(ji,jj,jk )
+ zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
+ zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
+
+ IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
+ ! spread correction amoung neigbourg wet cells (horizontal direction first)
+ ! as it is a rude correction corner and lateral cell have the same weight
+ !
+ z1_sum = 1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
+ !
+ ! lateral cells
+ IF (tmask(jip1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(jim1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(ji ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(ji ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+ !
+ ! corner cells
+ IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+ IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+ !
+ ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
+ ! spread correction amoung neigbourg wet cells (vertical direction)
+ CALL update_isfpts(zisfpts, jisf, ji , jj , jk+1, zdvol, zdsal, zdtem, 1.0_wp, 0)
+ ELSE
+ ! need to find where to put correction in later on
+ CALL update_isfpts(zisfpts, jisf, ji , jj , jk , zdvol, zdsal, zdtem, 1.0_wp, 1)
+ END IF
+ END IF
+ END_3D
!
! share data among all processes because for some point we need to find the closest wet point (could be on other process)
DO jproc=1,jpnij
@@ -686,19 +717,22 @@ CONTAINS
!==============================================================================
!
! mask
- risfcpl_cons_vol(:,:,: ) = risfcpl_cons_vol(:,:,: ) * tmask(:,:,:)
- risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
- risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
- !
- ! add lbclnk
- CALL lbc_lnk( 'isfcpl', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1.0_wp, risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1.0_wp, &
- & risfcpl_cons_vol(:,:,:) , 'T', 1.0_wp)
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ risfcpl_cons_vol(ji,jj,jk ) = risfcpl_cons_vol(ji,jj,jk ) * tmask(ji,jj,jk)
+ risfcpl_cons_tsc(ji,jj,jk,jp_sal) = risfcpl_cons_tsc(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+ risfcpl_cons_tsc(ji,jj,jk,jp_tem) = risfcpl_cons_tsc(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+ END_3D
!
! ssh correction (for dynspg_ts)
- DO jk = 1,jpk
- risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
- END DO
- risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ risfcpl_cons_ssh(ji,jj) = risfcpl_cons_ssh(ji,jj) + risfcpl_cons_vol(ji,jj,jk)
+ END_3D
+ DO_2D( 0, 0, 0, 0 )
+ risfcpl_cons_ssh(ji,jj) = risfcpl_cons_ssh(ji,jj) * r1_e1e2t(ji,jj)
+ END_2D
+ !
+ ! deallocate list of point receiving correction
+ DEALLOCATE( zisfpts )
!
END SUBROUTINE isfcpl_cons
!
@@ -736,7 +770,7 @@ CONTAINS
END IF
!
! update isfpts structure
- sisfpts(kpts) = isfcons(mig(ki,nn_hls), mjg(kj,nn_hls), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, &
+ sisfpts(kpts) = isfcons(mig(ki,0), mjg(kj,0), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, &
& glamt(ki,kj), gphit(ki,kj), ifind )
!
END SUBROUTINE update_isfpts
@@ -761,9 +795,8 @@ CONTAINS
iig = ki ; ijg = kj
IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
!
- ! fill the correction array
- DO jj = mj0(ijg,nn_hls),mj1(ijg,nn_hls)
- DO ji = mi0(iig,nn_hls),mi1(iig,nn_hls)
+ DO jj = mj0(ijg,0),mj1(ijg,0)
+ DO ji = mi0(iig,0),mi1(iig,0)
! correct the vol_flx and corresponding heat/salt flx in the closest cell
risfcpl_cons_vol(ji,jj,kk) = risfcpl_cons_vol(ji,jj,kk ) + pvolinc
risfcpl_cons_tsc(ji,jj,kk,jp_sal) = risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
diff --git a/src/OCE/ISF/isfdiags.F90 b/src/OCE/ISF/isfdiags.F90
index b9f977364..6f9ec1e04 100644
--- a/src/OCE/ISF/isfdiags.F90
+++ b/src/OCE/ISF/isfdiags.F90
@@ -13,6 +13,7 @@ MODULE isfdiags
!!----------------------------------------------------------------------
USE in_out_manager ! I/O manager
+ USE par_oce ! ocean space and time domain
USE dom_oce
USE isf_oce ! ice shelf variable
USE iom !
@@ -34,7 +35,7 @@ MODULE isfdiags
CONTAINS
- SUBROUTINE isf_diags_flx(Kmm, ktop, kbot, phtbl, pfrac, cdisf, pqfwf, pqoce, pqlat, pqhc)
+ SUBROUTINE isf_diags_flx( Kmm, ktop, kbot, phtbl, pfrac, cdisf, pfwf, pqoce, pqlat, pqhc )
!!---------------------------------------------------------------------
!! *** ROUTINE isf_diags_flx ***
!!
@@ -42,39 +43,37 @@ CONTAINS
!! from isf to oce fwf, latent heat, heat content fluxes
!!
!!----------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- !!-------------------------- IN -------------------------------------
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
- INTEGER , DIMENSION(jpi,jpj), INTENT(in) :: ktop , kbot ! top and bottom level of the tbl
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac ! thickness of the tbl and fraction of last cell affected by the tbl
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqfwf, pqoce, pqlat, pqhc ! 2d var to map in 3d
- CHARACTER(LEN=3), INTENT(in) :: cdisf ! parametrisation or interactive melt
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ INTEGER , DIMENSION(A2D(0)), INTENT(in) :: ktop , kbot ! top and bottom level of the tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: phtbl, pfrac ! thickness of the tbl and fraction of last cell affected by the tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pfwf, pqoce, pqlat, pqhc ! 2d var to map in 3d
+ CHARACTER(LEN=3), INTENT(in) :: cdisf ! parametrisation or interactive melt
!!---------------------------------------------------------------------
- CHARACTER(LEN=256) :: cvarqfwf , cvarqoce , cvarqlat , cvarqhc
- CHARACTER(LEN=256) :: cvarqfwf3d, cvarqoce3d, cvarqlat3d, cvarqhc3d
+ CHARACTER(LEN=256) :: cvarfwf , cvarqoce , cvarqlat , cvarqhc
+ CHARACTER(LEN=256) :: cvarfwf3d, cvarqoce3d, cvarqlat3d, cvarqhc3d
!!---------------------------------------------------------------------
!
! output melt
- cvarqfwf = 'fwfisf_'//cdisf ; cvarqfwf3d = 'fwfisf3d_'//cdisf
+ cvarfwf = 'fwfisf_'//cdisf ; cvarfwf3d = 'fwfisf3d_'//cdisf
cvarqoce = 'qoceisf_'//cdisf ; cvarqoce3d = 'qoceisf3d_'//cdisf
cvarqlat = 'qlatisf_'//cdisf ; cvarqlat3d = 'qlatisf3d_'//cdisf
cvarqhc = 'qhcisf_'//cdisf ; cvarqhc3d = 'qhcisf3d_'//cdisf
!
! output 2d melt rate, latent heat and heat content flux from the injected water
- CALL iom_put( TRIM(cvarqfwf), pqfwf(:,:) ) ! mass flux ( > 0 from isf to oce)
+ CALL iom_put( TRIM(cvarfwf) , pfwf(:,:) ) ! mass flux ( > 0 from isf to oce)
CALL iom_put( TRIM(cvarqoce), pqoce(:,:) ) ! oce to ice flux ( > 0 from isf to oce)
CALL iom_put( TRIM(cvarqlat), pqlat(:,:) ) ! latent heat flux ( > 0 from isf to oce)
CALL iom_put( TRIM(cvarqhc) , pqhc (:,:) ) ! heat content flux ( > 0 from isf to oce)
!
! output 3d Diagnostics
- IF ( iom_use( TRIM(cvarqfwf3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqfwf3d) , pqfwf(:,:))
+ IF ( iom_use( TRIM(cvarfwf3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarfwf3d) , pfwf(:,:))
IF ( iom_use( TRIM(cvarqoce3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqoce3d) , pqoce(:,:))
IF ( iom_use( TRIM(cvarqlat3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqlat3d) , pqoce(:,:))
IF ( iom_use( TRIM(cvarqhc3d) ) ) CALL isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, TRIM(cvarqhc3d) , pqhc (:,:))
!
END SUBROUTINE
- SUBROUTINE isf_diags_2dto3d(Kmm, ktop, kbot, phtbl, pfrac, cdvar, pvar2d)
+ SUBROUTINE isf_diags_2dto3d( Kmm, ktop, kbot, phtbl, pfrac, cdvar, pvar2d )
!!---------------------------------------------------------------------
!! *** ROUTINE isf_diags_2dto3d ***
!!
@@ -82,25 +81,23 @@ CONTAINS
!! (ie uniformaly spread into the top boundary layer or parametrisation layer)
!!
!!----------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- !!-------------------------- IN -------------------------------------
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
- INTEGER , DIMENSION(jpi,jpj), INTENT(in) :: ktop , kbot ! top and bottom level of the tbl
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac ! thickness of the tbl and fraction of last cell affected by the tbl
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d ! 2d var to map in 3d
- CHARACTER(LEN=*), INTENT(in) :: cdvar
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ INTEGER , DIMENSION(A2D(0)), INTENT(in) :: ktop , kbot ! top and bottom level of the tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: phtbl, pfrac ! thickness of the tbl and fraction of last cell affected by the tbl
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pvar2d ! 2d var to map in 3d
+ CHARACTER(LEN=*), INTENT(in) :: cdvar
!!---------------------------------------------------------------------
- INTEGER :: ji, jj, jk ! loop indices
- INTEGER :: ikt, ikb ! top and bottom level of the tbl
- REAL(wp), DIMENSION(jpi,jpj) :: zvar2d !
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvar3d ! 3d var to output
+ INTEGER :: ji, jj, jk ! loop indices
+ INTEGER :: ikt, ikb ! top and bottom level of the tbl
+ REAL(wp), DIMENSION(A2D(0)) :: zvar2d !
+ REAL(wp), DIMENSION(A2D(0),jpk) :: zvar3d ! 3d var to output
!!---------------------------------------------------------------------
!
! compute 3d output
zvar2d(:,:) = pvar2d(:,:) / phtbl(:,:)
zvar3d(:,:,:) = 0._wp
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
ikt = ktop(ji,jj)
ikb = kbot(ji,jj)
DO jk = ikt, ikb - 1
@@ -109,7 +106,7 @@ CONTAINS
zvar3d(ji,jj,ikb) = zvar2d(ji,jj) * e3t(ji,jj,ikb,Kmm) * pfrac(ji,jj)
END_2D
!
- CALL iom_put( TRIM(cdvar) , zvar3d(:,:,:))
+ CALL iom_put( TRIM(cdvar) , zvar3d(:,:,:) )
!
END SUBROUTINE isf_diags_2dto3d
diff --git a/src/OCE/ISF/isfdynatf.F90 b/src/OCE/ISF/isfdynatf.F90
index 3be8abf8e..9ecc95081 100644
--- a/src/OCE/ISF/isfdynatf.F90
+++ b/src/OCE/ISF/isfdynatf.F90
@@ -7,12 +7,13 @@ MODULE isfdynatf
!!-------------------------------------------------------------------------
!!-------------------------------------------------------------------------
- !! isfnxt : apply correction needed for the ice shelf to ensure conservation
+ !! isfdynatf : apply correction needed for the ice shelf to ensure conservation
!!-------------------------------------------------------------------------
USE isf_oce
USE phycst , ONLY: r1_rho0 ! physical constant
+ USE par_oce ! ocean space and time domain
USE dom_oce ! time and space domain
USE oce, ONLY : ssh ! sea-surface height for qco substitution
@@ -26,7 +27,6 @@ MODULE isfdynatf
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
-
CONTAINS
SUBROUTINE isf_dynatf ( kt, Kmm, pe3t_f, pcoef )
@@ -36,56 +36,45 @@ CONTAINS
!! ** Purpose : compute the ice shelf volume filter correction for cavity, param, ice sheet coupling case
!!
!!-------------------------- OUT -------------------------------------
- INTEGER , INTENT(in ) :: kt ! ocean time step
- INTEGER , INTENT(in ) :: Kmm ! ocean time level index
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3t_f ! time filtered scale factor to be corrected
+ INTEGER , INTENT(in ) :: kt ! ocean time step
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3t_f ! time filtered scale factor to be corrected
!
- REAL(wp) , INTENT(in ) :: pcoef ! rn_atfp * rn_Dt * r1_rho0
+ REAL(wp) , INTENT(in ) :: pcoef ! rn_atfp * rn_Dt * r1_rho0
!!--------------------------------------------------------------------
- INTEGER :: jk ! loop index
+ INTEGER :: ji, jj, jk ! loop index
+ REAL(wp) :: ztmp
!!--------------------------------------------------------------------
!
! ice shelf cavity
- IF ( ln_isfcav_mlt ) CALL isf_dynatf_mlt(Kmm, pe3t_f, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, pcoef)
+ IF( ln_isfcav_mlt ) THEN
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ ztmp = pcoef * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) * r1_rho0
+ !
+ DO jk = 1, jpkm1
+ pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) + tmask(ji,jj,jk) * ztmp * e3t(ji,jj,jk,Kmm)
+ END DO
+ END_2D
+ ENDIF
!
! ice shelf parametrised
- IF ( ln_isfpar_mlt ) CALL isf_dynatf_mlt(Kmm, pe3t_f, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, pcoef)
+ IF( ln_isfpar_mlt ) THEN
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ ztmp = pcoef * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) ) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) * r1_rho0
+ !
+ DO jk = 1, jpkm1
+ pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) + tmask(ji,jj,jk) * ztmp * e3t(ji,jj,jk,Kmm)
+ END DO
+ END_2D
+ ENDIF
!
- IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
- DO jk = 1, jpkm1
- pe3t_f(:,:,jk) = pe3t_f(:,:,jk) - pcoef * risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
- END DO
+ ! if coupled
+ IF( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
+ DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - pcoef * risfcpl_vol(ji,jj,jk) * r1_e1e2t(ji,jj)
+ END_3D
END IF
!
END SUBROUTINE isf_dynatf
- SUBROUTINE isf_dynatf_mlt ( Kmm, pe3t_f, ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, pcoef )
- !!--------------------------------------------------------------------
- !! *** ROUTINE isf_dynatf_mlt ***
- !!
- !! ** Purpose : compute the ice shelf volume filter correction for cavity or param
- !!
- !!-------------------------- IN -------------------------------------
- INTEGER , INTENT(in ) :: Kmm ! ocean time level index
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3t_f ! time-filtered scale factor to be corrected
- INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot ! top and bottom level of tbl
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfrac, phtbl ! fraction of bottom cell included in tbl, tbl thickness
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf , pfwf_b ! now/before fwf
- REAL(wp), INTENT(in ) :: pcoef ! rn_atfp * rn_Dt * r1_rho0
- !!----------------------------------------------------------------------
- INTEGER :: ji,jj,jk
- REAL(wp), DIMENSION(jpi,jpj) :: zfwfinc
- !!----------------------------------------------------------------------
- !
- ! compute fwf conservation correction
- zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / ( ht(:,:) + 1._wp - ssmask(:,:) ) * r1_rho0
- !
- ! add the increment
- DO jk = 1, jpkm1
- pe3t_f(:,:,jk) = pe3t_f(:,:,jk) + tmask(:,:,jk) * zfwfinc(:,:) &
- & * e3t(:,:,jk,Kmm)
- END DO
- !
- END SUBROUTINE isf_dynatf_mlt
-
END MODULE isfdynatf
diff --git a/src/OCE/ISF/isfhdiv.F90 b/src/OCE/ISF/isfhdiv.F90
index 58275063b..2f74bc536 100644
--- a/src/OCE/ISF/isfhdiv.F90
+++ b/src/OCE/ISF/isfhdiv.F90
@@ -15,6 +15,7 @@ MODULE isfhdiv
USE isf_oce ! ice shelf
+ USE par_oce ! ocean space and time domain
USE dom_oce ! time and space domain
USE phycst , ONLY: r1_rho0 ! physical constant
USE in_out_manager !
@@ -24,10 +25,10 @@ MODULE isfhdiv
PRIVATE
PUBLIC isf_hdiv
+
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
-
CONTAINS
SUBROUTINE isf_hdiv( kt, Kmm, phdiv )
@@ -39,25 +40,29 @@ CONTAINS
!! increment)
!!
!!----------------------------------------------------------------------
+ INTEGER, INTENT(in) :: kt
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdiv ! horizontal divergence
!!----------------------------------------------------------------------
- INTEGER, INTENT(in) :: kt
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
!
IF ( ln_isf ) THEN
!
#if defined key_RK3
! ice shelf cavity contribution (RK3)
- IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, phdiv)
+ IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, &
+ & fwfisf_cav, phdiv)
!
! ice shelf parametrisation contribution (RK3)
- IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, phdiv)
+ IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, &
+ fwfisf_par, phdiv)
#else
! ice shelf cavity contribution (MLF)
- IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, phdiv, fwfisf_cav_b)
+ IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, &
+ & fwfisf_cav, phdiv, fwfisf_cav_b)
!
! ice shelf parametrisation contribution (MLF)
- IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, phdiv, fwfisf_par_b)
+ IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, &
+ fwfisf_par, phdiv, fwfisf_par_b)
#endif
!
! ice sheet coupling contribution
@@ -81,7 +86,7 @@ CONTAINS
END SUBROUTINE isf_hdiv
- SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, phdiv, pfwf_b)
+ SUBROUTINE isf_hdiv_mlt( ktop, kbot, phtbl, pfrac, pfwf, phdiv, pfwf_b )
!!----------------------------------------------------------------------
!! *** SUBROUTINE sbc_isf_div ***
!!
@@ -92,45 +97,42 @@ CONTAINS
!!
!! ** Action : phdivn increased by the ice shelf outflow
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: phdiv
- !!----------------------------------------------------------------------
- INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfrac, phtbl
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
- REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in ) :: pfwf_b
+ INTEGER , DIMENSION(A2D(0)) , INTENT(in ) :: ktop , kbot
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: pfrac, phtbl
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: phdiv
+ REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in ) :: pfwf_b
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikt, ikb
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zhdiv
+ REAL(wp) :: zhdiv
!!----------------------------------------------------------------------
!
!== fwf distributed over several levels ==!
!
- ! compute integrated divergence correction
- DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+ ! update divergence at each level affected by ice shelf top boundary layer
+ DO_2D_OVR( 0, 0, 0, 0 )
+ ! compute integrated divergence correction
#if defined key_RK3
- zhdiv(ji,jj) = pfwf(ji,jj) * r1_rho0 / phtbl(ji,jj)
+ zhdiv = pfwf(ji,jj) * r1_rho0 / phtbl(ji,jj)
#else
- zhdiv(ji,jj) = 0.5_wp * ( pfwf(ji,jj) + pfwf_b(ji,jj) ) * r1_rho0 / phtbl(ji,jj)
+ zhdiv = 0.5_wp * ( pfwf(ji,jj) + pfwf_b(ji,jj) ) * r1_rho0 / phtbl(ji,jj)
#endif
- END_2D
- !
- ! update divergence at each level affected by ice shelf top boundary layer
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+ !
ikt = ktop(ji,jj)
ikb = kbot(ji,jj)
! level fully include in the ice shelf boundary layer
DO jk = ikt, ikb - 1
- phdiv(ji,jj,jk) = phdiv(ji,jj,jk) - zhdiv(ji,jj)
+ phdiv(ji,jj,jk) = phdiv(ji,jj,jk) - zhdiv
END DO
! level partially include in ice shelf boundary layer
- phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) - zhdiv(ji,jj) * pfrac(ji,jj)
+ phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) - zhdiv * pfrac(ji,jj)
END_2D
!
END SUBROUTINE isf_hdiv_mlt
- SUBROUTINE isf_hdiv_cpl(Kmm, pqvol, phdiv)
+ SUBROUTINE isf_hdiv_cpl( Kmm, pqvol, phdiv )
!!----------------------------------------------------------------------
!! *** SUBROUTINE isf_hdiv_cpl ***
!!
@@ -143,17 +145,15 @@ CONTAINS
!! ** Action : phdivn increased by the ice shelf outflow
!!
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
- !!----------------------------------------------------------------------
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pqvol
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pqvol
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
!!----------------------------------------------------------------------
- INTEGER :: ji, jj, jk
+ INTEGER :: ji, jj, jk
!!----------------------------------------------------------------------
!
- DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpk )
- phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + pqvol(ji,jj,jk) * r1_e1e2t(ji,jj) &
- & / e3t(ji,jj,jk,Kmm)
+ DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
+ phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + pqvol(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
END SUBROUTINE isf_hdiv_cpl
diff --git a/src/OCE/ISF/isfload.F90 b/src/OCE/ISF/isfload.F90
index a069cdf5d..ab1f9a3dc 100644
--- a/src/OCE/ISF/isfload.F90
+++ b/src/OCE/ISF/isfload.F90
@@ -12,6 +12,7 @@ MODULE isfload
USE isf_oce, ONLY: cn_isfload, rn_isfload_T, rn_isfload_S ! ice shelf variables
+ USE par_oce ! ocean space and time domain
USE dom_oce ! vertical scale factor
USE eosbn2 , ONLY: eos ! eos routine
@@ -81,7 +82,8 @@ CONTAINS
!!----------------------------------------------------------------------
!
! !- assume water displaced by the ice shelf is at T=rn_isfload_T and S=rn_isfload_S (rude)
- zts_top(:,:,jp_tem) = rn_isfload_T ; zts_top(:,:,jp_sal) = rn_isfload_S
+ zts_top(:,:,jp_tem) = rn_isfload_T
+ zts_top(:,:,jp_sal) = rn_isfload_S
!
DO jk = 1, jpk !- compute density of the water displaced by the ice shelf
#if defined key_qco && key_isf
diff --git a/src/OCE/ISF/isfpar.F90 b/src/OCE/ISF/isfpar.F90
index 31f363687..bc4787f73 100644
--- a/src/OCE/ISF/isfpar.F90
+++ b/src/OCE/ISF/isfpar.F90
@@ -16,14 +16,14 @@ MODULE isfpar
!!----------------------------------------------------------------------
USE isf_oce ! ice shelf
!
- USE isfrst , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
- USE isftbl , ONLY: isf_tbl_ktop, isf_tbl_lvl ! ice shelf top boundary layer properties subroutine
+ USE isfrst , ONLY: isfrst_read ! ice shelf restart read/write subroutine
+ USE isftbl , ONLY: isf_tbl_ktop ! ice shelf top boundary layer properties subroutine
USE isfparmlt, ONLY: isfpar_mlt ! ice shelf melt formulation subroutine
USE isfdiags , ONLY: isf_diags_flx ! ice shelf diags subroutine
USE isfutils , ONLY: debug, read_2dcstdta ! ice shelf debug subroutine
!
USE dom_oce , ONLY: bathy ! ocean space and time domain
- USE par_oce , ONLY: jpi,jpj ! ocean space and time domain
+ USE par_oce ! ocean space and time domain
USE phycst , ONLY: r1_rho0_rcp ! physical constants
!
USE in_out_manager ! I/O manager
@@ -44,7 +44,7 @@ MODULE isfpar
!!----------------------------------------------------------------------
CONTAINS
- SUBROUTINE isf_par( kt, Kmm, ptsc, pqfwf )
+ SUBROUTINE isf_par( kt, Kmm, ptsc, pfwf )
!!---------------------------------------------------------------------
!! *** ROUTINE isf_par ***
!!
@@ -60,28 +60,26 @@ CONTAINS
!! ** Convention : all fluxes are from isf to oce
!!
!!---------------------------------------------------------------------
- !!-------------------------- OUT --------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pqfwf
- REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc
- !!-------------------------- IN --------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(inout) :: pfwf
+ REAL(wp), DIMENSION(A2D(0),jpts), INTENT(inout) :: ptsc
!!---------------------------------------------------------------------
INTEGER :: ji, jj
- REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
+ REAL(wp), DIMENSION(A2D(0)) :: zqoce, zqhc, zqlat, zqh
!!---------------------------------------------------------------------
!
! compute heat content, latent heat and melt fluxes (2d)
- CALL isfpar_mlt( kt, Kmm, zqhc, zqoce, pqfwf )
+ CALL isfpar_mlt( kt, Kmm, zqhc, zqoce, pfwf )
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
! compute heat and water flux (from isf to oce)
- pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_par(ji,jj)
+ pfwf(ji,jj) = pfwf(ji,jj) * mskisf_par(ji,jj)
zqoce(ji,jj) = zqoce(ji,jj) * mskisf_par(ji,jj)
zqhc (ji,jj) = zqhc(ji,jj) * mskisf_par(ji,jj)
!
! compute latent heat flux (from isf to oce)
- zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2)
+ zqlat(ji,jj) = - pfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2)
!
! total heat flux (from isf to oce)
zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
@@ -91,18 +89,13 @@ CONTAINS
END_2D
!
! output fluxes
- CALL isf_diags_flx( Kmm, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc)
+ CALL isf_diags_flx( Kmm, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pfwf, zqoce, zqlat, zqhc )
!
-#if ! defined key_RK3
- ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
- IF (lrst_oce) CALL isfrst_write(kt, 'par', ptsc, pqfwf)
- !
-#endif
IF ( ln_isfdebug ) THEN
IF(lwp) WRITE(numout,*)
CALL debug('isf_par: ptsc T',ptsc(:,:,1))
CALL debug('isf_par: ptsc S',ptsc(:,:,2))
- CALL debug('isf_par: pqfwf fwf',pqfwf(:,:))
+ CALL debug('isf_par: pfwf fwf',pfwf(:,:))
IF(lwp) WRITE(numout,*)
END IF
!
@@ -115,34 +108,41 @@ CONTAINS
!! ** Purpose : initialisation of the variable needed for the parametrisation of ice shelf melt
!!
!!----------------------------------------------------------------------
- INTEGER :: ierr
- REAL(wp), DIMENSION(jpi,jpj) :: ztblmax, ztblmin
+ INTEGER :: ierr
+ INTEGER :: ji, jj ! dummy loop indices
+ REAL(wp), DIMENSION(A2D(0)) :: ztblmax, ztblmin
!!----------------------------------------------------------------------
!
- ! allocation
+ !==============
+ ! 0: allocation
+ !==============
CALL isf_alloc_par()
!
- ! initialisation
- misfkt_par(:,:) = 1 ; misfkb_par(:,:) = 1
- rhisf_tbl_par(:,:) = 1e-20 ; rfrac_tbl_par(:,:) = 0.0_wp
+ !==================
+ ! 1: initialisation
+ !==================
+ DO_2D( 0, 0, 0, 0 )
+ misfkt_par (ji,jj) = 1
+ misfkb_par (ji,jj) = 1
+ rhisf_tbl_par(ji,jj) = 1e-20
+ rfrac_tbl_par(ji,jj) = 0.0_wp
+ END_2D
!
! define isf tbl tickness, top and bottom indice
- CALL read_2dcstdta(TRIM(sn_isfpar_zmax%clname), TRIM(sn_isfpar_zmax%clvar), ztblmax)
- CALL read_2dcstdta(TRIM(sn_isfpar_zmin%clname), TRIM(sn_isfpar_zmin%clvar), ztblmin)
- !
- ! mask ice shelf parametrisation location
- ztblmax(:,:) = ztblmax(:,:) * ssmask(:,:)
- ztblmin(:,:) = ztblmin(:,:) * ssmask(:,:)
+ CALL read_2dcstdta( TRIM(sn_isfpar_zmax%clname), TRIM(sn_isfpar_zmax%clvar), ztblmax )
+ CALL read_2dcstdta( TRIM(sn_isfpar_zmin%clname), TRIM(sn_isfpar_zmin%clvar), ztblmin )
!
- ! if param used under an ice shelf overwrite ztblmin by the ice shelf draft
- WHERE ( risfdep > 0._wp .AND. ztblmin > 0._wp )
- ztblmin(:,:) = risfdep(:,:)
- END WHERE
- !
- ! ensure ztblmax <= bathy
- WHERE ( ztblmax(:,:) > bathy(:,:) )
- ztblmax(:,:) = bathy(:,:)
- END WHERE
+ DO_2D( 0, 0, 0, 0 )
+ ! mask ice shelf parametrisation location
+ ztblmax(ji,jj) = ztblmax(ji,jj) * ssmask(ji,jj)
+ ztblmin(ji,jj) = ztblmin(ji,jj) * ssmask(ji,jj)
+ !
+ ! if param used under an ice shelf overwrite ztblmin by the ice shelf draft
+ IF( risfdep(ji,jj) > 0._wp .AND. ztblmin(ji,jj) > 0._wp ) ztblmin(ji,jj) = risfdep(ji,jj)
+ !
+ ! ensure ztblmax <= bathy
+ ztblmax(ji,jj) = MIN( ztblmax(ji,jj), bathy(ji,jj) )
+ END_2D
!
! compute ktop and update ztblmin to gdepw_0(misfkt_par)
CALL isf_tbl_ktop(ztblmin, misfkt_par) ! out: misfkt_par
@@ -158,16 +158,22 @@ CONTAINS
END WHERE
!
#if ! defined key_RK3
+ !================
+ ! 2: read restart
+ !================
! MLF: read par variable from restart
- IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
+ IF ( ln_rstart ) CALL isfrst_read( 'par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b )
#endif
!
+ !==========================================
+ ! 3: specific allocation and initialisation (depending of scheme choice)
+ !==========================================
SELECT CASE ( TRIM(cn_isfpar_mlt) )
!
CASE ( 'spe' )
!
ALLOCATE( sf_isfpar_fwf(1), STAT=ierr )
- ALLOCATE( sf_isfpar_fwf(1)%fnow(jpi,jpj,1), sf_isfpar_fwf(1)%fdta(jpi,jpj,1,2) )
+ ALLOCATE( sf_isfpar_fwf(1)%fnow(A2D(0),1), sf_isfpar_fwf(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_isfdir, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )
!
IF(lwp) WRITE(numout,*)
@@ -185,7 +191,7 @@ CONTAINS
CASE ( 'oasis' )
!
ALLOCATE( sf_isfpar_fwf(1), STAT=ierr )
- ALLOCATE( sf_isfpar_fwf(1)%fnow(jpi,jpj,1), sf_isfpar_fwf(1)%fdta(jpi,jpj,1,2) )
+ ALLOCATE( sf_isfpar_fwf(1)%fnow(A2D(0),1), sf_isfpar_fwf(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_isfdir, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )
!
IF(lwp) WRITE(numout,*)
diff --git a/src/OCE/ISF/isfparmlt.F90 b/src/OCE/ISF/isfparmlt.F90
index 237e3bee3..fd7a65777 100644
--- a/src/OCE/ISF/isfparmlt.F90
+++ b/src/OCE/ISF/isfparmlt.F90
@@ -11,6 +11,7 @@ MODULE isfparmlt
USE isftbl , ONLY: isf_tbl ! ice shelf depth average
USE isfutils,ONLY: debug ! debug subroutine
+ USE par_oce ! ocean space and time domain
USE dom_oce ! ocean space and time domain
USE oce , ONLY: ts ! ocean dynamics and tracers
USE phycst , ONLY: rcp, rho0 ! physical constants
@@ -31,6 +32,7 @@ MODULE isfparmlt
!! * Substitutions
# include "domzgr_substitute.h90"
+# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
@@ -38,11 +40,8 @@ MODULE isfparmlt
!!----------------------------------------------------------------------
CONTAINS
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
- SUBROUTINE isfpar_mlt( kt, Kmm, pqhc, pqoce, pqfwf )
+ SUBROUTINE isfpar_mlt( kt, Kmm, pqhc, pqoce, pfwf )
!!---------------------------------------------------------------------
!! *** ROUTINE isfpar_mlt ***
!!
@@ -53,21 +52,26 @@ CONTAINS
!! 1 : Specified melt flux
!! 2 : Beckmann & Goose parameterization
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf, pqoce, pqhc ! fresh water, ice-ocean heat and heat content fluxes
- !!-------------------------- IN -------------------------------------
+ !!---------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pfwf, pqoce, pqhc ! fresh water, ice-ocean heat and heat content fluxes
!!---------------------------------------------------------------------
!
! Choose among the available ice shelf parametrisation
SELECT CASE ( cn_isfpar_mlt )
CASE ( 'spe' ) ! specified runoff in depth (Mathiot et al., 2017 in preparation)
- CALL isfpar_mlt_spe(kt, Kmm, pqhc, pqoce, pqfwf)
+ !
+ CALL isfpar_mlt_spe( kt, Kmm, pqhc, pqoce, pfwf )
+ !
CASE ( 'bg03' ) ! Beckmann and Goosse parametrisation
- CALL isfpar_mlt_bg03(kt, Kmm, pqhc, pqoce, pqfwf)
+ !
+ CALL isfpar_mlt_bg03( kt, Kmm, pqhc, pqoce, pfwf )
+ !
CASE ( 'oasis' )
- CALL isfpar_mlt_oasis( kt, Kmm, pqhc, pqoce, pqfwf)
+ !
+ CALL isfpar_mlt_oasis( kt, Kmm, pqhc, pqoce, pfwf )
+ !
CASE DEFAULT
CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)')
END SELECT
@@ -76,32 +80,28 @@ CONTAINS
IF(lwp) WRITE(numout,*) ''
CALL debug( 'isfpar_mlt qhc :', pqhc (:,:) )
CALL debug( 'isfpar_mlt qoce :', pqoce(:,:) )
- CALL debug( 'isfpar_mlt qfwf :', pqfwf(:,:) )
+ CALL debug( 'isfpar_mlt fwf :', pfwf(:,:) )
IF(lwp) WRITE(numout,*) ''
END IF
!
END SUBROUTINE isfpar_mlt
-! -------------------------------------------------------------------------------------------------------
-! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
-! -------------------------------------------------------------------------------------------------------
- SUBROUTINE isfpar_mlt_spe(kt, Kmm, pqhc, pqoce, pqfwf)
+ SUBROUTINE isfpar_mlt_spe( kt, Kmm, pqhc, pqoce, pfwf )
!!---------------------------------------------------------------------
!! *** ROUTINE isfpar_mlt_spe ***
!!
!! ** Purpose : prescribed ice shelf melting in case ice shelf cavities are closed.
!! data read into a forcing files.
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf, pqoce ! fresh water and ice-ocean heat fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER, INTENT(in) :: kt
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
!!--------------------------------------------------------------------
- INTEGER :: jk
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz
+ INTEGER, INTENT(in) :: kt
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pqhc, pfwf, pqoce ! fresh water and ice-ocean heat fluxes
+ !!--------------------------------------------------------------------
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztfrz3d, ztmp
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz
!!--------------------------------------------------------------------
!
! 0. ------------Read specified fwf from isf to oce
@@ -109,20 +109,26 @@ CONTAINS
!
! compute ptfrz
! 1. ------------Mean freezing point
- DO jk = 1,jpk
- CALL eos_fzp(ts(:,:,jk,jp_sal,Kmm), ztfrz3d(:,:,jk), gdept(:,:,jk,Kmm))
- END DO
- CALL isf_tbl(Kmm, ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ ztmp(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
+ END_3D
+ CALL eos_fzp( ts(A2D(0),:,jp_sal,Kmm), ztfrz3d(:,:,:), ztmp )
+ !
+ CALL isf_tbl( Kmm, ztfrz3d, 'T', misfkt_par, rhisf_tbl_par, & ! <<== in
+ & ztfrz, & ! ==>> out
+ & misfkb_par, rfrac_tbl_par ) ! <<== in (optional)
!
- pqfwf(:,:) = sf_isfpar_fwf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean/ice shelf flux assume to be equal to latent heat flux ( > 0 from isf to oce)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
+ DO_2D( 0, 0, 0, 0 )
+ pfwf(ji,jj) = sf_isfpar_fwf(1)%fnow(ji,jj,1) ! fresh water flux from the isf (fwfisf <0 mean melting) ( > 0 from isf to oce)
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocean/ice shelf flux assume to be equal to latent heat flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ END_2D
!
CALL iom_put('isftfrz_par', ztfrz(:,:) * mskisf_par(:,:) )
!
END SUBROUTINE isfpar_mlt_spe
- SUBROUTINE isfpar_mlt_bg03(kt, Kmm, pqhc, pqoce, pqfwf)
+ SUBROUTINE isfpar_mlt_bg03(kt, Kmm, pqhc, pqoce, pfwf)
!!---------------------------------------------------------------------
!! *** ROUTINE isfpar_mlt_bg03 ***
!!
@@ -137,45 +143,48 @@ CONTAINS
!! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
!! interaction for climate models", Ocean Modelling 5(2003) 157-170.
!!----------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf, pqoce ! fresh water and ice-ocean heat fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER, INTENT(in) :: kt
- INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ INTEGER, INTENT(in) :: kt
+ INTEGER, INTENT(in) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pqhc, pfwf, pqoce ! fresh water and ice-ocean heat fluxes
!!--------------------------------------------------------------------
- INTEGER :: jk
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d ! freezing point
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing point
- REAL(wp), DIMENSION(jpi,jpj) :: ztavg ! temperature avg
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztfrz3d, ztmp
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! freezing point
+ REAL(wp), DIMENSION(A2D(0)) :: ztavg ! temperature avg
!!----------------------------------------------------------------------
!
! 0. ------------Mean freezing point
- DO jk = 1,jpk
- CALL eos_fzp(ts(:,:,jk,jp_sal,Kmm), ztfrz3d(:,:,jk), gdept(:,:,jk,Kmm))
- END DO
- CALL isf_tbl(Kmm, ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ ztmp(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
+ END_3D
+ CALL eos_fzp( ts(A2D(0),:,jp_sal,Kmm), ztfrz3d(:,:,:), ztmp )
+
+ CALL isf_tbl( Kmm, ztfrz3d, 'T', misfkt_par, rhisf_tbl_par, & ! <<== in
+ & ztfrz, & ! ==>> out
+ & misfkb_par, rfrac_tbl_par ) ! <<== in (optional)
!
! 1. ------------Mean temperature
- CALL isf_tbl(Kmm, ts(:,:,:,jp_tem,Kmm), ztavg, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+ CALL isf_tbl( Kmm, ts(:,:,:,jp_tem,Kmm), 'T', misfkt_par, rhisf_tbl_par, & ! <<== in
+ & ztavg, & ! ==>> out
+ & misfkb_par, rfrac_tbl_par ) ! <<== in (optional)
!
! 2. ------------Net heat flux and fresh water flux due to the ice shelf
- pqfwf(:,:) = rho0 * rcp * rn_isfpar_bg03_gt0 * risfLeff(:,:) * e1t(:,:) * (ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:) / rLfusisf ! ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean/ice shelf flux assume to be equal to latent heat flux ( > 0 from isf to oce)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
+ DO_2D( 0, 0, 0, 0 )
+ pfwf(ji,jj) = rho0 * rcp * rn_isfpar_bg03_gt0 * risfLeff(ji,jj) * e1t(ji,jj) * ( ztavg(ji,jj) - ztfrz(ji,jj) ) &
+ & * r1_e1e2t(ji,jj) / rLfusisf ! ( > 0 from isf to oce)
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocean/ice shelf flux assume to be equal to latent heat flux ( > 0 from isf to oce)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ END_2D
!
! 3. ------------BG03 output
- ! output ttbl
- CALL iom_put('ttbl_par', ztavg(:,:) * mskisf_par(:,:) )
- !
- ! output thermal driving
- CALL iom_put('isfthermald_par',( ztavg(:,:) - ztfrz(:,:) ) * mskisf_par(:,:))
- !
- ! output freezing point used to define the thermal driving and heat content fluxes
- CALL iom_put('isftfrz_par', ztfrz(:,:) * mskisf_par(:,:) )
+ CALL iom_put('ttbl_par', ztavg(:,:) * mskisf_par(:,:) ) ! ttbl
+ CALL iom_put('isfthermald_par',( ztavg(:,:) - ztfrz(:,:) ) * mskisf_par(:,:) ) ! thermal driving
+ CALL iom_put('isftfrz_par', ztfrz(:,:) * mskisf_par(:,:) ) ! freezing point
+
!
END SUBROUTINE isfpar_mlt_bg03
- SUBROUTINE isfpar_mlt_oasis(kt, Kmm, pqhc , pqoce, pqfwf )
+ SUBROUTINE isfpar_mlt_oasis( kt, Kmm, pqhc , pqoce, pfwf )
!!----------------------------------------------------------------------
!! *** ROUTINE isfpar_mlt_oasis ***
!!
@@ -186,27 +195,29 @@ CONTAINS
!! - scale fwf and compute heat fluxes
!!
!!---------------------------------------------------------------------
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat content, latent heat and fwf fluxes
- !!-------------------------- IN -------------------------------------
- INTEGER , INTENT(in ) :: kt ! current time step
- INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ INTEGER , INTENT(in ) :: kt ! current time step
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pqhc, pqoce, pfwf ! heat content, latent heat and fwf fluxes
!!--------------------------------------------------------------------
- INTEGER :: jk ! loop index
- REAL(wp) :: zfwf_fld, zfwf_oasis ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
- REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! tbl freezing temperature
- REAL(wp), DIMENSION(jpi,jpj) :: zfwf ! 2d fwf map after scaling
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zfwf_fld, zfwf_oasis ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
+ REAL(wp), DIMENSION(A2D(0)) :: ztfrz ! tbl freezing temperature
+ REAL(wp), DIMENSION(A2D(0)) :: zfwf ! 2d fwf map after scaling
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ztfrz3d, ztmp
!!--------------------------------------------------------------------
!
! 0. ------------Read specified runoff
CALL fld_read ( kt, 1, sf_isfpar_fwf )
!
! 1. ------------Mean freezing point (needed for heat content flux)
- DO jk = 1,jpk
- CALL eos_fzp(ts(:,:,jk,jp_sal,Kmm), ztfrz3d(:,:,jk), gdept(:,:,jk,Kmm))
- END DO
- CALL isf_tbl(Kmm, ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ ztmp(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
+ END_3D
+ CALL eos_fzp( ts(A2D(0),:,jp_sal,Kmm), ztfrz3d(:,:,:), ztmp )
+ !
+ CALL isf_tbl( Kmm, ztfrz3d, 'T', misfkt_par, rhisf_tbl_par, & ! <<== in
+ & ztfrz, & ! ==>> out
+ & misfkb_par, rfrac_tbl_par ) ! <<== in (optional)
!
! 2. ------------Scale isf melt pattern with total amount from oasis
! ice shelf 2d map of fwf from isf to oce
@@ -214,22 +225,24 @@ CONTAINS
!
! compute glob sum from input file
! (PM) should we consider delay sum as in fwb ? (it will offset by 1 time step if I understood well)
- zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
+ zfwf_fld = glob_sum('isfcav_mlt', e1e2t(A2D(0)) * zfwf(:,:))
!
! compute glob sum from atm->oce ice shelf fwf
! (PM) should we consider delay sum as in fwb ?
- zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
+ zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(A2D(0)) * fwfisf_oasis(:,:))
!
! scale fwf
zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
!
! 3. -----------Define fwf and qoce
! ocean heat flux is assume to be equal to the latent heat
- pqfwf(:,:) = zfwf(:,:) ! fwf ( > 0 from isf to oce)
- pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( > 0 from isf to oce) (assumed to be the latent heat flux)
- pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
+ DO_2D( 0, 0, 0, 0 )
+ pfwf(ji,jj) = zfwf(ji,jj) ! fwf ( > 0 from isf to oce)
+ pqoce(ji,jj) = - pfwf(ji,jj) * rLfusisf ! ocean heat flux ( > 0 from isf to oce) (assumed to be the latent heat flux)
+ pqhc (ji,jj) = pfwf(ji,jj) * ztfrz(ji,jj) * rcp ! heat content flux ( > 0 from isf to oce)
+ END_2D
!
- CALL iom_put('isftfrz_par', ztfrz )
+ CALL iom_put('isftfrz_par', ztfrz(:,:) * mskisf_par(:,:) )
!
END SUBROUTINE isfpar_mlt_oasis
diff --git a/src/OCE/ISF/isfrst.F90 b/src/OCE/ISF/isfrst.F90
index 696605ff1..a6302a54a 100644
--- a/src/OCE/ISF/isfrst.F90
+++ b/src/OCE/ISF/isfrst.F90
@@ -10,7 +10,7 @@ MODULE isfrst
!! isfrst : read/write iceshelf variables in/from restart
!!----------------------------------------------------------------------
!
- USE par_oce, ONLY: jpi,jpj,jpk,jpts ! time and space domain
+ USE par_oce ! time and space domain
!
USE in_out_manager ! I/O manager
USE iom ! I/O library
@@ -21,6 +21,8 @@ MODULE isfrst
PUBLIC isfrst_read, isfrst_write ! iceshelf restart read and write
+ !! * Substitutions
+# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
@@ -33,13 +35,12 @@ CONTAINS
!!
!! isfrst_read : read iceshelf variables from restart
!!
- !!-------------------------- OUT --------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pfwf_b
- REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT( out) :: ptsc_b
- !!-------------------------- IN --------------------------------------
- CHARACTER(LEN=3) , INTENT(in ) :: cdisf
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
- REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: ptsc
+ !!----------------------------------------------------------------------
+ CHARACTER(LEN=3) , INTENT(in ) :: cdisf
+ REAL(wp), DIMENSION(A2D(0),jpts), INTENT(in ) :: ptsc
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
+ REAL(wp), DIMENSION(A2D(0),jpts), INTENT( out) :: ptsc_b
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pfwf_b
!!----------------------------------------------------------------------
CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
!!----------------------------------------------------------------------
@@ -68,11 +69,11 @@ CONTAINS
!!
!! isfrst_write : write iceshelf variables in restart
!!
- !!-------------------------- IN --------------------------------------
- INTEGER , INTENT(in ) :: kt
- CHARACTER(LEN=3) , INTENT(in ) :: cdisf
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
- REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: ptsc
+ !!---------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt
+ CHARACTER(LEN=3) , INTENT(in ) :: cdisf
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
+ REAL(wp), DIMENSION(A2D(0),jpts), INTENT(in ) :: ptsc
!!---------------------------------------------------------------------
CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
!!---------------------------------------------------------------------
diff --git a/src/OCE/ISF/isfstp.F90 b/src/OCE/ISF/isfstp.F90
index edc5f92aa..f742fee6b 100644
--- a/src/OCE/ISF/isfstp.F90
+++ b/src/OCE/ISF/isfstp.F90
@@ -13,20 +13,23 @@ MODULE isfstp
!! isfstp : compute iceshelf melt and heat flux
!!----------------------------------------------------------------------
USE isf_oce ! isf variables
+ USE isfrst , ONLY: isfrst_write ! ice shelf restart read/write subroutine
USE isfload, ONLY: isf_load ! ice shelf load
USE isftbl , ONLY: isf_tbl_lvl ! ice shelf boundary layer
USE isfpar , ONLY: isf_par, isf_par_init ! ice shelf parametrisation
USE isfcav , ONLY: isf_cav, isf_cav_init ! ice shelf cavity
USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
+ USE oce , ONLY: ssh
+ USE par_oce ! ocean space and time domain
USE dom_oce ! ocean space and time domain
- USE oce , ONLY: ssh ! sea surface height
USE domvvl, ONLY: ln_vvl_zstar ! zstar logical
USE zdfdrg, ONLY: r_Cdmin_top, r_ke0_top ! vertical physics: top/bottom drag coef.
!
USE lib_mpp, ONLY: ctl_stop, ctl_nam
USE fldread, ONLY: FLD, FLD_N
USE in_out_manager ! I/O manager
+ USE lbclnk
USE timing
IMPLICIT NONE
@@ -35,6 +38,7 @@ MODULE isfstp
PUBLIC isf_stp, isf_init, isf_nam ! routine called in sbcmod and divhor
!! * Substitutions
+# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -60,42 +64,49 @@ CONTAINS
INTEGER, INTENT(in) :: kt ! ocean time step
INTEGER, INTENT(in) :: Kmm ! ocean time level index
!
- INTEGER :: jk ! loop index
-#if defined key_qco
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! 3D workspace
-#endif
+ INTEGER :: ji, jj, jk, ikt ! loop index
+ REAL(wp), DIMENSION(A2D(0)) :: zhtmp ! temporary array for thickness
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ze3t ! 3D workspace for key_qco
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('isf')
!
+ ! temporary arrays for key_qco
+ DO_2D( 0 ,0, 0, 0 )
+ zhtmp(ji,jj) = ht(ji,jj)
+ DO jk = 1, jpk
+ ze3t(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
+ ENDDO
+ END_2D
+ !
!=======================================================================
! 1.: compute melt and associated heat fluxes in the ice shelf cavities
!=======================================================================
!
IF ( ln_isfcav_mlt ) THEN
!
+ ! --- before time step --- !
#if ! defined key_RK3
- ! MLF : need risf_cav_tsc_b update
- ! 1.1: before time step
- IF ( kt /= nit000 ) THEN
- risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
- fwfisf_cav_b(:,:) = fwfisf_cav(:,:)
+ IF ( kt /= nit000 ) THEN ! MLF : need risf_cav_tsc_b update
+ DO_2D( 0, 0, 0, 0 )
+ risf_cav_tsc_b(ji,jj,:) = risf_cav_tsc(ji,jj,:)
+ fwfisf_cav_b (ji,jj) = fwfisf_cav (ji,jj)
+ END_2D
END IF
#endif
!
- ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
- rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
-#if defined key_qco
- DO jk = 1, jpk
- ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
- END DO
- CALL isf_tbl_lvl( ht(:,:), ze3t , misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
-#else
- CALL isf_tbl_lvl( ht(:,:), e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
-#endif
+ ! --- deepest level (misfkb), thickness (rhisf) & fraction of deepest cell affected by tbl (rfrac) --- !
+ DO_2D( 0 ,0, 0, 0 )
+ ! limit the tbl to water depth and to the top level thickness
+ ikt = misfkt_cav(ji,jj) ! tbl top indices
+ rhisf_tbl_cav(ji,jj) = MAX( MIN( rn_htbl * mskisf_cav(ji,jj), zhtmp(ji,jj) ), ze3t(ji,jj,ikt) )
+ END_2D
+
+ CALL isf_tbl_lvl( ze3t, misfkt_cav, rhisf_tbl_cav, & ! <<== in
+ & misfkb_cav, rfrac_tbl_cav ) ! ==>> out
!
- ! 1.3: compute ice shelf melt
- CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav )
+ ! --- ice shelf melt (fwfisf) and temperature trend (risf) --- !
+ CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav(A2D(0)) ) ! <<==>> inout
!
END IF
!
@@ -105,37 +116,59 @@ CONTAINS
!
IF ( ln_isfpar_mlt ) THEN
!
+ ! --- before time step --- !
#if ! defined key_RK3
- ! MLF : need risf_par_tsc_b update
- ! 2.1: before time step
- IF ( kt /= nit000 ) THEN
- risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
- fwfisf_par_b (:,:) = fwfisf_par (:,:)
+ IF ( kt /= nit000 ) THEN ! MLF : need risf_par_tsc_b update
+ DO_2D( 0, 0, 0, 0 )
+ risf_par_tsc_b(ji,jj,:) = risf_par_tsc(ji,jj,:)
+ fwfisf_par_b (ji,jj) = fwfisf_par (ji,jj)
+ END_2D
END IF
#endif
!
- ! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+ ! --- deepest level (misfkb), thickness (rhisf) & fraction of deepest cell affected by tbl (rfrac) --- !
! by simplicity, we assume the top level where param applied do not change with time (done in init part)
- rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
-#if defined key_qco
- DO jk = 1, jpk
- ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
- END DO
- CALL isf_tbl_lvl( ht(:,:), ze3t , misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
-#else
- CALL isf_tbl_lvl( ht(:,:), e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
-#endif
+ ! limit the tbl to water depth and to the top level thickness
+ DO_2D( 0 ,0, 0, 0 )
+ ikt = misfkt_cav(ji,jj) ! tbl top indices
+ rhisf_tbl_par(ji,jj) = MAX( MIN( rhisf0_tbl_par(ji,jj), zhtmp(ji,jj) ), ze3t(ji,jj,ikt) )
+ END_2D
+
+ CALL isf_tbl_lvl( ze3t, misfkt_par, rhisf_tbl_par, & ! <<== in
+ & misfkb_par, rfrac_tbl_par ) ! ==>> out
!
- ! 2.3: compute ice shelf melt
- CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par )
+ ! --- ice shelf melt (fwfisf) and temperature trend (risf) --- !
+ CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par(A2D(0)) ) ! <<==>> inout
!
END IF
!
- !==================================================================================
- ! 3.: output specific restart variable in case of coupling with an ice sheet model
- !==================================================================================
!
- IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt, Kmm)
+ !clem: these lbc are needed since we calculate everything in the interior now
+ IF( ln_isfcpl ) THEN
+ CALL lbc_lnk( 'isf_stp', fwfisf_par , 'T', 1.0_wp, fwfisf_cav , 'T', 1.0_wp, &
+#if ! defined key_RK3
+ & fwfisf_par_b, 'T', 1.0_wp, fwfisf_cav_b, 'T', 1.0_wp, &
+#endif
+ & risfcpl_ssh, 'T', 1.0_wp, risfcpl_cons_ssh, 'T', 1.0_wp ) ! needed in dynspg_ts, stp2d
+ CALL lbc_lnk( 'isf_stp', risfcpl_vol, 'T', 1.0_wp ) ! needed in dynspg_ts, stp2d, sshwzv, dynatf
+ ELSE
+ CALL lbc_lnk( 'isf_stp', fwfisf_par , 'T', 1.0_wp, fwfisf_cav , 'T', 1.0_wp, &
+#if ! defined key_RK3
+ & fwfisf_par_b, 'T', 1.0_wp, fwfisf_cav_b, 'T', 1.0_wp &
+ & )
+#endif
+ ENDIF
+ !
+ !==================
+ ! 3.: write restart
+ !==================
+#if ! defined key_RK3
+ ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+ IF( ln_isfcav_mlt .AND. lrst_oce ) CALL isfrst_write( kt, 'cav', risf_cav_tsc , fwfisf_cav )
+ ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+ IF( ln_isfpar_mlt .AND. lrst_oce ) CALL isfrst_write( kt, 'par', risf_par_tsc , fwfisf_par )
+#endif
+ IF( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write( kt, Kmm )
!
IF( ln_timing ) CALL timing_stop('isf')
!
@@ -164,10 +197,19 @@ CONTAINS
CALL isf_nam() ! Read namelist
!
CALL isf_alloc() ! Allocate public array
+ !
+ ! initalisation of fwf and tsc array to 0
+ risfload (:,:) = 0._wp
+ fwfisf_oasis(:,:) = 0._wp ; fwfisf_par (:,:) = 0._wp ; fwfisf_cav(:,:) = 0._wp
+ risf_cav_tsc(:,:,:) = 0._wp ; risf_par_tsc(:,:,:) = 0._wp
+#if ! defined key_RK3
+ fwfisf_par_b (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp
+ risf_cav_tsc_b(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp
+#endif
!
CALL isf_ctl() ! check option compatibility
!
- IF( ln_isfcav ) CALL isf_load( Kmm, risfload ) ! compute ice shelf load
+ IF( ln_isfcav ) CALL isf_load( Kmm, risfload ) ! compute ice shelf load
!
! terminate routine now if no ice shelf melt formulation specify
IF( ln_isf ) THEN
@@ -200,11 +242,9 @@ CONTAINS
WRITE(numout,*)
!
IF ( ln_isf ) THEN
-#if key_qco
-# if ! defined key_isf
+#if defined key_qco && ! defined key_isf
CALL ctl_stop( 'STOP', 'isf_ctl: ice shelf requires both ln_isf=T AND key_isf activated' )
-# endif
-#endif
+#endif
WRITE(numout,*) ' Add debug print in isf module ln_isfdebug = ', ln_isfdebug
WRITE(numout,*)
WRITE(numout,*) ' melt inside the cavity ln_isfcav_mlt = ', ln_isfcav_mlt
diff --git a/src/OCE/ISF/isftbl.F90 b/src/OCE/ISF/isftbl.F90
index c6b3adc3b..8820b171e 100644
--- a/src/OCE/ISF/isftbl.F90
+++ b/src/OCE/ISF/isftbl.F90
@@ -15,20 +15,20 @@ MODULE isftbl
USE isf_oce ! ice shelf variables
+ USE par_oce ! ocean space and time domain
USE dom_oce ! vertical scale factor and depth
IMPLICIT NONE
PRIVATE
- PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
+ PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
-
CONTAINS
- SUBROUTINE isf_tbl( Kmm, pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
+ SUBROUTINE isf_tbl( Kmm, pvarin, cd_ptin, ktop, phtbl, pvarout, kbot, pfrac )
!!--------------------------------------------------------------------
!! *** SUBROUTINE isf_tbl ***
!!
@@ -39,77 +39,69 @@ CONTAINS
!! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
!! https://doi.org/10.1029/2007JC004368 , 2008
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pvarout ! 2d average of pvarin
- !!-------------------------- IN -------------------------------------
- INTEGER , INTENT(in ) :: Kmm ! ocean time level index
- CHARACTER(len=1) , INTENT(in ) :: cd_ptin ! point of variable in/out
- REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pvarin ! 3d variable to average over the tbl
- INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl ! tbl thickness
- !!-------------------------- IN OPTIONAL -----------------------------
- INTEGER, DIMENSION(jpi,jpj), OPTIONAL, INTENT(in ) :: kbot ! bottom level
- REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in ) :: pfrac ! fraction of bottom cell affected by tbl
!!--------------------------------------------------------------------
- INTEGER :: ji, jj ! loop index
- INTEGER , DIMENSION(jpi,jpj) :: ikbot ! bottom level of the tbl
- REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
- REAL(wp), DIMENSION(jpi,jpj) :: zhtbl ! thickness of the tbl
- REAL(wp), DIMENSION(jpi,jpj) :: zfrac ! thickness of the tbl
- INTEGER :: jk ! loop index
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t,ze3u,ze3v ! e3
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ REAL(wp), DIMENSION(A2D(0),jpk) , INTENT(in ) :: pvarin ! 3d variable to average over the tbl
+ CHARACTER(len=1) , INTENT(in ) :: cd_ptin ! point of variable in/out
+ INTEGER, DIMENSION(A2D(0)) , INTENT(in ) :: ktop ! top level
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: phtbl ! tbl thickness
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(out) :: pvarout ! 2d average of pvarin
+ INTEGER, DIMENSION(A2D(0)), OPTIONAL, INTENT(in ) :: kbot ! bottom level
+ REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in ) :: pfrac ! fraction of bottom cell affected by tbl
+ !!--------------------------------------------------------------------
+ INTEGER :: ji, jj, jk, ikt ! loop index
+ REAL(wp), DIMENSION(A2D(0)) :: zhtbl ! temporary array for thickness
+ INTEGER , DIMENSION(A2D(0)) :: ikbot ! bottom level of the tbl
+ REAL(wp), DIMENSION(A2D(0)) :: zfrac ! thickness of the tbl
+ REAL(wp), DIMENSION(A2D(0),jpk) :: ze3 ! e3
!!--------------------------------------------------------------------
!
SELECT CASE ( cd_ptin )
CASE ( 'U' )
!
- ! copy phtbl (phtbl is INTENT in as we don't want to change it)
- zhtbl = phtbl
+ DO_3D( 0 ,0, 0, 0, 1, jpk )
+ ze3(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
+ END_3D
!
- DO jk = 1, jpk
- ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
- END DO
! compute tbl lvl and thickness
- CALL isf_tbl_lvl( hu(:,:,Kmm), ze3u, ktop, ikbot, zhtbl, zfrac )
+ DO_2D( 0 ,0, 0, 0 )
+ ikt = ktop(ji,jj) ! tbl top indices
+ zhtbl(ji,jj) = MAX( MIN( phtbl(ji,jj), hu(ji,jj,Kmm) ), ze3(ji,jj,ikt) )
+ END_2D
+ CALL isf_tbl_lvl( ze3, ktop, zhtbl, & ! <<== in
+ & ikbot, zfrac ) ! ==>> out
!
! compute tbl property at U point
- CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, ze3u, pvarin, zvarout )
+ CALL isf_tbl_avg( miku(A2D(0)), ikbot, zhtbl, zfrac, ze3, pvarin, & ! <<== in
+ & pvarout ) ! ==>> out
!
- ! compute tbl property at T point
- pvarout(1,:) = 0._wp
- DO_2D( nn_hls-1, nn_hls, nn_hls, nn_hls )
- pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
- END_2D
- ! lbclnk not needed as a final communication is done after the computation of fwf
- !
CASE ( 'V' )
!
- ! copy phtbl (phtbl is INTENT in as we don't want to change it)
- zhtbl = phtbl
+ DO_3D( 0 ,0, 0, 0, 1, jpk )
+ ze3(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
+ END_3D
!
- DO jk = 1, jpk
- ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
- END DO
! compute tbl lvl and thickness
- CALL isf_tbl_lvl( hv(:,:,Kmm), ze3v, ktop, ikbot, zhtbl, zfrac )
+ DO_2D( 0 ,0, 0, 0 )
+ ikt = ktop(ji,jj) ! tbl top indices
+ zhtbl(ji,jj) = MAX( MIN( phtbl(ji,jj), hv(ji,jj,Kmm) ), ze3(ji,jj,ikt) )
+ END_2D
+ CALL isf_tbl_lvl( ze3, ktop, zhtbl, & ! <<== in
+ & ikbot, zfrac ) ! ==>> out
!
! compute tbl property at V point
- CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, ze3v, pvarin, zvarout )
- !
- ! pvarout is an averaging of wet point
- pvarout(:,1) = 0._wp
- DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls )
- pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
- END_2D
- ! lbclnk not needed as a final communication is done after the computation of fwf
+ CALL isf_tbl_avg( mikv(A2D(0)), ikbot, zhtbl, zfrac, ze3, pvarin, & ! <<== in
+ & pvarout ) ! ==>> out
!
CASE ( 'T' )
+ !
+ DO_3D( 0 ,0, 0, 0, 1, jpk )
+ ze3(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
+ END_3D
!
! compute tbl property at T point
- DO jk = 1, jpk
- ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
- END DO
- CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, ze3t, pvarin, pvarout )
+ CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, ze3, pvarin, & ! <<== in
+ & pvarout ) ! ==>> out
!
END SELECT
!
@@ -124,20 +116,19 @@ CONTAINS
!! ** Method : Depth average is made between the top level ktop and the bottom level kbot
!! over a thickness phtbl. The bottom level is partially counted (pfrac).
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pvarout ! tbl property averaged over phtbl between level ktop and kbot
- !!-------------------------- IN -------------------------------------
- INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop, kbot ! top and bottom level of the top boundary layer
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pvarin ! tbl property to average between ktop, kbot over phtbl
!!--------------------------------------------------------------------
- INTEGER :: ji,jj,jk ! loop indices
- INTEGER :: ikt, ikb ! top and bottom levels
+ INTEGER, DIMENSION(A2D(0)) , INTENT(in ) :: ktop, kbot ! top and bottom level of the top boundary layer
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
+ REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in ) :: pe3 ! vertical scale factor
+ REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in ) :: pvarin ! tbl property to average between ktop, kbot over phtbl
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(out) :: pvarout ! tbl property averaged over phtbl between level ktop and kbot
+ !!--------------------------------------------------------------------
+ INTEGER :: ji, jj, jk ! loop indices
+ INTEGER :: ikt, ikb ! top and bottom levels
!!--------------------------------------------------------------------
!
! compute tbl top.bottom level and thickness
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
!
! tbl top/bottom indices initialisation
ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
@@ -152,7 +143,7 @@ CONTAINS
END SUBROUTINE isf_tbl_avg
- SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
+ SUBROUTINE isf_tbl_lvl( pe3, ktop, phtbl, kbot, pfrac )
!!--------------------------------------------------------------------
!! *** ROUTINE isf_tbl_lvl ***
!!
@@ -160,96 +151,48 @@ CONTAINS
!! - thickness of the top boundary layer
!! - fraction of the bottom level affected by the tbl
!!
- !!-------------------------- OUT --------------------------------------
- INTEGER, DIMENSION(jpi,jpj) , INTENT( out) :: kbot ! bottom level of the top boundary layer
- REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pfrac ! fraction of bottom level in the tbl
- !!-------------------------- IN --------------------------------------
- INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level of the top boundary layer
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phw ! water column thickness
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
- !!-------------------------- INOUT ------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: phtbl ! top boundary layer thickness
+ !!--------------------------------------------------------------------
+ REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in ) :: pe3 ! vertical scale factor
+ INTEGER, DIMENSION(A2D(0)) , INTENT(in ) :: ktop ! top level of the top boundary layer
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: phtbl ! top boundary layer thickness
+ INTEGER, DIMENSION(A2D(0)) , INTENT(out) :: kbot ! bottom level of the top boundary layer
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(out) :: pfrac ! fraction of bottom level in the tbl
!!---------------------------------------------------------------------
- INTEGER :: ji,jj,jk
- INTEGER :: ikt, ikb
+ INTEGER :: ji, jj, jk
+ INTEGER :: ikt, ikb
!!---------------------------------------------------------------------
!
- ! get htbl
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- !
- ! tbl top/bottom indices initialisation
- ikt = ktop(ji,jj)
- !
- ! limit the tbl to water thickness.
- phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) )
+ DO_2D( 0, 0, 0, 0 )
!
- ! thickness of boundary layer must be at least the top level thickness
- phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) )
+ ikt = ktop(ji,jj) ! tbl top indices initialisation
!
- END_2D
- !
- ! get ktbl
- CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
- !
- ! get pfrac
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ ! --- get kbot --- !
+ ! ! determine the deepest level influenced by the boundary layer
+ ikb = ikt+1
+ DO WHILE( SUM( pe3(ji,jj,ikt:ikb-1) ) < phtbl(ji,jj ) ) ; ikb = ikb + 1 ; END DO
+ kbot(ji,jj) = ikb - 1
!
- ! tbl top/bottom indices initialisation
- ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
+ ikb = kbot(ji,jj) ! tbl bottom indices initialisation
!
- ! proportion of the bottom cell included in ice shelf boundary layer
+ ! --- get pfrac --- !
+ ! ! proportion of the bottom cell included in ice shelf boundary layer
pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb)
!
END_2D
!
END SUBROUTINE isf_tbl_lvl
!
- SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
- !!--------------------------------------------------------------------
- !! *** ROUTINE isf_tbl_bot ***
- !!
- !! ** Purpose : compute bottom level of the isf top boundary layer
- !!
- !!-------------------------- OUT -------------------------------------
- INTEGER, DIMENSION(jpi,jpj) , INTENT( out) :: kbot ! bottom level of the top boundary layer
- !!-------------------------- IN -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl ! top boundary layer thickness
- INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level of the top boundary layer
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
- !!--------------------------------------------------------------------
- INTEGER :: ji, jj
- INTEGER :: ikt, ikb
- !!--------------------------------------------------------------------
- !
- ! phtbl need to be bounded by water column thickness before
- ! test: if htbl = water column thickness, should return mbathy
- ! test: if htbl = 0 should return ktop (phtbl cap to pe3t(ji,jj,1))
- !
- ! get ktbl
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- !
- ! determine the deepest level influenced by the boundary layer
- ikt = ktop(ji,jj)
- ikb = ikt
- DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ; ikb = ikb + 1 ; END DO
- kbot(ji,jj) = ikb - 1
- !
- END_2D
- !
- END SUBROUTINE isf_tbl_kbot
- !
- SUBROUTINE isf_tbl_ktop(pdep, ktop)
+ SUBROUTINE isf_tbl_ktop( pdep, ktop )
!!--------------------------------------------------------------------
!! *** ROUTINE isf_tbl_top ***
!!
!! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
!!
- !!-------------------------- OUT -------------------------------------
- INTEGER, DIMENSION(jpi,jpj), INTENT( out) :: ktop ! top level affected by the ice shelf parametrisation
- !!-------------------------- IN -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdep ! top depth of the parametrisation influence
!!--------------------------------------------------------------------
- INTEGER :: ji,jj
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pdep ! top depth of the parametrisation influence
+ INTEGER, DIMENSION(A2D(0)), INTENT( out) :: ktop ! top level affected by the ice shelf parametrisation
+ !!--------------------------------------------------------------------
+ INTEGER :: ji, jj
INTEGER :: ikt
!!--------------------------------------------------------------------
!
@@ -260,7 +203,7 @@ CONTAINS
! test: this routine run on isfdraft should return mikt
! test: this routine run with pdep = 0 should return 1
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
! comput ktop
ikt = 2
DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ; ikt = ikt + 1 ; END DO
diff --git a/src/OCE/ISF/isfutils.F90 b/src/OCE/ISF/isfutils.F90
index 63c489c59..af103b46b 100644
--- a/src/OCE/ISF/isfutils.F90
+++ b/src/OCE/ISF/isfutils.F90
@@ -13,9 +13,9 @@ MODULE isfutils
USE iom , ONLY: iom_open, iom_get, iom_close, jpdom_global ! read input file
USE lib_fortran , ONLY: glob_sum, glob_min, glob_max ! compute global value
- USE par_oce , ONLY: jpi,jpj,jpk, jpnij, Nis0, Nie0, Njs0, Nje0 ! domain size
+ USE par_oce ! domain size
USE dom_oce , ONLY: narea ! local domain
- USE in_out_manager, ONLY: i8, wp, lwp, numout ! miscelenious
+ USE in_out_manager ! miscelenious
USE lib_mpp
IMPLICIT NONE
@@ -28,6 +28,8 @@ MODULE isfutils
PUBLIC read_2dcstdta, debug
+ !! * Substitutions
+# include "do_loop_substitute.h90"
CONTAINS
SUBROUTINE read_2dcstdta(cdfile, cdvar, pvar)
@@ -36,17 +38,16 @@ CONTAINS
!!
!! ** Purpose : read input file
!!
- !!-------------------------- OUT -------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pvar ! output variable
- !!-------------------------- IN -------------------------------------
- CHARACTER(len=*) , INTENT(in ) :: cdfile ! input file name
- CHARACTER(len=*) , INTENT(in ) :: cdvar ! variable name
+ !!--------------------------------------------------------------------
+ CHARACTER(len=*) , INTENT(in ) :: cdfile ! input file name
+ CHARACTER(len=*) , INTENT(in ) :: cdvar ! variable name
+ REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pvar ! output variable
!!--------------------------------------------------------------------
INTEGER :: inum
!!--------------------------------------------------------------------
CALL iom_open( TRIM(cdfile), inum )
- CALL iom_get( inum, jpdom_global, TRIM(cdvar), pvar)
+ CALL iom_get( inum, jpdom_global, TRIM(cdvar), pvar )
CALL iom_close(inum)
END SUBROUTINE read_2dcstdta
@@ -57,16 +58,16 @@ CONTAINS
!!
!! ** Purpose : add debug print for 2d variables
!!
- !!-------------------------- IN -------------------------------------
- CHARACTER(LEN=*) , INTENT(in ) :: cdtxt
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvar
!!--------------------------------------------------------------------
- REAL(wp) :: zmin, zmax, zsum
- INTEGER(i8) :: imodd, ip
- INTEGER :: imods
- INTEGER :: isums, idums
- INTEGER :: ji,jj,jk
- INTEGER, DIMENSION(jpnij) :: itmps
+ CHARACTER(LEN=*) , INTENT(in) :: cdtxt
+ REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pvar
+ !!--------------------------------------------------------------------
+ REAL(wp) :: zmin, zmax, zsum
+ INTEGER(i8) :: imodd, ip
+ INTEGER :: imods
+ INTEGER :: isums, idums
+ INTEGER :: ji, jj, jk
+ INTEGER, DIMENSION(jpnij) :: itmps
!!--------------------------------------------------------------------
!
! global min/max/sum to check data range and NaN
@@ -110,16 +111,16 @@ CONTAINS
!!
!! ** Purpose : add debug print for 3d variables
!!
- !!-------------------------- IN -------------------------------------
- CHARACTER(LEN=*) , INTENT(in ) :: cdtxt
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pvar
!!--------------------------------------------------------------------
- REAL(wp) :: zmin, zmax, zsum
- INTEGER(i8) :: imodd, ip
- INTEGER :: imods
- INTEGER :: isums, idums
- INTEGER :: ji,jj,jk
- INTEGER, DIMENSION(jpnij) :: itmps
+ CHARACTER(LEN=*) , INTENT(in) :: cdtxt
+ REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in) :: pvar
+ !!--------------------------------------------------------------------
+ REAL(wp) :: zmin, zmax, zsum
+ INTEGER(i8) :: imodd, ip
+ INTEGER :: imods
+ INTEGER :: isums, idums
+ INTEGER :: ji, jj, jk
+ INTEGER, DIMENSION(jpnij) :: itmps
!!--------------------------------------------------------------------
!
! global min/max/sum to check data range and NaN
diff --git a/src/OCE/SBC/sbc_oce.F90 b/src/OCE/SBC/sbc_oce.F90
index 56117e4fd..5df541eb5 100644
--- a/src/OCE/SBC/sbc_oce.F90
+++ b/src/OCE/SBC/sbc_oce.F90
@@ -188,8 +188,7 @@ CONTAINS
ALLOCATE( emp(jpi,jpj) , emp_b(jpi,jpj) , &
& STAT=ierr(2) )
!
- ALLOCATE( rnf (jpi,jpj) , rnf_b (jpi,jpj) , &
- & fwficb (jpi,jpj) , STAT=ierr(3) )
+ ALLOCATE( rnf(jpi,jpj), rnf_b(jpi,jpj), STAT=ierr(3) )
!
ALLOCATE( fr_i(jpi,jpj) , &
& ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , &
@@ -202,7 +201,7 @@ CONTAINS
& qns_tot(A2D(0)) , qsr_tot(A2D(0)) , qsr_hc(A2D(0),jpk) , qsr_hc_b(A2D(0),jpk) , STAT=ierr(5) )
!
ALLOCATE( sbc_tsc(A2D(0),jpts) , sbc_tsc_b(A2D(0),jpts) , &
- & sfx (A2D(0)) , sfx_b(A2D(0)) , emp_tot(A2D(0)), fmmflx(A2D(0)) ,&
+ & sfx (A2D(0)) , sfx_b(A2D(0)) , emp_tot(A2D(0)), fmmflx(A2D(0)), fwficb(A2D(0)), &
& wndm(A2D(0)) , taum (A2D(0)) , STAT=ierr(6) )
!
ALLOCATE( tprecip(A2D(0)) , sprecip(A2D(0)) , &
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index f681b98be..8b8c72576 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -1434,13 +1434,13 @@ CONTAINS
IF( srcv(jpr_cal)%laction ) zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
IF( srcv(jpr_icb)%laction ) THEN
- fwficb(A2D(0)) = frcv(jpr_icb)%z3(:,:,1)
- rnf (A2D(0)) = rnf(A2D(0)) + fwficb(A2D(0)) ! iceberg added to runfofs
+ fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+ rnf(A2D(0)) = rnf(A2D(0)) + fwficb(:,:) ! iceberg added to runfofs
ENDIF
!
! ice shelf fwf
IF( srcv(jpr_isf)%laction ) THEN
- fwfisf_oasis(A2D(0)) = frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf to the ocean ( > 0 = melting )
+ fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf to the ocean ( > 0 = melting )
END IF
IF( ln_mixcpl ) THEN ; emp(A2D(0)) = emp(A2D(0)) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
@@ -1748,11 +1748,11 @@ CONTAINS
zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
ENDIF
IF( srcv(jpr_icb)%laction ) THEN ! iceberg added to runoffs
- fwficb(A2D(0)) = frcv(jpr_icb)%z3(:,:,1)
- rnf (A2D(0)) = rnf(A2D(0)) + fwficb(A2D(0))
+ fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+ rnf(A2D(0)) = rnf(A2D(0)) + fwficb(:,:)
ENDIF
IF( srcv(jpr_isf)%laction ) THEN ! iceshelf (fwfisf > 0 mean melting)
- fwfisf_oasis(A2D(0)) = frcv(jpr_isf)%z3(:,:,1)
+ fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)
ENDIF
IF( ln_mixcpl ) THEN
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 82d0dd6e9..4561f6a45 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -51,7 +51,6 @@ MODULE sbcmod
USE sbcfwb ! surface boundary condition: freshwater budget
USE icbstp ! Icebergs
USE icb_oce , ONLY : ln_passive_mode ! iceberg interaction mode
- USE isf_oce , ONLY : ln_isf, l_isfoasis, fwfisf_oasis
USE traqsr ! active tracers: light penetration
USE sbcwave ! Wave module
USE bdy_oce , ONLY: ln_bdy
@@ -466,8 +465,7 @@ CONTAINS
IF( ll_sas .OR. ll_opa ) CALL lbc_lnk( 'sbcmod', sst_m, 'T', 1.0_wp, sss_m, 'T', 1.0_wp, ssh_m, 'T', 1.0_wp, &
& frq_m, 'T', 1.0_wp, e3t_m, 'T', 1.0_wp, fr_i , 'T', 1.0_wp )
!clem : these calls are needed for sbccpl => it needs an IF statement but it's complicated
- IF( ln_isf .AND. l_isfoasis ) CALL lbc_lnk( 'sbcmod', fwfisf_oasis, 'T', 1.0_wp )
- IF( ln_rnf .AND. l_rnfcpl ) CALL lbc_lnk( 'sbcmod', rnf, 'T', 1.0_wp, fwficb , 'T', 1.0_wp )
+ IF( ln_rnf .AND. l_rnfcpl ) CALL lbc_lnk( 'sbcmod', rnf, 'T', 1.0_wp )
!
IF( ln_icebergs ) THEN ! save pure stresses (with no ice-ocean stress) for use by icebergs
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
@@ -542,7 +540,7 @@ CONTAINS
! (in addition to snwice_mass)
IF( ln_rnf ) THEN
CALL lbc_lnk( 'sbcmod', utau, 'T', -1.0_wp, vtau , 'T', -1.0_wp, emp, 'T', 1.0_wp, &
- & rnf , 'T', 1.0_wp, fwficb, 'T', 1.0_wp ) ! fwficb is used on the halos in pisces (only)
+ & rnf , 'T', 1.0_wp )
ELSE
CALL lbc_lnk( 'sbcmod', utau, 'T', -1.0_wp, vtau , 'T', -1.0_wp, emp, 'T', 1.0_wp )
ENDIF
diff --git a/src/OCE/SBC/sbcrnf.F90 b/src/OCE/SBC/sbcrnf.F90
index 367e3c7cf..2562d0948 100644
--- a/src/OCE/SBC/sbcrnf.F90
+++ b/src/OCE/SBC/sbcrnf.F90
@@ -83,8 +83,8 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** ROUTINE sbc_rnf_alloc ***
!!----------------------------------------------------------------------
- ALLOCATE( rnfmsk(jpi,jpj) , rnfmsk_z(jpk) , &
- & h_rnf (jpi,jpj) , nk_rnf (jpi,jpj) , &
+ ALLOCATE( rnfmsk(jpi,jpj) , rnfmsk_z(jpk) , & ! needed over the whole domain by muscl (traadv_muscl)
+ & h_rnf (A2D(0)) , nk_rnf (A2D(0)) , &
& rnf_tsc_b(A2D(0),jpts) , rnf_tsc (A2D(0),jpts) , STAT=sbc_rnf_alloc )
!
CALL mpp_sum ( 'sbcrnf', sbc_rnf_alloc )
@@ -127,13 +127,13 @@ CONTAINS
IF( .NOT. l_rnfcpl ) THEN
rnf(A2D(0)) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * smask0(:,:) ! updated runoff value at time step kt
IF( ln_rnf_icb ) THEN
- fwficb(A2D(0)) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * smask0(:,:) ! updated runoff value at time step kt
- rnf(A2D(0)) = rnf(A2D(0)) + fwficb(A2D(0))
- qns(:,:) = qns(:,:) - fwficb(A2D(0)) * rLfus
+ fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * smask0(:,:) ! updated runoff value at time step kt
+ rnf(A2D(0)) = rnf(A2D(0)) + fwficb(:,:)
+ qns(:,:) = qns(:,:) - fwficb(:,:) * rLfus
!!qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus
!!qns_oce(:,:) = qns_oce(:,:) - fwficb(:,:) * rLfus
- CALL iom_put( 'iceberg_cea' , fwficb(A2D(0)) ) ! output iceberg flux
- CALL iom_put( 'hflx_icb_cea' , -fwficb(A2D(0)) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
+ CALL iom_put( 'iceberg_cea' , fwficb(:,:) ) ! output iceberg flux
+ CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
ENDIF
ENDIF
!
@@ -209,7 +209,7 @@ CONTAINS
!
IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN !== runoff distributed over several levels ==!
IF( ln_linssh ) THEN !* constant volume case : just apply the runoff input flow
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+ DO_2D_OVR( 0, 0, 0, 0 )
DO jk = 1, nk_rnf(ji,jj)
#if defined key_RK3
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
@@ -219,14 +219,10 @@ CONTAINS
END DO
END_2D
ELSE !* variable volume case
- DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) ! update the depth over which runoffs are distributed
+ DO_2D_OVR( 0, 0, 0, 0 ) ! update the depth over which runoffs are distributed
h_rnf(ji,jj) = 0._wp
- DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
- h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! to the bottom of the relevant grid box
- END DO
- END_2D
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! apply the runoff input flow
DO jk = 1, nk_rnf(ji,jj)
+ h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! recalculates h_rnf to be the depth in metres to the bottom of the relevant grid box
#if defined key_RK3
phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
#else
@@ -236,10 +232,8 @@ CONTAINS
END_2D
ENDIF
ELSE !== runoff put only at the surface ==!
- DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D_OVR( 0, 0, 0, 0 )
h_rnf (ji,jj) = e3t(ji,jj,1,Kmm) ! update h_rnf to be depth of top box
- END_2D
- DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
#if defined key_RK3
phdivn(ji,jj,1) = phdivn(ji,jj,1) - rnf(ji,jj) * r1_rho0 / e3t(ji,jj,1,Kmm) ! RK3: rnf forcing at n+1/2
#else
@@ -268,7 +262,7 @@ CONTAINS
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: nbrec ! temporary integer
REAL(wp) :: zacoef
- REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl
+ REAL(wp), DIMENSION(A2D(0),2) :: zrnfcl
!!
NAMELIST/namsbc_rnf/ cn_dir , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, &
& sn_rnf, sn_cnf , sn_i_rnf, sn_s_rnf , sn_t_rnf , sn_dep_rnf, &
@@ -337,7 +331,7 @@ CONTAINS
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(A2D(0),1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
fwficb(ji,jj) = 0._wp
END_2D
ENDIF
@@ -379,7 +373,7 @@ CONTAINS
CALL iom_get ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf, kfill = jpfillcopy ) ! read the river mouth. no 0 on halos!
CALL iom_close( inum ) ! close file
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the number of level over which river runoffs are applied
+ DO_2D( 0, 0, 0, 0 ) ! set the number of level over which river runoffs are applied
IF( h_rnf(ji,jj) > 0._wp ) THEN
jk = 2
DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
@@ -393,7 +387,7 @@ CONTAINS
ENDIF
END_2D
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth
+ DO_2D( 0, 0, 0, 0 ) ! set the associated depth
h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj)
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
@@ -423,14 +417,14 @@ CONTAINS
ELSEWHERE ; h_rnf(:,:) = 1._wp
ENDWHERE
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! take in account min depth of ocean rn_hmin
+ DO_2D( 0, 0, 0, 0 ) ! take in account min depth of ocean rn_hmin
IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
jk = mbkt(ji,jj)
h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) )
ENDIF
END_2D
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! number of levels on which runoffs are distributed
+ DO_2D( 0, 0, 0, 0 ) ! number of levels on which runoffs are distributed
IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
jk = 2
DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
@@ -441,7 +435,7 @@ CONTAINS
ENDIF
END_2D
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth
+ DO_2D( 0, 0, 0, 0 ) ! set the associated depth
h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj)
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
@@ -455,7 +449,7 @@ CONTAINS
CALL iom_close ( inum )
ENDIF
ELSE ! runoffs applied at the surface
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
nk_rnf(ji,jj) = 1
h_rnf (ji,jj) = e3t(ji,jj,1,Kmm)
END_2D
@@ -543,9 +537,9 @@ CONTAINS
ENDIF
!
! horizontal mask (read in NetCDF file)
- CALL iom_open ( cl_rnfile, inum ) ! open file
- CALL iom_get ( inum, jpdom_global, sn_cnf%clvar, rnfmsk ) ! read the river mouth array
- CALL iom_close( inum ) ! close file
+ CALL iom_open ( cl_rnfile, inum ) ! open file
+ CALL iom_get ( inum, jpdom_global, sn_cnf%clvar, rnfmsk(A2D(0)) ) ! read the river mouth array
+ CALL iom_close( inum ) ! close file
!
IF( l_clo_rnf ) CALL clo_rnf( rnfmsk ) ! closed sea inflow set as river mouth
!
@@ -556,6 +550,8 @@ CONTAINS
rnfmsk_z(4) = 0.25 ! **********
rnfmsk_z(5) = 0.125
!
+ CALL lbc_lnk( 'rnf_mouth', rnfmsk, 'T', 1._wp ) ! needed by muscl scheme (traadv_muscl)
+ !
END SUBROUTINE rnf_mouth
!!======================================================================
diff --git a/src/OCE/TRA/eosbn2.F90 b/src/OCE/TRA/eosbn2.F90
index 5a04a3ca0..e675be909 100644
--- a/src/OCE/TRA/eosbn2.F90
+++ b/src/OCE/TRA/eosbn2.F90
@@ -61,7 +61,7 @@ MODULE eosbn2
END INTERFACE
!
INTERFACE eos_fzp
- MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
+ MODULE PROCEDURE eos_fzp_3d, eos_fzp_2d, eos_fzp_0d
END INTERFACE
!
PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules
@@ -223,12 +223,20 @@ CONTAINS
!! TEOS-10 Manual, 2010
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in ) :: pts ! T-S
- INTEGER , INTENT(in ) :: Knn ! time-level
- REAL(wp), DIMENSION(:,:,: ), INTENT( out) :: prd ! in situ density
+ INTEGER , INTENT(in ) :: Knn ! time-level
+ REAL(wp), DIMENSION(:,:,: ) , INTENT( out) :: prd ! in situ density
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-insitu')
@@ -237,7 +245,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Knn) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem,Knn) * r1_T0 ! temperature
@@ -273,7 +281,7 @@ CONTAINS
!
CASE( np_seos ) !== simplified EOS ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - rn_T0
zs = pts (ji,jj,jk,jp_sal,Knn) - rn_S0
zh = gdept(ji,jj,jk,Knn)
@@ -340,14 +348,26 @@ CONTAINS
!! TEOS-10 Manual, 2010
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrd, ktdep
- REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
- ! ! 2 : salinity [psu]
- REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
- REAL(wp), DIMENSION(A2D_T(ktdep),JPK ), INTENT(in ) :: pdep ! depth [m]
+!!$ REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+!!$ ! ! 2 : salinity [psu]
+!!$ REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
+!!$ REAL(wp), DIMENSION(A2D_T(ktdep),JPK ), INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m]
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-insitu')
@@ -356,7 +376,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
!
zh = pdep(ji,jj,jk) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
@@ -392,7 +412,7 @@ CONTAINS
!
CASE( np_seos ) !== simplified EOS ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - rn_T0
zs = pts (ji,jj,jk,jp_sal) - rn_S0
zh = pdep (ji,jj,jk)
@@ -440,17 +460,30 @@ CONTAINS
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrd, ktrhop, ktdep
- REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
- ! ! 2 : salinity [psu]
- REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
- REAL(wp), DIMENSION(A2D_T(ktrhop),JPK ), INTENT( out) :: prhop ! potential density (surface referenced)
- REAL(wp), DIMENSION(A2D_T(ktdep) ,JPK ), INTENT(in ) :: pdep ! depth [m]
+!!$ REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+!!$ ! ! 2 : salinity [psu]
+!!$ REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
+!!$ REAL(wp), DIMENSION(A2D_T(ktrhop),JPK ), INTENT( out) :: prhop ! potential density (surface referenced)
+!!$ REAL(wp), DIMENSION(A2D_T(ktdep) ,JPK ), INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
+ REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prhop ! potential density (surface referenced)
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m]
!
INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
INTEGER :: jdof
REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-pot')
@@ -469,7 +502,7 @@ CONTAINS
zsign(jsmp+1) = -1._wp
END DO
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
!
! compute density (2*nn_sto_eos) times:
! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
@@ -519,7 +552,7 @@ CONTAINS
DEALLOCATE(zn0_sto,zn_sto,zsign)
! Non-stochastic equation of state
ELSE
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
!
zh = pdep(ji,jj,jk) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
@@ -557,7 +590,7 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - rn_T0
zs = pts (ji,jj,jk,jp_sal) - rn_S0
zh = pdep (ji,jj,jk)
@@ -605,14 +638,26 @@ CONTAINS
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktdep, ktrd
- REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
- ! ! 2 : salinity [psu]
- REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
- REAL(wp), DIMENSION(A2D_T(ktrd) ), INTENT( out) :: prd ! in situ density
+!!$ REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+!!$ ! ! 2 : salinity [psu]
+!!$ REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
+!!$ REAL(wp), DIMENSION(A2D_T(ktrd) ), INTENT( out) :: prd ! in situ density
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:) , INTENT( out) :: prd ! in situ density
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos2d')
@@ -623,7 +668,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
!
zh = pdep(ji,jj) * r1_Z0 ! depth
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
@@ -658,7 +703,7 @@ CONTAINS
!
CASE( np_seos ) !== simplified EOS ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
!
zt = pts (ji,jj,jp_tem) - rn_T0
zs = pts (ji,jj,jp_sal) - rn_S0
@@ -705,15 +750,26 @@ CONTAINS
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrhop
- REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
- ! ! 2 : salinity [psu]
- REAL(wp), DIMENSION(A2D_T(ktrhop) ), INTENT( out) :: prhop ! potential density (surface referenced)
+!!$ REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+!!$ ! ! 2 : salinity [psu]
+!!$ REAL(wp), DIMENSION(A2D_T(ktrhop) ), INTENT( out) :: prhop ! potential density (surface referenced)
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(:,:) , INTENT( out) :: prhop ! potential density (surface referenced)
!
INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
INTEGER :: jdof
REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-pot')
@@ -722,7 +778,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
!
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
@@ -743,7 +799,7 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
zt = pts (ji,jj,jp_tem) - rn_T0
zs = pts (ji,jj,jp_sal) - rn_S0
ztm = tmask(ji,jj,1)
@@ -787,12 +843,22 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER , INTENT(in ) :: ktts, ktab
- REAL(wp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
- REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
+!!$ REAL(wp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
+!!$ REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('rab_3d')
@@ -801,7 +867,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
@@ -854,7 +920,7 @@ CONTAINS
!
CASE( np_seos ) !== simplified EOS ==!
!
- DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - rn_T0 ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
@@ -903,13 +969,24 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER , INTENT(in ) :: ktts, ktdep, ktab
- REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
- REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
- REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
+!!$ REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
+!!$ REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
+!!$ REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(pts,1) ! 1st dimension
+ ipj = SIZE(pts,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('rab_2d')
@@ -920,7 +997,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
!
zh = pdep(ji,jj) * r1_Z0 ! depth
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
@@ -973,7 +1050,7 @@ CONTAINS
!
CASE( np_seos ) !== simplified EOS ==!
!
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
!
zt = pts (ji,jj,jp_tem) - rn_T0 ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
@@ -1219,6 +1296,81 @@ CONTAINS
!
END FUNCTION eos_pt_from_ct
+ SUBROUTINE eos_fzp_3d( psal, ptf, pdep )
+ !!
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: psal ! salinity [psu]
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ), OPTIONAL :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:,:), INTENT(out ) :: ptf ! freezing temperature [Celsius]
+ !!
+ CALL eos_fzp_3d_t( psal, ptf, is_tile(ptf), pdep )
+ END SUBROUTINE eos_fzp_3d
+
+ SUBROUTINE eos_fzp_3d_t( psal, ptf, kttf, pdep )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_fzp ***
+ !!
+ !! ** Purpose : Compute the freezing point temperature [Celsius]
+ !!
+ !! ** Method : UNESCO freezing point (ptf) in Celsius is given by
+ !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
+ !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
+ !!
+ !! Reference : UNESCO tech. papers in the marine science no. 28. 1978
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kttf
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: psal ! salinity [psu]
+ REAL(wp), DIMENSION(:,:,:), INTENT(in ), OPTIONAL :: pdep ! depth [m]
+ REAL(wp), DIMENSION(:,:,:), INTENT(out ) :: ptf ! freezing temperature [Celsius]
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt, zs, z1_S0 ! local scalars
+ INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
+ !!----------------------------------------------------------------------
+ !
+ ipi = SIZE(psal,1) ! 1st dimension
+ ipj = SIZE(psal,2) ! 2nd dimension
+ !
+ iisht = ( jpi - ipi ) / 2
+ ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
+ !
+ SELECT CASE ( neos )
+ !
+ CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
+ !
+ z1_S0 = 1._wp / 35.16504_wp
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpk )
+ zs= SQRT( ABS( psal(ji-iisht,jj-ijsht,jk) ) * z1_S0 ) ! square root salinity
+ ptf(ji-iisht,jj-ijsht,jk) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+ & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
+ ptf(ji-iisht,jj-ijsht,jk) = ptf(ji-iisht,jj-ijsht,jk) * psal(ji-iisht,jj-ijsht,jk)
+ END_3D
+ !
+ IF( PRESENT( pdep ) ) THEN
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpk )
+ ptf(ji-iisht,jj-ijsht,jk) = ptf(ji-iisht,jj-ijsht,jk) - 7.53e-4 * pdep(ji-iisht,jj-ijsht,jk)
+ END_3D
+ ENDIF
+ !
+ CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
+ !
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpk )
+ ptf(ji-iisht,jj-ijsht,jk) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(ji-iisht,jj-ijsht,jk) ) &
+ & - 2.154996e-4_wp * psal(ji-iisht,jj-ijsht,jk) ) * psal(ji-iisht,jj-ijsht,jk)
+ END_3D
+ !
+ IF( PRESENT( pdep ) ) THEN
+ DO_3D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht, 1, jpk )
+ ptf(ji-iisht,jj-ijsht,jk) = ptf(ji-iisht,jj-ijsht,jk) - 7.53e-4 * pdep(ji-iisht,jj-ijsht,jk)
+ END_3D
+ ENDIF
+ !
+ CASE DEFAULT
+ WRITE(ctmp1,*) ' bad flag value for neos = ', neos
+ CALL ctl_stop( 'eos_fzp_3d:', ctmp1 )
+ !
+ END SELECT
+ !
+ END SUBROUTINE eos_fzp_3d_t
SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
!!
diff --git a/src/OCE/TRA/traisf.F90 b/src/OCE/TRA/traisf.F90
index 610f9064f..d97faba3b 100644
--- a/src/OCE/TRA/traisf.F90
+++ b/src/OCE/TRA/traisf.F90
@@ -60,16 +60,20 @@ CONTAINS
!
#if defined key_RK3
! cavity case (RK3)
- IF ( ln_isfcav_mlt ) CALL isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, pts(:,:,:,:,Krhs))
+ IF ( ln_isfcav_mlt ) CALL isf_mlt( misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, &
+ & risf_cav_tsc, pts(A2D(0),:,:,Krhs) )
!
! parametrisation case (RK3)
- IF ( ln_isfpar_mlt ) CALL isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, pts(:,:,:,:,Krhs))
+ IF ( ln_isfpar_mlt ) CALL isf_mlt( misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, &
+ & risf_par_tsc, pts(A2D(0),:,:,Krhs) )
#else
! cavity case (MLF)
- IF ( ln_isfcav_mlt ) CALL isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, pts(:,:,:,:,Krhs), risf_cav_tsc_b)
+ IF ( ln_isfcav_mlt ) CALL isf_mlt( misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, &
+ & risf_cav_tsc, pts(A2D(0),:,:,Krhs), risf_cav_tsc_b )
!
! parametrisation case (MLF)
- IF ( ln_isfpar_mlt ) CALL isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, pts(:,:,:,:,Krhs), risf_par_tsc_b)
+ IF ( ln_isfpar_mlt ) CALL isf_mlt( misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, &
+ & risf_par_tsc, pts(A2D(0),:,:,Krhs), risf_par_tsc_b )
#endif
!
! ice sheet coupling case
@@ -81,11 +85,11 @@ CONTAINS
! half of it at nit000+1 (leap frog time step).
! in accordance to this, the heat content flux due to injected water need to be added in the temperature and salt trend
! at time step nit000 and nit000+1
- IF ( kt == nit000 ) CALL isf_cpl(Kmm, risfcpl_tsc , pts(:,:,:,:,Krhs))
- IF ( kt == nit000+1) CALL isf_cpl(Kmm, risfcpl_tsc*0.5_wp, pts(:,:,:,:,Krhs))
+ IF ( kt == nit000 ) CALL isf_cpl( Kmm, risfcpl_tsc , pts(A2D(0),:,:,Krhs) )
+ IF ( kt == nit000+1) CALL isf_cpl( Kmm, risfcpl_tsc*0.5_wp, pts(A2D(0),:,:,Krhs) )
!
! ensure 0 trend due to unconservation of the ice shelf coupling
- IF ( ln_isfcpl_cons ) CALL isf_cpl(Kmm, risfcpl_cons_tsc, pts(:,:,:,:,Krhs))
+ IF ( ln_isfcpl_cons ) CALL isf_cpl( Kmm, risfcpl_cons_tsc, pts(A2D(0),:,:,Krhs) )
!
END IF
!
@@ -113,39 +117,34 @@ CONTAINS
!! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend
!!
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) , INTENT(inout) :: pts
- INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl, pfrac
- REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: ptsc
- REAL(wp), DIMENSION(:,:,:) , OPTIONAL, INTENT(in ) :: ptsc_b
+ REAL(wp), DIMENSION(A2D(0),jpk,jpts) , INTENT(inout) :: pts
+ INTEGER , DIMENSION(A2D(0)) , INTENT(in ) :: ktop , kbot
+ REAL(wp), DIMENSION(A2D(0)) , INTENT(in ) :: phtbl, pfrac
+ REAL(wp), DIMENSION(A2D(0),jpts) , INTENT(in ) :: ptsc
+ REAL(wp), DIMENSION(A2D(0),jpts) , OPTIONAL, INTENT(in ) :: ptsc_b
!!
- INTEGER :: ji,jj,jk ! dummy loop index
- INTEGER :: ikt, ikb ! top and bottom level of the tbl
- REAL(wp), DIMENSION(A2D(nn_hls)) :: ztc ! total ice shelf tracer trend
+ INTEGER :: ji, jj, jk ! dummy loop index
+ INTEGER :: ikt, ikb ! top and bottom level of the tbl
+ REAL(wp) :: ztc ! total ice shelf tracer trend
!!----------------------------------------------------------------------
!
- ! compute 2d total trend due to isf
+ ! update pts(:,:,:,:,Krhs)
DO_2D( 0, 0, 0, 0 )
+ !
#if defined key_RK3
- ztc(ji,jj) = ptsc(ji,jj,jp_tem) / phtbl(ji,jj)
+ ztc = ptsc(ji,jj,jp_tem) / phtbl(ji,jj)
#else
- ztc(ji,jj) = 0.5_wp * ( ptsc(ji,jj,jp_tem) + ptsc_b(ji,jj,jp_tem) ) / phtbl(ji,jj)
+ ztc = 0.5_wp * ( ptsc(ji,jj,jp_tem) + ptsc_b(ji,jj,jp_tem) ) / phtbl(ji,jj)
#endif
- END_2D
- !
- ! update pts(:,:,:,:,Krhs)
- DO_2D( 0, 0, 0, 0 )
- !
+ ! level fully include in the ice shelf boundary layer
ikt = ktop(ji,jj)
ikb = kbot(ji,jj)
- !
- ! level fully include in the ice shelf boundary layer
DO jk = ikt, ikb - 1
- pts(ji,jj,jk,jp_tem) = pts(ji,jj,jk,jp_tem) + ztc(ji,jj)
+ pts(ji,jj,jk,jp_tem) = pts(ji,jj,jk,jp_tem) + ztc
END DO
!
! level partially include in ice shelf boundary layer
- pts(ji,jj,ikb,jp_tem) = pts(ji,jj,ikb,jp_tem) + ztc(ji,jj) * pfrac(ji,jj)
+ pts(ji,jj,ikb,jp_tem) = pts(ji,jj,ikb,jp_tem) + ztc * pfrac(ji,jj)
!
END_2D
!
@@ -159,9 +158,9 @@ CONTAINS
!! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend
!!
!!----------------------------------------------------------------------
- INTEGER , INTENT(in ) :: Kmm ! ocean time-level index
- REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: ptsc
- REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: ptsa
+ INTEGER , INTENT(in ) :: Kmm ! ocean time-level index
+ REAL(wp), DIMENSION(A2D(0),jpk,jpts), INTENT(in ) :: ptsc
+ REAL(wp), DIMENSION(A2D(0),jpk,jpts), INTENT(inout) :: ptsa
!!
INTEGER :: ji, jj, jk ! dummy loop index
!!----------------------------------------------------------------------
diff --git a/src/OCE/step.F90 b/src/OCE/step.F90
index d7036efb3..10629ea20 100644
--- a/src/OCE/step.F90
+++ b/src/OCE/step.F90
@@ -156,7 +156,9 @@ CONTAINS
IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries
IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice)
-
+!!$ IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
+ !clem: problem with isf and cpl: sbcfwb needs isf but isf needs fwf from sbccpl
+
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update stochastic parameters and random T/S fluctuations
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 2afa912fc..c8119ef2f 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -164,6 +164,8 @@ CONTAINS
IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries
IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice)
+!!$ IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
+ !clem: problem with isf and cpl: sbcfwb needs isf but isf needs fwf from sbccpl
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update stochastic parameters and random T/S fluctuations
diff --git a/src/OCE/stprk3.F90 b/src/OCE/stprk3.F90
index 337b55f62..51ad96fca 100644
--- a/src/OCE/stprk3.F90
+++ b/src/OCE/stprk3.F90
@@ -127,6 +127,8 @@ CONTAINS
IF( ln_bdy ) CALL bdy_dta ( kstp, Nbb ) ! update dynamic & tracer data at open boundaries
IF( ln_isf ) CALL isf_stp ( kstp, Nbb ) ! update iceshelf geometry
CALL sbc ( kstp, Nbb, Nbb ) ! Sea Boundary Condition (including sea-ice)
+!!$ IF( ln_isf ) CALL isf_stp ( kstp, Nbb ) ! update iceshelf geometry
+ !clem: problem with isf and cpl: sbcfwb needs isf but isf needs fwf from sbccpl
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update stochastic parameters and random T/S fluctuations
diff --git a/src/OFF/dtadyn.F90 b/src/OFF/dtadyn.F90
index 3a68e3128..a5649b3b1 100644
--- a/src/OFF/dtadyn.F90
+++ b/src/OFF/dtadyn.F90
@@ -136,15 +136,15 @@ CONTAINS
!
IF( l_ldfslp .AND. .NOT.ln_c1d ) CALL dta_dyn_slp( kt, Kbb, Kmm ) ! Computation of slopes
!
- ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
- ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
- wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange
- fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)
- fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction
- qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation
- emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P
+ ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask (:,:,:) ! temperature
+ ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask (:,:,:) ! salinity
+ wndm(:,:) = sf_dyn(jf_wnd)%fnow(A2D(0),1) * smask0(:,:) ! wind speed - needed for gas exchange
+ fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(A2D(0),1) * smask0(:,:) ! downward salt flux (v3.5+)
+ fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * ssmask(:,:) ! Sea-ice fraction
+ qsr (:,:) = sf_dyn(jf_qsr)%fnow(A2D(0),1) * smask0(:,:) ! solar radiation
+ emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * ssmask(:,:) ! E-P
IF( ln_dynrnf ) THEN
- rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P
+ rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * ssmask(:,:) ! E-P
IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_rnf( Kmm )
ENDIF
!
@@ -182,8 +182,8 @@ CONTAINS
rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl
CALL zdf_mxl( kt, Kmm ) ! In any case, we need mxl
!
- hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht
- avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient
+ hmld(:,:) = sf_dyn(jf_mld)%fnow(A2D(0),1) * smask0(:,:) ! mixed layer depht
+ avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient
avs(:,:,:) = avt(:,:,:)
!
IF( ln_trabbl .AND. .NOT.ln_c1d ) THEN ! diffusive Bottom boundary layer param
@@ -506,7 +506,7 @@ CONTAINS
CALL iom_close( inum ) ! close file
!
nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
IF( h_rnf(ji,jj) > 0._wp ) THEN
jk = 2
DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
@@ -521,15 +521,17 @@ CONTAINS
END_2D
!
! set the associated depth
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj)
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
END DO
END_2D
ELSE ! runoffs applied at the surface
- nk_rnf(:,:) = 1
- h_rnf (:,:) = e3t(:,:,1,Kmm)
+ DO_2D( 0, 0, 0, 0 )
+ nk_rnf(ji,jj) = 1
+ h_rnf (ji,jj) = e3t(ji,jj,1,Kmm)
+ END_2D
ENDIF
nkrnf_max = MAXVAL( nk_rnf(:,:) )
hrnf_max = MAXVAL( h_rnf(:,:) )
@@ -557,7 +559,7 @@ CONTAINS
!!----------------------------------------------------------------------
!
! update the depth over which runoffs are distributed
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! to the bottom of the relevant grid box
diff --git a/tests/ISOMIP+/MY_SRC/isf_oce.F90 b/tests/ISOMIP+/MY_SRC/isf_oce.F90
index 01aab242c..e013cddb9 100644
--- a/tests/ISOMIP+/MY_SRC/isf_oce.F90
+++ b/tests/ISOMIP+/MY_SRC/isf_oce.F90
@@ -13,8 +13,7 @@ MODULE isf_oce
!! isf : define and allocate ice shelf variables
!!----------------------------------------------------------------------
- USE par_oce , ONLY: jpi, jpj, jpk
- USE in_out_manager, ONLY: wp, jpts ! I/O manager
+ USE par_oce
USE lib_mpp , ONLY: ctl_stop, mpp_sum ! MPP library
USE fldread ! read input fields
@@ -22,7 +21,7 @@ MODULE isf_oce
PRIVATE
- PUBLIC isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl, isf_dealloc_cpl
+ PUBLIC isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl
!
!-------------------------------------------------------
! 0 : namelist parameter
@@ -65,7 +64,7 @@ MODULE isf_oce
!
REAL(wp), PARAMETER, PUBLIC :: rLfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg]
REAL(wp), PARAMETER, PUBLIC :: rcpisf = 2000.0_wp !: specific heat of ice shelf [J/kg/K]
- REAL(wp), PARAMETER, PUBLIC :: rkappa = 0._wp !: ISOMIP: no heat diffusivity [m2/s]
+ REAL(wp), PARAMETER, PUBLIC :: rkappa = 0._wp !: ISOMIP: no heat diffusivity through the ice-shelf [m2/s]
REAL(wp), PARAMETER, PUBLIC :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3]
REAL(wp), PARAMETER, PUBLIC :: rtsurf = -20.0 !: surface temperature [C]
!
@@ -79,33 +78,35 @@ MODULE isf_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_oasis
!
! 2.2 -------- ice shelf cavity melt namelist parameter -------------
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
- TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
+ TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
!
- REAL(wp) , PUBLIC :: risf_lamb1, risf_lamb2, risf_lamb3 ! freezing point linearization coeficient
+ REAL(wp) , PUBLIC :: risf_lamb1, risf_lamb2, risf_lamb3 !: freezing point linearization coeficient
!
! 2.3 -------- ice shelf param. melt namelist parameter -------------
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_par !:
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_par , misfkb_par !:
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_par, rfrac_tbl_par !:
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_par , fwfisf_par_b !: before and now net fwf from the ice shelf [kg/m2/s]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
- TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfpar_fwf !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_par !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_par , misfkb_par !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_par, rfrac_tbl_par !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_par , fwfisf_par_b !: before and now net fwf from the ice shelf [kg/m2/s]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
+ TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfpar_fwf !:
!
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf0_tbl_par !: thickness of tbl (initial value) [m]
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !:
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf0_tbl_par !: thickness of tbl (initial value) [m]
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !:
!
! 2.4 -------- coupling namelist parameter -------------
INTEGER , PUBLIC :: nstp_iscpl !:
REAL(wp), PUBLIC :: rdt_iscpl !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfcpl_ssh, risfcpl_cons_ssh, risfcpl_cons_ssh_b !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risfcpl_vol, risfcpl_cons_vol, risfcpl_cons_vol_b !:
- REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: risfcpl_tsc, risfcpl_cons_tsc, risfcpl_cons_tsc_b !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfcpl_ssh, risfcpl_cons_ssh !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risfcpl_vol, risfcpl_cons_vol !:
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: risfcpl_tsc, risfcpl_cons_tsc !:
!
+ !! * Substitutions
+# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
@@ -125,20 +126,12 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0 ! set to zero if no array to be allocated
- !
- ALLOCATE(risfLeff(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE(misfkt_par(jpi,jpj), misfkb_par(jpi,jpj), STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( rfrac_tbl_par(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( rhisf_tbl_par(jpi,jpj), rhisf0_tbl_par(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( mskisf_par(jpi,jpj), STAT=ialloc)
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( misfkt_par (A2D(0)) , misfkb_par (A2D(0)) , rfrac_tbl_par(A2D(0)) , &
+ & rhisf_tbl_par(A2D(0)) , rhisf0_tbl_par(A2D(0)) , &
+ & risfLeff (A2D(0)) , mskisf_par (A2D(0)) , STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
@@ -159,14 +152,10 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0 ! set to zero if no array to be allocated
- !
- ALLOCATE(misfkt_cav(jpi,jpj), misfkb_cav(jpi,jpj), STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( rfrac_tbl_cav(jpi,jpj), STAT=ialloc)
- ierr = ierr + ialloc
- !
- ALLOCATE( rhisf_tbl_cav(jpi,jpj), STAT=ialloc)
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( misfkt_cav(A2D(0)), misfkb_cav(A2D(0)), rfrac_tbl_cav(A2D(0)) , rhisf_tbl_cav(A2D(0)) , STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
@@ -185,46 +174,25 @@ CONTAINS
INTEGER :: ierr, ialloc
!!----------------------------------------------------------------------
ierr = 0
- !
- ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts) , risfcpl_vol(jpi,jpj,jpk) , STAT=ialloc )
+ ! ----------------- !
+ ! == FULL ARRAYS == !
+ ! ----------------- !
+ ALLOCATE( risfcpl_ssh (jpi,jpj) , risfcpl_vol (jpi,jpj,jpk) , &
+ & risfcpl_cons_ssh(jpi,jpj) , risfcpl_cons_vol(jpi,jpj,jpk) , STAT=ialloc )
+ ierr = ierr + ialloc
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( risfcpl_tsc (A2D(0),jpk,jpts) , &
+ & risfcpl_cons_tsc(A2D(0),jpk,jpts) , STAT=ialloc )
ierr = ierr + ialloc
- !
- risfcpl_tsc(:,:,:,:) = 0._wp ; risfcpl_vol(:,:,:) = 0._wp ; risfcpl_ssh(:,:) = 0._wp
-
- IF ( ln_isfcpl_cons ) THEN
- ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- risfcpl_cons_tsc(:,:,:,:) = 0._wp ; risfcpl_cons_vol(:,:,:) = 0._wp ; risfcpl_cons_ssh(:,:) = 0._wp
- !
- END IF
!
CALL mpp_sum ( 'isf', ierr )
IF( ierr /= 0 ) CALL ctl_stop('STOP','isfcpl: failed to allocate arrays.')
!
END SUBROUTINE isf_alloc_cpl
-
- SUBROUTINE isf_dealloc_cpl()
- !!---------------------------------------------------------------------
- !! *** ROUTINE isf_dealloc_cpl ***
- !!
- !! ** Purpose : de-allocate useless public 3d array used for ice sheet coupling
- !!
- !!----------------------------------------------------------------------
- INTEGER :: ierr, ialloc
- !!----------------------------------------------------------------------
- ierr = 0
- !
- DEALLOCATE( risfcpl_ssh , risfcpl_tsc , risfcpl_vol , STAT=ialloc )
- ierr = ierr + ialloc
- !
- CALL mpp_sum ( 'isf', ierr )
- IF( ierr /= 0 ) CALL ctl_stop('STOP','isfcpl: failed to deallocate arrays.')
- !
- END SUBROUTINE isf_dealloc_cpl
-
-
+
SUBROUTINE isf_alloc()
!!---------------------------------------------------------------------
!! *** ROUTINE isf_alloc ***
@@ -237,52 +205,29 @@ CONTAINS
!
ierr = 0 ! set to zero if no array to be allocated
!
- ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , &
- & fwfisf_oasis(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
+ ! ----------------- !
+ ! == FULL ARRAYS == !
+ ! ----------------- !
+ ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , risfload(jpi,jpj) , &
#if ! defined key_RK3
- ! MLF : need to allocate before arrays
- ALLOCATE( fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
- !
- ALLOCATE( risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
- ierr = ierr + ialloc
+ & fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , & ! MLF : need to allocate before arrays
#endif
- !
- ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
+ & STAT=ialloc )
ierr = ierr + ialloc
!
- ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc )
+ ! -------------------- !
+ ! == REDUCED ARRAYS == !
+ ! -------------------- !
+ ALLOCATE( fwfisf_oasis(A2D(0)) , risf_par_tsc (A2D(0),jpts) , risf_cav_tsc (A2D(0),jpts) , mskisf_cav(A2D(0)) , &
+#if ! defined key_RK3
+ & risf_par_tsc_b(A2D(0),jpts) , risf_cav_tsc_b(A2D(0),jpts) , & ! MLF : need to allocate before arrays
+#endif
+ & STAT=ialloc )
ierr = ierr + ialloc
!
CALL mpp_sum ( 'isf', ierr )
IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
!
- ! initalisation of fwf and tsc array to 0
- risfload (:,:) = 0._wp
- fwfisf_oasis(:,:) = 0._wp
-#if defined key_RK3
- fwfisf_par (:,:) = 0._wp
- fwfisf_cav (:,:) = 0._wp
- risf_cav_tsc(:,:,:) = 0._wp
- risf_par_tsc(:,:,:) = 0._wp
-#else
- fwfisf_par (:,:) = 0._wp ; fwfisf_par_b (:,:) = 0._wp
- fwfisf_cav (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp
- risf_cav_tsc(:,:,:) = 0._wp ; risf_cav_tsc_b(:,:,:) = 0._wp
- risf_par_tsc(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp
-#endif
- !
END SUBROUTINE isf_alloc
!!======================================================================
--
GitLab