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

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

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

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