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

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

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

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

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