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

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90 @ 6515

Last change on this file since 6515 was 6515, checked in by clem, 8 years ago

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