New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limadv.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 2760

Last change on this file since 2760 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 23.2 KB
Line 
1MODULE limadv 
2   !!======================================================================
3   !!                       ***  MODULE limadv   ***
4   !! LIM sea-ice model : sea-ice advection
5   !!======================================================================
6   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!--------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                     LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_adv_x  : advection of sea ice on x axis
15   !!   lim_adv_y  : advection of sea ice on y axis
16   !!----------------------------------------------------------------------
17   USE dom_oce          ! ocean domain
18   USE dom_ice          ! LIM-3 domain
19   USE ice              ! LIM-3 variables
20   USE lbclnk           ! lateral boundary condition - MPP exchanges
21   USE in_out_manager   ! I/O manager
22   USE prtctl           ! Print control
23   USE lib_mpp          ! MPP library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   lim_adv_x   ! called by lim_trp
29   PUBLIC   lim_adv_y   ! called by lim_trp
30
31   REAL(wp)  ::   epsi20 = 1.e-20_wp   ! constant values
32   REAL(wp)  ::   rzero  = 0._wp       !    -       -
33   REAL(wp)  ::   rone   = 1._wp       !    -       -
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
45      &                  psx, psxx, psy , psyy, psxy )
46      !!---------------------------------------------------------------------
47      !!                **  routine lim_adv_x  **
48      !! 
49      !! ** purpose :   Computes and adds the advection trend to sea-ice
50      !!              variable on i-axis
51      !!
52      !! ** method  :   Uses Prather second order scheme that advects tracers
53      !!              but also theirquadratic forms. The method preserves
54      !!              tracer structures by conserving second order moments.
55      !!
56      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
57      !!--------------------------------------------------------------------
58      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
59      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
60      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
61      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
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   ! local scalars
73      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
74      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
75      !---------------------------------------------------------------------
76
77      IF( wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
78         CALL ctl_stop('lim_adv_x: requested workspace arrays unavailable')   ;   RETURN
79      ENDIF
80
81      ! Limitation of moments.                                           
82
83      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
84
85      DO jj = 1, jpj
86         DO ji = 1, jpi
87            zslpmax = MAX( rzero, ps0(ji,jj) )
88            zs1max  = 1.5 * zslpmax
89            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
90            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
91               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
92            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
93
94            ps0 (ji,jj) = zslpmax 
95            psx (ji,jj) = zs1new      * zin0
96            psxx(ji,jj) = zs2new      * zin0
97            psy (ji,jj) = psy (ji,jj) * zin0
98            psyy(ji,jj) = psyy(ji,jj) * zin0
99            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
100         END DO
101      END DO
102
103      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
104      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
105
106      !  Calculate fluxes and moments between boxes i<-->i+1             
107      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
108         DO ji = 1, jpi
109            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
110            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
111            zalfq        =  zalf * zalf
112            zalf1        =  1.0 - zalf
113            zalf1q       =  zalf1 * zalf1
114            !
115            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
116            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
117            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
118            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
119            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
120            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
121            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
122
123            !  Readjust moments remaining in the box.
124            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
125            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
126            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
127            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
128            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
129            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
130            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
131         END DO
132      END DO
133
134      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
135         DO ji = 1, fs_jpim1
136            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
137            zalg  (ji,jj) = zalf
138            zalfq         = zalf * zalf
139            zalf1         = 1.0 - zalf
140            zalg1 (ji,jj) = zalf1
141            zalf1q        = zalf1 * zalf1
142            zalg1q(ji,jj) = zalf1q
143            !
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  *   psxx(ji+1,jj) * zalfq
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            !
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      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) )    &
227          CALL ctl_stop('lim_adv_x : failed to release workspace arrays')
228      !
229   END SUBROUTINE lim_adv_x
230
231
232   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
233      &                  psx, psxx, psy , psyy, psxy )
234      !!---------------------------------------------------------------------
235      !!                **  routine lim_adv_y  **
236      !!           
237      !! ** purpose :   Computes and adds the advection trend to sea-ice
238      !!              variable on y 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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
247      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
248      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
249      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
250      !
251      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
252      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
253      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
254      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
255      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
256      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
257      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
258      !!
259      INTEGER  ::   ji, jj                               ! dummy loop indices
260      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
261      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
262      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
263      !---------------------------------------------------------------------
264
265      IF( wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
266         CALL ctl_stop('lim_adv_y : requested workspace arrays unavailable')   ;   RETURN
267      ENDIF
268
269      ! Limitation of moments.
270
271      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
272
273      DO jj = 1, jpj
274         DO ji = 1, jpi
275            zslpmax = MAX( rzero, ps0(ji,jj) )
276            zs1max  = 1.5 * zslpmax
277            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
278            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
279               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
280            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
281            !
282            ps0 (ji,jj) = zslpmax 
283            psx (ji,jj) = psx (ji,jj) * zin0
284            psxx(ji,jj) = psxx(ji,jj) * zin0
285            psy (ji,jj) = zs1new * zin0
286            psyy(ji,jj) = zs2new * zin0
287            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
288         END DO
289      END DO
290
291      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
292      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
293
294      !  Calculate fluxes and moments between boxes j<-->j+1             
295      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
296         DO ji = 1, jpi
297            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
298            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
299            zalfq        =  zalf * zalf
300            zalf1        =  1.0 - zalf
301            zalf1q       =  zalf1 * zalf1
302            !
303            zfm (ji,jj)  =  zalf  * psm(ji,jj)
304            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
305            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
306            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
307            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
308            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
309            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
310            !
311            !  Readjust moments remaining in the box.
312            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
313            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
314            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
315            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
316            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
317            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
318            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
319         END DO
320      END DO
321      !
322      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
323         DO ji = 1, jpi
324            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
325            zalg  (ji,jj) = zalf
326            zalfq         = zalf * zalf
327            zalf1         = 1.0 - zalf
328            zalg1 (ji,jj) = zalf1
329            zalf1q        = zalf1 * zalf1
330            zalg1q(ji,jj) = zalf1q
331            !
332            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
333            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
334            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
335            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
336            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
337            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
338            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
339         END DO
340      END DO
341
342      !  Readjust moments remaining in the box.
343      DO jj = 2, jpj
344         DO ji = 1, jpi
345            zbt  =         zbet(ji,jj-1)
346            zbt1 = ( 1.0 - zbet(ji,jj-1) )
347            !
348            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
349            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
350            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
351            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
352            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
353            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
354            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
355         END DO
356      END DO
357
358      !   Put the temporary moments into appropriate neighboring boxes.   
359      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
360         DO ji = 1, jpi
361            zbt  =         zbet(ji,jj-1)
362            zbt1 = ( 1.0 - zbet(ji,jj-1) )
363            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
364            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
365            zalf1       = 1.0 - zalf
366            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
367            !
368            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
369            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
370               &                                               + zbt1 * psy(ji,jj) 
371            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
372               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
373               &                                               + zbt1 * psyy(ji,jj)
374            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
375               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
376               &                                                + zbt1 * psxy(ji,jj)
377            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
378            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
379         END DO
380      END DO
381
382      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
383         DO ji = 1, jpi
384            zbt  =         zbet(ji,jj)
385            zbt1 = ( 1.0 - zbet(ji,jj) )
386            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
387            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
388            zalf1       = 1.0 - zalf
389            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
390            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
391            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
392            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
393               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
394               &                                         + ( zalf1 - zalf ) * ztemp )                                )
395            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
396               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
397            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
398            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
399         END DO
400      END DO
401
402      !-- Lateral boundary conditions
403      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
404      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
405      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
406      CALL lbc_lnk( psxy, 'T',  1. )
407
408      IF(ln_ctl) THEN
409         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
410         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
411         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
412         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
413      ENDIF
414      !
415      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) )    &
416         CALL ctl_stop('lim_adv_y: failed to release workspace arrays')
417      !
418   END SUBROUTINE lim_adv_y
419
420#else
421   !!----------------------------------------------------------------------
422   !!   Default option            Dummy module         NO LIM sea-ice model
423   !!----------------------------------------------------------------------
424#endif
425
426   !!======================================================================
427END MODULE limadv
Note: See TracBrowser for help on using the repository browser.