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

Last change on this file since 8376 was 8376, checked in by clem, 5 years ago

to make the reference configuration orca2aghulas work

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