source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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