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

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

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90 @ 13745

Last change on this file since 13745 was 13630, checked in by mocavero, 4 years ago

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