source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 9295

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

Remove svn keywords

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