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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90 @ 13554

Last change on this file since 13554 was 13554, checked in by clem, 3 years ago

4.0-HEAD: 1st step to drastically reduce the number of communications in Parther advection scheme (SI3). It changes slightly run.stat because some loops are not written in the same order but outputs from a 1year long creg025 simulation are identical

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