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 @ 3750

Last change on this file since 3750 was 3625, checked in by acc, 11 years ago

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