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/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 2804

Last change on this file since 2804 was 2804, checked in by rblod, 13 years ago

dev_r2802_LOCEAN10_agrif_lim: first implementation see ticket #848
  • 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 lib_mpp            ! MPP library
25   USE prtctl             ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_adv_x_2   ! called by lim_trp
31   PUBLIC   lim_adv_y_2   ! called by lim_trp
32
33   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
34   REAL(wp)  ::   rzero  = 0.e0     !    -       -
35   REAL(wp)  ::   rone   = 1.e0     !    -       -
36   
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
48      &                    psx, psxx, psy , psyy, psxy )
49      !!---------------------------------------------------------------------
50      !!                **  routine lim_adv_x_2  **
51      !! 
52      !! ** purpose :   Computes and adds the advection trend to sea-ice
53      !!              variable on i-axis
54      !!
55      !! ** method  :   Uses Prather second order scheme that advects tracers
56      !!              but also theirquadratic forms. The method preserves
57      !!              tracer structures by conserving second order moments.
58      !!
59      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
60      !!--------------------------------------------------------------------
61      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
62      USE wrk_nemo, ONLY:   zf0  => wrk_2d_18 , zfx   => wrk_2d_19 , zfy    => wrk_2d_20 , zbet => wrk_2d_21   ! 2D workspace
63      USE wrk_nemo, ONLY:   zfm  => wrk_2d_22 , zfxx  => wrk_2d_23 , zfyy   => wrk_2d_24 , zfxy => wrk_2d_25   !  -      -
64      USE wrk_nemo, ONLY:   zalg => wrk_2d_26 , zalg1 => wrk_2d_27 , zalg1q => wrk_2d_28                       !  -      -
65      !
66      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
67      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
68      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
69      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
70      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
71      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
72      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
73      !
74      INTEGER  ::   ji, jj                               ! dummy loop indices
75      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
76      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
77      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
78      !---------------------------------------------------------------------
79
80      IF( wrk_in_use(2, 18,19,20,21,22,23,24,25,26,27,28) ) THEN
81         CALL ctl_stop( 'lim_adv_x_2 : requested workspace arrays unavailable.' )   ;   RETURN
82      ENDIF
83
84      ! Limitation of moments.                                           
85
86      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
87
88      DO jj = 1, jpj
89         DO ji = 1, jpi
90            zslpmax = MAX( rzero, ps0(ji,jj) )
91            zs1max  = 1.5 * zslpmax
92            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
93            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
94               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
95            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
96            !
97            ps0 (ji,jj) = zslpmax 
98            psx (ji,jj) = zs1new      * zin0
99            psxx(ji,jj) = zs2new      * zin0
100            psy (ji,jj) = psy (ji,jj) * zin0
101            psyy(ji,jj) = psyy(ji,jj) * zin0
102            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
103         END DO
104      END DO
105
106      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
107      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
108
109      !  Calculate fluxes and moments between boxes i<-->i+1             
110      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
111         DO ji = 1, jpi
112            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
113            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
114            zalfq        =  zalf * zalf
115            zalf1        =  1.0 - zalf
116            zalf1q       =  zalf1 * zalf1
117            !
118            zfm (ji,jj)  =  zalf  * psm(ji,jj)
119            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) )
120            zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
121            zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj)
122            zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) )
123            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
124            zfyy(ji,jj)  =  zalf  * psyy(ji,jj)
125            !
126            !  Readjust moments remaining in the box.
127            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
128            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
129            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
130            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
131            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
132            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
133            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
134         END DO
135      END DO
136
137      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
138         DO ji = 1, fs_jpim1
139            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
140            zalg  (ji,jj) = zalf
141            zalfq         = zalf * zalf
142            zalf1         = 1.0 - zalf
143            zalg1 (ji,jj) = zalf1
144            zalf1q        = zalf1 * zalf1
145            zalg1q(ji,jj) = zalf1q
146            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj)
147            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
148            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
149            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj)
150            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) )
151            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)
152            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj)
153         END DO
154      END DO
155
156      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
157         DO ji = fs_2, fs_jpim1
158            zbt  =       zbet(ji-1,jj)
159            zbt1 = 1.0 - zbet(ji-1,jj)
160            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
161            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
162            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
163            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
164            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
165            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
166            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
167         END DO
168      END DO
169
170      !   Put the temporary moments into appropriate neighboring boxes.   
171      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
172         DO ji = fs_2, fs_jpim1
173            zbt  =       zbet(ji-1,jj)
174            zbt1 = 1.0 - zbet(ji-1,jj)
175            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
176            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
177            zalf1       = 1.0 - zalf
178            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
179            !
180            ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj)
181            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
182            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
183               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
184               &        + zbt1 * psxx(ji,jj)
185            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
186               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
187               &         + zbt1 * psxy(ji,jj)
188            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
189            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
190         END DO
191      END DO
192
193      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
194         DO ji = fs_2, fs_jpim1
195            zbt  =       zbet(ji,jj)
196            zbt1 = 1.0 - zbet(ji,jj)
197            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
198            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
199            zalf1       = 1.0 - zalf
200            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
201            !
202            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
203            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
204            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)   &
205               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )        &
206               &                                      + ( zalf1 - zalf ) * ztemp )                                 )
207            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)          &
208               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
209            psy(ji,jj)  = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
210            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
211         END DO
212      END DO
213
214      !-- Lateral boundary conditions
215      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
216      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
217      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
218      CALL lbc_lnk( psxy, 'T',  1. )
219
220      IF(ln_ctl) THEN
221         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
222         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
223         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
224         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
225      ENDIF
226      !
227      IF( wrk_not_released(2, 18,19,20,21,22,23,24,25,26,27,28) )   &
228          CALL ctl_stop( 'lim_adv_x_2 : failed to release workspace arrays.' )
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_18 , zfx   => wrk_2d_19 , zfy    => wrk_2d_20 , zbet => wrk_2d_21   ! 2D workspace
249      USE wrk_nemo, ONLY:   zfm  => wrk_2d_22 , zfxx  => wrk_2d_23 , zfyy   => wrk_2d_24 , zfxy => wrk_2d_25   !  -      -
250      USE wrk_nemo, ONLY:   zalg => wrk_2d_26 , zalg1 => wrk_2d_27 , zalg1q => wrk_2d_28                       !  -      -
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, 18,19,20,21,22,23,24,25,26,27,28) ) THEN
267         CALL ctl_stop( 'lim_adv_y_2 : requested workspace arrays unavailable.' )   ;   RETURN
268      ENDIF
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, 18,19,20,21,22,23,24,25,26,27,28) )   &
417          CALL ctl_stop( 'lim_adv_y_2 : failed to release workspace arrays.' )
418      !
419   END SUBROUTINE lim_adv_y_2
420
421#else
422   !!----------------------------------------------------------------------
423   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
424   !!----------------------------------------------------------------------
425#endif
426
427   !!======================================================================
428END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.