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.
iceadv_prather.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90 @ 8409

Last change on this file since 8409 was 8409, checked in by clem, 7 years ago

change calls in icestp.F90 for advection
File size: 22.1 KB
Line 
1MODULE iceadv_prather 
2   !!======================================================================
3   !!                       ***  MODULE iceadv_prather   ***
4   !! LIM sea-ice model : 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'                                     LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_adv_x  : advection of sea ice on x axis
15   !!   ice_adv_y  : advection of sea ice on y axis
16   !!----------------------------------------------------------------------
17   USE dom_oce          ! ocean domain
18   USE ice              ! LIM-3 variables
19   !
20   USE lbclnk           ! lateral boundary condition - MPP exchanges
21   USE in_out_manager   ! I/O manager
22   USE prtctl           ! Print control
23   USE lib_mpp          ! MPP library
24   USE lib_fortran      ! to use key_nosignedzero
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_adv_x   ! called by iceadv
30   PUBLIC   ice_adv_y   ! called by iceadv
31
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
36   !! $Id: iceadv.F90 6746 2016-06-27 17:20:57Z clem $
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE ice_adv_x( pdf, put , pcrh, psm , ps0 ,   &
42      &                  psx, psxx, psy , psyy, psxy )
43      !!---------------------------------------------------------------------
44      !!                **  routine ice_adv_x  **
45      !! 
46      !! ** purpose :   Computes and adds the advection trend to sea-ice
47      !!              variable on i-axis
48      !!
49      !! ** method  :   Uses Prather second order scheme that advects tracers
50      !!              but also theirquadratic forms. The method preserves
51      !!              tracer structures by conserving second order moments.
52      !!
53      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
54      !!--------------------------------------------------------------------
55      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
56      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)
57      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
60      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
62      !!
63      INTEGER  ::   ji, jj                               ! dummy loop indices
64      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
65      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
66      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
67      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
68      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
69      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
70      !---------------------------------------------------------------------
71
72      ! Limitation of moments.                                           
73
74      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
75
76      DO jj = 1, jpj
77         DO ji = 1, jpi
78            zslpmax = MAX( 0._wp, ps0(ji,jj) )
79            zs1max  = 1.5 * zslpmax
80            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
81            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
82               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
83            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
84
85            ps0 (ji,jj) = zslpmax 
86            psx (ji,jj) = zs1new      * rswitch
87            psxx(ji,jj) = zs2new      * rswitch
88            psy (ji,jj) = psy (ji,jj) * rswitch
89            psyy(ji,jj) = psyy(ji,jj) * rswitch
90            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
91         END DO
92      END DO
93
94      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
95      psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
96
97      !  Calculate fluxes and moments between boxes i<-->i+1             
98      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
99         DO ji = 1, jpi
100            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
101            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
102            zalfq        =  zalf * zalf
103            zalf1        =  1.0 - zalf
104            zalf1q       =  zalf1 * zalf1
105            !
106            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
107            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
108            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
109            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
110            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
111            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
112            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
113
114            !  Readjust moments remaining in the box.
115            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
116            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
117            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
118            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
119            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
120            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
121            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
122         END DO
123      END DO
124
125      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
126         DO ji = 1, fs_jpim1
127            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
128            zalg  (ji,jj) = zalf
129            zalfq         = zalf * zalf
130            zalf1         = 1.0 - zalf
131            zalg1 (ji,jj) = zalf1
132            zalf1q        = zalf1 * zalf1
133            zalg1q(ji,jj) = zalf1q
134            !
135            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
136            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
137            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
138            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
139            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
140            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
141            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
142         END DO
143      END DO
144
145      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
146         DO ji = fs_2, fs_jpim1
147            zbt  =       zbet(ji-1,jj)
148            zbt1 = 1.0 - zbet(ji-1,jj)
149            !
150            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
151            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
152            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
153            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
154            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
155            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
156            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
157         END DO
158      END DO
159
160      !   Put the temporary moments into appropriate neighboring boxes.   
161      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
162         DO ji = fs_2, fs_jpim1
163            zbt  =       zbet(ji-1,jj)
164            zbt1 = 1.0 - zbet(ji-1,jj)
165            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
166            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
167            zalf1       = 1.0 - zalf
168            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
169            !
170            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
171            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
172            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
173               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
174               &                                                + zbt1 * psxx(ji,jj)
175            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
176               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
177               &                                                + zbt1 * psxy(ji,jj)
178            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
179            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
180         END DO
181      END DO
182
183      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
184         DO ji = fs_2, fs_jpim1
185            zbt  =       zbet(ji,jj)
186            zbt1 = 1.0 - zbet(ji,jj)
187            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
188            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
189            zalf1       = 1.0 - zalf
190            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
191            !
192            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
193            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
194            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
195               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
196               &                                      + ( zalf1 - zalf ) * ztemp ) )
197            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
198               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
199            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
200            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
201         END DO
202      END DO
203
204      !-- Lateral boundary conditions
205      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
206         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
207         &              , psxx, 'T',  1., psyy, 'T',  1.   &
208         &              , psxy, 'T',  1. )
209
210      IF(ln_ctl) THEN
211         CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
212         CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
213         CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
214         CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_x: psxy :')
215      ENDIF
216      !
217   END SUBROUTINE ice_adv_x
218
219
220   SUBROUTINE ice_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
221      &                  psx, psxx, psy , psyy, psxy )
222      !!---------------------------------------------------------------------
223      !!                **  routine ice_adv_y  **
224      !!           
225      !! ** purpose :   Computes and adds the advection trend to sea-ice
226      !!              variable on y axis
227      !!
228      !! ** method  :   Uses Prather second order scheme that advects tracers
229      !!              but also their quadratic forms. The method preserves
230      !!              tracer structures by conserving second order moments.
231      !!
232      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
233      !!---------------------------------------------------------------------
234      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
235      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)
236      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
237      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
238      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
239      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
240      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
241      !!
242      INTEGER  ::   ji, jj                               ! dummy loop indices
243      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
244      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
245      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
246      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
247      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
248      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
249      !---------------------------------------------------------------------
250
251      ! Limitation of moments.
252
253      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
254
255      DO jj = 1, jpj
256         DO ji = 1, jpi
257            zslpmax = MAX( 0._wp, ps0(ji,jj) )
258            zs1max  = 1.5 * zslpmax
259            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
260            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
261               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
262            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
263            !
264            ps0 (ji,jj) = zslpmax 
265            psx (ji,jj) = psx (ji,jj) * rswitch
266            psxx(ji,jj) = psxx(ji,jj) * rswitch
267            psy (ji,jj) = zs1new * rswitch
268            psyy(ji,jj) = zs2new * rswitch
269            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
270         END DO
271      END DO
272
273      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
274      psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
275
276      !  Calculate fluxes and moments between boxes j<-->j+1             
277      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
278         DO ji = 1, jpi
279            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
280            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
281            zalfq        =  zalf * zalf
282            zalf1        =  1.0 - zalf
283            zalf1q       =  zalf1 * zalf1
284            !
285            zfm (ji,jj)  =  zalf  * psm(ji,jj)
286            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
287            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
288            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
289            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
290            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
291            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
292            !
293            !  Readjust moments remaining in the box.
294            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
295            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
296            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
297            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
298            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
299            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
300            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
301         END DO
302      END DO
303      !
304      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
305         DO ji = 1, jpi
306            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
307            zalg  (ji,jj) = zalf
308            zalfq         = zalf * zalf
309            zalf1         = 1.0 - zalf
310            zalg1 (ji,jj) = zalf1
311            zalf1q        = zalf1 * zalf1
312            zalg1q(ji,jj) = zalf1q
313            !
314            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
315            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
316            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
317            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
318            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
319            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
320            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
321         END DO
322      END DO
323
324      !  Readjust moments remaining in the box.
325      DO jj = 2, jpj
326         DO ji = 1, jpi
327            zbt  =         zbet(ji,jj-1)
328            zbt1 = ( 1.0 - zbet(ji,jj-1) )
329            !
330            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
331            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
332            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
333            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
334            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
335            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
336            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
337         END DO
338      END DO
339
340      !   Put the temporary moments into appropriate neighboring boxes.   
341      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
342         DO ji = 1, jpi
343            zbt  =         zbet(ji,jj-1)
344            zbt1 = ( 1.0 - zbet(ji,jj-1) )
345            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
346            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
347            zalf1       = 1.0 - zalf
348            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
349            !
350            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
351            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
352               &                                               + zbt1 * psy(ji,jj) 
353            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
354               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
355               &                                               + zbt1 * psyy(ji,jj)
356            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
357               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
358               &                                                + zbt1 * psxy(ji,jj)
359            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
360            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
361         END DO
362      END DO
363
364      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
365         DO ji = 1, jpi
366            zbt  =         zbet(ji,jj)
367            zbt1 = ( 1.0 - zbet(ji,jj) )
368            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
369            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
370            zalf1       = 1.0 - zalf
371            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
372            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
373            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
374            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
375               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
376               &                                         + ( zalf1 - zalf ) * ztemp )                                )
377            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
378               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
379            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
380            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
381         END DO
382      END DO
383
384      !-- Lateral boundary conditions
385      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
386         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
387         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
388         &              , psxy, 'T',  1. )
389
390      IF(ln_ctl) THEN
391         CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
392         CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
393         CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
394         CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_y: psxy :')
395      ENDIF
396      !
397   END SUBROUTINE ice_adv_y
398
399#else
400   !!----------------------------------------------------------------------
401   !!   Default option            Dummy module         NO LIM sea-ice model
402   !!----------------------------------------------------------------------
403#endif
404
405   !!======================================================================
406END MODULE iceadv_prather
Note: See TracBrowser for help on using the repository browser.