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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 11 years ago

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