New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_pra.F90 in NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icedyn_adv_pra.F90 @ 12197

Last change on this file since 12197 was 12197, checked in by clem, 4 years ago

correct ticket #2349
  • Property svn:keywords set to Id
File size: 59.7 KB
Line 
1MODULE icedyn_adv_pra 
2   !!======================================================================
3   !!                       ***  MODULE icedyn_adv_pra   ***
4   !!   sea-ice : advection => Prather scheme
5   !!======================================================================
6   !! History :       !  2008-03  (M. Vancoppenolle) original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!--------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_adv_pra : advection of sea ice using Prather scheme
14   !!   adv_x, adv_y    : Prather scheme applied in i- and j-direction, resp.
15   !!   adv_pra_init    : initialisation of the Prather scheme
16   !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constant
19   USE dom_oce        ! ocean domain
20   USE ice            ! sea-ice variables
21   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28   USE lbclnk         ! lateral boundary conditions (or mpp links)
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_dyn_adv_pra   ! called by icedyn_adv
34   PUBLIC   adv_pra_init      ! called by icedyn_adv
35
36   ! Moments for advection
37   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   ! ice thickness
38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! ice concentration
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
57      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
58      !!----------------------------------------------------------------------
59      !!                **  routine ice_dyn_adv_pra  **
60      !! 
61      !! ** purpose :   Computes and adds the advection trend to sea-ice
62      !!
63      !! ** method  :   Uses Prather second order scheme that advects tracers
64      !!                but also their quadratic forms. The method preserves
65      !!                tracer structures by conserving second order moments.
66      !!
67      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
68      !!----------------------------------------------------------------------
69      INTEGER                     , INTENT(in   ) ::   kt         ! time step
70      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
71      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
72      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
73      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
74      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
75      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
76      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
79      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
80      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
82      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
83      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
84      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
85      !
86      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices
87      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
88      REAL(wp) ::   zdt                     !   -      -
89      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
90      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2
91      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx
92      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max
93      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
95      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
97      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
98      !!----------------------------------------------------------------------
99      !
100      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
101      !
102      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
103      DO jl = 1, jpl
104         DO jj = 2, jpjm1
105            DO ji = fs_2, fs_jpim1
106               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
107                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
108                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
109                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
110               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
111                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
112                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
113                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
114               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
115                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
116                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
117                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
118            END DO
119         END DO
120      END DO
121      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
122      !
123      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
124      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
125      !              this should not affect too much the stability
126      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
127      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
128     
129      ! non-blocking global communication send zcflnow and receive zcflprv
130      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
131
132      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
133      ELSE                         ;   icycle = 1
134      ENDIF
135      zdt = rdt_ice / REAL(icycle)
136     
137      ! --- transport --- !
138      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
139      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
140
141      DO jt = 1, icycle
142
143         ! record at_i before advection (for open water)
144         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
145         
146         ! --- transported fields --- !                                       
147         DO jl = 1, jpl
148            zarea(:,:,jl) = e1e2t(:,:)
149            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
150            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
151            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
152            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
153            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
154            DO jk = 1, nlay_s
155               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
156            END DO
157            DO jk = 1, nlay_i
158               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
159            END DO
160            IF ( ln_pnd_H12 ) THEN
161               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
162               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
163            ENDIF
164         END DO
165         !
166         !                                                                  !--------------------------------------------!
167         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
168            !                                                               !--------------------------------------------!
169            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
170            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
171            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
172            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
173            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
174            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
175            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
176            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
177            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
178            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
179            !
180            DO jk = 1, nlay_s                                                                           !--- snow heat content
181               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
182                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
183               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
184                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
185            END DO
186            DO jk = 1, nlay_i                                                                           !--- ice heat content
187               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
188                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
189               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
190                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
191            END DO
192            !
193            IF ( ln_pnd_H12 ) THEN
194               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
195               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
196               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
197               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
198            ENDIF
199            !                                                               !--------------------------------------------!
200         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==!
201            !                                                               !--------------------------------------------!
202            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
203            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
204            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
205            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
206            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
207            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
208            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
209            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
210            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
211            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
212            DO jk = 1, nlay_s                                                                           !--- snow heat content
213               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
214                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
215               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
216                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
217            END DO
218            DO jk = 1, nlay_i                                                                           !--- ice heat content
219               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
220                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
221               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
222                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
223            END DO
224            IF ( ln_pnd_H12 ) THEN
225               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
226               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
227               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
228               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
229            ENDIF
230            !
231         ENDIF
232
233         ! --- Recover the properties from their contents --- !
234         DO jl = 1, jpl
235            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
236            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
237            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
238            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
239            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
240            DO jk = 1, nlay_s
241               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
242            END DO
243            DO jk = 1, nlay_i
244               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
245            END DO
246            IF ( ln_pnd_H12 ) THEN
247               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
248               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
249            ENDIF
250         END DO
251         !
252         ! derive open water from ice concentration
253         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
254         DO jj = 2, jpjm1
255            DO ji = fs_2, fs_jpim1
256               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
257                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
258            END DO
259         END DO
260         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. )
261         !
262         ! --- Ensure non-negative fields --- !
263         !     Remove negative values (conservation is ensured)
264         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
265         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
266         !
267         ! --- Make sure ice thickness is not too big --- !
268         !     (because ice thickness can be too large where ice concentration is very small)
269         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
270         !
271         ! --- Ensure snow load is not too big --- !
272         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
273         !
274      END DO
275      !
276      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
277      !
278   END SUBROUTINE ice_dyn_adv_pra
279   
280   
281   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
282      &              psx, psxx, psy , psyy, psxy )
283      !!----------------------------------------------------------------------
284      !!                **  routine adv_x  **
285      !! 
286      !! ** purpose :   Computes and adds the advection trend to sea-ice
287      !!                variable on x axis
288      !!----------------------------------------------------------------------
289      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
290      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
291      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
292      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
293      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
294      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
295      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
296      !!
297      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
298      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
299      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
300      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
301      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
302      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
303      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
304      !-----------------------------------------------------------------------
305      !
306      jcat = SIZE( ps0 , 3 )   ! size of input arrays
307      !
308      DO jl = 1, jcat   ! loop on categories
309         !
310         ! Limitation of moments.                                           
311         DO jj = 2, jpjm1
312            DO ji = 1, jpi
313               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
314               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
315               !
316               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
317               zs1max  = 1.5 * zslpmax
318               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
319               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
320                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
321               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
322
323               ps0 (ji,jj,jl) = zslpmax 
324               psx (ji,jj,jl) = zs1new         * rswitch
325               psxx(ji,jj,jl) = zs2new         * rswitch
326               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
327               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
328               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
329            END DO
330         END DO
331
332         !  Calculate fluxes and moments between boxes i<-->i+1             
333         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0
334            DO ji = 1, jpi
335               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
336               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
337               zalfq        =  zalf * zalf
338               zalf1        =  1.0 - zalf
339               zalf1q       =  zalf1 * zalf1
340               !
341               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
342               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
343               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
344               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
345               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
346               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
347               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
348
349               !  Readjust moments remaining in the box.
350               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
351               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
352               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
353               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
354               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
355               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
356               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
357            END DO
358         END DO
359
360         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
361            DO ji = 1, fs_jpim1
362               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
363               zalg  (ji,jj) = zalf
364               zalfq         = zalf * zalf
365               zalf1         = 1.0 - zalf
366               zalg1 (ji,jj) = zalf1
367               zalf1q        = zalf1 * zalf1
368               zalg1q(ji,jj) = zalf1q
369               !
370               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
371               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
372                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
373               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
374               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
375               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
376               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
377               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
378            END DO
379         END DO
380
381         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
382            DO ji = fs_2, fs_jpim1
383               zbt  =       zbet(ji-1,jj)
384               zbt1 = 1.0 - zbet(ji-1,jj)
385               !
386               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
387               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
388               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
389               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
390               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
391               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
392               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
393            END DO
394         END DO
395
396         !   Put the temporary moments into appropriate neighboring boxes.   
397         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
398            DO ji = fs_2, fs_jpim1
399               zbt  =       zbet(ji-1,jj)
400               zbt1 = 1.0 - zbet(ji-1,jj)
401               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
402               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
403               zalf1         = 1.0 - zalf
404               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
405               !
406               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
407               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
408               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
409                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
410                  &            + zbt1 * psxx(ji,jj,jl)
411               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
412                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
413                  &            + zbt1 * psxy(ji,jj,jl)
414               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
415               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
416            END DO
417         END DO
418
419         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0.
420            DO ji = fs_2, fs_jpim1
421               zbt  =       zbet(ji,jj)
422               zbt1 = 1.0 - zbet(ji,jj)
423               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
424               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
425               zalf1         = 1.0 - zalf
426               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
427               !
428               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
429               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
430               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
431                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
432                  &                                           + ( zalf1 - zalf ) * ztemp ) )
433               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
434                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
435               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
436               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
437            END DO
438         END DO
439
440      END DO
441
442      !-- Lateral boundary conditions
443      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
444         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
445         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
446      !
447   END SUBROUTINE adv_x
448
449
450   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
451      &              psx, psxx, psy , psyy, psxy )
452      !!---------------------------------------------------------------------
453      !!                **  routine adv_y  **
454      !!           
455      !! ** purpose :   Computes and adds the advection trend to sea-ice
456      !!                variable on y axis
457      !!---------------------------------------------------------------------
458      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
459      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
460      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
461      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
462      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
463      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
464      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
465      !!
466      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
467      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
468      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
469      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
470      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
471      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
472      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
473      !---------------------------------------------------------------------
474      !
475      jcat = SIZE( ps0 , 3 )   ! size of input arrays
476      !     
477      DO jl = 1, jcat   ! loop on categories
478         !
479         ! Limitation of moments.
480         DO jj = 1, jpj
481            DO ji = fs_2, fs_jpim1
482               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
483               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
484               !
485               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
486               zs1max  = 1.5 * zslpmax
487               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
488               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
489                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
490               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
491               !
492               ps0 (ji,jj,jl) = zslpmax 
493               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
494               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
495               psy (ji,jj,jl) = zs1new         * rswitch
496               psyy(ji,jj,jl) = zs2new         * rswitch
497               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
498            END DO
499         END DO
500 
501         !  Calculate fluxes and moments between boxes j<-->j+1             
502         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
503            DO ji = fs_2, fs_jpim1
504               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
505               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
506               zalfq        =  zalf * zalf
507               zalf1        =  1.0 - zalf
508               zalf1q       =  zalf1 * zalf1
509               !
510               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
511               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
512               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
513               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
514               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
515               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
516               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
517               !
518               !  Readjust moments remaining in the box.
519               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
520               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
521               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
522               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
523               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
524               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
525               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
526            END DO
527         END DO
528         !
529         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
530            DO ji = fs_2, fs_jpim1
531               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
532               zalg  (ji,jj) = zalf
533               zalfq         = zalf * zalf
534               zalf1         = 1.0 - zalf
535               zalg1 (ji,jj) = zalf1
536               zalf1q        = zalf1 * zalf1
537               zalg1q(ji,jj) = zalf1q
538               !
539               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl)
540               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
541                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
542               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
543               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
544               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
545               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
546               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
547            END DO
548         END DO
549
550         !  Readjust moments remaining in the box.
551         DO jj = 2, jpjm1
552            DO ji = fs_2, fs_jpim1
553               zbt  =         zbet(ji,jj-1)
554               zbt1 = ( 1.0 - zbet(ji,jj-1) )
555               !
556               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
557               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
558               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
559               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
560               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
561               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
562               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
563            END DO
564         END DO
565
566         !   Put the temporary moments into appropriate neighboring boxes.   
567         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
568            DO ji = fs_2, fs_jpim1
569               zbt  =       zbet(ji,jj-1)
570               zbt1 = 1.0 - zbet(ji,jj-1)
571               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
572               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
573               zalf1         = 1.0 - zalf
574               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
575               !
576               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
577               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
578                  &             + zbt1 * psy(ji,jj,jl) 
579               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
580                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
581                  &             + zbt1 * psyy(ji,jj,jl)
582               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
583                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
584                  &             + zbt1 * psxy(ji,jj,jl)
585               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
586               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
587            END DO
588         END DO
589
590         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
591            DO ji = fs_2, fs_jpim1
592               zbt  =       zbet(ji,jj)
593               zbt1 = 1.0 - zbet(ji,jj)
594               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
595               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
596               zalf1         = 1.0 - zalf
597               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
598               !
599               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
600               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
601               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
602                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
603                  &                                            + ( zalf1 - zalf ) * ztemp ) )
604               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
605                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
606               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
607               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
608            END DO
609         END DO
610
611      END DO
612
613      !-- Lateral boundary conditions
614      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
615         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
616         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
617      !
618   END SUBROUTINE adv_y
619
620
621   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s )
622      !!-------------------------------------------------------------------
623      !!                  ***  ROUTINE Hbig  ***
624      !!
625      !! ** Purpose : Thickness correction in case advection scheme creates
626      !!              abnormally tick ice or snow
627      !!
628      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
629      !!                 (before advection) and reduce it by adapting ice concentration
630      !!              2- check whether snow thickness is larger than the surrounding 9-points
631      !!                 (before advection) and reduce it by sending the excess in the ocean
632      !!
633      !! ** input   : Max thickness of the surrounding 9-points
634      !!-------------------------------------------------------------------
635      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step
636      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
637      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip
638      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
639      !
640      INTEGER  ::   ji, jj, jl         ! dummy loop indices
641      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra
642      !!-------------------------------------------------------------------
643      !
644      z1_dt = 1._wp / pdt
645      !
646      DO jl = 1, jpl
647
648         DO jj = 1, jpj
649            DO ji = 1, jpi
650               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
651                  !
652                  !                               ! -- check h_ip -- !
653                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
654                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
655                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
656                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
657                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
658                     ENDIF
659                  ENDIF
660                  !
661                  !                               ! -- check h_i -- !
662                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
663                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
664                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
665                     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)
666                  ENDIF
667                  !
668                  !                               ! -- check h_s -- !
669                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
670                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
671                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
672                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
673                     !
674                     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
675                     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
676                     !
677                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
678                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
679                  ENDIF           
680                  !                 
681               ENDIF
682            END DO
683         END DO
684      END DO 
685      !
686   END SUBROUTINE Hbig
687
688
689   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
690      !!-------------------------------------------------------------------
691      !!                  ***  ROUTINE Hsnow  ***
692      !!
693      !! ** Purpose : 1- Check snow load after advection
694      !!              2- Correct pond concentration to avoid a_ip > a_i
695      !!
696      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
697      !!              then put the snow excess in the ocean
698      !!
699      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards
700      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
701      !!              make the snow very thick (if concentration decreases drastically)
702      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
703      !!-------------------------------------------------------------------
704      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
705      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
706      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
707      !
708      INTEGER  ::   ji, jj, jl   ! dummy loop indices
709      REAL(wp) ::   z1_dt, zvs_excess, zfra
710      !!-------------------------------------------------------------------
711      !
712      z1_dt = 1._wp / pdt
713      !
714      ! -- check snow load -- !
715      DO jl = 1, jpl
716         DO jj = 1, jpj
717            DO ji = 1, jpi
718               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
719                  !
720                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
721                  !
722                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
723                     ! put snow excess in the ocean
724                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
725                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
726                     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
727                     ! correct snow volume and heat content
728                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
729                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
730                  ENDIF
731                  !
732               ENDIF
733            END DO
734         END DO
735      END DO
736      !
737      !-- correct pond concentration to avoid a_ip > a_i -- !
738      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
739      !
740   END SUBROUTINE Hsnow
741
742
743   SUBROUTINE adv_pra_init
744      !!-------------------------------------------------------------------
745      !!                  ***  ROUTINE adv_pra_init  ***
746      !!
747      !! ** Purpose :   allocate and initialize arrays for Prather advection
748      !!-------------------------------------------------------------------
749      INTEGER ::   ierr
750      !!-------------------------------------------------------------------
751      !
752      !                             !* allocate prather fields
753      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
754         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
755         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
756         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
757         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
758         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
759         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
760         !
761         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
762         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
763         !
764         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
765         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
766         &      STAT = ierr )
767      !
768      CALL mpp_sum( 'icedyn_adv_pra', ierr )
769      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
770      !
771      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
772      !
773   END SUBROUTINE adv_pra_init
774
775
776   SUBROUTINE adv_pra_rst( cdrw, kt )
777      !!---------------------------------------------------------------------
778      !!                   ***  ROUTINE adv_pra_rst  ***
779      !!                     
780      !! ** Purpose :   Read or write file in restart file
781      !!
782      !! ** Method  :   use of IOM library
783      !!----------------------------------------------------------------------
784      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
785      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
786      !
787      INTEGER ::   jk, jl   ! dummy loop indices
788      INTEGER ::   iter     ! local integer
789      INTEGER ::   id1      ! local integer
790      CHARACTER(len=25) ::   znam
791      CHARACTER(len=2)  ::   zchar, zchar1
792      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
793      !!----------------------------------------------------------------------
794      !
795      !                                      !==========================!
796      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
797         !                                   !==========================!
798         !
799         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
800         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
801         ENDIF
802         !
803         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
804            !
805            !                                                        ! ice thickness
806            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
807            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
808            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
809            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
810            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
811            !                                                        ! snow thickness
812            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
813            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
814            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
815            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
816            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
817            !                                                        ! ice concentration
818            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
819            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
820            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
821            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
822            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
823            !                                                        ! ice salinity
824            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
825            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
826            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
827            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
828            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
829            !                                                        ! ice age
830            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
831            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
832            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
833            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
834            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
835            !                                                        ! snow layers heat content
836            DO jk = 1, nlay_s
837               WRITE(zchar1,'(I2.2)') jk
838               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
839               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
840               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
841               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
842               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
843            END DO
844            !                                                        ! ice layers heat content
845            DO jk = 1, nlay_i
846               WRITE(zchar1,'(I2.2)') jk
847               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
848               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
849               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
850               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
851               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
852            END DO
853            !
854            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
855               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
856               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
857               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
858               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
859               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
860               !                                                     ! melt pond volume
861               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
862               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
863               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
864               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
865               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
866            ENDIF
867            !
868         ELSE                                   !**  start rheology from rest  **!
869            !
870            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
871            !
872            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
873            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
874            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
875            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
876            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
877            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
878            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
879            IF( ln_pnd_H12 ) THEN
880               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
881               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
882            ENDIF
883         ENDIF
884         !
885         !                                   !=====================================!
886      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
887         !                                   !=====================================!
888         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
889         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
890         !
891         !
892         ! In case Prather scheme is used for advection, write second order moments
893         ! ------------------------------------------------------------------------
894         !
895         !                                                           ! ice thickness
896         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
897         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
898         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
899         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
900         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
901         !                                                           ! snow thickness
902         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
903         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
904         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
905         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
906         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
907         !                                                           ! ice concentration
908         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
909         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
910         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
911         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
912         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
913         !                                                           ! ice salinity
914         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
915         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
916         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
917         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
918         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
919         !                                                           ! ice age
920         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
921         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
922         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
923         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
924         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
925         !                                                           ! snow layers heat content
926         DO jk = 1, nlay_s
927            WRITE(zchar1,'(I2.2)') jk
928            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
929            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
930            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
931            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
932            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
933         END DO
934         !                                                           ! ice layers heat content
935         DO jk = 1, nlay_i
936            WRITE(zchar1,'(I2.2)') jk
937            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
938            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
939            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
940            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
941            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
942         END DO
943         !
944         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
945            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
946            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
947            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
948            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
949            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
950            !                                                        ! melt pond volume
951            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
952            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
953            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
954            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
955            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
956         ENDIF
957         !
958      ENDIF
959      !
960   END SUBROUTINE adv_pra_rst
961
962#else
963   !!----------------------------------------------------------------------
964   !!   Default option            Dummy module        NO SI3 sea-ice model
965   !!----------------------------------------------------------------------
966#endif
967
968   !!======================================================================
969END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.