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

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

adapt ice advection and rheology to nn_hls=2. Number of mpi communications are reduced. I also changed lbc_lnk routine to be able to do lbc on 30 variables at once.

  • Property svn:keywords set to Id
File size: 72.1 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            &                          , 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            &                          , 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            IF( ln_pnd_lids ) THEN
283               CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
284                  &                          , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
285                  &                          , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
286                  &                          , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  &
287                  &                          , 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            ELSE
290               CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
291                  &                          , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
292                  &                          , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
293                  &                          , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )
294            ENDIF
295         ENDIF
296
297         ! --- Recover the properties from their contents --- !
298         DO jl = 1, jpl
299            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
300            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
301            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
302            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
303            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
304            DO jk = 1, nlay_s
305               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
306            END DO
307            DO jk = 1, nlay_i
308               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
309            END DO
310            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
311               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
312               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
313               IF ( ln_pnd_lids ) THEN
314                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
315               ENDIF
316            ENDIF
317         END DO
318         !
319         ! derive open water from ice concentration
320         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
321         DO_2D( 0, 0, 0, 0 )
322            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
323               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
324         END_2D
325         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp )
326         !
327         ! --- diagnostics --- !
328         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
329            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
330            &                                        - zdiag_adv_mass(:,:) ) * z1_dt
331         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
332            &                                        - zdiag_adv_salt(:,:) ) * z1_dt
333         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
334            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
335            &                                        - zdiag_adv_heat(:,:) ) * z1_dt
336         !
337         ! --- Ensure non-negative fields --- !
338         !     Remove negative values (conservation is ensured)
339         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
340         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 )
341         !
342         ! --- Make sure ice thickness is not too big --- !
343         !     (because ice thickness can be too large where ice concentration is very small)
344         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, &
345            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
346         !
347         ! --- Ensure snow load is not too big --- !
348         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
349         !
350      END DO
351      !
352      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
353      !
354   END SUBROUTINE ice_dyn_adv_pra
355
356
357   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
358      &              psx, psxx, psy , psyy, psxy )
359      !!----------------------------------------------------------------------
360      !!                **  routine adv_x  **
361      !!
362      !! ** purpose :   Computes and adds the advection trend to sea-ice
363      !!                variable on x axis
364      !!----------------------------------------------------------------------
365      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
366      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
367      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
368      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
369      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
370      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
371      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
372      !!
373      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
374      INTEGER  ::   jj0                                  ! dummy loop indices
375      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
376      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
377      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
378      REAL(wp) ::   zpsm, zps0
379      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
380      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
381      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
382      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
383      !-----------------------------------------------------------------------
384      ! in order to avoid lbc_lnk (communications):
385      !    jj loop must be 1:jpj   if adv_x is called first
386      !                and 2:jpj-1 if adv_x is called second
387      jj0 = NINT(pcrh)
388      !
389      jcat = SIZE( ps0 , 3 )   ! size of input arrays
390      !
391      DO jl = 1, jcat   ! loop on categories
392         !
393         ! Limitation of moments.
394         DO jj = Njs0 - jj0, Nje0 + jj0
395
396            DO ji = Nis0 - 1, Nie0 + 1
397
398               zpsm  = psm (ji,jj,jl) ! optimization
399               zps0  = ps0 (ji,jj,jl)
400               zpsx  = psx (ji,jj,jl)
401               zpsxx = psxx(ji,jj,jl)
402               zpsy  = psy (ji,jj,jl)
403               zpsyy = psyy(ji,jj,jl)
404               zpsxy = psxy(ji,jj,jl)
405
406               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)
407               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
408               !
409               zslpmax = MAX( 0._wp, zps0 )
410               zs1max  = 1.5 * zslpmax
411               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) )
412               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
413               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
414
415               zps0  = zslpmax
416               zpsx  = zs1new  * rswitch
417               zpsxx = zs2new  * rswitch
418               zpsy  = zpsy    * rswitch
419               zpsyy = zpsyy   * rswitch
420               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
421
422               !  Calculate fluxes and moments between boxes i<-->i+1
423               !                                !  Flux from i to i+1 WHEN u GT 0
424               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
425               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
426               zalfq        =  zalf * zalf
427               zalf1        =  1.0 - zalf
428               zalf1q       =  zalf1 * zalf1
429               !
430               zfm (ji,jj)  =  zalf  *   zpsm
431               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
432               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx )
433               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq
434               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy )
435               zfxy(ji,jj)  =  zalfq *   zpsxy
436               zfyy(ji,jj)  =  zalf  *   zpsyy
437
438               !                                !  Readjust moments remaining in the box.
439               zpsm  =  zpsm  - zfm(ji,jj)
440               zps0  =  zps0  - zf0(ji,jj)
441               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
442               zpsxx =  zalf1  * zalf1q * zpsxx
443               zpsy  =  zpsy  - zfy (ji,jj)
444               zpsyy =  zpsyy - zfyy(ji,jj)
445               zpsxy =  zalf1q * zpsxy
446               !
447               psm (ji,jj,jl) = zpsm ! optimization
448               ps0 (ji,jj,jl) = zps0
449               psx (ji,jj,jl) = zpsx
450               psxx(ji,jj,jl) = zpsxx
451               psy (ji,jj,jl) = zpsy
452               psyy(ji,jj,jl) = zpsyy
453               psxy(ji,jj,jl) = zpsxy
454               !
455            END DO
456
457            DO ji = Nis0 - 1, Nie0
458               !                                !  Flux from i+1 to i when u LT 0.
459               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
460               zalg  (ji,jj) = zalf
461               zalfq         = zalf * zalf
462               zalf1         = 1.0 - zalf
463               zalg1 (ji,jj) = zalf1
464               zalf1q        = zalf1 * zalf1
465               zalg1q(ji,jj) = zalf1q
466               !
467               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
468               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
469                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
470               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
471               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
472               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
473               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
474               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
475            END DO
476
477            DO ji = Nis0, Nie0
478               !
479               zpsm  = psm (ji,jj,jl) ! optimization
480               zps0  = ps0 (ji,jj,jl)
481               zpsx  = psx (ji,jj,jl)
482               zpsxx = psxx(ji,jj,jl)
483               zpsy  = psy (ji,jj,jl)
484               zpsyy = psyy(ji,jj,jl)
485               zpsxy = psxy(ji,jj,jl)
486               !                                !  Readjust moments remaining in the box.
487               zbt  =       zbet(ji-1,jj)
488               zbt1 = 1.0 - zbet(ji-1,jj)
489               !
490               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
491               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
492               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
493               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
494               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) )
495               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
496               zpsxy = zalg1q(ji-1,jj) * zpsxy
497
498               !   Put the temporary moments into appropriate neighboring boxes.
499               !                                !   Flux from i to i+1 IF u GT 0.
500               zbt   =       zbet(ji-1,jj)
501               zbt1  = 1.0 - zbet(ji-1,jj)
502               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
503               zalf  = zbt * zfm(ji-1,jj) / zpsm
504               zalf1 = 1.0 - zalf
505               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
506               !
507               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
508               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
509               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            &
510                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
511                  &            + zbt1 * zpsxx
512               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            &
513                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  &
514                  &            + zbt1 * zpsxy
515               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy
516               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
517
518               !                                !  Flux from i+1 to i IF u LT 0.
519               zbt   =       zbet(ji,jj)
520               zbt1  = 1.0 - zbet(ji,jj)
521               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
522               zalf  = zbt1 * zfm(ji,jj) / zpsm
523               zalf1 = 1.0 - zalf
524               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
525               !
526               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) )
527               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
528               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
529                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    &
530                  &                         + ( zalf1 - zalf ) * ztemp ) )
531               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
532                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
533               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) )
534               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
535               !
536               psm (ji,jj,jl) = zpsm  ! optimization
537               ps0 (ji,jj,jl) = zps0
538               psx (ji,jj,jl) = zpsx
539               psxx(ji,jj,jl) = zpsxx
540               psy (ji,jj,jl) = zpsy
541               psyy(ji,jj,jl) = zpsyy
542               psxy(ji,jj,jl) = zpsxy
543            END DO
544            !
545         END DO
546         !
547      END DO
548      !
549   END SUBROUTINE adv_x
550
551
552   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
553      &              psx, psxx, psy , psyy, psxy )
554      !!---------------------------------------------------------------------
555      !!                **  routine adv_y  **
556      !!
557      !! ** purpose :   Computes and adds the advection trend to sea-ice
558      !!                variable on y axis
559      !!---------------------------------------------------------------------
560      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
561      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
562      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
563      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
564      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
565      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments
566      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
567      !!
568      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
569      INTEGER  ::   ji0                                  ! dummy loop indices
570      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
571      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
572      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
573      REAL(wp) ::   zpsm, zps0
574      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
575      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
576      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
577      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
578      !---------------------------------------------------------------------
579      ! in order to avoid lbc_lnk (communications):
580      !    ji loop must be 1:jpi   if adv_y is called first
581      !                and 2:jpi-1 if adv_y is called second
582      ji0 = NINT(pcrh)
583      !
584      jcat = SIZE( ps0 , 3 )   ! size of input arrays
585      !
586      DO jl = 1, jcat   ! loop on categories
587         !
588         ! Limitation of moments.
589         DO_2D( ji0, ji0, 1, 1 )
590            !
591            zpsm  = psm (ji,jj,jl) ! optimization
592            zps0  = ps0 (ji,jj,jl)
593            zpsx  = psx (ji,jj,jl)
594            zpsxx = psxx(ji,jj,jl)
595            zpsy  = psy (ji,jj,jl)
596            zpsyy = psyy(ji,jj,jl)
597            zpsxy = psxy(ji,jj,jl)
598            !
599            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
600            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  )
601            !
602            zslpmax = MAX( 0._wp, zps0 )
603            zs1max  = 1.5 * zslpmax
604            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) )
605            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
606            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
607            !
608            zps0  = zslpmax
609            zpsx  = zpsx  * rswitch
610            zpsxx = zpsxx * rswitch
611            zpsy  = zs1new         * rswitch
612            zpsyy = zs2new         * rswitch
613            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
614
615            !  Calculate fluxes and moments between boxes j<-->j+1
616            !                                !  Flux from j to j+1 WHEN v GT 0
617            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
618            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
619            zalfq        =  zalf * zalf
620            zalf1        =  1.0 - zalf
621            zalf1q       =  zalf1 * zalf1
622            !
623            zfm (ji,jj)  =  zalf  * zpsm
624            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )
625            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy )
626            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy
627            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy )
628            zfxy(ji,jj)  =  zalfq * zpsxy
629            zfxx(ji,jj)  =  zalf  * zpsxx
630            !
631            !                                !  Readjust moments remaining in the box.
632            zpsm   =  zpsm  - zfm(ji,jj)
633            zps0   =  zps0  - zf0(ji,jj)
634            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy )
635            zpsyy  =  zalf1 * zalf1q * zpsyy
636            zpsx   =  zpsx  - zfx(ji,jj)
637            zpsxx  =  zpsxx - zfxx(ji,jj)
638            zpsxy  =  zalf1q * zpsxy
639            !
640            psm (ji,jj,jl) = zpsm ! optimization
641            ps0 (ji,jj,jl) = zps0
642            psx (ji,jj,jl) = zpsx
643            psxx(ji,jj,jl) = zpsxx
644            psy (ji,jj,jl) = zpsy
645            psyy(ji,jj,jl) = zpsyy
646            psxy(ji,jj,jl) = zpsxy
647         END_2D
648         !
649         DO_2D( ji0, ji0, 1, 0 )
650            !                                !  Flux from j+1 to j when v LT 0.
651            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
652            zalg  (ji,jj) = zalf
653            zalfq         = zalf * zalf
654            zalf1         = 1.0 - zalf
655            zalg1 (ji,jj) = zalf1
656            zalf1q        = zalf1 * zalf1
657            zalg1q(ji,jj) = zalf1q
658            !
659            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl)
660            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
661               &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
662            zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
663            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
664            zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
665            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
666            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
667         END_2D
668
669         DO_2D( ji0, ji0, 0, 0 )
670            !                                !  Readjust moments remaining in the box.
671            zbt  =         zbet(ji,jj-1)
672            zbt1 = ( 1.0 - zbet(ji,jj-1) )
673            !
674            zpsm  = psm (ji,jj,jl) ! optimization
675            zps0  = ps0 (ji,jj,jl)
676            zpsx  = psx (ji,jj,jl)
677            zpsxx = psxx(ji,jj,jl)
678            zpsy  = psy (ji,jj,jl)
679            zpsyy = psyy(ji,jj,jl)
680            zpsxy = psxy(ji,jj,jl)
681            !
682            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
683            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
684            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
685            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
686            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) )
687            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
688            zpsxy = zalg1q(ji,jj-1) * zpsxy
689
690            !   Put the temporary moments into appropriate neighboring boxes.
691            !                                !   Flux from j to j+1 IF v GT 0.
692            zbt   =       zbet(ji,jj-1)
693            zbt1  = 1.0 - zbet(ji,jj-1)
694            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm
695            zalf  = zbt * zfm(ji,jj-1) / zpsm
696            zalf1 = 1.0 - zalf
697            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
698            !
699            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
700            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  &
701               &             + zbt1 * zpsy
702            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           &
703               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
704               &             + zbt1 * zpsyy
705            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             &
706               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  &
707               &             + zbt1 * zpsxy
708            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx
709            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
710
711            !                                !  Flux from j+1 to j IF v LT 0.
712            zbt   =       zbet(ji,jj)
713            zbt1  = 1.0 - zbet(ji,jj)
714            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
715            zalf  = zbt1 * zfm(ji,jj) / zpsm
716            zalf1 = 1.0 - zalf
717            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
718            !
719            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) )
720            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
721            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
722               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     &
723               &                         + ( zalf1 - zalf ) * ztemp ) )
724            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
725               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
726            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) )
727            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
728            !
729            psm (ji,jj,jl) = zpsm ! optimization
730            ps0 (ji,jj,jl) = zps0
731            psx (ji,jj,jl) = zpsx
732            psxx(ji,jj,jl) = zpsxx
733            psy (ji,jj,jl) = zpsy
734            psyy(ji,jj,jl) = zpsyy
735            psxy(ji,jj,jl) = zpsxy
736         END_2D
737         !
738      END DO
739      !
740   END SUBROUTINE adv_y
741
742
743   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
744      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
745      !!-------------------------------------------------------------------
746      !!                  ***  ROUTINE Hbig  ***
747      !!
748      !! ** Purpose : Thickness correction in case advection scheme creates
749      !!              abnormally tick ice or snow
750      !!
751      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
752      !!                 (before advection) and reduce it by adapting ice concentration
753      !!              2- check whether snow thickness is larger than the surrounding 9-points
754      !!                 (before advection) and reduce it by sending the excess in the ocean
755      !!
756      !! ** input   : Max thickness of the surrounding 9-points
757      !!-------------------------------------------------------------------
758      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step
759      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts
760      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max
761      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max
762      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
763      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
764      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
765      !
766      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
767      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
768      !!-------------------------------------------------------------------
769      !
770      z1_dt = 1._wp / pdt
771      !
772      DO jl = 1, jpl
773         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
774            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
775               !
776               !                               ! -- check h_ip -- !
777               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
778               IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
779                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
780                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
781                     pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
782                  ENDIF
783               ENDIF
784               !
785               !                               ! -- check h_i -- !
786               ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
787               zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
788               IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
789                  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)
790               ENDIF
791               !
792               !                               ! -- check h_s -- !
793               ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
794               zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
795               IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
796                  zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
797                  !
798                  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
799                  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
800                  !
801                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
802                  pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
803               ENDIF
804               !
805               !                               ! -- check s_i -- !
806               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
807               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
808               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
809                  zfra = psi_max(ji,jj,jl) / zsi
810                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
811                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
812               ENDIF
813               !
814            ENDIF
815         END_2D
816      END DO
817      !
818      !                                           ! -- check e_i/v_i -- !
819      DO jl = 1, jpl
820         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
821            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
822               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
823               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
824               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
825                  zfra = pei_max(ji,jj,jk,jl) / zei
826                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
827                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
828               ENDIF
829            ENDIF
830         END_3D
831      END DO
832      !                                           ! -- check e_s/v_s -- !
833      DO jl = 1, jpl
834         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
835            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
836               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
837               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
838               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
839                  zfra = pes_max(ji,jj,jk,jl) / zes
840                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
841                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
842               ENDIF
843            ENDIF
844         END_3D
845      END DO
846      !
847   END SUBROUTINE Hbig
848
849
850   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
851      !!-------------------------------------------------------------------
852      !!                  ***  ROUTINE Hsnow  ***
853      !!
854      !! ** Purpose : 1- Check snow load after advection
855      !!              2- Correct pond concentration to avoid a_ip > a_i
856      !!
857      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
858      !!              then put the snow excess in the ocean
859      !!
860      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards
861      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
862      !!              make the snow very thick (if concentration decreases drastically)
863      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
864      !!-------------------------------------------------------------------
865      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
866      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
867      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
868      !
869      INTEGER  ::   ji, jj, jl   ! dummy loop indices
870      REAL(wp) ::   z1_dt, zvs_excess, zfra
871      !!-------------------------------------------------------------------
872      !
873      z1_dt = 1._wp / pdt
874      !
875      ! -- check snow load -- !
876      DO jl = 1, jpl
877         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
878            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
879               !
880               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
881               !
882               IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
883                  ! put snow excess in the ocean
884                  zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
885                  wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
886                  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
887                  ! correct snow volume and heat content
888                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
889                  pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
890               ENDIF
891               !
892            ENDIF
893         END_2D
894      END DO
895      !
896      !-- correct pond concentration to avoid a_ip > a_i -- !
897      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
898      !
899   END SUBROUTINE Hsnow
900
901
902   SUBROUTINE adv_pra_init
903      !!-------------------------------------------------------------------
904      !!                  ***  ROUTINE adv_pra_init  ***
905      !!
906      !! ** Purpose :   allocate and initialize arrays for Prather advection
907      !!-------------------------------------------------------------------
908      INTEGER ::   ierr
909      !!-------------------------------------------------------------------
910      !
911      !                             !* allocate prather fields
912      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
913         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
914         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
915         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
916         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
917         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
918         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
919         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   &
920         !
921         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
922         &      syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl)                             , &
923         !
924         &      sxe  (jpi,jpj,nlay_i,jpl) , sye  (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
925         &      syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl)                             , &
926         &      STAT = ierr )
927      !
928      CALL mpp_sum( 'icedyn_adv_pra', ierr )
929      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
930      !
931      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
932      !
933   END SUBROUTINE adv_pra_init
934
935
936   SUBROUTINE adv_pra_rst( cdrw, kt )
937      !!---------------------------------------------------------------------
938      !!                   ***  ROUTINE adv_pra_rst  ***
939      !!
940      !! ** Purpose :   Read or write file in restart file
941      !!
942      !! ** Method  :   use of IOM library
943      !!----------------------------------------------------------------------
944      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
945      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
946      !
947      INTEGER ::   jk, jl   ! dummy loop indices
948      INTEGER ::   iter     ! local integer
949      INTEGER ::   id1      ! local integer
950      CHARACTER(len=25) ::   znam
951      CHARACTER(len=2)  ::   zchar, zchar1
952      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
953      !!----------------------------------------------------------------------
954      !
955      !                                      !==========================!
956      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
957         !                                   !==========================!
958         !
959         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
960         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
961         ENDIF
962         !
963         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
964            !
965            !                                                        ! ice thickness
966            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp )
967            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp )
968            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice )
969            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice )
970            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice )
971            !                                                        ! snow thickness
972            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp )
973            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp )
974            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  )
975            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  )
976            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  )
977            !                                                        ! ice concentration
978            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp )
979            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp )
980            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   )
981            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   )
982            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   )
983            !                                                        ! ice salinity
984            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp )
985            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp )
986            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal )
987            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal )
988            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal )
989            !                                                        ! ice age
990            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp )
991            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp )
992            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage )
993            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage )
994            CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage )
995            !                                                        ! snow layers heat content
996            DO jk = 1, nlay_s
997               WRITE(zchar1,'(I2.2)') jk
998               znam = 'sxc0'//'_l'//zchar1
999               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
1000               znam = 'syc0'//'_l'//zchar1
1001               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
1002               znam = 'sxxc0'//'_l'//zchar1
1003               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
1004               znam = 'syyc0'//'_l'//zchar1
1005               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
1006               znam = 'sxyc0'//'_l'//zchar1
1007               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
1008            END DO
1009            !                                                        ! ice layers heat content
1010            DO jk = 1, nlay_i
1011               WRITE(zchar1,'(I2.2)') jk
1012               znam = 'sxe'//'_l'//zchar1
1013               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
1014               znam = 'sye'//'_l'//zchar1
1015               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:)
1016               znam = 'sxxe'//'_l'//zchar1
1017               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
1018               znam = 'syye'//'_l'//zchar1
1019               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
1020               znam = 'sxye'//'_l'//zchar1
1021               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
1022            END DO
1023            !
1024            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction
1025               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
1026                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
1027                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
1028                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
1029                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
1030                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
1031                  !                                                     ! melt pond volume
1032                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
1033                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
1034                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
1035                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
1036                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
1037               ELSE
1038                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction
1039                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume
1040               ENDIF
1041                  !
1042               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
1043                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
1044                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
1045                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
1046                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
1047                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
1048                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl )
1049                  ELSE
1050                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume
1051                  ENDIF
1052               ENDIF
1053            ENDIF
1054            !
1055         ELSE                                   !**  start rheology from rest  **!
1056            !
1057            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
1058            !
1059            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
1060            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
1061            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
1062            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
1063            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
1064            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
1065            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
1066            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
1067               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
1068               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
1069               IF ( ln_pnd_lids ) THEN
1070                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
1071               ENDIF
1072            ENDIF
1073         ENDIF
1074         !
1075         !                                   !=====================================!
1076      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
1077         !                                   !=====================================!
1078         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
1079         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1080         !
1081         !
1082         ! In case Prather scheme is used for advection, write second order moments
1083         ! ------------------------------------------------------------------------
1084         !
1085         !                                                           ! ice thickness
1086         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
1087         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
1088         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
1089         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
1090         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
1091         !                                                           ! snow thickness
1092         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
1093         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
1094         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
1095         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
1096         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
1097         !                                                           ! ice concentration
1098         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
1099         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
1100         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
1101         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
1102         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
1103         !                                                           ! ice salinity
1104         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
1105         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
1106         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
1107         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
1108         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
1109         !                                                           ! ice age
1110         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
1111         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
1112         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
1113         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
1114         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
1115         !                                                           ! snow layers heat content
1116         DO jk = 1, nlay_s
1117            WRITE(zchar1,'(I2.2)') jk
1118            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)
1119            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1120            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)
1121            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1122            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)
1123            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1124            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)
1125            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1126            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)
1127            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1128         END DO
1129         !                                                           ! ice layers heat content
1130         DO jk = 1, nlay_i
1131            WRITE(zchar1,'(I2.2)') jk
1132            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)
1133            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1134            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)
1135            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1136            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)
1137            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1138            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)
1139            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1140            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)
1141            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1142         END DO
1143         !
1144         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction
1145            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
1146            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
1147            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
1148            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
1149            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
1150            !                                                        ! melt pond volume
1151            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
1152            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
1153            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
1154            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
1155            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
1156            !
1157            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
1158               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
1159               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
1160               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
1161               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
1162               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
1163            ENDIF
1164         ENDIF
1165         !
1166      ENDIF
1167      !
1168   END SUBROUTINE adv_pra_rst
1169
1170   SUBROUTINE icemax3D( pice , pmax )
1171      !!---------------------------------------------------------------------
1172      !!                   ***  ROUTINE icemax3D ***
1173      !! ** Purpose :  compute the max of the 9 points around
1174      !!----------------------------------------------------------------------
1175      REAL(wp), DIMENSION(:,:,:), INTENT(in ) ::   pice   ! input
1176      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pmax   ! output
1177      !
1178      REAL(wp), DIMENSION(Nis0:Nie0) ::   zmax1, zmax2
1179      REAL(wp)                       ::   zmax3
1180      INTEGER  ::   ji, jj, jl   ! dummy loop indices
1181      !!----------------------------------------------------------------------
1182      ! basic version: get the max of epsi20 + 9 neighbours
1183!!$      DO jl = 1, jpl
1184!!$         DO_2D( 0, 0, 0, 0 )
1185!!$            pmax(ji,jj,jl) = MAX( epsi20, pice(ji-1,jj-1,jl), pice(ji,jj-1,jl), pice(ji+1,jj-1,jl),   &
1186!!$               &                          pice(ji-1,jj  ,jl), pice(ji,jj  ,jl), pice(ji+1,jj  ,jl),   &
1187!!$               &                          pice(ji-1,jj+1,jl), pice(ji,jj+1,jl), pice(ji+1,jj+1,jl) )
1188!!$         END_2D
1189!!$      END DO
1190      ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
1191      DO jl = 1, jpl
1192         DO ji = Nis0, Nie0
1193            zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jl), pice(ji-1,Njs0-1,jl), pice(ji+1,Njs0-1,jl) )
1194            zmax2(ji) = MAX( epsi20, pice(ji,Njs0  ,jl), pice(ji-1,Njs0  ,jl), pice(ji+1,Njs0  ,jl) )
1195         END DO
1196         DO_2D( 0, 0, 0, 0 )
1197            zmax3 = MAX( epsi20, pice(ji,jj+1,jl), pice(ji-1,jj+1,jl), pice(ji+1,jj+1,jl) )
1198            pmax(ji,jj,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
1199            zmax1(ji) = zmax2(ji)
1200            zmax2(ji) = zmax3
1201         END_2D
1202      END DO
1203   END SUBROUTINE icemax3D
1204
1205   SUBROUTINE icemax4D( pice , pmax )
1206      !!---------------------------------------------------------------------
1207      !!                   ***  ROUTINE icemax4D ***
1208      !! ** Purpose :  compute the max of the 9 points around
1209      !!----------------------------------------------------------------------
1210      REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) ::   pice   ! input
1211      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) ::   pmax   ! output
1212      !
1213      REAL(wp), DIMENSION(Nis0:Nie0) ::   zmax1, zmax2
1214      REAL(wp)                       ::   zmax3
1215      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
1216      !!----------------------------------------------------------------------
1217      jlay = SIZE( pice , 3 )   ! size of input arrays
1218      ! basic version: get the max of epsi20 + 9 neighbours
1219!!$      DO jl = 1, jpl
1220!!$         DO jk = 1, jlay
1221!!$            DO_2D( 0, 0, 0, 0 )
1222!!$               pmax(ji,jj,jk,jl) = MAX( epsi20, pice(ji-1,jj-1,jk,jl), pice(ji,jj-1,jk,jl), pice(ji+1,jj-1,jk,jl),   &
1223!!$                  &                             pice(ji-1,jj  ,jk,jl), pice(ji,jj  ,jk,jl), pice(ji+1,jj  ,jk,jl),   &
1224!!$                  &                             pice(ji-1,jj+1,jk,jl), pice(ji,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
1225!!$            END_2D
1226!!$         END DO
1227!!$      END DO
1228      ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
1229      DO jl = 1, jpl
1230         DO jk = 1, jlay
1231            DO ji = Nis0, Nie0
1232               zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jk,jl), pice(ji-1,Njs0-1,jk,jl), pice(ji+1,Njs0-1,jk,jl) )
1233               zmax2(ji) = MAX( epsi20, pice(ji,Njs0  ,jk,jl), pice(ji-1,Njs0  ,jk,jl), pice(ji+1,Njs0  ,jk,jl) )
1234            END DO
1235            DO_2D( 0, 0, 0, 0 )
1236               zmax3 = MAX( epsi20, pice(ji,jj+1,jk,jl), pice(ji-1,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
1237               pmax(ji,jj,jk,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
1238               zmax1(ji) = zmax2(ji)
1239               zmax2(ji) = zmax3
1240            END_2D
1241         END DO
1242      END DO
1243   END SUBROUTINE icemax4D
1244
1245#else
1246   !!----------------------------------------------------------------------
1247   !!   Default option            Dummy module        NO SI3 sea-ice model
1248   !!----------------------------------------------------------------------
1249#endif
1250
1251   !!======================================================================
1252END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.