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 branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90 @ 8884

Last change on this file since 8884 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864
File size: 53.8 KB
Line 
1MODULE icedyn_adv_pra 
2   !!======================================================================
3   !!                       ***  MODULE icedyn_adv_pra   ***
4   !!   sea-ice : advection => Prather scheme
5   !!======================================================================
6   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b.c.
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!--------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                       ESIM sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_dyn_adv_pra : advection of sea ice using Prather scheme
15   !!   adv_x, adv_y    : Prather scheme applied in i- and j-direction, resp.
16   !!   adv_pra_init    : initialisation of the Prather scheme
17   !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero
18   !!----------------------------------------------------------------------
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   !
23   USE in_out_manager ! I/O manager
24   USE iom            ! I/O manager library
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28   USE prtctl         ! Print control
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     ! lead fraction
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow thermal content
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
52   !! $Id: icedyn_adv_pra.F90 6746 2016-06-27 17:20:57Z clem $
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  &
58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, 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(inout) ::   pato_i     ! open water area
74      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
75      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
76      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
77      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
78      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
79      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
80      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
81      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
82      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
83      !
84      INTEGER  ::   jk, jl, jt              ! dummy loop indices
85      INTEGER  ::   initad                  ! number of sub-timestep for the advection
86      REAL(wp) ::   zcfl , zusnit           !   -      -
87      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea
88      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw
89      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
90      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
91      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
92      !!----------------------------------------------------------------------
93      !
94      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
95      !
96      ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           &
97         &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   &
98         &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
99         &      z0ei (jpi,jpj,nlay_i,jpl)                                                           )
100      !
101      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !       
102      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
103      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
104      IF( lk_mpp )   CALL mpp_max( zcfl )
105     
106      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
107      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
108      ENDIF
109     
110      zarea(:,:) = e1e2t(:,:)
111      !-------------------------
112      ! transported fields                                       
113      !-------------------------
114      z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area
115      DO jl = 1, jpl
116         z0snw(:,:,jl) = pv_s (:,:,  jl) * e1e2t(:,:)     ! Snow volume
117         z0ice(:,:,jl) = pv_i (:,:,  jl) * e1e2t(:,:)     ! Ice  volume
118         z0ai (:,:,jl) = pa_i (:,:,  jl) * e1e2t(:,:)     ! Ice area
119         z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content
120         z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content
121         z0es (:,:,jl) = pe_s (:,:,1,jl) * e1e2t(:,:)     ! Snow heat content
122         DO jk = 1, nlay_i
123            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
124         END DO
125         IF ( ln_pnd_H12 ) THEN
126            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
127            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
128         ENDIF
129      END DO
130
131      !                                                    !--------------------------------------------!
132      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
133         !                                                 !--------------------------------------------!
134         DO jt = 1, initad
135            CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
136               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
137            CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
138               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
139            DO jl = 1, jpl
140               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
141                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
142               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
143                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
144               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
145                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
146               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
147                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
148               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
149                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
150               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
151                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
152               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
153                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
154               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
155                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
156               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
157                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
158               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
159                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
160               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
161                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
162               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
163                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
164               DO jk = 1, nlay_i                                                               !--- ice heat contents ---
165                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
166                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
167                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
168                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
169                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
170                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
171               END DO
172               IF ( ln_pnd_H12 ) THEN
173                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
174                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
175                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
176                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
177                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
178                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
179                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
180                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
181               ENDIF
182            END DO
183         END DO
184      !                                                    !--------------------------------------------!
185      ELSE                                                 !== even ice time step:  adv_y then adv_x  ==!
186         !                                                 !--------------------------------------------!
187         DO jt = 1, initad
188            CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
189               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
190            CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
191               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
192            DO jl = 1, jpl
193               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
194                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
195               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
196                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
197               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
198                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
199               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
200                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
201               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
202                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
203               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
204                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
205               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
206                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
207               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
208                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
209               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
210                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
211               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
212                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
213               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
214                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
215               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
216                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
217               DO jk = 1, nlay_i                                                             !--- ice heat contents ---
218                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
219                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
220                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
221                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
222                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
223                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
224               END DO
225               IF ( ln_pnd_H12 ) THEN
226                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
227                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
228                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
229                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
230                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
231                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
232                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
233                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
234               ENDIF
235            END DO
236         END DO
237      ENDIF
238
239      !-------------------------------------------
240      ! Recover the properties from their contents
241      !-------------------------------------------
242      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
243      DO jl = 1, jpl
244         pv_i (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
245         pv_s (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
246         psv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
247         poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
248         pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
249         pe_s (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
250         DO jk = 1, nlay_i
251            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
252         END DO
253         ! MV MP 2016
254         IF ( ln_pnd_H12 ) THEN
255            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
256            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
257         ENDIF
258         ! END MV MP 2016
259      END DO
260      !
261      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei )
262      !
263      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
264      !
265   END SUBROUTINE ice_dyn_adv_pra
266   
267   
268   SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   &
269      &              psx, psxx, psy , psyy, psxy )
270      !!----------------------------------------------------------------------
271      !!                **  routine adv_x  **
272      !! 
273      !! ** purpose :   Computes and adds the advection trend to sea-ice
274      !!                variable on x axis
275      !!----------------------------------------------------------------------
276      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
277      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
279      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
280      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
281      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
282      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
283      !!
284      INTEGER  ::   ji, jj                               ! dummy loop indices
285      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
286      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
287      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
288      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
289      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
290      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
291      !-----------------------------------------------------------------------
292
293      ! Limitation of moments.                                           
294
295      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
296
297      DO jj = 1, jpj
298         DO ji = 1, jpi
299            zslpmax = MAX( 0._wp, ps0(ji,jj) )
300            zs1max  = 1.5 * zslpmax
301            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
302            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
303               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
304            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
305
306            ps0 (ji,jj) = zslpmax 
307            psx (ji,jj) = zs1new      * rswitch
308            psxx(ji,jj) = zs2new      * rswitch
309            psy (ji,jj) = psy (ji,jj) * rswitch
310            psyy(ji,jj) = psyy(ji,jj) * rswitch
311            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
312         END DO
313      END DO
314
315      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
316      psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
317
318      !  Calculate fluxes and moments between boxes i<-->i+1             
319      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
320         DO ji = 1, jpi
321            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
322            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
323            zalfq        =  zalf * zalf
324            zalf1        =  1.0 - zalf
325            zalf1q       =  zalf1 * zalf1
326            !
327            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
328            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
329            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
330            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
331            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
332            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
333            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
334
335            !  Readjust moments remaining in the box.
336            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
337            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
338            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
339            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
340            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
341            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
342            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
343         END DO
344      END DO
345
346      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
347         DO ji = 1, fs_jpim1
348            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
349            zalg  (ji,jj) = zalf
350            zalfq         = zalf * zalf
351            zalf1         = 1.0 - zalf
352            zalg1 (ji,jj) = zalf1
353            zalf1q        = zalf1 * zalf1
354            zalg1q(ji,jj) = zalf1q
355            !
356            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
357            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
358            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
359            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
360            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
361            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
362            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
363         END DO
364      END DO
365
366      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
367         DO ji = fs_2, fs_jpim1
368            zbt  =       zbet(ji-1,jj)
369            zbt1 = 1.0 - zbet(ji-1,jj)
370            !
371            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
372            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
373            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
374            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
375            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
376            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
377            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
378         END DO
379      END DO
380
381      !   Put the temporary moments into appropriate neighboring boxes.   
382      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
383         DO ji = fs_2, fs_jpim1
384            zbt  =       zbet(ji-1,jj)
385            zbt1 = 1.0 - zbet(ji-1,jj)
386            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
387            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
388            zalf1       = 1.0 - zalf
389            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
390            !
391            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
392            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
393            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
394               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
395               &                                                + zbt1 * psxx(ji,jj)
396            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
397               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
398               &                                                + zbt1 * psxy(ji,jj)
399            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
400            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
401         END DO
402      END DO
403
404      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
405         DO ji = fs_2, fs_jpim1
406            zbt  =       zbet(ji,jj)
407            zbt1 = 1.0 - zbet(ji,jj)
408            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
409            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
410            zalf1       = 1.0 - zalf
411            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
412            !
413            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
414            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
415            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
416               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
417               &                                      + ( zalf1 - zalf ) * ztemp ) )
418            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
419               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
420            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
421            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
422         END DO
423      END DO
424
425      !-- Lateral boundary conditions
426      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
427         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
428         &              , psxx, 'T',  1., psyy, 'T',  1.   &
429         &              , psxy, 'T',  1. )
430
431      IF(ln_ctl) THEN
432         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
433         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
434         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
435         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :')
436      ENDIF
437      !
438   END SUBROUTINE adv_x
439
440
441   SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
442      &              psx, psxx, psy , psyy, psxy )
443      !!---------------------------------------------------------------------
444      !!                **  routine adv_y  **
445      !!           
446      !! ** purpose :   Computes and adds the advection trend to sea-ice
447      !!                variable on y axis
448      !!---------------------------------------------------------------------
449      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
450      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
451      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
452      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
453      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
454      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
455      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
456      !!
457      INTEGER  ::   ji, jj                               ! dummy loop indices
458      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
459      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
460      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
461      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
462      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
463      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
464      !---------------------------------------------------------------------
465
466      ! Limitation of moments.
467
468      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
469
470      DO jj = 1, jpj
471         DO ji = 1, jpi
472            zslpmax = MAX( 0._wp, ps0(ji,jj) )
473            zs1max  = 1.5 * zslpmax
474            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
475            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
476               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
477            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
478            !
479            ps0 (ji,jj) = zslpmax 
480            psx (ji,jj) = psx (ji,jj) * rswitch
481            psxx(ji,jj) = psxx(ji,jj) * rswitch
482            psy (ji,jj) = zs1new * rswitch
483            psyy(ji,jj) = zs2new * rswitch
484            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
485         END DO
486      END DO
487
488      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
489      psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
490
491      !  Calculate fluxes and moments between boxes j<-->j+1             
492      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
493         DO ji = 1, jpi
494            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
495            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
496            zalfq        =  zalf * zalf
497            zalf1        =  1.0 - zalf
498            zalf1q       =  zalf1 * zalf1
499            !
500            zfm (ji,jj)  =  zalf  * psm(ji,jj)
501            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
502            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
503            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
504            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
505            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
506            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
507            !
508            !  Readjust moments remaining in the box.
509            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
510            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
511            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
512            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
513            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
514            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
515            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
516         END DO
517      END DO
518      !
519      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
520         DO ji = 1, jpi
521            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
522            zalg  (ji,jj) = zalf
523            zalfq         = zalf * zalf
524            zalf1         = 1.0 - zalf
525            zalg1 (ji,jj) = zalf1
526            zalf1q        = zalf1 * zalf1
527            zalg1q(ji,jj) = zalf1q
528            !
529            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
530            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
531            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
532            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
533            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
534            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
535            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
536         END DO
537      END DO
538
539      !  Readjust moments remaining in the box.
540      DO jj = 2, jpj
541         DO ji = 1, jpi
542            zbt  =         zbet(ji,jj-1)
543            zbt1 = ( 1.0 - zbet(ji,jj-1) )
544            !
545            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
546            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
547            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
548            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
549            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
550            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
551            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
552         END DO
553      END DO
554
555      !   Put the temporary moments into appropriate neighboring boxes.   
556      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
557         DO ji = 1, jpi
558            zbt  =         zbet(ji,jj-1)
559            zbt1 = ( 1.0 - zbet(ji,jj-1) )
560            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
561            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
562            zalf1       = 1.0 - zalf
563            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
564            !
565            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
566            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
567               &                                               + zbt1 * psy(ji,jj) 
568            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
569               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
570               &                                               + zbt1 * psyy(ji,jj)
571            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
572               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
573               &                                                + zbt1 * psxy(ji,jj)
574            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
575            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
576         END DO
577      END DO
578
579      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
580         DO ji = 1, jpi
581            zbt  =         zbet(ji,jj)
582            zbt1 = ( 1.0 - zbet(ji,jj) )
583            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
584            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
585            zalf1       = 1.0 - zalf
586            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
587            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
588            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
589            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
590               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
591               &                                         + ( zalf1 - zalf ) * ztemp )                                )
592            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
593               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
594            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
595            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
596         END DO
597      END DO
598
599      !-- Lateral boundary conditions
600      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
601         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
602         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
603         &              , psxy, 'T',  1. )
604
605      IF(ln_ctl) THEN
606         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
607         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
608         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
609         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :')
610      ENDIF
611      !
612   END SUBROUTINE adv_y
613
614
615   SUBROUTINE adv_pra_init
616      !!-------------------------------------------------------------------
617      !!                  ***  ROUTINE adv_pra_init  ***
618      !!
619      !! ** Purpose :   allocate and initialize arrays for Prather advection
620      !!-------------------------------------------------------------------
621      INTEGER ::   ierr
622      !!-------------------------------------------------------------------
623      !
624      !                             !* allocate prather fields
625      ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   &
626         &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
627         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
628         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
629         &      sxc0 (jpi,jpj,jpl) , syc0 (jpi,jpj,jpl) , sxxc0 (jpi,jpj,jpl) , syyc0 (jpi,jpj,jpl) , sxyc0 (jpi,jpj,jpl) ,   &
630         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
631         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
632         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
633         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
634         &      sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) , &
635         &      syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl)                            , &
636         &      STAT = ierr )
637      !
638      IF( lk_mpp    )   CALL mpp_sum( ierr )
639      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
640      !
641      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
642      !
643   END SUBROUTINE adv_pra_init
644
645
646   SUBROUTINE adv_pra_rst( cdrw, kt )
647      !!---------------------------------------------------------------------
648      !!                   ***  ROUTINE adv_pra_rst  ***
649      !!                     
650      !! ** Purpose :   Read or write RHG file in restart file
651      !!
652      !! ** Method  :   use of IOM library
653      !!----------------------------------------------------------------------
654      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
655      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
656      !
657      INTEGER ::   jk, jl   ! dummy loop indices
658      INTEGER ::   iter     ! local integer
659      INTEGER ::   id1      ! local integer
660      CHARACTER(len=25) ::   znam
661      CHARACTER(len=2)  ::   zchar, zchar1
662      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
663      !!----------------------------------------------------------------------
664      !
665      !                                      !==========================!
666      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
667         !                                   !==========================!
668         !
669         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0
670         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
671         ENDIF
672         !
673         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
674            !
675            !                                                        ! ice thickness
676            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
677            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
678            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
679            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
680            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
681            !                                                        ! snow thickness
682            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
683            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
684            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
685            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
686            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
687            !                                                        ! lead fraction
688            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
689            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
690            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
691            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
692            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
693            !                                                        ! snow thermal content
694            CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   )
695            CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   )
696            CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  )
697            CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  )
698            CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  )
699            !                                                        ! ice salinity
700            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
701            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
702            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
703            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
704            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
705            !                                                        ! ice age
706            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
707            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
708            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
709            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
710            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
711            !                                                        ! open water in sea ice
712            CALL iom_get( numrir, jpdom_autoglo, 'sxopw ', sxopw  )
713            CALL iom_get( numrir, jpdom_autoglo, 'syopw ', syopw  )
714            CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )
715            CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
716            CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
717            !                                                        ! ice layers heat content
718            DO jk = 1, nlay_i 
719               WRITE(zchar1,'(I2.2)') jk
720               znam = 'sxe'//'_il'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
721               znam = 'sye'//'_il'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
722               znam = 'sxxe'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
723               znam = 'syye'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
724               znam = 'sxye'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
725            END DO
726            !
727            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
728               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
729               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
730               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
731               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
732               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
733               !                                                     ! melt pond volume
734               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
735               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
736               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
737               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
738               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
739            ENDIF
740            !
741         ELSE                                   !**  start rheology from rest  **!
742            !
743            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
744            !
745            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
746            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
747            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction
748            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow thermal content
749            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
750            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
751            sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
752            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
753            IF( ln_pnd_H12 ) THEN
754               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
755               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
756            ENDIF
757         ENDIF
758         !
759         !                                   !=====================================!
760      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
761         !                                   !=====================================!
762         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
763         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
764         !
765         !
766         ! In case Prather scheme is used for advection, write second order moments
767         ! ------------------------------------------------------------------------
768         !
769         !                                                           ! ice thickness
770         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
771         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
772         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
773         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
774         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
775         !                                                           ! snow thickness
776         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
777         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
778         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
779         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
780         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
781         !                                                           ! lead fraction
782         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
783         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
784         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
785         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
786         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
787         !                                                           ! snow thermal content
788         CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0   )
789         CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0   )
790         CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0  )
791         CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0  )
792         CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0  )
793         !                                                           ! ice salinity
794         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
795         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
796         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
797         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
798         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
799         !                                                           ! ice age
800         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
801         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
802         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
803         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
804         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
805         !                                                           ! open water in sea ice
806         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ', sxopw  )
807         CALL iom_rstput( iter, nitrst, numriw, 'syopw ', syopw  )
808         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )
809         CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
810         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
811         !                                                           ! ice layers heat content
812         DO jk = 1, nlay_i 
813            WRITE(zchar1,'(I2.2)') jk
814            znam = 'sxe'//'_il'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
815            znam = 'sye'//'_il'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
816            znam = 'sxxe'//'_il'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
817            znam = 'syye'//'_il'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
818            znam = 'sxye'//'_il'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
819         END DO
820         !
821         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
822            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
823            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
824            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
825            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
826            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
827            !                                                        ! melt pond volume
828            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
829            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
830            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
831            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
832            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
833         ENDIF
834         !
835      ENDIF
836      !
837   END SUBROUTINE adv_pra_rst
838
839#else
840   !!----------------------------------------------------------------------
841   !!   Default option            Dummy module        NO ESIM sea-ice model
842   !!----------------------------------------------------------------------
843#endif
844
845   !!======================================================================
846END MODULE icedyn_adv_pra
Note: See TracBrowser for help on using the repository browser.