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 branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 4921

Last change on this file since 4921 was 4921, checked in by timgraham, 9 years ago

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