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

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

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

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