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 trunk/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 3564

Last change on this file since 3564 was 3558, checked in by rblod, 11 years ago

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