source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_pra.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 8 months ago

Add:

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