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 @ 106

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

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