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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90 @ 8373

Last change on this file since 8373 was 8373, checked in by clem, 7 years ago

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