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 @ 12720

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

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