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

Last change on this file since 13637 was 13637, checked in by mathiot, 4 years ago

fix ticket #2553 in trunk
  • Property svn:keywords set to Id
File size: 70.0 KB
Line 
1MODULE icedyn_adv_pra 
2   !!======================================================================
3   !!                       ***  MODULE icedyn_adv_pra   ***
4   !!   sea-ice : advection => Prather scheme
5   !!======================================================================
6   !! History :       !  2008-03  (M. Vancoppenolle) original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!--------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_adv_pra : advection of sea ice using Prather scheme
14   !!   adv_x, adv_y    : Prather scheme applied in i- and j-direction, resp.
15   !!   adv_pra_init    : initialisation of the Prather scheme
16   !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constant
19   USE dom_oce        ! ocean domain
20   USE ice            ! sea-ice variables
21   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28   USE lbclnk         ! lateral boundary conditions (or mpp links)
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_dyn_adv_pra   ! called by icedyn_adv
34   PUBLIC   adv_pra_init      ! called by icedyn_adv
35
36   ! Moments for advection
37   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   ! ice thickness
38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! ice concentration
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume
47
48   !! * Substitutions
49#  include "do_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
59      !!----------------------------------------------------------------------
60      !!                **  routine ice_dyn_adv_pra  **
61      !! 
62      !! ** purpose :   Computes and adds the advection trend to sea-ice
63      !!
64      !! ** method  :   Uses Prather second order scheme that advects tracers
65      !!                but also their quadratic forms. The method preserves
66      !!                tracer structures by conserving second order moments.
67      !!
68      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
69      !!----------------------------------------------------------------------
70      INTEGER                     , INTENT(in   ) ::   kt         ! time step
71      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
72      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
73      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
74      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
75      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
76      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
79      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
80      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
82      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
83      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness
85      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
86      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
87      !
88      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
89      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
90      REAL(wp) ::   zdt, 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  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
992               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
993               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
994               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
995               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
996            END DO
997            !                                                        ! ice layers heat content
998            DO jk = 1, nlay_i
999               WRITE(zchar1,'(I2.2)') jk
1000               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
1001               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:)
1002               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
1003               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
1004               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
1005            END DO
1006            !
1007            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction
1008               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
1009                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
1010                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
1011                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
1012                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
1013                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
1014                  !                                                     ! melt pond volume
1015                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
1016                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
1017                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
1018                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
1019                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
1020               ELSE
1021                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction
1022                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume
1023               ENDIF
1024                  !
1025               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
1026                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
1027                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
1028                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
1029                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
1030                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
1031                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl )
1032                  ELSE
1033                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume
1034                  ENDIF
1035               ENDIF
1036            ENDIF
1037            !
1038         ELSE                                   !**  start rheology from rest  **!
1039            !
1040            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
1041            !
1042            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
1043            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
1044            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
1045            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
1046            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
1047            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
1048            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
1049            IF( ln_pnd_LEV ) THEN
1050               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
1051               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
1052               IF ( ln_pnd_lids ) THEN
1053                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
1054               ENDIF
1055            ENDIF
1056         ENDIF
1057         !
1058         !                                   !=====================================!
1059      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
1060         !                                   !=====================================!
1061         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
1062         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1063         !
1064         !
1065         ! In case Prather scheme is used for advection, write second order moments
1066         ! ------------------------------------------------------------------------
1067         !
1068         !                                                           ! ice thickness
1069         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
1070         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
1071         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
1072         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
1073         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
1074         !                                                           ! snow thickness
1075         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
1076         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
1077         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
1078         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
1079         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
1080         !                                                           ! ice concentration
1081         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
1082         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
1083         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
1084         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
1085         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
1086         !                                                           ! ice salinity
1087         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
1088         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
1089         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
1090         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
1091         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
1092         !                                                           ! ice age
1093         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
1094         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
1095         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
1096         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
1097         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
1098         !                                                           ! snow layers heat content
1099         DO jk = 1, nlay_s
1100            WRITE(zchar1,'(I2.2)') jk
1101            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1102            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1103            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1104            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1105            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1106         END DO
1107         !                                                           ! ice layers heat content
1108         DO jk = 1, nlay_i
1109            WRITE(zchar1,'(I2.2)') jk
1110            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1111            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1112            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1113            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1114            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1115         END DO
1116         !
1117         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction
1118            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
1119            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
1120            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
1121            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
1122            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
1123            !                                                        ! melt pond volume
1124            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
1125            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
1126            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
1127            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
1128            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
1129            !
1130            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
1131               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
1132               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
1133               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
1134               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
1135               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
1136            ENDIF
1137         ENDIF
1138         !
1139      ENDIF
1140      !
1141   END SUBROUTINE adv_pra_rst
1142
1143   SUBROUTINE icemax3D( pice , pmax )
1144      !!---------------------------------------------------------------------
1145      !!                   ***  ROUTINE icemax3D ***                     
1146      !! ** Purpose :  compute the max of the 9 points around
1147      !!----------------------------------------------------------------------
1148      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input
1149      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output
1150      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
1151      INTEGER  ::   ji, jj, jl   ! dummy loop indices
1152      !!----------------------------------------------------------------------
1153      DO jl = 1, jpl
1154         DO jj = Njs0-1, Nje0+1   
1155            DO ji = Nis0, Nie0
1156               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) )
1157            END DO
1158         END DO
1159         DO jj = Njs0, Nje0   
1160            DO ji = Nis0, Nie0
1161               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
1162            END DO
1163         END DO
1164      END DO
1165   END SUBROUTINE icemax3D
1166
1167   SUBROUTINE icemax4D( pice , pmax )
1168      !!---------------------------------------------------------------------
1169      !!                   ***  ROUTINE icemax4D ***                     
1170      !! ** Purpose :  compute the max of the 9 points around
1171      !!----------------------------------------------------------------------
1172      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input
1173      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output
1174      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
1175      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
1176      !!----------------------------------------------------------------------
1177      jlay = SIZE( pice , 3 )   ! size of input arrays
1178      DO jl = 1, jpl
1179         DO jk = 1, jlay
1180            DO jj = Njs0-1, Nje0+1   
1181               DO ji = Nis0, Nie0
1182                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )
1183               END DO
1184            END DO
1185            DO jj = Njs0, Nje0   
1186               DO ji = Nis0, Nie0
1187                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
1188               END DO
1189            END DO
1190         END DO
1191      END DO
1192   END SUBROUTINE icemax4D
1193
1194#else
1195   !!----------------------------------------------------------------------
1196   !!   Default option            Dummy module        NO SI3 sea-ice model
1197   !!----------------------------------------------------------------------
1198#endif
1199
1200   !!======================================================================
1201END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.