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

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

source: NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/ICE/icedyn_adv_pra.F90 @ 11692

Last change on this file since 11692 was 11692, checked in by francesca, 4 years ago

Update branch to integrate the development starting from the current v4.01 ready trunk

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