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

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

First optimisation of LIM3, limrhg optimisation induces computation change
File size: 21.8 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   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limadv.F90,v 1.4 2005/03/27 18:34:41 opalod Exp $
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.