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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_pra.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

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