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/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90 @ 12744

Last change on this file since 12744 was 12744, checked in by clem, 4 years ago

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