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.
iceadv_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/iceadv_umx.F90 @ 8409

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

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