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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90 @ 7512

Last change on this file since 7512 was 7037, checked in by mocavero, 8 years ago

ORCA2_LIM_PISCES hybrid version update
  • Property svn:keywords set to Id
File size: 23.9 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 wrk_nemo       ! work arrays
26   USE prtctl         ! Print control
27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   lim_adv_x_2   ! called by lim_trp
33   PUBLIC   lim_adv_y_2   ! called by lim_trp
34
35   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
36   REAL(wp)  ::   rzero  = 0.e0     !    -       -
37   REAL(wp)  ::   rone   = 1.e0     !    -       -
38   
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE lim_adv_x_2( pdf, put , pcrh, psm , ps0 ,   &
50      &                    psx, psxx, psy , psyy, psxy )
51      !!---------------------------------------------------------------------
52      !!                **  routine lim_adv_x_2  **
53      !! 
54      !! ** purpose :   Computes and adds the advection trend to sea-ice
55      !!              variable on i-axis
56      !!
57      !! ** method  :   Uses Prather second order scheme that advects tracers
58      !!              but also theirquadratic forms. The method preserves
59      !!              tracer structures by conserving second order moments.
60      !!
61      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
62      !!--------------------------------------------------------------------
63      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
64      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
65      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
66      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
67      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
68      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
69      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
70      !
71      INTEGER  ::   ji, jj                               ! dummy loop indices
72      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
73      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
74      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
75      REAL(wp), DIMENSION(:,:), POINTER ::   zf0, zfx , zfy , zbet   ! 2D workspace
76      REAL(wp), DIMENSION(:,:), POINTER ::   zfm, zfxx, zfyy, zfxy   !  -      -
77      REAL(wp), DIMENSION(:,:), POINTER ::   zalg, zalg1, zalg1q     !  -      -
78      !---------------------------------------------------------------------
79
80      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
81      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
82
83      ! Limitation of moments.                                           
84
85      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
86
87!$OMP PARALLEL
88!$OMP DO schedule(static) private(jj,ji,zslpmax,zs1max,zs1new,zs2new,zin0)
89      DO jj = 1, jpj
90         DO ji = 1, jpi
91            zslpmax = MAX( rzero, ps0(ji,jj) )
92            zs1max  = 1.5 * zslpmax
93            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
94            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
95               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
96            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
97            !
98            ps0 (ji,jj) = zslpmax 
99            psx (ji,jj) = zs1new      * zin0
100            psxx(ji,jj) = zs2new      * zin0
101            psy (ji,jj) = psy (ji,jj) * zin0
102            psyy(ji,jj) = psyy(ji,jj) * zin0
103            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
104         END DO
105      END DO
106!$OMP END DO NOWAIT
107
108      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
109!$OMP DO schedule(static) private(jj,ji)
110      DO jj = 1, jpj
111         DO ji = 1, jpi
112            psm (ji,jj)  = MAX( pcrh * area(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj) , epsi20 )
113         END DO
114      END DO
115
116      !  Calculate fluxes and moments between boxes i<-->i+1             
117!$OMP DO schedule(static) private(jj,ji,zalf,zalfq,zalf1,zalf1q)
118      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
119         DO ji = 1, jpi
120            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) )
121            zalf         =  MAX( rzero, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
122            zalfq        =  zalf * zalf
123            zalf1        =  1.0 - zalf
124            zalf1q       =  zalf1 * zalf1
125            !
126            zfm (ji,jj)  =  zalf  * psm(ji,jj)
127            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) )
128            zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
129            zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj)
130            zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) )
131            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
132            zfyy(ji,jj)  =  zalf  * psyy(ji,jj)
133            !
134            !  Readjust moments remaining in the box.
135            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
136            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
137            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
138            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
139            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
140            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
141            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
142         END DO
143      END DO
144!$OMP END DO
145
146!$OMP DO schedule(static) private(jj,ji,zalf,zalfq,zalf1,zalf1q)
147      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
148         DO ji = 1, fs_jpim1
149            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
150            zalg  (ji,jj) = zalf
151            zalfq         = zalf * zalf
152            zalf1         = 1.0 - zalf
153            zalg1 (ji,jj) = zalf1
154            zalf1q        = zalf1 * zalf1
155            zalg1q(ji,jj) = zalf1q
156            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj)
157            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
158            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
159            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj)
160            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) )
161            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)
162            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj)
163         END DO
164      END DO
165
166!$OMP DO schedule(static) private(jj,ji,zbt,zbt1)
167      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
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) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
172            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
173            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
174            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
175            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
176            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
177            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
178         END DO
179      END DO
180
181      !   Put the temporary moments into appropriate neighboring boxes.   
182!$OMP DO schedule(static) private(jj,ji,zbt,zbt1,zalf,zalf1,ztemp)
183      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
184         DO ji = fs_2, fs_jpim1
185            zbt  =       zbet(ji-1,jj)
186            zbt1 = 1.0 - zbet(ji-1,jj)
187            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
188            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
189            zalf1       = 1.0 - zalf
190            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
191            !
192            ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj)
193            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
194            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
195               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
196               &        + zbt1 * psxx(ji,jj)
197            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
198               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
199               &         + zbt1 * psxy(ji,jj)
200            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
201            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
202         END DO
203      END DO
204
205!$OMP DO schedule(static) private(jj,ji,zbt,zbt1,zalf,zalf1,ztemp)
206      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
207         DO ji = fs_2, fs_jpim1
208            zbt  =       zbet(ji,jj)
209            zbt1 = 1.0 - zbet(ji,jj)
210            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
211            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
212            zalf1       = 1.0 - zalf
213            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
214            !
215            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
216            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
217            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)   &
218               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )        &
219               &                                      + ( zalf1 - zalf ) * ztemp )                                 )
220            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)          &
221               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
222            psy(ji,jj)  = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
223            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
224         END DO
225      END DO
226!$OMP END DO NOWAIT
227!$OMP END PARALLEL
228
229      !-- Lateral boundary conditions
230      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
231      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
232      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
233      CALL lbc_lnk( psxy, 'T',  1. )
234
235      IF(ln_ctl) THEN
236         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
237         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
238         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
239         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
240      ENDIF
241      !
242      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
243      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
244      !
245   END SUBROUTINE lim_adv_x_2
246
247
248   SUBROUTINE lim_adv_y_2( pdf, pvt , pcrh, psm , ps0 ,   &
249      &                    psx, psxx, psy , psyy, psxy )
250      !!---------------------------------------------------------------------
251      !!                **  routine lim_adv_y_2  **
252      !!           
253      !! ** purpose :   Computes and adds the advection trend to sea-ice
254      !!              variable on j-axis
255      !!
256      !! ** method  :   Uses Prather second order scheme that advects tracers
257      !!              but also their quadratic forms. The method preserves
258      !!              tracer structures by conserving second order moments.
259      !!
260      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
261      !!---------------------------------------------------------------------
262      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
263      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
264      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
265      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
266      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
267      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
268      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
269      !!
270      INTEGER  ::   ji, jj                               ! dummy loop indices
271      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
272      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
273      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
274      REAL(wp), DIMENSION(:,:), POINTER ::   zf0, zfx , zfy , zbet   ! 2D workspace
275      REAL(wp), DIMENSION(:,:), POINTER ::   zfm, zfxx, zfyy, zfxy   !  -      -
276      REAL(wp), DIMENSION(:,:), POINTER ::   zalg, zalg1, zalg1q     !  -      -
277      !---------------------------------------------------------------------
278
279      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
280      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
281
282      ! Limitation of moments.
283
284      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
285
286!$OMP PARALLEL
287!$OMP DO schedule(static) private(jj,ji,zslpmax,zs1max,zs1new,zs2new,zin0)
288      DO jj = 1, jpj
289         DO ji = 1, jpi
290            zslpmax = MAX( rzero, ps0(ji,jj) )
291            zs1max  = 1.5 * zslpmax
292            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
293            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
294               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
295            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
296            !
297            ps0 (ji,jj) = zslpmax 
298            psx (ji,jj) = psx (ji,jj) * zin0
299            psxx(ji,jj) = psxx(ji,jj) * zin0
300            psy (ji,jj) = zs1new * zin0
301            psyy(ji,jj) = zs2new * zin0
302            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
303         END DO
304      END DO
305!$OMP END DO NOWAIT
306
307      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
308!$OMP DO schedule(static) private(jj,ji,zslpmax,zs1max,zs1new,zs2new,zin0)
309      DO jj = 1, jpj
310         DO ji = 1, jpi
311            psm(ji,jj)  = MAX(  pcrh * area(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj) , epsi20  )
312         END DO
313      END DO
314
315      !  Calculate fluxes and moments between boxes j<-->j+1             
316!$OMP DO schedule(static) private(jj,ji,zalf,zalfq,zalf1,zalf1q)
317      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
318         DO ji = 1, jpi
319            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
320            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
321            zalfq        =  zalf * zalf
322            zalf1        =  1.0 - zalf
323            zalf1q       =  zalf1 * zalf1
324            zfm (ji,jj)  =  zalf  * psm(ji,jj)
325            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
326            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
327            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
328            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
329            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
330            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
331            !
332            !  Readjust moments remaining in the box.
333            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
334            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
335            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
336            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
337            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
338            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
339            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
340         END DO
341      END DO
342      !
343!$OMP DO schedule(static) private(jj,ji,zalf,zalfq,zalf1,zalf1q)
344      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
345         DO ji = 1, jpi
346            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
347            zalg  (ji,jj) = zalf
348            zalfq         = zalf * zalf
349            zalf1         = 1.0 - zalf
350            zalg1 (ji,jj) = zalf1
351            zalf1q        = zalf1 * zalf1
352            zalg1q(ji,jj) = zalf1q
353            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1)
354            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
355            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
356            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1)
357            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) )
358            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1)
359            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1)
360         END DO
361      END DO
362 
363      !  Readjust moments remaining in the box.
364!$OMP DO schedule(static) private(jj,ji,zbt,zbt1)
365      DO jj = 2, jpj
366         DO ji = 1, jpi
367            zbt  =         zbet(ji,jj-1)
368            zbt1 = ( 1.0 - zbet(ji,jj-1) )
369            !
370            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
371            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
372            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
373            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj)
374            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
375            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
376            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
377         END DO
378      END DO
379
380      !   Put the temporary moments into appropriate neighboring boxes.   
381!$OMP DO schedule(static) private(jj,ji,zbt,zbt1,zalf,zalf1,ztemp)
382      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
383         DO ji = 1, jpi
384            zbt  =         zbet(ji,jj-1)
385            zbt1 = ( 1.0 - zbet(ji,jj-1) )
386            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
387            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
388            zalf1       = 1.0 - zalf
389            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
390            !
391            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj)
392            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
393               &        + zbt1 * psy(ji,jj) 
394            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
395               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
396               &        + zbt1 * psyy(ji,jj)
397            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              &
398               &                 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   &
399               &        + zbt1 * psxy(ji,jj)
400            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
401            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
402         END DO
403      END DO
404      !
405!$OMP DO schedule(static) private(jj,ji,zbt,zbt1,zalf,zalf1,ztemp)
406      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
407         DO ji = 1, jpi
408            zbt  =         zbet(ji,jj)
409            zbt1 = ( 1.0 - zbet(ji,jj) )
410            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
411            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
412            zalf1       = 1.0 - zalf
413            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
414            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
415            psy(ji,jj)  = zbt  * psy(ji,jj)  &
416               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
417            psyy(ji,jj) = zbt  * psyy(ji,jj)                                                                         &
418               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)                           &
419               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) )
420            psxy(ji,jj) = zbt  * psxy(ji,jj)                                  &
421               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   &
422               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) )
423            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
424            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
425         END DO
426      END DO
427!$OMP END DO NOWAIT
428!$OMP END PARALLEL
429
430      !-- Lateral boundary conditions
431      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
432      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
433      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
434      CALL lbc_lnk( psxy, 'T',  1. )
435
436      IF(ln_ctl) THEN
437         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
438         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
439         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
440         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
441      ENDIF
442      !
443      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
444      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
445      !
446   END SUBROUTINE lim_adv_y_2
447
448#else
449   !!----------------------------------------------------------------------
450   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
451   !!----------------------------------------------------------------------
452#endif
453
454   !!======================================================================
455END MODULE limadv_2
Note: See TracBrowser for help on using the repository browser.