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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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