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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90 @ 5429

Last change on this file since 5429 was 5429, checked in by smasson, 9 years ago

merge dev_r5302_CNRS18_HPC_scalability into the trunk, see #1523
  • Property svn:keywords set to Id
File size: 22.6 KB
RevLine 
[825]1MODULE limadv 
2   !!======================================================================
3   !!                       ***  MODULE limadv   ***
4   !! LIM sea-ice model : sea-ice advection
5   !!======================================================================
[1530]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.
[2715]8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
[1530]9   !!--------------------------------------------------------------------
[825]10#if defined key_lim3
11   !!----------------------------------------------------------------------
[834]12   !!   'key_lim3'                                     LIM3 sea-ice model
[825]13   !!----------------------------------------------------------------------
14   !!   lim_adv_x  : advection of sea ice on x axis
15   !!   lim_adv_y  : advection of sea ice on y axis
16   !!----------------------------------------------------------------------
[4161]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
[825]26
27   IMPLICIT NONE
28   PRIVATE
29
[1530]30   PUBLIC   lim_adv_x   ! called by lim_trp
31   PUBLIC   lim_adv_y   ! called by lim_trp
[825]32
[868]33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
[825]35   !!----------------------------------------------------------------------
[4161]36   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]37   !! $Id$
[2715]38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[825]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
[1530]48      !!              variable on i-axis
[825]49      !!
[1530]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.
[825]55      !!--------------------------------------------------------------------
[1530]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
[5123]65      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
[2715]66      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
67      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
[3294]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         !  -      -
[825]71      !---------------------------------------------------------------------
72
[3294]73      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
74      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
[2715]75
[825]76      ! Limitation of moments.                                           
77
[1530]78      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
[825]79
80      DO jj = 1, jpj
81         DO ji = 1, jpi
[4688]82            zslpmax = MAX( 0._wp, ps0(ji,jj) )
[825]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) )  )
[5123]87            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
[825]88
89            ps0 (ji,jj) = zslpmax 
[5123]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
[825]95         END DO
96      END DO
97
98      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
[5123]99      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
[825]100
101      !  Calculate fluxes and moments between boxes i<-->i+1             
[868]102      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
[825]103         DO ji = 1, jpi
[4688]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)
[825]106            zalfq        =  zalf * zalf
107            zalf1        =  1.0 - zalf
108            zalf1q       =  zalf1 * zalf1
[1530]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)
[825]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
[868]129      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
130         DO ji = 1, fs_jpim1
[4688]131            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
[825]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
[1530]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)
[825]146         END DO
147      END DO
148
149      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
[1529]150         DO ji = fs_2, fs_jpim1
[825]151            zbt  =       zbet(ji-1,jj)
152            zbt1 = 1.0 - zbet(ji-1,jj)
[1530]153            !
[825]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.
[868]166         DO ji = fs_2, fs_jpim1
[825]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)
[1530]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)                               &
[825]177               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
[1530]178               &                                                + zbt1 * psxx(ji,jj)
[825]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) ) )   &
[1530]181               &                                                + zbt1 * psxy(ji,jj)
[825]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.
[868]188         DO ji = fs_2, fs_jpim1
[825]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
[1530]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) )  )
[825]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      !-- Lateral boundary conditions
[5429]209      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
210         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
211         &              , psxx, 'T',  1., psyy, 'T',  1.   &
212         &              , psxy, 'T',  1. )
[825]213
[921]214      IF(ln_ctl) THEN
[825]215         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
216         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
217         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
218         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
[921]219      ENDIF
[1530]220      !
[3294]221      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
222      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
[2715]223      !
[825]224   END SUBROUTINE lim_adv_x
225
226
227   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
228      &                  psx, psxx, psy , psyy, psxy )
229      !!---------------------------------------------------------------------
230      !!                **  routine lim_adv_y  **
231      !!           
232      !! ** purpose :   Computes and adds the advection trend to sea-ice
[1530]233      !!              variable on y axis
[825]234      !!
235      !! ** method  :   Uses Prather second order scheme that advects tracers
[1530]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
[825]248      !!
[1530]249      INTEGER  ::   ji, jj                               ! dummy loop indices
[5123]250      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
[1530]251      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
252      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
[3294]253      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
254      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
255      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
[825]256      !---------------------------------------------------------------------
257
[3294]258      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
259      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
[2715]260
[825]261      ! Limitation of moments.
262
263      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
264
[921]265      DO jj = 1, jpj
266         DO ji = 1, jpi
[4688]267            zslpmax = MAX( 0._wp, ps0(ji,jj) )
[921]268            zs1max  = 1.5 * zslpmax
269            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
270            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
271               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
[5123]272            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
[1530]273            !
[921]274            ps0 (ji,jj) = zslpmax 
[5123]275            psx (ji,jj) = psx (ji,jj) * rswitch
276            psxx(ji,jj) = psxx(ji,jj) * rswitch
277            psy (ji,jj) = zs1new * rswitch
278            psyy(ji,jj) = zs2new * rswitch
279            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
[921]280         END DO
281      END DO
[825]282
[921]283      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
[5123]284      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
[825]285
[921]286      !  Calculate fluxes and moments between boxes j<-->j+1             
[1530]287      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
[921]288         DO ji = 1, jpi
[4688]289            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
290            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
[921]291            zalfq        =  zalf * zalf
292            zalf1        =  1.0 - zalf
293            zalf1q       =  zalf1 * zalf1
[1530]294            !
[921]295            zfm (ji,jj)  =  zalf  * psm(ji,jj)
296            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
297            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
298            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
299            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
300            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
301            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
[1530]302            !
[921]303            !  Readjust moments remaining in the box.
304            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
305            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
306            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
307            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
308            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
309            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
310            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
311         END DO
312      END DO
[1530]313      !
[921]314      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
315         DO ji = 1, jpi
[4688]316            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
[921]317            zalg  (ji,jj) = zalf
318            zalfq         = zalf * zalf
319            zalf1         = 1.0 - zalf
320            zalg1 (ji,jj) = zalf1
321            zalf1q        = zalf1 * zalf1
322            zalg1q(ji,jj) = zalf1q
[1530]323            !
324            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
325            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
326            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
327            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
328            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
329            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
330            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
[921]331         END DO
332      END DO
[825]333
[921]334      !  Readjust moments remaining in the box.
335      DO jj = 2, jpj
336         DO ji = 1, jpi
337            zbt  =         zbet(ji,jj-1)
338            zbt1 = ( 1.0 - zbet(ji,jj-1) )
[1530]339            !
[921]340            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
341            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
342            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
[1530]343            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
[921]344            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
345            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
346            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
347         END DO
348      END DO
[825]349
[921]350      !   Put the temporary moments into appropriate neighboring boxes.   
351      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
352         DO ji = 1, jpi
353            zbt  =         zbet(ji,jj-1)
354            zbt1 = ( 1.0 - zbet(ji,jj-1) )
355            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
356            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
357            zalf1       = 1.0 - zalf
358            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
[1530]359            !
360            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
[921]361            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
[1530]362               &                                               + zbt1 * psy(ji,jj) 
[921]363            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
364               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
[1530]365               &                                               + zbt1 * psyy(ji,jj)
366            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
367               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
368               &                                                + zbt1 * psxy(ji,jj)
[921]369            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
370            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
371         END DO
372      END DO
[825]373
[921]374      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
375         DO ji = 1, jpi
376            zbt  =         zbet(ji,jj)
377            zbt1 = ( 1.0 - zbet(ji,jj) )
378            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
379            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
380            zalf1       = 1.0 - zalf
[1530]381            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
382            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
383            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
384            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
385               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
386               &                                         + ( zalf1 - zalf ) * ztemp )                                )
387            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
388               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
389            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
390            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
[921]391         END DO
392      END DO
393
[825]394      !-- Lateral boundary conditions
[5429]395      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
396         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
397         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
398         &              , psxy, 'T',  1. )
[825]399
[921]400      IF(ln_ctl) THEN
[825]401         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
402         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
403         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
404         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
[921]405      ENDIF
[1530]406      !
[3294]407      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
408      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
[2715]409      !
[825]410   END SUBROUTINE lim_adv_y
411
412#else
413   !!----------------------------------------------------------------------
414   !!   Default option            Dummy module         NO LIM sea-ice model
415   !!----------------------------------------------------------------------
416#endif
417
[1530]418   !!======================================================================
[825]419END MODULE limadv
Note: See TracBrowser for help on using the repository browser.