source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limadv_prather.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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