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 @ 6746

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

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