New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limadv_2.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
  • Property svn:keywords set to Id
File size: 23.0 KB
RevLine 
[821]1MODULE limadv_2 
[3]2   !!======================================================================
[821]3   !!                       ***  MODULE limadv_2   ***
4   !! LIM 2.0 sea-ice model : sea-ice advection
[3]5   !!======================================================================
[1530]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   !!--------------------------------------------------------------------
[821]12#if defined key_lim2
[3]13   !!----------------------------------------------------------------------
[821]14   !!   'key_lim2'                                    LIM 2.0 sea-ice model
[88]15   !!----------------------------------------------------------------------
[821]16   !!   lim_adv_x_2  : advection of sea ice on x axis
17   !!   lim_adv_y_2  : advection of sea ice on y axis
[3]18   !!----------------------------------------------------------------------
19   USE dom_oce
[821]20   USE dom_ice_2
[1465]21   USE ice_2
[3]22   USE lbclnk
[1530]23   USE in_out_manager     ! I/O manager
[2636]24   USE lib_mpp            ! MPP library
[1530]25   USE prtctl             ! Print control
[3]26
27   IMPLICIT NONE
28   PRIVATE
29
[1530]30   PUBLIC   lim_adv_x_2   ! called by lim_trp
31   PUBLIC   lim_adv_y_2   ! called by lim_trp
[3]32
[1530]33   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
34   REAL(wp)  ::   rzero  = 0.e0     !    -       -
35   REAL(wp)  ::   rone   = 1.e0     !    -       -
36   
[1529]37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
[3]39   !!----------------------------------------------------------------------
[2528]40   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
[1156]41   !! $Id$
[2528]42   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]43   !!----------------------------------------------------------------------
44
45CONTAINS
46
[821]47   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
48      &                    psx, psxx, psy , psyy, psxy )
[3]49      !!---------------------------------------------------------------------
[821]50      !!                **  routine lim_adv_x_2  **
[3]51      !! 
52      !! ** purpose :   Computes and adds the advection trend to sea-ice
[1530]53      !!              variable on i-axis
[3]54      !!
[1530]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.
[3]60      !!--------------------------------------------------------------------
[2633]61      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
[2613]62      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
63      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
64      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
65      !
[1530]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
[2613]73      !
[1530]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          !    -         -
[3]78      !---------------------------------------------------------------------
79
[2633]80      IF( wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
[2613]81         CALL ctl_stop( 'lim_adv_x_2 : requested workspace arrays unavailable.' )   ;   RETURN
[2636]82      ENDIF
[2590]83
[3]84      ! Limitation of moments.                                           
85
[1530]86      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
[3]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
[1530]96            !
[3]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             
[1529]110      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
[3]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
[1530]117            !
[3]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)
[1530]125            !
[3]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
[1529]137      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
138         DO ji = 1, fs_jpim1
[3]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.
[1529]157         DO ji = fs_2, fs_jpim1
[3]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.
[1529]172         DO ji = fs_2, fs_jpim1
[3]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)
[1530]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)                               &
[3]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.
[1529]194         DO ji = fs_2, fs_jpim1
[3]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)
[1530]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) )
[3]211         END DO
212      END DO
213
214      !-- Lateral boundary conditions
[1530]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. )
[1510]218      CALL lbc_lnk( psxy, 'T',  1. )
[3]219
[1529]220      IF(ln_ctl) THEN
[258]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 :')
[3]225      ENDIF
[1530]226      !
[2636]227      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) )   &
228          CALL ctl_stop( 'lim_adv_x_2 : failed to release workspace arrays.' )
[2590]229      !
[821]230   END SUBROUTINE lim_adv_x_2
[3]231
232
[821]233   SUBROUTINE lim_adv_y_2( pdf, pvt , pcrh, psm , ps0 ,   &
234      &                    psx, psxx, psy , psyy, psxy )
[3]235      !!---------------------------------------------------------------------
[821]236      !!                **  routine lim_adv_y_2  **
[3]237      !!           
238      !! ** purpose :   Computes and adds the advection trend to sea-ice
[1530]239      !!              variable on j-axis
[3]240      !!
241      !! ** method  :   Uses Prather second order scheme that advects tracers
[1530]242      !!              but also their quadratic forms. The method preserves
243      !!              tracer structures by conserving second order moments.
[3]244      !!
[1530]245      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
[3]246      !!---------------------------------------------------------------------
[2633]247      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
[2613]248      USE wrk_nemo, ONLY:   zf0  => wrk_2d_11 , zfx   => wrk_2d_12 , zfy    => wrk_2d_13 , zbet => wrk_2d_14   ! 2D workspace
249      USE wrk_nemo, ONLY:   zfm  => wrk_2d_15 , zfxx  => wrk_2d_16 , zfyy   => wrk_2d_17 , zfxy => wrk_2d_18   !  -      -
250      USE wrk_nemo, ONLY:   zalg => wrk_2d_19 , zalg1 => wrk_2d_20 , zalg1q => wrk_2d_21                       !  -      -
[2590]251      !!
[1530]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          !    -         -
[3]264      !---------------------------------------------------------------------
265
[2633]266      IF(wrk_in_use(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
[2613]267         CALL ctl_stop( 'lim_adv_y_2 : requested workspace arrays unavailable.' )   ;   RETURN
[2590]268      END IF
269
[3]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
[1530]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
[3]291
[1530]292      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
293      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
[3]294
[1530]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
[3]340 
[1530]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
[3]356
[1530]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
[3]402
403      !-- Lateral boundary conditions
[1530]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. )
[1510]407      CALL lbc_lnk( psxy, 'T',  1. )
[3]408
[258]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 :')
[3]414      ENDIF
[1530]415      !
[2633]416      IF( wrk_not_released(2, 11,12,13,14,15,16,17,18,19,20,21) ) THEN
[2613]417        CALL ctl_stop( 'lim_adv_y_2 : failed to release workspace arrays.' )
[2590]418      END IF
419      !
[821]420   END SUBROUTINE lim_adv_y_2
[3]421
422#else
[88]423   !!----------------------------------------------------------------------
[821]424   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
[88]425   !!----------------------------------------------------------------------
[3]426#endif
427
[1530]428   !!======================================================================
[821]429END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.