source: NEMO/trunk/src/ICE/icedyn_adv_pra.F90 @ 11812

Last change on this file since 11812 was 11812, checked in by smasson, 12 months ago

trunk: add missing USE phycst in icedyn_adv_pra, needed when not using key_iomput

  • Property svn:keywords set to Id
File size: 53.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
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  &
57      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
58      !!----------------------------------------------------------------------
59      !!                **  routine ice_dyn_adv_pra  **
60      !! 
61      !! ** purpose :   Computes and adds the advection trend to sea-ice
62      !!
63      !! ** method  :   Uses Prather second order scheme that advects tracers
64      !!                but also their quadratic forms. The method preserves
65      !!                tracer structures by conserving second order moments.
66      !!
67      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
68      !!----------------------------------------------------------------------
69      INTEGER                     , INTENT(in   ) ::   kt         ! time step
70      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
71      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
72      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
73      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
74      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
75      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
76      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
79      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
80      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
81      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
82      !
83      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices
84      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
85      REAL(wp) ::   zdt                     !   -      -
86      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
87      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2
88      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx
89      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
90      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
91      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
92      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
93      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
94      !!----------------------------------------------------------------------
95      !
96      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
97      !
98      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
99      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
100      !              this should not affect too much the stability
101      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
102      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
103     
104      ! non-blocking global communication send zcflnow and receive zcflprv
105      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
106
107      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
108      ELSE                         ;   icycle = 1
109      ENDIF
110      zdt = rdt_ice / REAL(icycle)
111     
112      ! --- transport --- !
113      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
114      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
115
116      DO jt = 1, icycle
117
118         ! record at_i before advection (for open water)
119         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
120         
121         ! --- transported fields --- !                                       
122         DO jl = 1, jpl
123            zarea(:,:,jl) = e1e2t(:,:)
124            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
125            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
126            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
127            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
128            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
129            DO jk = 1, nlay_s
130               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
131            END DO
132            DO jk = 1, nlay_i
133               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
134            END DO
135            IF ( ln_pnd_H12 ) THEN
136               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
137               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
138            ENDIF
139         END DO
140         !
141         !                                                                  !--------------------------------------------!
142         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
143            !                                                               !--------------------------------------------!
144            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
145            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
146            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
147            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
148            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
149            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
150            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
151            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
152            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
153            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
154            !
155            DO jk = 1, nlay_s                                                                           !--- snow heat content
156               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
157                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
158               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
159                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
160            END DO
161            DO jk = 1, nlay_i                                                                           !--- ice heat content
162               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
163                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
164               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
165                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
166            END DO
167            !
168            IF ( ln_pnd_H12 ) THEN
169               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
170               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
171               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
172               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
173            ENDIF
174            !                                                               !--------------------------------------------!
175         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==!
176            !                                                               !--------------------------------------------!
177            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
178            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
179            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
180            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
181            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
182            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
183            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
184            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
185            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
186            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
187            DO jk = 1, nlay_s                                                                           !--- snow heat content
188               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
189                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
190               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
191                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
192            END DO
193            DO jk = 1, nlay_i                                                                           !--- ice heat content
194               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
195                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
196               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
197                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
198            END DO
199            IF ( ln_pnd_H12 ) THEN
200               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
201               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
202               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
203               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
204            ENDIF
205            !
206         ENDIF
207
208         ! --- Recover the properties from their contents --- !
209         DO jl = 1, jpl
210            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
211            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
212            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
213            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
214            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
215            DO jk = 1, nlay_s
216               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
217            END DO
218            DO jk = 1, nlay_i
219               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
220            END DO
221            IF ( ln_pnd_H12 ) THEN
222               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
223               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
224            ENDIF
225         END DO
226         !
227         ! derive open water from ice concentration
228         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
229         DO jj = 2, jpjm1
230            DO ji = fs_2, fs_jpim1
231               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
232                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
233            END DO
234         END DO
235         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. )
236         !
237         ! --- Ensure non-negative fields --- !
238         !     Remove negative values (conservation is ensured)
239         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
240         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
241         !
242         ! --- Ensure snow load is not too big --- !
243         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
244         !
245      END DO
246      !
247      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
248      !
249   END SUBROUTINE ice_dyn_adv_pra
250   
251   
252   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
253      &              psx, psxx, psy , psyy, psxy )
254      !!----------------------------------------------------------------------
255      !!                **  routine adv_x  **
256      !! 
257      !! ** purpose :   Computes and adds the advection trend to sea-ice
258      !!                variable on x axis
259      !!----------------------------------------------------------------------
260      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
261      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
262      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
263      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
264      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
265      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
266      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
267      !!
268      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
269      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
270      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
271      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
272      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
273      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
274      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
275      !-----------------------------------------------------------------------
276      !
277      jcat = SIZE( ps0 , 3 )   ! size of input arrays
278      !
279      DO jl = 1, jcat   ! loop on categories
280         !
281         ! Limitation of moments.                                           
282         DO jj = 2, jpjm1
283            DO ji = 1, jpi
284               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
285               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
286               !
287               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
288               zs1max  = 1.5 * zslpmax
289               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
290               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
291                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
292               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
293
294               ps0 (ji,jj,jl) = zslpmax 
295               psx (ji,jj,jl) = zs1new         * rswitch
296               psxx(ji,jj,jl) = zs2new         * rswitch
297               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
298               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
299               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
300            END DO
301         END DO
302
303         !  Calculate fluxes and moments between boxes i<-->i+1             
304         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0
305            DO ji = 1, jpi
306               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
307               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
308               zalfq        =  zalf * zalf
309               zalf1        =  1.0 - zalf
310               zalf1q       =  zalf1 * zalf1
311               !
312               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
313               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
314               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
315               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
316               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
317               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
318               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
319
320               !  Readjust moments remaining in the box.
321               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
322               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
323               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
324               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
325               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
326               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
327               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
328            END DO
329         END DO
330
331         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
332            DO ji = 1, fs_jpim1
333               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
334               zalg  (ji,jj) = zalf
335               zalfq         = zalf * zalf
336               zalf1         = 1.0 - zalf
337               zalg1 (ji,jj) = zalf1
338               zalf1q        = zalf1 * zalf1
339               zalg1q(ji,jj) = zalf1q
340               !
341               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
342               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
343                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
344               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
345               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
346               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
347               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
348               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
349            END DO
350         END DO
351
352         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
353            DO ji = fs_2, fs_jpim1
354               zbt  =       zbet(ji-1,jj)
355               zbt1 = 1.0 - zbet(ji-1,jj)
356               !
357               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
358               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
359               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
360               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
361               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
362               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
363               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
364            END DO
365         END DO
366
367         !   Put the temporary moments into appropriate neighboring boxes.   
368         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
369            DO ji = fs_2, fs_jpim1
370               zbt  =       zbet(ji-1,jj)
371               zbt1 = 1.0 - zbet(ji-1,jj)
372               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
373               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
374               zalf1         = 1.0 - zalf
375               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
376               !
377               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
378               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
379               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
380                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
381                  &            + zbt1 * psxx(ji,jj,jl)
382               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
383                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
384                  &            + zbt1 * psxy(ji,jj,jl)
385               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
386               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
387            END DO
388         END DO
389
390         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0.
391            DO ji = fs_2, fs_jpim1
392               zbt  =       zbet(ji,jj)
393               zbt1 = 1.0 - zbet(ji,jj)
394               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
395               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
396               zalf1         = 1.0 - zalf
397               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
398               !
399               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
400               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
401               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
402                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
403                  &                                           + ( zalf1 - zalf ) * ztemp ) )
404               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
405                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
406               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
407               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
408            END DO
409         END DO
410
411      END DO
412
413      !-- Lateral boundary conditions
414      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
415         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
416         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
417      !
418   END SUBROUTINE adv_x
419
420
421   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
422      &              psx, psxx, psy , psyy, psxy )
423      !!---------------------------------------------------------------------
424      !!                **  routine adv_y  **
425      !!           
426      !! ** purpose :   Computes and adds the advection trend to sea-ice
427      !!                variable on y axis
428      !!---------------------------------------------------------------------
429      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
430      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
431      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
432      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
433      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
434      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
435      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
436      !!
437      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
438      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
439      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
440      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
441      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
442      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
443      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
444      !---------------------------------------------------------------------
445      !
446      jcat = SIZE( ps0 , 3 )   ! size of input arrays
447      !     
448      DO jl = 1, jcat   ! loop on categories
449         !
450         ! Limitation of moments.
451         DO jj = 1, jpj
452            DO ji = fs_2, fs_jpim1
453               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
454               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
455               !
456               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
457               zs1max  = 1.5 * zslpmax
458               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
459               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
460                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
461               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
462               !
463               ps0 (ji,jj,jl) = zslpmax 
464               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
465               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
466               psy (ji,jj,jl) = zs1new         * rswitch
467               psyy(ji,jj,jl) = zs2new         * rswitch
468               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
469            END DO
470         END DO
471 
472         !  Calculate fluxes and moments between boxes j<-->j+1             
473         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
474            DO ji = fs_2, fs_jpim1
475               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
476               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
477               zalfq        =  zalf * zalf
478               zalf1        =  1.0 - zalf
479               zalf1q       =  zalf1 * zalf1
480               !
481               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
482               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
483               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
484               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
485               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
486               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
487               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
488               !
489               !  Readjust moments remaining in the box.
490               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
491               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
492               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
493               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
494               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
495               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
496               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
497            END DO
498         END DO
499         !
500         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
501            DO ji = fs_2, fs_jpim1
502               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
503               zalg  (ji,jj) = zalf
504               zalfq         = zalf * zalf
505               zalf1         = 1.0 - zalf
506               zalg1 (ji,jj) = zalf1
507               zalf1q        = zalf1 * zalf1
508               zalg1q(ji,jj) = zalf1q
509               !
510               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl)
511               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
512                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
513               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
514               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
515               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
516               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
517               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
518            END DO
519         END DO
520
521         !  Readjust moments remaining in the box.
522         DO jj = 2, jpjm1
523            DO ji = fs_2, fs_jpim1
524               zbt  =         zbet(ji,jj-1)
525               zbt1 = ( 1.0 - zbet(ji,jj-1) )
526               !
527               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
528               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
529               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
530               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
531               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
532               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
533               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
534            END DO
535         END DO
536
537         !   Put the temporary moments into appropriate neighboring boxes.   
538         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
539            DO ji = fs_2, fs_jpim1
540               zbt  =       zbet(ji,jj-1)
541               zbt1 = 1.0 - zbet(ji,jj-1)
542               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
543               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
544               zalf1         = 1.0 - zalf
545               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
546               !
547               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
548               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
549                  &             + zbt1 * psy(ji,jj,jl) 
550               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
551                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
552                  &             + zbt1 * psyy(ji,jj,jl)
553               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
554                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
555                  &             + zbt1 * psxy(ji,jj,jl)
556               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
557               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
558            END DO
559         END DO
560
561         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
562            DO ji = fs_2, fs_jpim1
563               zbt  =       zbet(ji,jj)
564               zbt1 = 1.0 - zbet(ji,jj)
565               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
566               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
567               zalf1         = 1.0 - zalf
568               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
569               !
570               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
571               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
572               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
573                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
574                  &                                            + ( zalf1 - zalf ) * ztemp ) )
575               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
576                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
577               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
578               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
579            END DO
580         END DO
581
582      END DO
583
584      !-- Lateral boundary conditions
585      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
586         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
587         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
588      !
589   END SUBROUTINE adv_y
590
591
592   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
593      !!-------------------------------------------------------------------
594      !!                  ***  ROUTINE Hsnow  ***
595      !!
596      !! ** Purpose : 1- Check snow load after advection
597      !!              2- Correct pond concentration to avoid a_ip > a_i
598      !!
599      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
600      !!              then put the snow excess in the ocean
601      !!
602      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards
603      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
604      !!              make the snow very thick (if concentration decreases drastically)
605      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
606      !!-------------------------------------------------------------------
607      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
608      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
609      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
610      !
611      INTEGER  ::   ji, jj, jl   ! dummy loop indices
612      REAL(wp) ::   z1_dt, zvs_excess, zfra
613      !!-------------------------------------------------------------------
614      !
615      z1_dt = 1._wp / pdt
616      !
617      ! -- check snow load -- !
618      DO jl = 1, jpl
619         DO jj = 1, jpj
620            DO ji = 1, jpi
621               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
622                  !
623                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
624                  !
625                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
626                     ! put snow excess in the ocean
627                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
628                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
629                     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
630                     ! correct snow volume and heat content
631                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
632                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
633                  ENDIF
634                  !
635               ENDIF
636            END DO
637         END DO
638      END DO
639      !
640      !-- correct pond concentration to avoid a_ip > a_i -- !
641      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
642      !
643   END SUBROUTINE Hsnow
644
645
646   SUBROUTINE adv_pra_init
647      !!-------------------------------------------------------------------
648      !!                  ***  ROUTINE adv_pra_init  ***
649      !!
650      !! ** Purpose :   allocate and initialize arrays for Prather advection
651      !!-------------------------------------------------------------------
652      INTEGER ::   ierr
653      !!-------------------------------------------------------------------
654      !
655      !                             !* allocate prather fields
656      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
657         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
658         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
659         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
660         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
661         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
662         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
663         !
664         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
665         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
666         !
667         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
668         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
669         &      STAT = ierr )
670      !
671      CALL mpp_sum( 'icedyn_adv_pra', ierr )
672      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
673      !
674      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
675      !
676   END SUBROUTINE adv_pra_init
677
678
679   SUBROUTINE adv_pra_rst( cdrw, kt )
680      !!---------------------------------------------------------------------
681      !!                   ***  ROUTINE adv_pra_rst  ***
682      !!                     
683      !! ** Purpose :   Read or write file in restart file
684      !!
685      !! ** Method  :   use of IOM library
686      !!----------------------------------------------------------------------
687      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
688      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
689      !
690      INTEGER ::   jk, jl   ! dummy loop indices
691      INTEGER ::   iter     ! local integer
692      INTEGER ::   id1      ! local integer
693      CHARACTER(len=25) ::   znam
694      CHARACTER(len=2)  ::   zchar, zchar1
695      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
696      !!----------------------------------------------------------------------
697      !
698      !                                      !==========================!
699      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
700         !                                   !==========================!
701         !
702         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
703         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
704         ENDIF
705         !
706         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
707            !
708            !                                                        ! ice thickness
709            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
710            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
711            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
712            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
713            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
714            !                                                        ! snow thickness
715            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
716            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
717            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
718            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
719            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
720            !                                                        ! ice concentration
721            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
722            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
723            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
724            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
725            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
726            !                                                        ! ice salinity
727            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
728            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
729            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
730            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
731            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
732            !                                                        ! ice age
733            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
734            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
735            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
736            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
737            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
738            !                                                        ! snow layers heat content
739            DO jk = 1, nlay_s
740               WRITE(zchar1,'(I2.2)') jk
741               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
742               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
743               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
744               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
745               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
746            END DO
747            !                                                        ! ice layers heat content
748            DO jk = 1, nlay_i
749               WRITE(zchar1,'(I2.2)') jk
750               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
751               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
752               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
753               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
754               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
755            END DO
756            !
757            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
758               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
759               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
760               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
761               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
762               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
763               !                                                     ! melt pond volume
764               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
765               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
766               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
767               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
768               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
769            ENDIF
770            !
771         ELSE                                   !**  start rheology from rest  **!
772            !
773            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
774            !
775            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
776            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
777            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
778            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
779            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
780            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
781            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
782            IF( ln_pnd_H12 ) THEN
783               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
784               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
785            ENDIF
786         ENDIF
787         !
788         !                                   !=====================================!
789      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
790         !                                   !=====================================!
791         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
792         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
793         !
794         !
795         ! In case Prather scheme is used for advection, write second order moments
796         ! ------------------------------------------------------------------------
797         !
798         !                                                           ! ice thickness
799         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
800         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
801         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
802         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
803         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
804         !                                                           ! snow thickness
805         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
806         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
807         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
808         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
809         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
810         !                                                           ! ice concentration
811         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
812         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
813         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
814         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
815         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
816         !                                                           ! ice salinity
817         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
818         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
819         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
820         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
821         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
822         !                                                           ! ice age
823         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
824         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
825         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
826         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
827         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
828         !                                                           ! snow layers heat content
829         DO jk = 1, nlay_s
830            WRITE(zchar1,'(I2.2)') jk
831            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
832            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
833            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
834            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
835            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
836         END DO
837         !                                                           ! ice layers heat content
838         DO jk = 1, nlay_i
839            WRITE(zchar1,'(I2.2)') jk
840            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
841            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
842            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
843            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
844            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
845         END DO
846         !
847         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
848            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
849            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
850            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
851            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
852            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
853            !                                                        ! melt pond volume
854            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
855            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
856            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
857            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
858            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
859         ENDIF
860         !
861      ENDIF
862      !
863   END SUBROUTINE adv_pra_rst
864
865#else
866   !!----------------------------------------------------------------------
867   !!   Default option            Dummy module        NO SI3 sea-ice model
868   !!----------------------------------------------------------------------
869#endif
870
871   !!======================================================================
872END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.