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_umx.F90 in branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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