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

source: branches/2015/dev_r5302_CNRS18_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 5372

Last change on this file since 5372 was 5372, checked in by mcastril, 9 years ago

ticket #1523 Message Packing
  • Property svn:keywords set to Id
File size: 23.1 KB
Line 
1MODULE limadv 
2   !!======================================================================
3   !!                       ***  MODULE limadv   ***
4   !! LIM sea-ice model : sea-ice advection
5   !!======================================================================
6   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!--------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                     LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_adv_x  : advection of sea ice on x axis
15   !!   lim_adv_y  : advection of sea ice on y axis
16   !!----------------------------------------------------------------------
17   USE dom_oce          ! ocean domain
18   USE dom_ice          ! LIM-3 domain
19   USE ice              ! LIM-3 variables
20   USE lbclnk           ! lateral boundary condition - MPP exchanges
21   USE in_out_manager   ! I/O manager
22   USE prtctl           ! Print control
23   USE lib_mpp          ! MPP library
24   USE wrk_nemo         ! work arrays
25   USE lib_fortran      ! to use key_nosignedzero
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_adv_x   ! called by lim_trp
31   PUBLIC   lim_adv_y   ! called by lim_trp
32
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
43      &                  psx, psxx, psy , psyy, psxy )
44      !!---------------------------------------------------------------------
45      !!                **  routine lim_adv_x  **
46      !! 
47      !! ** purpose :   Computes and adds the advection trend to sea-ice
48      !!              variable on i-axis
49      !!
50      !! ** method  :   Uses Prather second order scheme that advects tracers
51      !!              but also theirquadratic forms. The method preserves
52      !!              tracer structures by conserving second order moments.
53      !!
54      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
55      !!--------------------------------------------------------------------
56      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
57      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
60      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
62      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
63      !!
64      INTEGER  ::   ji, jj                               ! dummy loop indices
65      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
66      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
67      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
68      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
69      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
70      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q         !  -      -
71      !---------------------------------------------------------------------
72
73      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
74      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
75
76      ! Limitation of moments.                                           
77
78      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
79
80      DO jj = 1, jpj
81         DO ji = 1, jpi
82            zslpmax = MAX( 0._wp, ps0(ji,jj) )
83            zs1max  = 1.5 * zslpmax
84            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
85            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
86               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
87            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
88
89            ps0 (ji,jj) = zslpmax 
90            psx (ji,jj) = zs1new      * rswitch
91            psxx(ji,jj) = zs2new      * rswitch
92            psy (ji,jj) = psy (ji,jj) * rswitch
93            psyy(ji,jj) = psyy(ji,jj) * rswitch
94            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
95         END DO
96      END DO
97
98      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
99      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
100
101      !  Calculate fluxes and moments between boxes i<-->i+1             
102      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
103         DO ji = 1, jpi
104            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
105            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
106            zalfq        =  zalf * zalf
107            zalf1        =  1.0 - zalf
108            zalf1q       =  zalf1 * zalf1
109            !
110            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
111            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
112            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
113            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
114            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
115            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
116            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
117
118            !  Readjust moments remaining in the box.
119            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
120            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
121            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
122            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
123            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
124            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
125            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
126         END DO
127      END DO
128
129      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
130         DO ji = 1, fs_jpim1
131            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
132            zalg  (ji,jj) = zalf
133            zalfq         = zalf * zalf
134            zalf1         = 1.0 - zalf
135            zalg1 (ji,jj) = zalf1
136            zalf1q        = zalf1 * zalf1
137            zalg1q(ji,jj) = zalf1q
138            !
139            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
140            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
141            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
142            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
143            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
144            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
145            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
146         END DO
147      END DO
148
149      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
150         DO ji = fs_2, fs_jpim1
151            zbt  =       zbet(ji-1,jj)
152            zbt1 = 1.0 - zbet(ji-1,jj)
153            !
154            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
155            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
156            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
157            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
158            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
159            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
160            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
161         END DO
162      END DO
163
164      !   Put the temporary moments into appropriate neighboring boxes.   
165      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
166         DO ji = fs_2, fs_jpim1
167            zbt  =       zbet(ji-1,jj)
168            zbt1 = 1.0 - zbet(ji-1,jj)
169            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
170            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
171            zalf1       = 1.0 - zalf
172            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
173            !
174            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
175            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
176            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
177               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
178               &                                                + zbt1 * psxx(ji,jj)
179            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
180               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
181               &                                                + zbt1 * psxy(ji,jj)
182            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
183            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
184         END DO
185      END DO
186
187      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
188         DO ji = fs_2, fs_jpim1
189            zbt  =       zbet(ji,jj)
190            zbt1 = 1.0 - zbet(ji,jj)
191            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
192            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
193            zalf1       = 1.0 - zalf
194            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
195            !
196            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
197            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
198            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
199               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
200               &                                      + ( zalf1 - zalf ) * ztemp ) )
201            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
202               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
203            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
204            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
205         END DO
206      END DO
207
208
209#if defined key_multisend
210      CALL lbc_lnk_multi(  psm , 'T',  1. ,  ps0 , 'T',  1. ,  psx , 'T', -1. ,  psy , 'T', -1. ,  psxx, 'T',  1.  ,  psyy, 'T',  1. ,  psxy, 'T',  1. )
211#else
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#endif
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      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
227      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
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      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
247      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
248      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
249      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
250      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
251      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
252      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
253      !!
254      INTEGER  ::   ji, jj                               ! dummy loop indices
255      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
256      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
257      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
258      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
259      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
260      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
261      !---------------------------------------------------------------------
262
263      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
264      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
265
266      ! Limitation of moments.
267
268      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
269
270      DO jj = 1, jpj
271         DO ji = 1, jpi
272            zslpmax = MAX( 0._wp, ps0(ji,jj) )
273            zs1max  = 1.5 * zslpmax
274            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
275            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
276               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
277            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
278            !
279            ps0 (ji,jj) = zslpmax 
280            psx (ji,jj) = psx (ji,jj) * rswitch
281            psxx(ji,jj) = psxx(ji,jj) * rswitch
282            psy (ji,jj) = zs1new * rswitch
283            psyy(ji,jj) = zs2new * rswitch
284            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
285         END DO
286      END DO
287
288      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
289      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
290
291      !  Calculate fluxes and moments between boxes j<-->j+1             
292      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
293         DO ji = 1, jpi
294            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
295            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
296            zalfq        =  zalf * zalf
297            zalf1        =  1.0 - zalf
298            zalf1q       =  zalf1 * zalf1
299            !
300            zfm (ji,jj)  =  zalf  * psm(ji,jj)
301            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
302            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
303            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
304            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
305            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
306            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
307            !
308            !  Readjust moments remaining in the box.
309            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
310            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
311            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
312            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
313            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
314            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
315            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
316         END DO
317      END DO
318      !
319      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
320         DO ji = 1, jpi
321            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
322            zalg  (ji,jj) = zalf
323            zalfq         = zalf * zalf
324            zalf1         = 1.0 - zalf
325            zalg1 (ji,jj) = zalf1
326            zalf1q        = zalf1 * zalf1
327            zalg1q(ji,jj) = zalf1q
328            !
329            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
330            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
331            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
332            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
333            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
334            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
335            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
336         END DO
337      END DO
338
339      !  Readjust moments remaining in the box.
340      DO jj = 2, jpj
341         DO ji = 1, jpi
342            zbt  =         zbet(ji,jj-1)
343            zbt1 = ( 1.0 - zbet(ji,jj-1) )
344            !
345            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
346            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
347            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
348            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
349            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
350            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
351            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
352         END DO
353      END DO
354
355      !   Put the temporary moments into appropriate neighboring boxes.   
356      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
357         DO ji = 1, jpi
358            zbt  =         zbet(ji,jj-1)
359            zbt1 = ( 1.0 - zbet(ji,jj-1) )
360            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
361            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
362            zalf1       = 1.0 - zalf
363            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
364            !
365            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
366            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
367               &                                               + zbt1 * psy(ji,jj) 
368            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
369               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
370               &                                               + zbt1 * psyy(ji,jj)
371            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
372               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
373               &                                                + zbt1 * psxy(ji,jj)
374            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
375            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
376         END DO
377      END DO
378
379      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
380         DO ji = 1, jpi
381            zbt  =         zbet(ji,jj)
382            zbt1 = ( 1.0 - zbet(ji,jj) )
383            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
384            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
385            zalf1       = 1.0 - zalf
386            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
387            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
388            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
389            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
390               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
391               &                                         + ( zalf1 - zalf ) * ztemp )                                )
392            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
393               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
394            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
395            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
396         END DO
397      END DO
398
399#if defined key_multisend
400      CALL lbc_lnk_multi(  psm , 'T',  1. ,  ps0 , 'T',  1. ,  psx , 'T', -1. ,  psy , 'T', -1. ,  psxx, 'T',  1.  ,  psyy, 'T',  1. ,  psxy, 'T',  1. )
401#else
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#endif
408
409
410      IF(ln_ctl) THEN
411         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
412         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
413         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
414         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
415      ENDIF
416      !
417      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
418      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
419      !
420   END SUBROUTINE lim_adv_y
421
422#else
423   !!----------------------------------------------------------------------
424   !!   Default option            Dummy module         NO LIM sea-ice model
425   !!----------------------------------------------------------------------
426#endif
427
428   !!======================================================================
429END MODULE limadv
Note: See TracBrowser for help on using the repository browser.