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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_pra.F90 in NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90 @ 12949

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

make sure ice and snow temperatures do not overshoot after advection (both with Prather or UMx), as well as ice salinity. Include the salinity bounds (rn_simin, rn_simax) inside the icethd_sal routine. There is no more problems with super cold temperatures now. Sette tests passed.
  • Property svn:keywords set to Id
File size: 67.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 "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      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
365      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
366      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
367      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
368      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
369      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
370      !-----------------------------------------------------------------------
371      !
372      jcat = SIZE( ps0 , 3 )   ! size of input arrays
373      !
374      DO jl = 1, jcat   ! loop on categories
375         !
376         ! Limitation of moments.                                           
377         DO jj = 2, jpjm1
378            DO ji = 1, jpi
379               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
380               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
381               !
382               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
383               zs1max  = 1.5 * zslpmax
384               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
385               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
386                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
387               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
388
389               ps0 (ji,jj,jl) = zslpmax 
390               psx (ji,jj,jl) = zs1new         * rswitch
391               psxx(ji,jj,jl) = zs2new         * rswitch
392               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
393               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
394               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
395            END DO
396         END DO
397
398         !  Calculate fluxes and moments between boxes i<-->i+1             
399         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0
400            DO ji = 1, jpi
401               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
402               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
403               zalfq        =  zalf * zalf
404               zalf1        =  1.0 - zalf
405               zalf1q       =  zalf1 * zalf1
406               !
407               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
408               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
409               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
410               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
411               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
412               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
413               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
414
415               !  Readjust moments remaining in the box.
416               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
417               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
418               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
419               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
420               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
421               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
422               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
423            END DO
424         END DO
425
426         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
427            DO ji = 1, fs_jpim1
428               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
429               zalg  (ji,jj) = zalf
430               zalfq         = zalf * zalf
431               zalf1         = 1.0 - zalf
432               zalg1 (ji,jj) = zalf1
433               zalf1q        = zalf1 * zalf1
434               zalg1q(ji,jj) = zalf1q
435               !
436               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
437               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
438                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
439               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
440               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
441               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
442               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
443               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
444            END DO
445         END DO
446
447         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
448            DO ji = fs_2, fs_jpim1
449               zbt  =       zbet(ji-1,jj)
450               zbt1 = 1.0 - zbet(ji-1,jj)
451               !
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            END DO
460         END DO
461
462         !   Put the temporary moments into appropriate neighboring boxes.   
463         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
464            DO ji = fs_2, fs_jpim1
465               zbt  =       zbet(ji-1,jj)
466               zbt1 = 1.0 - zbet(ji-1,jj)
467               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
468               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
469               zalf1         = 1.0 - zalf
470               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
471               !
472               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
473               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
474               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
475                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
476                  &            + zbt1 * psxx(ji,jj,jl)
477               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
478                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
479                  &            + zbt1 * psxy(ji,jj,jl)
480               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
481               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
482            END DO
483         END DO
484
485         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0.
486            DO ji = fs_2, fs_jpim1
487               zbt  =       zbet(ji,jj)
488               zbt1 = 1.0 - zbet(ji,jj)
489               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
490               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
491               zalf1         = 1.0 - zalf
492               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
493               !
494               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
495               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
496               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
497                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
498                  &                                           + ( zalf1 - zalf ) * ztemp ) )
499               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
500                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
501               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
502               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
503            END DO
504         END DO
505
506      END DO
507
508      !-- Lateral boundary conditions
509      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
510         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
511         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
512      !
513   END SUBROUTINE adv_x
514
515
516   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
517      &              psx, psxx, psy , psyy, psxy )
518      !!---------------------------------------------------------------------
519      !!                **  routine adv_y  **
520      !!           
521      !! ** purpose :   Computes and adds the advection trend to sea-ice
522      !!                variable on y axis
523      !!---------------------------------------------------------------------
524      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
525      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
526      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
527      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
528      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
529      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
530      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
531      !!
532      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
533      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
534      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
535      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
536      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
537      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
538      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
539      !---------------------------------------------------------------------
540      !
541      jcat = SIZE( ps0 , 3 )   ! size of input arrays
542      !     
543      DO jl = 1, jcat   ! loop on categories
544         !
545         ! Limitation of moments.
546         DO jj = 1, jpj
547            DO ji = fs_2, fs_jpim1
548               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
549               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
550               !
551               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
552               zs1max  = 1.5 * zslpmax
553               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
554               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
555                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
556               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
557               !
558               ps0 (ji,jj,jl) = zslpmax 
559               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
560               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
561               psy (ji,jj,jl) = zs1new         * rswitch
562               psyy(ji,jj,jl) = zs2new         * rswitch
563               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
564            END DO
565         END DO
566 
567         !  Calculate fluxes and moments between boxes j<-->j+1             
568         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
569            DO ji = fs_2, fs_jpim1
570               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
571               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
572               zalfq        =  zalf * zalf
573               zalf1        =  1.0 - zalf
574               zalf1q       =  zalf1 * zalf1
575               !
576               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
577               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
578               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
579               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
580               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
581               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
582               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
583               !
584               !  Readjust moments remaining in the box.
585               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
586               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
587               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
588               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
589               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
590               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
591               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
592            END DO
593         END DO
594         !
595         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
596            DO ji = fs_2, fs_jpim1
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         !  Readjust moments remaining in the box.
617         DO jj = 2, jpjm1
618            DO ji = fs_2, fs_jpim1
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            END DO
630         END DO
631
632         !   Put the temporary moments into appropriate neighboring boxes.   
633         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
634            DO ji = fs_2, fs_jpim1
635               zbt  =       zbet(ji,jj-1)
636               zbt1 = 1.0 - zbet(ji,jj-1)
637               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
638               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
639               zalf1         = 1.0 - zalf
640               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
641               !
642               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
643               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
644                  &             + zbt1 * psy(ji,jj,jl) 
645               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
646                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
647                  &             + zbt1 * psyy(ji,jj,jl)
648               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
649                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
650                  &             + zbt1 * psxy(ji,jj,jl)
651               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
652               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
653            END DO
654         END DO
655
656         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
657            DO ji = fs_2, fs_jpim1
658               zbt  =       zbet(ji,jj)
659               zbt1 = 1.0 - zbet(ji,jj)
660               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
661               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
662               zalf1         = 1.0 - zalf
663               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
664               !
665               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
666               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
667               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
668                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
669                  &                                            + ( zalf1 - zalf ) * ztemp ) )
670               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
671                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
672               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
673               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
674            END DO
675         END DO
676
677      END DO
678
679      !-- Lateral boundary conditions
680      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
681         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
682         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
683      !
684   END SUBROUTINE adv_y
685
686
687   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
688      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
689      !!-------------------------------------------------------------------
690      !!                  ***  ROUTINE Hbig  ***
691      !!
692      !! ** Purpose : Thickness correction in case advection scheme creates
693      !!              abnormally tick ice or snow
694      !!
695      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
696      !!                 (before advection) and reduce it by adapting ice concentration
697      !!              2- check whether snow thickness is larger than the surrounding 9-points
698      !!                 (before advection) and reduce it by sending the excess in the ocean
699      !!
700      !! ** input   : Max thickness of the surrounding 9-points
701      !!-------------------------------------------------------------------
702      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step
703      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts
704      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max
705      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max
706      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
707      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
708      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
709      !
710      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
711      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
712      !!-------------------------------------------------------------------
713      !
714      z1_dt = 1._wp / pdt
715      !
716      DO jl = 1, jpl
717         DO jj = 1, jpj
718            DO ji = 1, jpi
719               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
720                  !
721                  !                               ! -- check h_ip -- !
722                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
723                  IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
724                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
725                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
726                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
727                     ENDIF
728                  ENDIF
729                  !
730                  !                               ! -- check h_i -- !
731                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
732                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
733                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
734                     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)
735                  ENDIF
736                  !
737                  !                               ! -- check h_s -- !
738                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
739                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
740                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
741                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
742                     !
743                     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
744                     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
745                     !
746                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
747                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
748                  ENDIF           
749                  !                 
750                  !                               ! -- check s_i -- !
751                  ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
752                  zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
753                  IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
754                     zfra = psi_max(ji,jj,jl) / zsi
755                     sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
756                     psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
757                  ENDIF
758                  !
759               ENDIF
760            END DO
761         END DO
762      END DO 
763      !
764      !                                           ! -- check e_i/v_i -- !
765      DO jl = 1, jpl
766         DO jk = 1, nlay_i
767            DO jj = 1, jpj
768               DO ji = 1, jpi
769                  IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
770                     ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
771                     zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
772                     IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
773                        zfra = pei_max(ji,jj,jk,jl) / zei
774                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
775                        pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
776                     ENDIF
777                  ENDIF
778               END DO
779            END DO
780         END DO
781      END DO
782      !                                           ! -- check e_s/v_s -- !
783      DO jl = 1, jpl
784         DO jk = 1, nlay_s
785            DO jj = 1, jpj
786               DO ji = 1, jpi
787                  IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
788                     ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
789                     zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
790                     IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
791                        zfra = pes_max(ji,jj,jk,jl) / zes
792                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
793                        pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
794                     ENDIF
795                  ENDIF
796               END DO
797            END DO
798         END DO
799      END DO
800      !
801   END SUBROUTINE Hbig
802
803
804   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
805      !!-------------------------------------------------------------------
806      !!                  ***  ROUTINE Hsnow  ***
807      !!
808      !! ** Purpose : 1- Check snow load after advection
809      !!              2- Correct pond concentration to avoid a_ip > a_i
810      !!
811      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
812      !!              then put the snow excess in the ocean
813      !!
814      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards
815      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
816      !!              make the snow very thick (if concentration decreases drastically)
817      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
818      !!-------------------------------------------------------------------
819      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
820      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
821      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
822      !
823      INTEGER  ::   ji, jj, jl   ! dummy loop indices
824      REAL(wp) ::   z1_dt, zvs_excess, zfra
825      !!-------------------------------------------------------------------
826      !
827      z1_dt = 1._wp / pdt
828      !
829      ! -- check snow load -- !
830      DO jl = 1, jpl
831         DO jj = 1, jpj
832            DO ji = 1, jpi
833               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
834                  !
835                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
836                  !
837                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
838                     ! put snow excess in the ocean
839                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
840                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
841                     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
842                     ! correct snow volume and heat content
843                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
844                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
845                  ENDIF
846                  !
847               ENDIF
848            END DO
849         END DO
850      END DO
851      !
852      !-- correct pond concentration to avoid a_ip > a_i -- !
853      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
854      !
855   END SUBROUTINE Hsnow
856
857
858   SUBROUTINE adv_pra_init
859      !!-------------------------------------------------------------------
860      !!                  ***  ROUTINE adv_pra_init  ***
861      !!
862      !! ** Purpose :   allocate and initialize arrays for Prather advection
863      !!-------------------------------------------------------------------
864      INTEGER ::   ierr
865      !!-------------------------------------------------------------------
866      !
867      !                             !* allocate prather fields
868      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
869         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
870         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
871         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
872         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
873         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
874         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
875         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   &
876         !
877         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
878         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
879         !
880         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
881         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
882         &      STAT = ierr )
883      !
884      CALL mpp_sum( 'icedyn_adv_pra', ierr )
885      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
886      !
887      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
888      !
889   END SUBROUTINE adv_pra_init
890
891
892   SUBROUTINE adv_pra_rst( cdrw, kt )
893      !!---------------------------------------------------------------------
894      !!                   ***  ROUTINE adv_pra_rst  ***
895      !!                     
896      !! ** Purpose :   Read or write file in restart file
897      !!
898      !! ** Method  :   use of IOM library
899      !!----------------------------------------------------------------------
900      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
901      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
902      !
903      INTEGER ::   jk, jl   ! dummy loop indices
904      INTEGER ::   iter     ! local integer
905      INTEGER ::   id1      ! local integer
906      CHARACTER(len=25) ::   znam
907      CHARACTER(len=2)  ::   zchar, zchar1
908      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
909      !!----------------------------------------------------------------------
910      !
911      !                                      !==========================!
912      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
913         !                                   !==========================!
914         !
915         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
916         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
917         ENDIF
918         !
919         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
920            !
921            !                                                        ! ice thickness
922            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
923            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
924            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
925            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
926            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
927            !                                                        ! snow thickness
928            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
929            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
930            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
931            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
932            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
933            !                                                        ! ice concentration
934            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
935            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
936            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
937            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
938            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
939            !                                                        ! ice salinity
940            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
941            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
942            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
943            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
944            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
945            !                                                        ! ice age
946            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
947            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
948            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
949            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
950            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
951            !                                                        ! snow layers heat content
952            DO jk = 1, nlay_s
953               WRITE(zchar1,'(I2.2)') jk
954               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
955               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
956               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
957               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
958               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
959            END DO
960            !                                                        ! ice layers heat content
961            DO jk = 1, nlay_i
962               WRITE(zchar1,'(I2.2)') jk
963               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
964               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
965               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
966               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
967               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
968            END DO
969            !
970            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction
971               IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
972                  CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
973                  CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
974                  CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
975                  CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
976                  CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
977                  !                                                     ! melt pond volume
978                  CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
979                  CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
980                  CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
981                  CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
982                  CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
983               ELSE
984                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction
985                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume
986               ENDIF
987                  !
988               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
989                  IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
990                     CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl  )
991                     CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl  )
992                     CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl )
993                     CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl )
994                     CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl )
995                  ELSE
996                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume
997                  ENDIF
998               ENDIF
999            ENDIF
1000            !
1001         ELSE                                   !**  start rheology from rest  **!
1002            !
1003            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
1004            !
1005            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
1006            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
1007            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
1008            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
1009            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
1010            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
1011            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
1012            IF( ln_pnd_LEV ) THEN
1013               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
1014               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
1015               IF ( ln_pnd_lids ) THEN
1016                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
1017               ENDIF
1018            ENDIF
1019         ENDIF
1020         !
1021         !                                   !=====================================!
1022      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
1023         !                                   !=====================================!
1024         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
1025         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1026         !
1027         !
1028         ! In case Prather scheme is used for advection, write second order moments
1029         ! ------------------------------------------------------------------------
1030         !
1031         !                                                           ! ice thickness
1032         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
1033         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
1034         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
1035         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
1036         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
1037         !                                                           ! snow thickness
1038         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
1039         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
1040         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
1041         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
1042         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
1043         !                                                           ! ice concentration
1044         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
1045         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
1046         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
1047         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
1048         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
1049         !                                                           ! ice salinity
1050         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
1051         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
1052         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
1053         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
1054         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
1055         !                                                           ! ice age
1056         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
1057         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
1058         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
1059         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
1060         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
1061         !                                                           ! snow layers heat content
1062         DO jk = 1, nlay_s
1063            WRITE(zchar1,'(I2.2)') jk
1064            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1065            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1066            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1067            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1068            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1069         END DO
1070         !                                                           ! ice layers heat content
1071         DO jk = 1, nlay_i
1072            WRITE(zchar1,'(I2.2)') jk
1073            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1074            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1075            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1076            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1077            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1078         END DO
1079         !
1080         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction
1081            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
1082            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
1083            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
1084            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
1085            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
1086            !                                                        ! melt pond volume
1087            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
1088            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
1089            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
1090            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
1091            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
1092            !
1093            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
1094               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
1095               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
1096               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
1097               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
1098               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
1099            ENDIF
1100         ENDIF
1101         !
1102      ENDIF
1103      !
1104   END SUBROUTINE adv_pra_rst
1105
1106#else
1107   !!----------------------------------------------------------------------
1108   !!   Default option            Dummy module        NO SI3 sea-ice model
1109   !!----------------------------------------------------------------------
1110#endif
1111
1112   !!======================================================================
1113END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.