MODULE icedyn_adv_pra
!!======================================================================
!! *** MODULE icedyn_adv_pra ***
!! sea-ice : advection => Prather scheme
!!======================================================================
!! History : ! 2008-03 (M. Vancoppenolle) original code
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!--------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_dyn_adv_pra : advection of sea ice using Prather scheme
!! adv_x, adv_y : Prather scheme applied in i- and j-direction, resp.
!! adv_pra_init : initialisation of the Prather scheme
!! adv_pra_rst : read/write Prather field in ice restart file, or initialized to zero
!!----------------------------------------------------------------------
USE phycst ! physical constant
USE dom_oce ! ocean domain
USE ice ! sea-ice variables
USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call
USE icevar ! sea-ice: operations
!
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE
PRIVATE
PUBLIC ice_dyn_adv_pra ! called by icedyn_adv
PUBLIC adv_pra_init ! called by icedyn_adv
! Moments for advection
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvl , syvl , sxxvl , syyvl , sxyvl ! melt pond lid volume
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn_adv_pra.F90 15049 2021-06-23 16:17:30Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
!!----------------------------------------------------------------------
!! ** routine ice_dyn_adv_pra **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!!
!! ** method : Uses Prather second order scheme that advects tracers
!! but also their quadratic forms. The method preserves
!! tracer structures by conserving second order moments.
!!
!! Reference: Prather, 1986, JGR, 91, D6. 6671-6681.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid thickness
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
!
INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zdt, z1_dt ! - -
REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0vl
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
!! diagnostics
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
!
! --- Record max of the surrounding 9-pts (for call Hbig) --- !
! thickness and salinity
WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:)
ELSEWHERE ; zs_i(:,:,:) = 0._wp
END WHERE
CALL icemax3D( ph_i , zhi_max )
CALL icemax3D( ph_s , zhs_max )
CALL icemax3D( ph_ip, zhip_max)
CALL icemax3D( zs_i , zsi_max )
CALL lbc_lnk( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
!
! enthalpies
DO jk = 1, nlay_i
WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:)
ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp
END WHERE
END DO
DO jk = 1, nlay_s
WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:)
ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp
END WHERE
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
! this should not affect too much the stability
zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) )
zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) )
! non-blocking global communication send zcflnow and receive zcflprv
CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
IF( zcflprv(1) > .5 ) THEN ; icycle = 2
ELSE ; icycle = 1
ENDIF
zdt = rDt_ice / REAL(icycle)
z1_dt = 1._wp / zdt
! --- transport --- !
zudy(:,:) = pu_ice(:,:) * e2u(:,:)
zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
! --- transported fields --- !
DO jl = 1, jpl
zarea(:,:,jl) = e1e2t(:,:)
z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume
z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume
z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area
z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content
z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content
DO jk = 1, nlay_s
z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
END DO
DO jk = 1, nlay_i
z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction
z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume
IF ( ln_pnd_lids ) THEN
z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:) ! Melt pond lid volume
ENDIF
ENDIF
END DO
!
! !--------------------------------------------!
IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
! !--------------------------------------------!
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
!
DO jk = 1, nlay_s !--- snow heat content
CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
END DO
DO jk = 1, nlay_i !--- ice heat content
CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
END DO
!
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
IF ( ln_pnd_lids ) THEN
CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
ENDIF
ENDIF
! !--------------------------------------------!
ELSE !== even ice time step: adv_y then adv_x ==!
! !--------------------------------------------!
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
DO jk = 1, nlay_s !--- snow heat content
CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
END DO
DO jk = 1, nlay_i !--- ice heat content
CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
IF ( ln_pnd_lids ) THEN
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
ENDIF
ENDIF
!
ENDIF
! --- Lateral boundary conditions --- !
! caution: for gradients (sx and sy) the sign changes
CALL lbc_lnk( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume
& , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp &
& , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume
& , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp &
& , z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity
& , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp &
& , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration
& , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp &
& , z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age
& , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy
& , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy
& , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp )
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
IF( ln_pnd_lids ) THEN
CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction
& , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp &
& , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume
& , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp &
& , z0vl , 'T', 1._wp, sxvl , 'T', -1._wp, syvl , 'T', -1._wp & ! melt pond lid volume
& , sxxvl, 'T', 1._wp, syyvl, 'T', 1._wp, sxyvl, 'T', 1._wp )
ELSE
CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction
& , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp &
& , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume
& , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp )
ENDIF
ENDIF
! --- Recover the properties from their contents --- !
DO jl = 1, jpl
pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
DO jk = 1, nlay_s
pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
END DO
DO jk = 1, nlay_i
pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
IF ( ln_pnd_lids ) THEN
pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
ENDIF
ENDIF
END DO
!
! derive open water from ice concentration
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
DO_2D( 0, 0, 0, 0 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
END_2D
CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp )
!
! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
& - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt
!
! --- Ensure non-negative fields --- !
! Remove negative values (conservation is ensured)
! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
!
! --- Make sure ice thickness is not too big --- !
! (because ice thickness can be too large where ice concentration is very small)
CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, &
& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
!
! --- Ensure snow load is not too big --- !
CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
!
END DO
!
IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file
!
END SUBROUTINE ice_dyn_adv_pra
SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , &
& psx, psxx, psy , psyy, psxy )
!!----------------------------------------------------------------------
!! ** routine adv_x **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!! variable on x axis
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! the time step
REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
!!
INTEGER :: ji, jj, jl, jcat ! dummy loop indices
INTEGER :: jj0 ! dummy loop indices
REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars
REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
REAL(wp) :: zpsm, zps0
REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
!-----------------------------------------------------------------------
! in order to avoid lbc_lnk (communications):
! jj loop must be 1:jpj if adv_x is called first
! and 2:jpj-1 if adv_x is called second
jj0 = NINT(pcrh)
!
jcat = SIZE( ps0 , 3 ) ! size of input arrays
!
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
DO jj = Njs0 - jj0, Nje0 + jj0
DO ji = Nis0 - 1, Nie0 + 1
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
!
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
zps0 = zslpmax
zpsx = zs1new * rswitch
zpsxx = zs2new * rswitch
zpsy = zpsy * rswitch
zpsyy = zpsyy * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes i<-->i+1
! ! Flux from i to i+1 WHEN u GT 0
zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
zfm (ji,jj) = zalf * zpsm
zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
zfxx(ji,jj) = zalf * zpsxx * zalfq
zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
zfxy(ji,jj) = zalfq * zpsxy
zfyy(ji,jj) = zalf * zpsyy
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
zpsxx = zalf1 * zalf1q * zpsxx
zpsy = zpsy - zfy (ji,jj)
zpsyy = zpsyy - zfyy(ji,jj)
zpsxy = zalf1q * zpsxy
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
!
END DO
DO ji = Nis0 - 1, Nie0
! ! Flux from i+1 to i when u LT 0.
zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
zalg (ji,jj) = zalf
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalg1 (ji,jj) = zalf1
zalf1q = zalf1 * zalf1
zalg1q(ji,jj) = zalf1q
!
zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
& - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq
zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
END DO
DO ji = Nis0, Nie0
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
! ! Readjust moments remaining in the box.
zbt = zbet(ji-1,jj)
zbt1 = 1.0 - zbet(ji-1,jj)
!
zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
zpsxy = zalg1q(ji-1,jj) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from i to i+1 IF u GT 0.
zbt = zbet(ji-1,jj)
zbt1 = 1.0 - zbet(ji-1,jj)
zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
zalf = zbt * zfm(ji-1,jj) / zpsm
zalf1 = 1.0 - zalf
ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
!
zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx &
& + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
& + zbt1 * zpsxx
zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy &
& + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) &
& + zbt1 * zpsxy
zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy
zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
! ! Flux from i+1 to i IF u LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
zalf = zbt1 * zfm(ji,jj) / zpsm
zalf1 = 1.0 - zalf
ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
!
zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
& + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) &
& + ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
& + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) )
zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
END DO
!
END DO
!
END DO
!
END SUBROUTINE adv_x
SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , &
& psx, psxx, psy , psyy, psxy )
!!---------------------------------------------------------------------
!! ** routine adv_y **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!! variable on y axis
!!---------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! time step
REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
!!
INTEGER :: ji, jj, jl, jcat ! dummy loop indices
INTEGER :: ji0 ! dummy loop indices
REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars
REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
REAL(wp) :: zpsm, zps0
REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
!---------------------------------------------------------------------
! in order to avoid lbc_lnk (communications):
! ji loop must be 1:jpi if adv_y is called first
! and 2:jpi-1 if adv_y is called second
ji0 = NINT(pcrh)
!
jcat = SIZE( ps0 , 3 ) ! size of input arrays
!
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
DO_2D( ji0, ji0, 1, 1 )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
!
! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
!
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) )
zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
!
zps0 = zslpmax
zpsx = zpsx * rswitch
zpsxx = zpsxx * rswitch
zpsy = zs1new * rswitch
zpsyy = zs2new * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes j<-->j+1
! ! Flux from j to j+1 WHEN v GT 0
zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
zfm (ji,jj) = zalf * zpsm
zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) )
zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy )
zfyy(ji,jj) = zalf * zalfq * zpsyy
zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
zfxy(ji,jj) = zalfq * zpsxy
zfxx(ji,jj) = zalf * zpsxx
!
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy )
zpsyy = zalf1 * zalf1q * zpsyy
zpsx = zpsx - zfx(ji,jj)
zpsxx = zpsxx - zfxx(ji,jj)
zpsxy = zalf1q * zpsxy
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
END_2D
!
DO_2D( ji0, ji0, 1, 0 )
! ! Flux from j+1 to j when v LT 0.
zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
zalg (ji,jj) = zalf
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalg1 (ji,jj) = zalf1
zalf1q = zalf1 * zalf1
zalg1q(ji,jj) = zalf1q
!
zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl)
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq
zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
END_2D
DO_2D( ji0, ji0, 0, 0 )
! ! Readjust moments remaining in the box.
zbt = zbet(ji,jj-1)
zbt1 = ( 1.0 - zbet(ji,jj-1) )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
!
zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
zpsxy = zalg1q(ji,jj-1) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from j to j+1 IF v GT 0.
zbt = zbet(ji,jj-1)
zbt1 = 1.0 - zbet(ji,jj-1)
zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm
zalf = zbt * zfm(ji,jj-1) / zpsm
zalf1 = 1.0 - zalf
ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
!
zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) &
& + zbt1 * zpsy
zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy &
& + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
& + zbt1 * zpsyy
zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy &
& + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) &
& + zbt1 * zpsxy
zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx
zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
! ! Flux from j+1 to j IF v LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
zalf = zbt1 * zfm(ji,jj) / zpsm
zalf1 = 1.0 - zalf
ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
!
zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
& + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) &
& + ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
& + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) )
zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
END_2D
!
END DO
!
END SUBROUTINE adv_y
SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
!!-------------------------------------------------------------------
!! *** ROUTINE Hbig ***
!!
!! ** Purpose : Thickness correction in case advection scheme creates
!! abnormally tick ice or snow
!!
!! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
!! (before advection) and reduce it by adapting ice concentration
!! 2- check whether snow thickness is larger than the surrounding 9-points
!! (before advection) and reduce it by sending the excess in the ocean
!!
!! ** input : Max thickness of the surrounding 9-points
!!-------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
!
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
!!-------------------------------------------------------------------
!
z1_dt = 1._wp / pdt
!
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
!
! ! -- check h_ip -- !
! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
ENDIF
ENDIF
!
! ! -- check h_i -- !
! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
ENDIF
!
! ! -- check h_s -- !
! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
!
wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
!
pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
ENDIF
!
! ! -- check s_i -- !
! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = psi_max(ji,jj,jl) / zsi
sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
ENDIF
!
ENDIF
END_2D
END DO
!
! ! -- check e_i/v_i -- !
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = pei_max(ji,jj,jk,jl) / zei
hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
ENDIF
ENDIF
END_3D
END DO
! ! -- check e_s/v_s -- !
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = pes_max(ji,jj,jk,jl) / zes
hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
ENDIF
ENDIF
END_3D
END DO
!
END SUBROUTINE Hbig
SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
!!-------------------------------------------------------------------
!! *** ROUTINE Hsnow ***
!!
!! ** Purpose : 1- Check snow load after advection
!! 2- Correct pond concentration to avoid a_ip > a_i
!!
!! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface
!! then put the snow excess in the ocean
!!
!! ** Notes : This correction is crucial because of the call to routine icecor afterwards
!! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
!! make the snow very thick (if concentration decreases drastically)
!! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
!!-------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
!
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_dt, zvs_excess, zfra
!!-------------------------------------------------------------------
!
z1_dt = 1._wp / pdt
!
! -- check snow load -- !
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
!
zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
!
IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface
! put snow excess in the ocean
zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
! correct snow volume and heat content
pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
ENDIF
!
ENDIF
END_2D
END DO
!
!-- correct pond concentration to avoid a_ip > a_i -- !
WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
!
END SUBROUTINE Hsnow
SUBROUTINE adv_pra_init
!!-------------------------------------------------------------------
!! *** ROUTINE adv_pra_init ***
!!
!! ** Purpose : allocate and initialize arrays for Prather advection
!!-------------------------------------------------------------------
INTEGER :: ierr
!!-------------------------------------------------------------------
!
! !* allocate prather fields
ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , &
& sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , &
& sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , &
& sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , &
& sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , &
& sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , &
& sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , &
& sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) , &
!
& sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
& syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl) , &
!
& sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
& syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl) , &
& STAT = ierr )
!
CALL mpp_sum( 'icedyn_adv_pra', ierr )
IF( ierr /= 0 ) CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
!
CALL adv_pra_rst( 'READ' ) !* read or initialize all required files
!
END SUBROUTINE adv_pra_init
SUBROUTINE adv_pra_rst( cdrw, kt )
!!---------------------------------------------------------------------
!! *** ROUTINE adv_pra_rst ***
!!
!! ** Purpose : Read or write file in restart file
!!
!! ** Method : use of IOM library
!!----------------------------------------------------------------------
CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
!
INTEGER :: jk, jl ! dummy loop indices
INTEGER :: iter ! local integer
INTEGER :: id1 ! local integer
CHARACTER(len=25) :: znam
CHARACTER(len=2) :: zchar, zchar1
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
! !==========================!
IF( TRIM(cdrw) == 'READ' ) THEN !== Read or initialize ==!
! !==========================!
!
IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0
ELSE ; id1 = 0 ! no restart: id1=0
ENDIF
!
IF( id1 > 0 ) THEN !** Read the restart file **!
!
! ! ice thickness
CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice )
CALL iom_get( numrir, jpdom_auto, 'syyice', syyice )
CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice )
! ! snow thickness
CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn )
CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn )
CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn )
! ! ice concentration
CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sya' , sya , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa )
CALL iom_get( numrir, jpdom_auto, 'syya' , syya )
CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya )
! ! ice salinity
CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal )
CALL iom_get( numrir, jpdom_auto, 'syysal', syysal )
CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal )
! ! ice age
CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage )
CALL iom_get( numrir, jpdom_auto, 'syyage', syyage )
CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage )
! ! snow layers heat content
DO jk = 1, nlay_s
WRITE(zchar1,'(I2.2)') jk
znam = 'sxc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)
znam = 'syc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; syc0 (:,:,jk,:) = z3d(:,:,:)
znam = 'sxxc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)
znam = 'syyc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)
znam = 'sxyc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)
END DO
! ! ice layers heat content
DO jk = 1, nlay_i
WRITE(zchar1,'(I2.2)') jk
znam = 'sxe'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxe (:,:,jk,:) = z3d(:,:,:)
znam = 'sye'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sye (:,:,jk,:) = z3d(:,:,:)
znam = 'sxxe'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:)
znam = 'syye'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:)
znam = 'sxye'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:)
END DO
!
IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction
IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
! ! melt pond volume
CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
ELSE
sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction
sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume
ENDIF
!
IF ( ln_pnd_lids ) THEN ! melt pond lid volume
IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl )
ELSE
sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume
ENDIF
ENDIF
ENDIF
!
ELSE !** start rheology from rest **!
!
IF(lwp) WRITE(numout,*) ' ==>> start from rest OR previous run without Prather, set moments to 0'
!
sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness
sxsn = 0._wp ; sysn = 0._wp ; sxxsn = 0._wp ; syysn = 0._wp ; sxysn = 0._wp ! snow thickness
sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration
sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity
sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age
sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content
sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content
IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction
sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume
IF ( ln_pnd_lids ) THEN
sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume
ENDIF
ENDIF
ENDIF
!
! !=====================================!
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN !== write in the ice restart file ==!
! !=====================================!
IF(lwp) WRITE(numout,*) '---- ice-adv-rst ----'
iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
!
!
! In case Prather scheme is used for advection, write second order moments
! ------------------------------------------------------------------------
!
! ! ice thickness
CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice )
CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice )
CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
! ! snow thickness
CALL iom_rstput( iter, nitrst, numriw, 'sxsn' , sxsn )
CALL iom_rstput( iter, nitrst, numriw, 'sysn' , sysn )
CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn )
CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn )
CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn )
! ! ice concentration
CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa )
CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya )
CALL iom_rstput( iter, nitrst, numriw, 'sxxa' , sxxa )
CALL iom_rstput( iter, nitrst, numriw, 'syya' , syya )
CALL iom_rstput( iter, nitrst, numriw, 'sxya' , sxya )
! ! ice salinity
CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal )
CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal )
CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
! ! ice age
CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage )
CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage )
CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
! ! snow layers heat content
DO jk = 1, nlay_s
WRITE(zchar1,'(I2.2)') jk
znam = 'sxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxc0 (:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'syc0'//'_l'//zchar1 ; z3d(:,:,:) = syc0 (:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'sxxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxxc0(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'syyc0'//'_l'//zchar1 ; z3d(:,:,:) = syyc0(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'sxyc0'//'_l'//zchar1 ; z3d(:,:,:) = sxyc0(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
END DO
! ! ice layers heat content
DO jk = 1, nlay_i
WRITE(zchar1,'(I2.2)') jk
znam = 'sxe'//'_l'//zchar1 ; z3d(:,:,:) = sxe (:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'sye'//'_l'//zchar1 ; z3d(:,:,:) = sye (:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'sxxe'//'_l'//zchar1 ; z3d(:,:,:) = sxxe(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'syye'//'_l'//zchar1 ; z3d(:,:,:) = syye(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
znam = 'sxye'//'_l'//zchar1 ; z3d(:,:,:) = sxye(:,:,jk,:)
CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
END DO
!
IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction
CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap )
CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap )
CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
! ! melt pond volume
CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp )
CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp )
CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
!
IF ( ln_pnd_lids ) THEN ! melt pond lid volume
CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl )
CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl )
CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
ENDIF
ENDIF
!
ENDIF
!
END SUBROUTINE adv_pra_rst
SUBROUTINE icemax3D( pice , pmax )
!!---------------------------------------------------------------------
!! *** ROUTINE icemax3D ***
!! ** Purpose : compute the max of the 9 points around
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pice ! input
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pmax ! output
!
REAL(wp), DIMENSION(Nis0:Nie0) :: zmax1, zmax2
REAL(wp) :: zmax3
INTEGER :: ji, jj, jl ! dummy loop indices
!!----------------------------------------------------------------------
! basic version: get the max of epsi20 + 9 neighbours
!!$ DO jl = 1, jpl
!!$ DO_2D( 0, 0, 0, 0 )
!!$ pmax(ji,jj,jl) = MAX( epsi20, pice(ji-1,jj-1,jl), pice(ji,jj-1,jl), pice(ji+1,jj-1,jl), &
!!$ & pice(ji-1,jj ,jl), pice(ji,jj ,jl), pice(ji+1,jj ,jl), &
!!$ & pice(ji-1,jj+1,jl), pice(ji,jj+1,jl), pice(ji+1,jj+1,jl) )
!!$ END_2D
!!$ END DO
! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
DO jl = 1, jpl
DO ji = Nis0, Nie0
zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jl), pice(ji-1,Njs0-1,jl), pice(ji+1,Njs0-1,jl) )
zmax2(ji) = MAX( epsi20, pice(ji,Njs0 ,jl), pice(ji-1,Njs0 ,jl), pice(ji+1,Njs0 ,jl) )
END DO
DO_2D( 0, 0, 0, 0 )
zmax3 = MAX( epsi20, pice(ji,jj+1,jl), pice(ji-1,jj+1,jl), pice(ji+1,jj+1,jl) )
pmax(ji,jj,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
zmax1(ji) = zmax2(ji)
zmax2(ji) = zmax3
END_2D
END DO
END SUBROUTINE icemax3D
SUBROUTINE icemax4D( pice , pmax )
!!---------------------------------------------------------------------
!! *** ROUTINE icemax4D ***
!! ** Purpose : compute the max of the 9 points around
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pice ! input
REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pmax ! output
!
REAL(wp), DIMENSION(Nis0:Nie0) :: zmax1, zmax2
REAL(wp) :: zmax3
INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices
!!----------------------------------------------------------------------
jlay = SIZE( pice , 3 ) ! size of input arrays
! basic version: get the max of epsi20 + 9 neighbours
!!$ DO jl = 1, jpl
!!$ DO jk = 1, jlay
!!$ DO_2D( 0, 0, 0, 0 )
!!$ pmax(ji,jj,jk,jl) = MAX( epsi20, pice(ji-1,jj-1,jk,jl), pice(ji,jj-1,jk,jl), pice(ji+1,jj-1,jk,jl), &
!!$ & pice(ji-1,jj ,jk,jl), pice(ji,jj ,jk,jl), pice(ji+1,jj ,jk,jl), &
!!$ & pice(ji-1,jj+1,jk,jl), pice(ji,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
!!$ END_2D
!!$ END DO
!!$ END DO
! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
DO jl = 1, jpl
DO jk = 1, jlay
DO ji = Nis0, Nie0
zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jk,jl), pice(ji-1,Njs0-1,jk,jl), pice(ji+1,Njs0-1,jk,jl) )
zmax2(ji) = MAX( epsi20, pice(ji,Njs0 ,jk,jl), pice(ji-1,Njs0 ,jk,jl), pice(ji+1,Njs0 ,jk,jl) )
END DO
DO_2D( 0, 0, 0, 0 )
zmax3 = MAX( epsi20, pice(ji,jj+1,jk,jl), pice(ji-1,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
pmax(ji,jj,jk,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
zmax1(ji) = zmax2(ji)
zmax2(ji) = zmax3
END_2D
END DO
END DO
END SUBROUTINE icemax4D
#else
!!----------------------------------------------------------------------
!! Default option Dummy module NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE icedyn_adv_pra