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

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90 @ 13589

Last change on this file since 13589 was 13589, checked in by clem, 3 years ago

4.0-HEAD: rewrite heat budget of sea ice to make it perfectly conservative by construction. Also, activating ln_icediachk now gives an ascii file (icedrift_diagnostics.ascii) containing mass, salt and heat global conservation issues (if any). In addition, 2D drift files can be outputed (icedrift_mass...) but since I did not want to change the xml files, one has to uncomment some lines in icectl.F90 and add the corresponding fields in field_def_oce.xml

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