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/2020/dev_12905_xios_restart/src/ICE – NEMO

source: NEMO/branches/2020/dev_12905_xios_restart/src/ICE/icedyn_adv_pra.F90 @ 13022

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

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