source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90 @ 8554

Last change on this file since 8554 was 8534, checked in by clem, 3 years ago

changes in style - part6 - pure cosmetics

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