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

Last change on this file since 13564 was 13564, checked in by clem, 4 years ago

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