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 trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limadv.F90 @ 1529

Last change on this file since 1529 was 1529, checked in by ctlod, 15 years ago

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