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.
limadv.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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