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

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

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90 @ 14288

Last change on this file since 14288 was 13877, checked in by mocavero, 4 years ago

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