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_2.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.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.0 KB
Line 
1MODULE limadv_2 
2   !!======================================================================
3   !!                       ***  MODULE limadv_2   ***
4   !! LIM 2.0 sea-ice model : sea-ice advection
5   !!======================================================================
6   !! History :  OPA  ! 2000-01 (LIM)  Original code
7   !!                 ! 2001-05 (G. Madec, R. Hordoir) Doctor norm
8   !!   NEMO     1.0  ! 2003-10 (C. Ethe) F90, module
9   !!             -   ! 2003-12 (R. Hordoir, G. Madec) mpp
10   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
11   !!--------------------------------------------------------------------
12#if defined key_lim2
13   !!----------------------------------------------------------------------
14   !!   'key_lim2'                                    LIM 2.0 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   lim_adv_x_2  : advection of sea ice on x axis
17   !!   lim_adv_y_2  : advection of sea ice on y axis
18   !!----------------------------------------------------------------------
19   USE dom_oce
20   USE dom_ice_2
21   USE ice_2
22   USE lbclnk
23   USE in_out_manager     ! I/O manager
24   USE prtctl             ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   lim_adv_x_2   ! called by lim_trp
30   PUBLIC   lim_adv_y_2   ! called by lim_trp
31
32   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
33   REAL(wp)  ::   rzero  = 0.e0     !    -       -
34   REAL(wp)  ::   rone   = 1.e0     !    -       -
35   
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
47      &                    psx, psxx, psy , psyy, psxy )
48      !!---------------------------------------------------------------------
49      !!                **  routine lim_adv_x_2  **
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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
61      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
62      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
63      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
64      !
65      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
66      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
67      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
68      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
69      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
70      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
71      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
72      !
73      INTEGER  ::   ji, jj                               ! dummy loop indices
74      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
75      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
76      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
77      !---------------------------------------------------------------------
78
79      IF( wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
80         CALL ctl_stop( 'lim_adv_x_2 : requested workspace arrays unavailable.' )   ;   RETURN
81      END IF
82
83      ! Limitation of moments.                                           
84
85      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
86
87      DO jj = 1, jpj
88         DO ji = 1, jpi
89            zslpmax = MAX( rzero, ps0(ji,jj) )
90            zs1max  = 1.5 * zslpmax
91            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
92            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
93               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
94            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
95            !
96            ps0 (ji,jj) = zslpmax 
97            psx (ji,jj) = zs1new      * zin0
98            psxx(ji,jj) = zs2new      * zin0
99            psy (ji,jj) = psy (ji,jj) * zin0
100            psyy(ji,jj) = psyy(ji,jj) * zin0
101            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
102         END DO
103      END DO
104
105      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
106      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
107
108      !  Calculate fluxes and moments between boxes i<-->i+1             
109      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
110         DO ji = 1, jpi
111            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
112            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
113            zalfq        =  zalf * zalf
114            zalf1        =  1.0 - zalf
115            zalf1q       =  zalf1 * zalf1
116            !
117            zfm (ji,jj)  =  zalf  * psm(ji,jj)
118            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) )
119            zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
120            zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj)
121            zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) )
122            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
123            zfyy(ji,jj)  =  zalf  * psyy(ji,jj)
124            !
125            !  Readjust moments remaining in the box.
126            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
127            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
128            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
129            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
130            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
131            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
132            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
133         END DO
134      END DO
135
136      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
137         DO ji = 1, fs_jpim1
138            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
139            zalg  (ji,jj) = zalf
140            zalfq         = zalf * zalf
141            zalf1         = 1.0 - zalf
142            zalg1 (ji,jj) = zalf1
143            zalf1q        = zalf1 * zalf1
144            zalg1q(ji,jj) = zalf1q
145            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj)
146            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
147            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
148            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj)
149            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) )
150            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)
151            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj)
152         END DO
153      END DO
154
155      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
156         DO ji = fs_2, fs_jpim1
157            zbt  =       zbet(ji-1,jj)
158            zbt1 = 1.0 - zbet(ji-1,jj)
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      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
227         CALL ctl_stop( 'lim_adv_x_2 : failed to release workspace arrays.' )
228      END IF
229      !
230   END SUBROUTINE lim_adv_x_2
231
232
233   SUBROUTINE lim_adv_y_2( pdf, pvt , pcrh, psm , ps0 ,   &
234      &                    psx, psxx, psy , psyy, psxy )
235      !!---------------------------------------------------------------------
236      !!                **  routine lim_adv_y_2  **
237      !!           
238      !! ** purpose :   Computes and adds the advection trend to sea-ice
239      !!              variable on j-axis
240      !!
241      !! ** method  :   Uses Prather second order scheme that advects tracers
242      !!              but also their quadratic forms. The method preserves
243      !!              tracer structures by conserving second order moments.
244      !!
245      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
246      !!---------------------------------------------------------------------
247      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
248      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
249      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
250      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
251      !!
252      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
253      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
254      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
255      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
256      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
257      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
258      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
259      !!
260      INTEGER  ::   ji, jj                               ! dummy loop indices
261      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
262      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
263      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
264      !---------------------------------------------------------------------
265
266      IF(wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
267         CALL ctl_stop( 'lim_adv_y_2 : requested workspace arrays unavailable.' )   ;   RETURN
268      END IF
269
270      ! Limitation of moments.
271
272      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
273
274      DO jj = 1, jpj
275         DO ji = 1, jpi
276            zslpmax = MAX( rzero, ps0(ji,jj) )
277            zs1max  = 1.5 * zslpmax
278            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
279            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
280               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
281            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
282            !
283            ps0 (ji,jj) = zslpmax 
284            psx (ji,jj) = psx (ji,jj) * zin0
285            psxx(ji,jj) = psxx(ji,jj) * zin0
286            psy (ji,jj) = zs1new * zin0
287            psyy(ji,jj) = zs2new * zin0
288            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
289         END DO
290      END DO
291
292      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
293      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
294
295      !  Calculate fluxes and moments between boxes j<-->j+1             
296      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
297         DO ji = 1, jpi
298            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
299            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
300            zalfq        =  zalf * zalf
301            zalf1        =  1.0 - zalf
302            zalf1q       =  zalf1 * zalf1
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            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1)
332            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
333            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
334            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1)
335            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) )
336            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)
337            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1)
338         END DO
339      END DO
340 
341      !  Readjust moments remaining in the box.
342      DO jj = 2, jpj
343         DO ji = 1, jpi
344            zbt  =         zbet(ji,jj-1)
345            zbt1 = ( 1.0 - zbet(ji,jj-1) )
346            !
347            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
348            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
349            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
350            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj)
351            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
352            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
353            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
354         END DO
355      END DO
356
357      !   Put the temporary moments into appropriate neighboring boxes.   
358      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
359         DO ji = 1, jpi
360            zbt  =         zbet(ji,jj-1)
361            zbt1 = ( 1.0 - zbet(ji,jj-1) )
362            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
363            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
364            zalf1       = 1.0 - zalf
365            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
366            !
367            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj)
368            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
369               &        + zbt1 * psy(ji,jj) 
370            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
371               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
372               &        + zbt1 * psyy(ji,jj)
373            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              &
374               &                 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   &
375               &        + zbt1 * psxy(ji,jj)
376            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
377            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
378         END DO
379      END DO
380      !
381      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
382         DO ji = 1, jpi
383            zbt  =         zbet(ji,jj)
384            zbt1 = ( 1.0 - zbet(ji,jj) )
385            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
386            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
387            zalf1       = 1.0 - zalf
388            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
389            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
390            psy(ji,jj)  = zbt  * psy(ji,jj)  &
391               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
392            psyy(ji,jj) = zbt  * psyy(ji,jj)                                                                         &
393               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)                           &
394               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
395            psxy(ji,jj) = zbt  * psxy(ji,jj)                                  &
396               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   &
397               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )
398            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
399            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
400         END DO
401      END DO
402
403      !-- Lateral boundary conditions
404      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
405      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
406      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
407      CALL lbc_lnk( psxy, 'T',  1. )
408
409      IF(ln_ctl) THEN
410         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
411         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
412         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
413         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
414      ENDIF
415      !
416      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
417        CALL ctl_stop( 'lim_adv_y_2 : failed to release workspace arrays.' )
418      END IF
419      !
420   END SUBROUTINE lim_adv_y_2
421
422#else
423   !!----------------------------------------------------------------------
424   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
425   !!----------------------------------------------------------------------
426#endif
427
428   !!======================================================================
429END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.