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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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