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

Last change on this file since 11612 was 11612, checked in by clem, 2 years ago

make Prather the default advection scheme for ice simulations and drastically reduce the number of communications in Prather routine (by a factor of 10)

  • Property svn:keywords set to Id
File size: 52.2 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     ! lead fraction
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      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
236      !
237   END SUBROUTINE ice_dyn_adv_pra
238   
239   
240   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
241      &              psx, psxx, psy , psyy, psxy )
242      !!----------------------------------------------------------------------
243      !!                **  routine adv_x  **
244      !! 
245      !! ** purpose :   Computes and adds the advection trend to sea-ice
246      !!                variable on x axis
247      !!----------------------------------------------------------------------
248      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
249      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
250      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
251      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
252      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
253      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
254      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
255      !!
256      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
257      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
258      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
259      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
260      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
261      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
262      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
263      !-----------------------------------------------------------------------
264      !
265      jcat = SIZE( ps0 , 3 )   ! size of input arrays
266      !
267      DO jl = 1, jcat   ! loop on categories
268         !
269         ! Limitation of moments.                                           
270         DO jj = 2, jpjm1
271            DO ji = 1, jpi
272               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
273               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
274               !
275               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
276               zs1max  = 1.5 * zslpmax
277               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
278               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
279                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
280               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
281
282               ps0 (ji,jj,jl) = zslpmax 
283               psx (ji,jj,jl) = zs1new         * rswitch
284               psxx(ji,jj,jl) = zs2new         * rswitch
285               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
286               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
287               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
288            END DO
289         END DO
290
291         !  Calculate fluxes and moments between boxes i<-->i+1             
292         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0
293            DO ji = 1, jpi
294               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
295               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl)
296               zalfq        =  zalf * zalf
297               zalf1        =  1.0 - zalf
298               zalf1q       =  zalf1 * zalf1
299               !
300               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
301               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
302               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
303               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
304               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
305               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
306               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
307
308               !  Readjust moments remaining in the box.
309               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
310               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
311               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
312               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
313               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
314               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
315               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
316            END DO
317         END DO
318
319         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
320            DO ji = 1, fs_jpim1
321               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl) 
322               zalg  (ji,jj) = zalf
323               zalfq         = zalf * zalf
324               zalf1         = 1.0 - zalf
325               zalg1 (ji,jj) = zalf1
326               zalf1q        = zalf1 * zalf1
327               zalg1q(ji,jj) = zalf1q
328               !
329               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
330               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
331                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
332               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
333               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
334               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
335               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
336               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
337            END DO
338         END DO
339
340         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
341            DO ji = fs_2, fs_jpim1
342               zbt  =       zbet(ji-1,jj)
343               zbt1 = 1.0 - zbet(ji-1,jj)
344               !
345               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
346               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
347               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
348               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
349               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
350               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
351               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
352            END DO
353         END DO
354
355         !   Put the temporary moments into appropriate neighboring boxes.   
356         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
357            DO ji = fs_2, fs_jpim1
358               zbt  =       zbet(ji-1,jj)
359               zbt1 = 1.0 - zbet(ji-1,jj)
360               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
361               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
362               zalf1         = 1.0 - zalf
363               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
364               !
365               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
366               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
367               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
368                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
369                  &            + zbt1 * psxx(ji,jj,jl)
370               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
371                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
372                  &            + zbt1 * psxy(ji,jj,jl)
373               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
374               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
375            END DO
376         END DO
377
378         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0.
379            DO ji = fs_2, fs_jpim1
380               zbt  =       zbet(ji,jj)
381               zbt1 = 1.0 - zbet(ji,jj)
382               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
383               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
384               zalf1         = 1.0 - zalf
385               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
386               !
387               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
388               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
389               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
390                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
391                  &                                           + ( zalf1 - zalf ) * ztemp ) )
392               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
393                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
394               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
395               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
396            END DO
397         END DO
398
399      END DO
400
401      !-- Lateral boundary conditions
402      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
403         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
404         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
405      !
406   END SUBROUTINE adv_x
407
408
409   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
410      &              psx, psxx, psy , psyy, psxy )
411      !!---------------------------------------------------------------------
412      !!                **  routine adv_y  **
413      !!           
414      !! ** purpose :   Computes and adds the advection trend to sea-ice
415      !!                variable on y axis
416      !!---------------------------------------------------------------------
417      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
418      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
419      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
420      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
421      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
422      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
423      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
424      !!
425      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
426      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
427      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
428      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
429      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
430      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
431      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
432      !---------------------------------------------------------------------
433      !
434      jcat = SIZE( ps0 , 3 )   ! size of input arrays
435      !     
436      DO jl = 1, jcat   ! loop on categories
437         !
438         ! Limitation of moments.
439         DO jj = 1, jpj
440            DO ji = fs_2, fs_jpim1
441               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
442               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
443               !
444               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
445               zs1max  = 1.5 * zslpmax
446               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
447               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
448                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
449               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
450               !
451               ps0 (ji,jj,jl) = zslpmax 
452               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
453               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
454               psy (ji,jj,jl) = zs1new         * rswitch
455               psyy(ji,jj,jl) = zs2new         * rswitch
456               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
457            END DO
458         END DO
459 
460         !  Calculate fluxes and moments between boxes j<-->j+1             
461         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
462            DO ji = fs_2, fs_jpim1
463               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
464               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl)
465               zalfq        =  zalf * zalf
466               zalf1        =  1.0 - zalf
467               zalf1q       =  zalf1 * zalf1
468               !
469               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
470               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
471               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
472               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
473               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
474               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
475               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
476               !
477               !  Readjust moments remaining in the box.
478               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
479               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
480               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
481               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
482               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
483               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
484               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
485            END DO
486         END DO
487         !
488         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
489            DO ji = fs_2, fs_jpim1
490               zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl) 
491               zalg  (ji,jj) = zalf
492               zalfq         = zalf * zalf
493               zalf1         = 1.0 - zalf
494               zalg1 (ji,jj) = zalf1
495               zalf1q        = zalf1 * zalf1
496               zalg1q(ji,jj) = zalf1q
497               !
498               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl)
499               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
500                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
501               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
502               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
503               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
504               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
505               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
506            END DO
507         END DO
508
509         !  Readjust moments remaining in the box.
510         DO jj = 2, jpjm1
511            DO ji = fs_2, fs_jpim1
512               zbt  =         zbet(ji,jj-1)
513               zbt1 = ( 1.0 - zbet(ji,jj-1) )
514               !
515               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
516               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
517               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
518               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
519               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
520               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
521               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
522            END DO
523         END DO
524
525         !   Put the temporary moments into appropriate neighboring boxes.   
526         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
527            DO ji = fs_2, fs_jpim1
528               zbt  =       zbet(ji,jj-1)
529               zbt1 = 1.0 - zbet(ji,jj-1)
530               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
531               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
532               zalf1         = 1.0 - zalf
533               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
534               !
535               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
536               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
537                  &             + zbt1 * psy(ji,jj,jl) 
538               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
539                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
540                  &             + zbt1 * psyy(ji,jj,jl)
541               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
542                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
543                  &             + zbt1 * psxy(ji,jj,jl)
544               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
545               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
546            END DO
547         END DO
548
549         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
550            DO ji = fs_2, fs_jpim1
551               zbt  =       zbet(ji,jj)
552               zbt1 = 1.0 - zbet(ji,jj)
553               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
554               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
555               zalf1         = 1.0 - zalf
556               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
557               !
558               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
559               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
560               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
561                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
562                  &                                            + ( zalf1 - zalf ) * ztemp ) )
563               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
564                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
565               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
566               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
567            END DO
568         END DO
569
570      END DO
571
572      !-- Lateral boundary conditions
573      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
574         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
575         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
576      !
577   END SUBROUTINE adv_y
578
579
580   SUBROUTINE adv_pra_init
581      !!-------------------------------------------------------------------
582      !!                  ***  ROUTINE adv_pra_init  ***
583      !!
584      !! ** Purpose :   allocate and initialize arrays for Prather advection
585      !!-------------------------------------------------------------------
586      INTEGER ::   ierr
587      !!-------------------------------------------------------------------
588      !
589      !                             !* allocate prather fields
590      ALLOCATE( sxopw(jpi,jpj,1)   , syopw(jpi,jpj,1)   , sxxopw(jpi,jpj,1)   , syyopw(jpi,jpj,1)   , sxyopw(jpi,jpj,1)   ,   &
591         &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
592         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
593         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
594         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
595         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
596         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
597         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
598         !
599         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
600         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
601         !
602         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
603         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
604         &      STAT = ierr )
605      !
606      CALL mpp_sum( 'icedyn_adv_pra', ierr )
607      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
608      !
609      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
610      !
611   END SUBROUTINE adv_pra_init
612
613
614   SUBROUTINE adv_pra_rst( cdrw, kt )
615      !!---------------------------------------------------------------------
616      !!                   ***  ROUTINE adv_pra_rst  ***
617      !!                     
618      !! ** Purpose :   Read or write file in restart file
619      !!
620      !! ** Method  :   use of IOM library
621      !!----------------------------------------------------------------------
622      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
623      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
624      !
625      INTEGER ::   jk, jl   ! dummy loop indices
626      INTEGER ::   iter     ! local integer
627      INTEGER ::   id1      ! local integer
628      CHARACTER(len=25) ::   znam
629      CHARACTER(len=2)  ::   zchar, zchar1
630      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
631      !!----------------------------------------------------------------------
632      !
633      !                                      !==========================!
634      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
635         !                                   !==========================!
636         !
637         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0
638         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
639         ENDIF
640         !
641         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
642            !
643            !                                                        ! ice thickness
644            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
645            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
646            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
647            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
648            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
649            !                                                        ! snow thickness
650            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
651            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
652            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
653            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
654            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
655            !                                                        ! lead fraction
656            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
657            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
658            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
659            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
660            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
661            !                                                        ! ice salinity
662            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
663            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
664            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
665            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
666            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
667            !                                                        ! ice age
668            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
669            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
670            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
671            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
672            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
673            !                                                        ! open water in sea ice
674            CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw  )
675            CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw  )
676            CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )
677            CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
678            CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
679            !                                                        ! snow layers heat content
680            DO jk = 1, nlay_s
681               WRITE(zchar1,'(I2.2)') jk
682               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
683               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
684               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
685               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
686               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
687            END DO
688            !                                                        ! ice layers heat content
689            DO jk = 1, nlay_i
690               WRITE(zchar1,'(I2.2)') jk
691               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
692               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
693               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
694               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
695               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
696            END DO
697            !
698            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
699               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
700               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
701               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
702               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
703               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
704               !                                                     ! melt pond volume
705               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
706               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
707               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
708               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
709               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
710            ENDIF
711            !
712         ELSE                                   !**  start rheology from rest  **!
713            !
714            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
715            !
716            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
717            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
718            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction
719            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
720            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
721            sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
722            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
723            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
724            IF( ln_pnd_H12 ) THEN
725               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
726               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
727            ENDIF
728         ENDIF
729         !
730         !                                   !=====================================!
731      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
732         !                                   !=====================================!
733         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
734         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
735         !
736         !
737         ! In case Prather scheme is used for advection, write second order moments
738         ! ------------------------------------------------------------------------
739         !
740         !                                                           ! ice thickness
741         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
742         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
743         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
744         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
745         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
746         !                                                           ! snow thickness
747         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
748         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
749         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
750         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
751         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
752         !                                                           ! lead fraction
753         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
754         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
755         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
756         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
757         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
758         !                                                           ! ice salinity
759         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
760         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
761         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
762         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
763         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
764         !                                                           ! ice age
765         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
766         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
767         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
768         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
769         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
770         !                                                           ! open water in sea ice
771         CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw  )
772         CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw  )
773         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )
774         CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
775         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
776         !                                                           ! snow layers heat content
777         DO jk = 1, nlay_s
778            WRITE(zchar1,'(I2.2)') jk
779            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
780            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
781            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
782            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
783            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
784         END DO
785         !                                                           ! ice layers heat content
786         DO jk = 1, nlay_i
787            WRITE(zchar1,'(I2.2)') jk
788            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
789            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
790            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
791            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
792            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
793         END DO
794         !
795         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
796            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
797            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
798            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
799            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
800            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
801            !                                                        ! melt pond volume
802            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
803            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
804            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
805            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
806            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
807         ENDIF
808         !
809      ENDIF
810      !
811   END SUBROUTINE adv_pra_rst
812
813#else
814   !!----------------------------------------------------------------------
815   !!   Default option            Dummy module        NO SI3 sea-ice model
816   !!----------------------------------------------------------------------
817#endif
818
819   !!======================================================================
820END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.