source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_pra.F90 @ 12379

Last change on this file since 12379 was 12379, checked in by dancopsey, 14 months ago

Add meltpond lid thickness as a new prognostic.

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