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_2.F90 in tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_2 – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_2/limadv_2.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 22.0 KB
Line 
1MODULE limadv_2 
2   !!======================================================================
3   !!                       ***  MODULE limadv_2   ***
4   !! LIM 2.0 sea-ice model : sea-ice advection
5   !!======================================================================
6   !! History :  OPA  ! 2000-01 (LIM)  Original code
7   !!                 ! 2001-05 (G. Madec, R. Hordoir) Doctor norm
8   !!   NEMO     1.0  ! 2003-10 (C. Ethe) F90, module
9   !!             -   ! 2003-12 (R. Hordoir, G. Madec) mpp
10   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
11   !!--------------------------------------------------------------------
12#if defined key_lim2
13   !!----------------------------------------------------------------------
14   !!   'key_lim2'                                    LIM 2.0 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   lim_adv_x_2  : advection of sea ice on x axis
17   !!   lim_adv_y_2  : advection of sea ice on y axis
18   !!----------------------------------------------------------------------
19   USE dom_oce
20   USE dom_ice_2
21   USE ice_2
22   USE lbclnk
23   USE in_out_manager     ! I/O manager
24   USE prtctl             ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   lim_adv_x_2   ! called by lim_trp
30   PUBLIC   lim_adv_y_2   ! called by lim_trp
31
32   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
33   REAL(wp)  ::   rzero  = 0.e0     !    -       -
34   REAL(wp)  ::   rone   = 1.e0     !    -       -
35   
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009)
40   !! $Id: limadv_2.F90 1530 2009-07-23 16:22:51Z ctlod $
41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
47      &                    psx, psxx, psy , psyy, psxy )
48      !!---------------------------------------------------------------------
49      !!                **  routine lim_adv_x_2  **
50      !! 
51      !! ** purpose :   Computes and adds the advection trend to sea-ice
52      !!              variable on i-axis
53      !!
54      !! ** method  :   Uses Prather second order scheme that advects tracers
55      !!              but also theirquadratic forms. The method preserves
56      !!              tracer structures by conserving second order moments.
57      !!
58      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
59      !!--------------------------------------------------------------------
60      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
61      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
62      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
63      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
64      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
65      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
66      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
67      !!
68      INTEGER  ::   ji, jj                               ! dummy loop indices
69      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
70      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
71      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
72      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
73      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
74      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
75      !---------------------------------------------------------------------
76
77      ! Limitation of moments.                                           
78
79      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
80
81      DO jj = 1, jpj
82         DO ji = 1, jpi
83            zslpmax = MAX( rzero, ps0(ji,jj) )
84            zs1max  = 1.5 * zslpmax
85            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
86            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
87               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
88            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
89            !
90            ps0 (ji,jj) = zslpmax 
91            psx (ji,jj) = zs1new      * zin0
92            psxx(ji,jj) = zs2new      * zin0
93            psy (ji,jj) = psy (ji,jj) * zin0
94            psyy(ji,jj) = psyy(ji,jj) * zin0
95            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
96         END DO
97      END DO
98
99      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
100      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
101
102      !  Calculate fluxes and moments between boxes i<-->i+1             
103      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
104         DO ji = 1, jpi
105            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
106            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
107            zalfq        =  zalf * zalf
108            zalf1        =  1.0 - zalf
109            zalf1q       =  zalf1 * zalf1
110            !
111            zfm (ji,jj)  =  zalf  * psm(ji,jj)
112            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) )
113            zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
114            zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj)
115            zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) )
116            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
117            zfyy(ji,jj)  =  zalf  * psyy(ji,jj)
118            !
119            !  Readjust moments remaining in the box.
120            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
121            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
122            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
123            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
124            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
125            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
126            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
127         END DO
128      END DO
129
130      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
131         DO ji = 1, fs_jpim1
132            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
133            zalg  (ji,jj) = zalf
134            zalfq         = zalf * zalf
135            zalf1         = 1.0 - zalf
136            zalg1 (ji,jj) = zalf1
137            zalf1q        = zalf1 * zalf1
138            zalg1q(ji,jj) = zalf1q
139            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj)
140            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
141            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
142            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj)
143            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) )
144            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)
145            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj)
146         END DO
147      END DO
148
149      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
150         DO ji = fs_2, fs_jpim1
151            zbt  =       zbet(ji-1,jj)
152            zbt1 = 1.0 - zbet(ji-1,jj)
153            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
154            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
155            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
156            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
157            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
158            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
159            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
160         END DO
161      END DO
162
163      !   Put the temporary moments into appropriate neighboring boxes.   
164      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
165         DO ji = fs_2, fs_jpim1
166            zbt  =       zbet(ji-1,jj)
167            zbt1 = 1.0 - zbet(ji-1,jj)
168            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
169            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
170            zalf1       = 1.0 - zalf
171            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
172            !
173            ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj)
174            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
175            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
176               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
177               &        + zbt1 * psxx(ji,jj)
178            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
179               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
180               &         + zbt1 * psxy(ji,jj)
181            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
182            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
183         END DO
184      END DO
185
186      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
187         DO ji = fs_2, fs_jpim1
188            zbt  =       zbet(ji,jj)
189            zbt1 = 1.0 - zbet(ji,jj)
190            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
191            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
192            zalf1       = 1.0 - zalf
193            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
194            !
195            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
196            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
197            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)   &
198               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )        &
199               &                                      + ( zalf1 - zalf ) * ztemp )                                 )
200            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)          &
201               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
202            psy(ji,jj)  = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
203            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
204         END DO
205      END DO
206
207      !-- Lateral boundary conditions
208      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
209      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
210      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
211      CALL lbc_lnk( psxy, 'T',  1. )
212
213      IF(ln_ctl) THEN
214         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
215         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
216         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
217         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
218      ENDIF
219      !
220   END SUBROUTINE lim_adv_x_2
221
222
223   SUBROUTINE lim_adv_y_2( pdf, pvt , pcrh, psm , ps0 ,   &
224      &                    psx, psxx, psy , psyy, psxy )
225      !!---------------------------------------------------------------------
226      !!                **  routine lim_adv_y_2  **
227      !!           
228      !! ** purpose :   Computes and adds the advection trend to sea-ice
229      !!              variable on j-axis
230      !!
231      !! ** method  :   Uses Prather second order scheme that advects tracers
232      !!              but also their quadratic forms. The method preserves
233      !!              tracer structures by conserving second order moments.
234      !!
235      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
236      !!---------------------------------------------------------------------
237      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
238      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
239      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
240      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
241      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
242      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
243      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
244      !!
245      INTEGER  ::   ji, jj                               ! dummy loop indices
246      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
247      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
248      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
249      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
250      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
251      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
252      !---------------------------------------------------------------------
253
254      ! Limitation of moments.
255
256      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
257
258      DO jj = 1, jpj
259         DO ji = 1, jpi
260            zslpmax = MAX( rzero, ps0(ji,jj) )
261            zs1max  = 1.5 * zslpmax
262            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
263            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
264               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
265            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
266            !
267            ps0 (ji,jj) = zslpmax 
268            psx (ji,jj) = psx (ji,jj) * zin0
269            psxx(ji,jj) = psxx(ji,jj) * zin0
270            psy (ji,jj) = zs1new * zin0
271            psyy(ji,jj) = zs2new * zin0
272            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
273         END DO
274      END DO
275
276      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
277      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
278
279      !  Calculate fluxes and moments between boxes j<-->j+1             
280      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
281         DO ji = 1, jpi
282            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
283            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
284            zalfq        =  zalf * zalf
285            zalf1        =  1.0 - zalf
286            zalf1q       =  zalf1 * zalf1
287            zfm (ji,jj)  =  zalf  * psm(ji,jj)
288            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
289            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
290            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
291            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
292            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
293            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
294            !
295            !  Readjust moments remaining in the box.
296            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
297            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
298            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
299            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
300            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
301            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
302            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
303         END DO
304      END DO
305      !
306      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
307         DO ji = 1, jpi
308            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
309            zalg  (ji,jj) = zalf
310            zalfq         = zalf * zalf
311            zalf1         = 1.0 - zalf
312            zalg1 (ji,jj) = zalf1
313            zalf1q        = zalf1 * zalf1
314            zalg1q(ji,jj) = zalf1q
315            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1)
316            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
317            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
318            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1)
319            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) )
320            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)
321            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1)
322         END DO
323      END DO
324 
325      !  Readjust moments remaining in the box.
326      DO jj = 2, jpj
327         DO ji = 1, jpi
328            zbt  =         zbet(ji,jj-1)
329            zbt1 = ( 1.0 - zbet(ji,jj-1) )
330            !
331            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
332            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
333            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
334            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj)
335            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
336            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
337            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
338         END DO
339      END DO
340
341      !   Put the temporary moments into appropriate neighboring boxes.   
342      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
343         DO ji = 1, jpi
344            zbt  =         zbet(ji,jj-1)
345            zbt1 = ( 1.0 - zbet(ji,jj-1) )
346            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
347            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
348            zalf1       = 1.0 - zalf
349            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
350            !
351            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj)
352            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
353               &        + zbt1 * psy(ji,jj) 
354            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
355               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
356               &        + zbt1 * psyy(ji,jj)
357            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              &
358               &                 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   &
359               &        + zbt1 * psxy(ji,jj)
360            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
361            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
362         END DO
363      END DO
364      !
365      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
366         DO ji = 1, jpi
367            zbt  =         zbet(ji,jj)
368            zbt1 = ( 1.0 - zbet(ji,jj) )
369            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
370            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
371            zalf1       = 1.0 - zalf
372            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
373            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
374            psy(ji,jj)  = zbt  * psy(ji,jj)  &
375               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
376            psyy(ji,jj) = zbt  * psyy(ji,jj)                                                                         &
377               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)                           &
378               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
379            psxy(ji,jj) = zbt  * psxy(ji,jj)                                  &
380               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   &
381               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )
382            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
383            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
384         END DO
385      END DO
386
387      !-- Lateral boundary conditions
388      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
389      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
390      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
391      CALL lbc_lnk( psxy, 'T',  1. )
392
393      IF(ln_ctl) THEN
394         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
395         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
396         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
397         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
398      ENDIF
399      !
400   END SUBROUTINE lim_adv_y_2
401
402#else
403   !!----------------------------------------------------------------------
404   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
405   !!----------------------------------------------------------------------
406#endif
407
408   !!======================================================================
409END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.