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

Last change on this file since 13613 was 13613, checked in by smasson, 4 years ago

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