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

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

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

Last change on this file since 13970 was 13970, checked in by andmirek, 3 years ago

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