New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_pra.F90 in NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/ICE/icedyn_adv_pra.F90 @ 11870

Last change on this file since 11870 was 11870, checked in by andmirek, 4 years ago

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