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/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

All wrk_alloc removed
  • Property svn:keywords set to Id
File size: 22.1 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   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_adv_x_2   ! called by lim_trp
32   PUBLIC   lim_adv_y_2   ! called by lim_trp
33
34   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
35   REAL(wp)  ::   rzero  = 0.e0     !    -       -
36   REAL(wp)  ::   rone   = 1.e0     !    -       -
37   
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
49      &                    psx, psxx, psy , psyy, psxy )
50      !!---------------------------------------------------------------------
51      !!                **  routine lim_adv_x_2  **
52      !! 
53      !! ** purpose :   Computes and adds the advection trend to sea-ice
54      !!              variable on i-axis
55      !!
56      !! ** method  :   Uses Prather second order scheme that advects tracers
57      !!              but also theirquadratic forms. The method preserves
58      !!              tracer structures by conserving second order moments.
59      !!
60      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
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   ! temporary scalars
72      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
73      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
74      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
75      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
76      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
77      !---------------------------------------------------------------------
78
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  * zalfq * psxx(ji,jj)
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            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  * zalfq * psxx(ji+1,jj)
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            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
157            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
158            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
159            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
160            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
161            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
162            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
163         END DO
164      END DO
165
166      !   Put the temporary moments into appropriate neighboring boxes.   
167      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
168         DO ji = fs_2, fs_jpim1
169            zbt  =       zbet(ji-1,jj)
170            zbt1 = 1.0 - zbet(ji-1,jj)
171            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
172            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
173            zalf1       = 1.0 - zalf
174            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
175            !
176            ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj)
177            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
178            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
179               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
180               &        + zbt1 * psxx(ji,jj)
181            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
182               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
183               &         + zbt1 * psxy(ji,jj)
184            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
185            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
186         END DO
187      END DO
188
189      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
190         DO ji = fs_2, fs_jpim1
191            zbt  =       zbet(ji,jj)
192            zbt1 = 1.0 - zbet(ji,jj)
193            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
194            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
195            zalf1       = 1.0 - zalf
196            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
197            !
198            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
199            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
200            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)   &
201               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )        &
202               &                                      + ( zalf1 - zalf ) * ztemp )                                 )
203            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)          &
204               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
205            psy(ji,jj)  = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
206            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
207         END DO
208      END DO
209
210      !-- Lateral boundary conditions
211      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
212      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
213      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
214      CALL lbc_lnk( psxy, 'T',  1. )
215
216      IF(ln_ctl) THEN
217         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
218         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
219         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
220         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
221      ENDIF
222      !
223      !
224   END SUBROUTINE lim_adv_x_2
225
226
227   SUBROUTINE lim_adv_y_2( pdf, pvt , pcrh, psm , ps0 ,   &
228      &                    psx, psxx, psy , psyy, psxy )
229      !!---------------------------------------------------------------------
230      !!                **  routine lim_adv_y_2  **
231      !!           
232      !! ** purpose :   Computes and adds the advection trend to sea-ice
233      !!              variable on j-axis
234      !!
235      !! ** method  :   Uses Prather second order scheme that advects tracers
236      !!              but also their quadratic forms. The method preserves
237      !!              tracer structures by conserving second order moments.
238      !!
239      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
240      !!---------------------------------------------------------------------
241      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
242      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
243      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
244      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
246      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
247      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
248      !!
249      INTEGER  ::   ji, jj                               ! dummy loop indices
250      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
251      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
252      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
253      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
254      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
255      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
256      !---------------------------------------------------------------------
257
258
259      ! Limitation of moments.
260
261      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
262
263      DO jj = 1, jpj
264         DO ji = 1, jpi
265            zslpmax = MAX( rzero, ps0(ji,jj) )
266            zs1max  = 1.5 * zslpmax
267            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
268            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
269               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
270            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
271            !
272            ps0 (ji,jj) = zslpmax 
273            psx (ji,jj) = psx (ji,jj) * zin0
274            psxx(ji,jj) = psxx(ji,jj) * zin0
275            psy (ji,jj) = zs1new * zin0
276            psyy(ji,jj) = zs2new * zin0
277            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
278         END DO
279      END DO
280
281      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
282      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
283
284      !  Calculate fluxes and moments between boxes j<-->j+1             
285      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
286         DO ji = 1, jpi
287            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
288            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
289            zalfq        =  zalf * zalf
290            zalf1        =  1.0 - zalf
291            zalf1q       =  zalf1 * zalf1
292            zfm (ji,jj)  =  zalf  * psm(ji,jj)
293            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
294            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
295            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
296            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
297            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
298            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
299            !
300            !  Readjust moments remaining in the box.
301            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
302            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
303            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
304            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
305            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
306            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
307            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
308         END DO
309      END DO
310      !
311      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
312         DO ji = 1, jpi
313            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
314            zalg  (ji,jj) = zalf
315            zalfq         = zalf * zalf
316            zalf1         = 1.0 - zalf
317            zalg1 (ji,jj) = zalf1
318            zalf1q        = zalf1 * zalf1
319            zalg1q(ji,jj) = zalf1q
320            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1)
321            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
322            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
323            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1)
324            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) )
325            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)
326            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1)
327         END DO
328      END DO
329 
330      !  Readjust moments remaining in the box.
331      DO jj = 2, jpj
332         DO ji = 1, jpi
333            zbt  =         zbet(ji,jj-1)
334            zbt1 = ( 1.0 - zbet(ji,jj-1) )
335            !
336            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
337            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
338            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
339            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj)
340            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
341            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
342            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
343         END DO
344      END DO
345
346      !   Put the temporary moments into appropriate neighboring boxes.   
347      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
348         DO ji = 1, jpi
349            zbt  =         zbet(ji,jj-1)
350            zbt1 = ( 1.0 - zbet(ji,jj-1) )
351            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
352            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
353            zalf1       = 1.0 - zalf
354            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
355            !
356            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj)
357            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
358               &        + zbt1 * psy(ji,jj) 
359            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
360               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
361               &        + zbt1 * psyy(ji,jj)
362            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              &
363               &                 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   &
364               &        + zbt1 * psxy(ji,jj)
365            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
366            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
367         END DO
368      END DO
369      !
370      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
371         DO ji = 1, jpi
372            zbt  =         zbet(ji,jj)
373            zbt1 = ( 1.0 - zbet(ji,jj) )
374            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
375            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
376            zalf1       = 1.0 - zalf
377            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
378            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
379            psy(ji,jj)  = zbt  * psy(ji,jj)  &
380               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
381            psyy(ji,jj) = zbt  * psyy(ji,jj)                                                                         &
382               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)                           &
383               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
384            psxy(ji,jj) = zbt  * psxy(ji,jj)                                  &
385               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   &
386               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )
387            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
388            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
389         END DO
390      END DO
391
392      !-- Lateral boundary conditions
393      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
394      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
395      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
396      CALL lbc_lnk( psxy, 'T',  1. )
397
398      IF(ln_ctl) THEN
399         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
400         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
401         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
402         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
403      ENDIF
404      !
405      !
406   END SUBROUTINE lim_adv_y_2
407
408#else
409   !!----------------------------------------------------------------------
410   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
411   !!----------------------------------------------------------------------
412#endif
413
414   !!======================================================================
415END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.