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/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_adv_pra.F90 @ 12937

Last change on this file since 12937 was 12937, checked in by dancopsey, 3 years ago

Merge in Clem's branch. It was originally here:

svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch

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