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/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 3764

Last change on this file since 3764 was 3764, checked in by smasson, 11 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012
  • Property svn:keywords set to Id
File size: 22.9 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 ice            ! LIM-3 variables
19   USE dom_ice        ! LIM-3 domain
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    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_adv_x   ! called by lim_trp
32   PUBLIC   lim_adv_y   ! called by lim_trp
33
34   REAL(wp)  ::   epsi20 = 1.e-20_wp   ! constant values
35   REAL(wp)  ::   rzero  = 0._wp       !    -       -
36   REAL(wp)  ::   rone   = 1._wp       !    -       -
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
48      &                  psx, psxx, psy , psyy, psxy )
49      !!---------------------------------------------------------------------
50      !!                **  routine lim_adv_x  **
51      !! 
52      !! ** purpose :   Computes and adds the advection trend to sea-ice
53      !!              variable on i-axis
54      !!
55      !! ** method  :   Uses Prather second order scheme that advects tracers
56      !!              but also theirquadratic forms. The method preserves
57      !!              tracer structures by conserving second order moments.
58      !!
59      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
60      !!--------------------------------------------------------------------
61      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
62      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
63      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
64      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
65      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
66      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
67      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
68      !!
69      INTEGER  ::   ji, jj                               ! dummy loop indices
70      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! local scalars
71      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
72      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
73      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
74      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
75      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q         !  -      -
76      !---------------------------------------------------------------------
77
78      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
79      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
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      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
227      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
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      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
247      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
248      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
249      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
250      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
251      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
252      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
253      !!
254      INTEGER  ::   ji, jj                               ! dummy loop indices
255      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
256      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
257      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
258      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
259      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
260      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
261      !---------------------------------------------------------------------
262
263      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
264      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
265
266      ! Limitation of moments.
267
268      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
269
270      DO jj = 1, jpj
271         DO ji = 1, jpi
272            zslpmax = MAX( rzero, ps0(ji,jj) )
273            zs1max  = 1.5 * zslpmax
274            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
275            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
276               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
277            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
278            !
279            ps0 (ji,jj) = zslpmax 
280            psx (ji,jj) = psx (ji,jj) * zin0
281            psxx(ji,jj) = psxx(ji,jj) * zin0
282            psy (ji,jj) = zs1new * zin0
283            psyy(ji,jj) = zs2new * zin0
284            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
285         END DO
286      END DO
287
288      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
289      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
290
291      !  Calculate fluxes and moments between boxes j<-->j+1             
292      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
293         DO ji = 1, jpi
294            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
295            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
296            zalfq        =  zalf * zalf
297            zalf1        =  1.0 - zalf
298            zalf1q       =  zalf1 * zalf1
299            !
300            zfm (ji,jj)  =  zalf  * psm(ji,jj)
301            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
302            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
303            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
304            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
305            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
306            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
307            !
308            !  Readjust moments remaining in the box.
309            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
310            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
311            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
312            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
313            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
314            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
315            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
316         END DO
317      END DO
318      !
319      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
320         DO ji = 1, jpi
321            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
322            zalg  (ji,jj) = zalf
323            zalfq         = zalf * zalf
324            zalf1         = 1.0 - zalf
325            zalg1 (ji,jj) = zalf1
326            zalf1q        = zalf1 * zalf1
327            zalg1q(ji,jj) = zalf1q
328            !
329            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
330            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
331            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
332            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
333            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
334            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
335            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
336         END DO
337      END DO
338
339      !  Readjust moments remaining in the box.
340      DO jj = 2, jpj
341         DO ji = 1, jpi
342            zbt  =         zbet(ji,jj-1)
343            zbt1 = ( 1.0 - zbet(ji,jj-1) )
344            !
345            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
346            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
347            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
348            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
349            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
350            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
351            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
352         END DO
353      END DO
354
355      !   Put the temporary moments into appropriate neighboring boxes.   
356      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
357         DO ji = 1, jpi
358            zbt  =         zbet(ji,jj-1)
359            zbt1 = ( 1.0 - zbet(ji,jj-1) )
360            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
361            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
362            zalf1       = 1.0 - zalf
363            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
364            !
365            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
366            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
367               &                                               + zbt1 * psy(ji,jj) 
368            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
369               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
370               &                                               + zbt1 * psyy(ji,jj)
371            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
372               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
373               &                                                + zbt1 * psxy(ji,jj)
374            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
375            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
376         END DO
377      END DO
378
379      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
380         DO ji = 1, jpi
381            zbt  =         zbet(ji,jj)
382            zbt1 = ( 1.0 - zbet(ji,jj) )
383            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
384            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
385            zalf1       = 1.0 - zalf
386            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
387            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
388            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
389            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
390               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
391               &                                         + ( zalf1 - zalf ) * ztemp )                                )
392            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
393               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
394            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
395            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
396         END DO
397      END DO
398
399      !-- Lateral boundary conditions
400      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
401      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
402      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
403      CALL lbc_lnk( psxy, 'T',  1. )
404
405      IF(ln_ctl) THEN
406         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
407         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
408         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
409         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
410      ENDIF
411      !
412      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
413      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
414      !
415   END SUBROUTINE lim_adv_y
416
417#else
418   !!----------------------------------------------------------------------
419   !!   Default option            Dummy module         NO LIM sea-ice model
420   !!----------------------------------------------------------------------
421#endif
422
423   !!======================================================================
424END MODULE limadv
Note: See TracBrowser for help on using the repository browser.