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

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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_pra.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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