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 branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 3522

Last change on this file since 3522 was 3294, checked in by rblod, 12 years ago

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