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 – NEMO

source: trunk/NEMO/LIM_SRC/limadv.F90 @ 12

Last change on this file since 12 was 3, checked in by opalod, 20 years ago

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