source: NEMO/trunk/src/ICE/icedyn_adv_pra.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 7 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 57.5 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            !
773            !                                                        ! ice thickness
774            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
775            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
776            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
777            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
778            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
779            !                                                        ! snow thickness
780            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
781            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
782            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
783            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
784            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
785            !                                                        ! ice concentration
786            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
787            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
788            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
789            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
790            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
791            !                                                        ! ice salinity
792            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
793            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
794            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
795            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
796            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
797            !                                                        ! ice age
798            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
799            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
800            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
801            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
802            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
803            !                                                        ! snow layers heat content
804            DO jk = 1, nlay_s
805               WRITE(zchar1,'(I2.2)') jk
806               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
807               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
808               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
809               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
810               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
811            END DO
812            !                                                        ! ice layers heat content
813            DO jk = 1, nlay_i
814               WRITE(zchar1,'(I2.2)') jk
815               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
816               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
817               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
818               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
819               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
820            END DO
821            !
822            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
823               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
824               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
825               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
826               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
827               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
828               !                                                     ! melt pond volume
829               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
830               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
831               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
832               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
833               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
834            ENDIF
835            !
836         ELSE                                   !**  start rheology from rest  **!
837            !
838            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
839            !
840            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
841            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
842            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
843            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
844            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
845            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
846            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
847            IF( ln_pnd_H12 ) THEN
848               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
849               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
850            ENDIF
851         ENDIF
852         !
853         !                                   !=====================================!
854      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
855         !                                   !=====================================!
856         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
857         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
858         !
859         !
860         ! In case Prather scheme is used for advection, write second order moments
861         ! ------------------------------------------------------------------------
862         !
863         !                                                           ! ice thickness
864         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
865         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
866         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
867         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
868         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
869         !                                                           ! snow thickness
870         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
871         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
872         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
873         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
874         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
875         !                                                           ! ice concentration
876         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
877         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
878         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
879         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
880         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
881         !                                                           ! ice salinity
882         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
883         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
884         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
885         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
886         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
887         !                                                           ! ice age
888         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
889         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
890         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
891         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
892         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
893         !                                                           ! snow layers heat content
894         DO jk = 1, nlay_s
895            WRITE(zchar1,'(I2.2)') jk
896            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
897            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
898            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
899            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
900            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
901         END DO
902         !                                                           ! ice layers heat content
903         DO jk = 1, nlay_i
904            WRITE(zchar1,'(I2.2)') jk
905            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
906            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
907            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
908            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
909            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
910         END DO
911         !
912         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
913            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
914            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
915            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
916            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
917            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
918            !                                                        ! melt pond volume
919            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
920            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
921            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
922            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
923            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
924         ENDIF
925         !
926      ENDIF
927      !
928   END SUBROUTINE adv_pra_rst
929
930#else
931   !!----------------------------------------------------------------------
932   !!   Default option            Dummy module        NO SI3 sea-ice model
933   !!----------------------------------------------------------------------
934#endif
935
936   !!======================================================================
937END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.