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/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90 @ 9499

Last change on this file since 9499 was 9499, checked in by davestorkey, 6 years ago

branches/UKMO/dev_merge_2017_CICE_interface : clear SVN keywords.

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