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

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

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icedyn_adv_pra.F90 @ 11573

Last change on this file since 11573 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

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