source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

File size: 29.2 KB
Line 
1MODULE limadv_umx
2   !!==============================================================================
3   !!                       ***  MODULE  limadv_umx  ***
4   !! LIM sea-ice model : sea-ice advection using the ULTIMATE-MACHO scheme
5   !!==============================================================================
6   !! History :  3.5  !  2014-11  (C. Rousset, G. Madec)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   lim_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
11   !!   ultimate      : compute a tracer value at velocity points using ULTIMATE scheme at various orders
12   !!   macho         :
13   !!   nonosc_2d     : compute monotonic tracer fluxes by a non-oscillatory algorithm
14   !!----------------------------------------------------------------------
15   USE phycst         ! physical constant
16   USE dom_oce        ! ocean domain
17   USE sbc_oce        ! ocean surface boundary condition
18   USE ice            ! ice variables
19   !
20   USE in_out_manager ! I/O manager
21   USE lbclnk         ! lateral boundary conditions -- MPP exchanges
22   USE lib_mpp        ! MPP library
23   USE wrk_nemo       ! work arrays
24   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
25   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_adv_umx    ! routine called by limtrp.F90
31     
32   REAL(wp) ::   z1_6   = 1._wp / 6._wp   ! =1/6
33   REAL(wp) ::   z1_120 = 1._wp / 120._wp ! =1/120
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id: limadv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE lim_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE lim_adv_umx  ***
47      !!
48      !! **  Purpose :   Compute the now trend due to total advection of
49      !!       tracers and add it to the general trend of tracer equations
50      !!
51      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
52      !!       corrected flux (monotonic correction)
53      !!       note: - this advection scheme needs a leap-frog time scheme
54      !!
55      !! ** Action : - pt  the after advective tracer
56      !!----------------------------------------------------------------------
57      INTEGER                     , INTENT(in   )           ::   kt         ! number of iteration
58      REAL(wp)                    , INTENT(in   )           ::   pdt        ! tracer time-step
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components => u*e2
60      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc        ! tracer content field
62      !
63      INTEGER  ::   ji, jj           ! dummy loop indices 
64      REAL(wp) ::   ztra             ! local scalar
65      REAL(wp) ::   zfp_ui, zfp_vj   !   -      -
66      REAL(wp) ::   zfm_ui, zfm_vj   !   -      -
67      REAL(wp), POINTER, DIMENSION(:,:) :: zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v
68      !!----------------------------------------------------------------------
69      !
70      IF( nn_timing == 1 )  CALL timing_start('lim_adv_umx')
71      !
72      CALL wrk_alloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v )
73      !
74      !
75      !  upstream advection with initial mass fluxes & intermediate update
76      ! --------------------------------------------------------------------
77      DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction
78         DO ji = 1, fs_jpim1   ! vector opt.
79            zfp_ui = puc(ji,jj) + ABS( puc(ji,jj) )
80            zfm_ui = puc(ji,jj) - ABS( puc(ji,jj) )
81            zfp_vj = pvc(ji,jj) + ABS( pvc(ji,jj) )
82            zfm_vj = pvc(ji,jj) - ABS( pvc(ji,jj) )
83            zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * ptc(ji,jj) + zfm_ui * ptc(ji+1,jj  ) )
84            zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * ptc(ji,jj) + zfm_vj * ptc(ji  ,jj+1) )
85         END DO
86      END DO
87     
88      DO jj = 2, jpjm1            ! total intermediate advective trends
89         DO ji = fs_2, fs_jpim1   ! vector opt.
90            ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
91               &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj)
92            !
93            ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ] 
94            zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme
95         END DO
96      END DO
97      CALL lbc_lnk( zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
98     
99      ! High order (_ho) fluxes
100      ! -----------------------
101      SELECT CASE( nn_limadv_ord )
102      CASE ( 20 )                          ! centered second order
103         DO jj = 2, jpjm1
104            DO ji = fs_2, fs_jpim1   ! vector opt.
105               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) )
106               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) )
107            END DO
108         END DO
109         !
110      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
111         CALL macho( kt, nn_limadv_ord, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
112         !
113         DO jj = 2, jpjm1
114            DO ji = fs_2, fs_jpim1   ! vector opt.
115               zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj)
116               zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj)
117            END DO
118         END DO
119         !
120      END SELECT
121         
122      ! antidiffusive flux : high order minus low order
123      ! --------------------------------------------------
124      DO jj = 2, jpjm1
125         DO ji = fs_2, fs_jpim1   ! vector opt.
126            zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)
127            zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj)
128         END DO
129      END DO
130      CALL lbc_lnk_multi( zfu_ho, 'U', -1., zfv_ho, 'V', -1. )         ! Lateral bondary conditions
131     
132      ! monotonicity algorithm
133      ! -------------------------
134      CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )
135     
136      ! final trend with corrected fluxes
137      ! ------------------------------------
138      DO jj = 2, jpjm1
139         DO ji = fs_2, fs_jpim1   ! vector opt. 
140            ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   &
141               &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj) 
142            ptc(ji,jj) = ptc(ji,jj) + pdt * ztra
143         END DO
144      END DO
145      CALL lbc_lnk( ptc(:,:) , 'T',  1. )
146      !
147      !
148      CALL wrk_dealloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v )
149      !
150      IF( nn_timing == 1 )  CALL timing_stop('lim_adv_umx')
151      !
152   END SUBROUTINE lim_adv_umx
153
154
155   SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
156      !!---------------------------------------------------------------------
157      !!                    ***  ROUTINE ultimate_x  ***
158      !!     
159      !! **  Purpose :   compute 
160      !!
161      !! **  Method  :   ... ???
162      !!                 TIM = transient interpolation Modeling
163      !!
164      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
165      !!----------------------------------------------------------------------
166      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
167      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
168      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
169      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields
170      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components
171      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
172      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points
173      !
174      INTEGER  ::   ji, jj    ! dummy loop indices
175      REAL(wp) ::   zc_box    !   -      -
176      REAL(wp), POINTER, DIMENSION(:,:) :: zzt
177      !!----------------------------------------------------------------------
178      !
179      IF( nn_timing == 1 )  CALL timing_start('macho')
180      !
181      CALL wrk_alloc( jpi,jpj,   zzt )
182      !
183      IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==!
184         !
185         !                                                           !--  ultimate interpolation of pt at u-point  --!
186         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )
187         !
188         !                                                           !--  advective form update in zzt  --!
189         DO jj = 2, jpjm1
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
192                  &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj)
193               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
194            END DO
195         END DO
196         CALL lbc_lnk( zzt, 'T', 1. )
197         !
198         !                                                           !--  ultimate interpolation of pt at v-point  --!
199         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )
200         !
201      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==!
202         !
203         !                                                           !--  ultimate interpolation of pt at v-point  --!
204         CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v )
205         !
206         !                                                           !--  advective form update in zzt  --!
207         DO jj = 2, jpjm1
208            DO ji = fs_2, fs_jpim1
209               zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
210                  &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj)
211               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
212            END DO
213         END DO
214         CALL lbc_lnk( zzt, 'T', 1. )
215         !
216         !                                                           !--  ultimate interpolation of pt at u-point  --!
217         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )
218         !     
219      ENDIF     
220      !
221      CALL wrk_dealloc( jpi,jpj,   zzt )
222      !
223      IF( nn_timing == 1 )  CALL timing_stop('macho')
224      !
225   END SUBROUTINE macho
226
227
228   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )
229      !!---------------------------------------------------------------------
230      !!                    ***  ROUTINE ultimate_x  ***
231      !!     
232      !! **  Purpose :   compute 
233      !!
234      !! **  Method  :   ... ???
235      !!                 TIM = transient interpolation Modeling
236      !!
237      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
238      !!----------------------------------------------------------------------
239      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
240      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
241      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component
242      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
243      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point
244      !
245      INTEGER  ::   ji, jj       ! dummy loop indices
246      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
247      REAL(wp), POINTER, DIMENSION(:,:) :: ztu1, ztu2, ztu3, ztu4
248      !!----------------------------------------------------------------------
249      !
250      IF( nn_timing == 1 )  CALL timing_start('ultimate_x')
251      !
252      CALL wrk_alloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 )
253      !
254      !                                                     !--  Laplacian in i-direction  --!
255      DO jj = 2, jpjm1         ! First derivative (gradient)
256         DO ji = 1, fs_jpim1
257            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
258         END DO
259         !                     ! Second derivative (Laplacian)
260         DO ji = fs_2, fs_jpim1
261            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
262         END DO
263      END DO
264      CALL lbc_lnk( ztu2, 'T', 1. )
265      !
266      !                                                     !--  BiLaplacian in i-direction  --!
267      DO jj = 2, jpjm1         ! Third derivative
268         DO ji = 1, fs_jpim1
269            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
270         END DO
271         !                     ! Fourth derivative
272         DO ji = fs_2, fs_jpim1
273            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
274         END DO
275      END DO
276      CALL lbc_lnk( ztu4, 'T', 1. )
277      !
278      !
279      SELECT CASE (k_order )
280      !
281      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
282         !       
283         DO jj = 1, jpj
284            DO ji = 1, fs_jpim1   ! vector opt.
285               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   &
286                  &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
287            END DO
288         END DO
289         !
290      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
291         !
292         DO jj = 1, jpj
293            DO ji = 1, fs_jpim1   ! vector opt.
294               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
295               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
296                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
297            END DO
298         END DO
299         CALL lbc_lnk( pt_u(:,:) , 'U',  1. )
300         
301      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
302         !
303         DO jj = 1, jpj
304            DO ji = 1, fs_jpim1   ! vector opt.
305               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
306               zdx2 = e1u(ji,jj) * e1u(ji,jj)
307!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
308               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
309                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
310                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
311                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
312            END DO
313         END DO
314         !
315      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
316         !
317         DO jj = 1, jpj
318            DO ji = 1, fs_jpim1   ! vector opt.
319               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
320               zdx2 = e1u(ji,jj) * e1u(ji,jj)
321!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
322               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
323                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
324                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
325                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
326            END DO
327         END DO
328         !
329      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
330         !
331         DO jj = 1, jpj
332            DO ji = 1, fs_jpim1   ! vector opt.
333               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
334               zdx2 = e1u(ji,jj) * e1u(ji,jj)
335!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
336               zdx4 = zdx2 * zdx2
337               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
338                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
339                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
340                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
341                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
342                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
343            END DO
344         END DO
345         !
346      END SELECT
347      !
348      CALL wrk_dealloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 )
349      !
350      IF( nn_timing == 1 )  CALL timing_stop('ultimate_x')
351      !
352   END SUBROUTINE ultimate_x
353   
354 
355   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )
356      !!---------------------------------------------------------------------
357      !!                    ***  ROUTINE ultimate_y  ***
358      !!     
359      !! **  Purpose :   compute 
360      !!
361      !! **  Method  :   ... ???
362      !!                 TIM = transient interpolation Modeling
363      !!
364      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
365      !!----------------------------------------------------------------------
366      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
367      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
368      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component
369      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
370      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point
371      !
372      INTEGER  ::   ji, jj       ! dummy loop indices
373      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
374      REAL(wp), POINTER, DIMENSION(:,:) :: ztv1, ztv2, ztv3, ztv4
375      !!----------------------------------------------------------------------
376      !
377      IF( nn_timing == 1 )  CALL timing_start('ultimate_y')
378      !
379      CALL wrk_alloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 )
380      !
381      !                                                     !--  Laplacian in j-direction  --!
382      DO jj = 1, jpjm1         ! First derivative (gradient)
383         DO ji = fs_2, fs_jpim1
384            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
385         END DO
386      END DO
387      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
388         DO ji = fs_2, fs_jpim1
389            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
390         END DO
391      END DO
392      CALL lbc_lnk( ztv2, 'T', 1. )
393      !
394      !                                                     !--  BiLaplacian in j-direction  --!
395      DO jj = 1, jpjm1         ! First derivative
396         DO ji = fs_2, fs_jpim1
397            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
398         END DO
399      END DO
400      DO jj = 2, jpjm1         ! Second derivative
401         DO ji = fs_2, fs_jpim1
402            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
403         END DO
404      END DO
405      CALL lbc_lnk( ztv4, 'T', 1. )
406      !
407      !
408      SELECT CASE (k_order )
409         !
410      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
411         !       
412         DO jj = 1, jpjm1
413            DO ji = 1, jpi
414               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  &
415                  &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
416            END DO
417         END DO
418         !
419      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
420         DO jj = 1, jpjm1
421            DO ji = 1, jpi
422               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
423               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
424                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
425            END DO
426         END DO
427         CALL lbc_lnk( pt_v(:,:) , 'V',  1. )
428         !
429      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
430         !
431         DO jj = 1, jpjm1
432            DO ji = 1, jpi
433               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
434               zdy2 = e2v(ji,jj) * e2v(ji,jj)
435!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
436               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
437                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
438                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
439                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
440            END DO
441         END DO
442         !
443      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
444         !
445         DO jj = 1, jpjm1
446            DO ji = 1, jpi
447               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
448               zdy2 = e2v(ji,jj) * e2v(ji,jj)
449!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
450               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
451                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
452                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
453                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
454            END DO
455         END DO
456         !
457      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
458         !
459         DO jj = 1, jpjm1
460            DO ji = 1, jpi
461               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
462               zdy2 = e2v(ji,jj) * e2v(ji,jj)
463!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
464               zdy4 = zdy2 * zdy2
465               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
466                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
467                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
468                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
469                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
470                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
471            END DO
472         END DO
473         !
474      END SELECT
475      !
476      CALL wrk_dealloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 )
477      !
478      IF( nn_timing == 1 )  CALL timing_stop('ultimate_y')
479      !
480   END SUBROUTINE ultimate_y
481   
482 
483   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )
484      !!---------------------------------------------------------------------
485      !!                    ***  ROUTINE nonosc  ***
486      !!     
487      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
488      !!       scheme and the before field by a nonoscillatory algorithm
489      !!
490      !! **  Method  :   ... ???
491      !!       warning : pbef and paft must be masked, but the boundaries
492      !!       conditions on the fluxes are not necessary zalezak (1979)
493      !!       drange (1995) multi-dimensional forward-in-time and upstream-
494      !!       in-space based differencing for fluid
495      !!----------------------------------------------------------------------
496      REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step
497      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field
498      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions
499      !
500      INTEGER  ::   ji, jj    ! dummy loop indices
501      INTEGER  ::   ikm1      ! local integer
502      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars
503      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
504      REAL(wp), POINTER, DIMENSION(:,:) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv
505      !!----------------------------------------------------------------------
506      !
507      IF( nn_timing == 1 )  CALL timing_start('nonosc_2d')
508      !
509      CALL wrk_alloc( jpi,jpj,   zbetup, zbetdo, zbup, zbdo, zmsk, zdiv )
510      !
511      zbig = 1.e+40_wp
512      zsml = 1.e-15_wp
513
514      ! clem test
515      DO jj = 2, jpjm1
516         DO ji = fs_2, fs_jpim1   ! vector opt. 
517            zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   &
518               &              + pbb(ji,jj) - pbb(ji  ,jj-1) ) 
519         END DO
520      END DO
521      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
522
523      ! Determine ice masks for before and after tracers
524      WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp
525      ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1)
526      END WHERE
527
528      ! Search local extrema
529      ! --------------------
530      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
531!      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   &
532!         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  )
533!      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   &
534!         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  )
535      zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
536         &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
537      zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
538         &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
539
540      z1_dt = 1._wp / pdt
541      DO jj = 2, jpjm1
542         DO ji = fs_2, fs_jpim1   ! vector opt.
543            !
544            zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood
545               &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    )
546            zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   &
547               &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    )
548               !
549            zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux
550               & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) )
551            zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   &
552               & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) )
553               !
554            zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
555            zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt
556            zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt
557         END DO
558      END DO
559      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
560
561      ! monotonic flux in the i & j direction (paa & pbb)
562      ! -------------------------------------
563      DO jj = 2, jpjm1
564         DO ji = fs_2, fs_jpim1   ! vector opt.
565            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
566            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
567            zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) )
568            !
569            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
570            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
571            zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) )
572            !
573            paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu )
574            pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv )
575            !
576         END DO
577      END DO
578      CALL lbc_lnk_multi( paa, 'U', -1., pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
579      !
580      CALL wrk_dealloc( jpi,jpj,   zbetup, zbetdo, zbup, zbdo, zmsk, zdiv )
581      !
582      IF( nn_timing == 1 )  CALL timing_stop('nonosc_2d')
583      !
584   END SUBROUTINE nonosc_2d
585
586   !!======================================================================
587END MODULE limadv_umx
Note: See TracBrowser for help on using the repository browser.