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

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

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