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

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

Last change on this file since 13497 was 13497, checked in by techene, 4 years ago

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