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_r11756_SI3restart_XIOS/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/ICE/icedyn_adv_pra.F90 @ 11837

Last change on this file since 11837 was 11837, checked in by andmirek, 4 years ago

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