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

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

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/ICE/icedyn_adv_pra.F90 @ 14345

Last change on this file since 14345 was 13945, checked in by mocavero, 4 years ago

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