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 branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90 @ 8752

Last change on this file since 8752 was 8752, checked in by dancopsey, 6 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

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