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

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210
  • Property svn:keywords set to Id
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_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_oce         ! ice variables
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 2.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         !i bug   DO ji = 1, jpim1
117         !i    DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
118         DO ji = 1, jpi
119            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
120            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
121            zalfq        =  zalf * zalf
122            zalf1        =  1.0 - zalf
123            zalf1q       =  zalf1 * zalf1
124            zfm (ji,jj)  =  zalf  * psm(ji,jj)
125            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) )
126            zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
127            zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj)
128            zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) )
129            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
130            zfyy(ji,jj)  =  zalf  * psyy(ji,jj)
131
132            !  Readjust moments remaining in the box.
133            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
134            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
135            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
136            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
137            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
138            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
139            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
140         END DO
141      END DO
142
143      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
144         !i    DO jj = 1, fs_jpjm1                   !  Flux from i+1 to i when u LT 0.
145         DO ji = 1, fs_jpim1
146            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
147            zalg  (ji,jj) = zalf
148            zalfq         = zalf * zalf
149            zalf1         = 1.0 - zalf
150            zalg1 (ji,jj) = zalf1
151            zalf1q        = zalf1 * zalf1
152            zalg1q(ji,jj) = zalf1q
153            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj)
154            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
155            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
156            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj)
157            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) )
158            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)
159            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj)
160         END DO
161      END DO
162
163      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
164         DO ji = fs_2, jpi
165            zbt  =       zbet(ji-1,jj)
166            zbt1 = 1.0 - zbet(ji-1,jj)
167            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
168            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
169            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
170            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
171            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
172            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
173            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
174         END DO
175      END DO
176
177      !   Put the temporary moments into appropriate neighboring boxes.   
178      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
179         DO ji = fs_2, fs_jpim1
180            zbt  =       zbet(ji-1,jj)
181            zbt1 = 1.0 - zbet(ji-1,jj)
182            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
183            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
184            zalf1       = 1.0 - zalf
185            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
186            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj)
187            psx(ji,jj)  = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
188            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)   &
189               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
190               &        + zbt1 * psxx(ji,jj)
191            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
192               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
193               &         + zbt1 * psxy(ji,jj)
194            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
195            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
196         END DO
197      END DO
198
199      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
200         DO ji = fs_2, fs_jpim1
201            zbt  =       zbet(ji,jj)
202            zbt1 = 1.0 - zbet(ji,jj)
203            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
204            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
205            zalf1       = 1.0 - zalf
206            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
207            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
208            psx(ji,jj)  = zbt  * psx(ji,jj)   &
209               &        + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
210            psxx(ji,jj) = zbt  * psxx(ji,jj)   &
211               &        + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
212               &                 + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
213            psxy(ji,jj) = zbt  * psxy(ji,jj)   &
214               &        + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
215               &                 + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
216            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
217            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
218         END DO
219      END DO
220
221      !-- Lateral boundary conditions
222      CALL lbc_lnk( psm , 'T', 1. )
223      CALL lbc_lnk( ps0 , 'T', 1. )
224      CALL lbc_lnk( psx , 'T', 1. )
225      CALL lbc_lnk( psxx, 'T', 1. )
226      CALL lbc_lnk( psy , 'T', 1. )
227      CALL lbc_lnk( psyy, 'T', 1. )
228      CALL lbc_lnk( psxy, 'T', 1. )
229
230      IF(ln_ctl) THEN
231         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
232         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
233         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
234         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
235      ENDIF
236
237   END SUBROUTINE lim_adv_x
238
239
240   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
241      &                  psx, psxx, psy , psyy, psxy )
242      !!---------------------------------------------------------------------
243      !!                **  routine lim_adv_y  **
244      !!           
245      !! ** purpose :   Computes and adds the advection trend to sea-ice
246      !!      variable on y axis
247      !!
248      !! ** method  :   Uses Prather second order scheme that advects tracers
249      !!      but also their quadratic forms. The method preserves tracer
250      !!      structures by conserving second order moments.
251      !!
252      !! History :
253      !!   1.0  !  00-01 (LIM)
254      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
255      !!   2.0  !  03-10 (C. Ethe) F90, module
256      !!        !  03-12 (R. Hordoir, G. Madec) mpp
257      !!---------------------------------------------------------------------
258      !! * Arguments
259      REAL(wp),                     INTENT(in)  :: &
260         pdf,        &  ! ???
261         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y
262      !              ! = 0. : lim_adv_x is called after  lim_adv_y
263      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: &
264         pvt            ! j-direction ice velocity at ocean V-point (m/s)
265      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: &
266         psm , ps0 , psx , psy,   &
267         psxx, psyy, psxy
268
269      !! * Local Variables
270      INTEGER  ::   ji, jj                ! dummy loop indices
271      REAL(wp) ::   &
272         zrdt, zslpmax, zin0, ztemp,  &   ! temporary scalars
273         zs1max, zs1new, zs2new,      &   !    "         "
274         zalf, zalfq, zalf1, zalf1q,  &   !    "         "
275         zbt , zbt1                       !
276      REAL(wp), DIMENSION(jpi,jpj)  :: &
277         zf0 , zfx , zfy ,             &  ! temporary workspace
278         zfxx, zfyy, zfxy,             &  !    "           "
279         zfm , zbet,                   &  !    "           "
280         zalg, zalg1, zalg1q              !    "           "
281      !---------------------------------------------------------------------
282
283      ! Limitation of moments.
284
285      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
286
287      DO jj = 1, jpj
288         DO ji = 1, jpi
289            zslpmax = MAX( rzero, ps0(ji,jj) )
290            zs1max  = 1.5 * zslpmax
291            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
292            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
293               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
294            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
295            ps0 (ji,jj) = zslpmax 
296            psx (ji,jj)  = psx (ji,jj) * zin0
297            psxx(ji,jj)  = psxx(ji,jj) * zin0
298            psy (ji,jj) = zs1new * zin0
299            psyy(ji,jj) = zs2new * zin0
300            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
301         END DO
302      END DO
303
304      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
305      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
306
307      !  Calculate fluxes and moments between boxes j<-->j+1             
308      !!bug  DO jj = 2, jpjm1
309      DO jj = 1, jpj
310         DO ji = 1, jpi
311            !!bug     DO ji = 1, jpim1
312            !  Flux from j to j+1 WHEN v GT 0   
313            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
314            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
315            zalfq        =  zalf * zalf
316            zalf1        =  1.0 - zalf
317            zalf1q       =  zalf1 * zalf1
318            zfm (ji,jj)  =  zalf  * psm(ji,jj)
319            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
320            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
321            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
322            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
323            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
324            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
325
326            !  Readjust moments remaining in the box.
327            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
328            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
329            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
330            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
331            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
332            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
333            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
334         END DO
335      END DO
336
337      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
338         DO ji = 1, jpi
339            !i     DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
340            !i        DO ji = 2, jpim1
341            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
342            zalg  (ji,jj) = zalf
343            zalfq         = zalf * zalf
344            zalf1         = 1.0 - zalf
345            zalg1 (ji,jj) = zalf1
346            zalf1q        = zalf1 * zalf1
347            zalg1q(ji,jj) = zalf1q
348            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1)
349            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
350            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
351            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1)
352            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) )
353            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)
354            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1)
355         END DO
356      END DO
357
358      !  Readjust moments remaining in the box.
359      DO jj = 2, jpj
360         DO ji = 1, jpi
361            zbt  =         zbet(ji,jj-1)
362            zbt1 = ( 1.0 - zbet(ji,jj-1) )
363            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
364            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
365            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
366            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj)
367            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
368            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
369            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
370         END DO
371      END DO
372
373      !   Put the temporary moments into appropriate neighboring boxes.   
374      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
375         DO ji = 1, jpi
376            zbt  =         zbet(ji,jj-1)
377            zbt1 = ( 1.0 - zbet(ji,jj-1) )
378            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
379            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
380            zalf1       = 1.0 - zalf
381            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
382            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj)
383
384            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
385               &        + zbt1 * psy(ji,jj) 
386
387            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
388               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
389               &        + zbt1 * psyy(ji,jj)
390
391            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              &
392               + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   &
393               + zbt1 * psxy(ji,jj)
394            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
395            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
396         END DO
397      END DO
398
399      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
400         DO ji = 1, jpi
401            zbt  =         zbet(ji,jj)
402            zbt1 = ( 1.0 - zbet(ji,jj) )
403            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
404            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
405            zalf1       = 1.0 - zalf
406            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
407            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
408            psy(ji,jj)  = zbt  * psy(ji,jj)  &
409               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
410            psyy(ji,jj) = zbt  * psyy(ji,jj)  &
411               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
412               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
413            psxy(ji,jj) = zbt  * psxy(ji,jj)   &
414               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   &
415               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )
416            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
417            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
418         END DO
419      END DO
420
421      !-- Lateral boundary conditions
422      CALL lbc_lnk( psm , 'T', 1. )
423      CALL lbc_lnk( ps0 , 'T', 1. )
424      CALL lbc_lnk( psx , 'T', 1. )
425      CALL lbc_lnk( psxx, 'T', 1. )
426      CALL lbc_lnk( psy , 'T', 1. )
427      CALL lbc_lnk( psyy, 'T', 1. )
428      CALL lbc_lnk( psxy, 'T', 1. )
429
430      IF(ln_ctl) THEN
431         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
432         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
433         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
434         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
435      ENDIF
436
437   END SUBROUTINE lim_adv_y
438
439#else
440   !!----------------------------------------------------------------------
441   !!   Default option            Dummy module         NO LIM sea-ice model
442   !!----------------------------------------------------------------------
443CONTAINS
444   SUBROUTINE lim_adv_x         ! Empty routine
445   END SUBROUTINE lim_adv_x
446   SUBROUTINE lim_adv_y           ! Empty routine
447   END SUBROUTINE lim_adv_y
448
449#endif
450
451END MODULE limadv
Note: See TracBrowser for help on using the repository browser.