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

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

changes in style - part4 - clarify ice advection routines
File size: 33.9 KB
Line 
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.6  !  2014-11  (C. Rousset, G. Madec)  Original code
7   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                    LIM 3.0 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
13   !!   ultimate_x(_y): compute a tracer value at velocity points using ULTIMATE scheme at various orders
14   !!   macho         : ???
15   !!   nonosc_2d     : compute monotonic tracer fluxes by a non-oscillatory algorithm
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constant
18   USE dom_oce        ! ocean domain
19   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition
20   USE ice            ! sea-ice variables
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     
33   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
34   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
35
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
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, pu_ice, pv_ice,  &
46      &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE ice_adv_umx  ***
49      !!
50      !! **  Purpose :   Compute the now trend due to total advection of
51      !!                 tracers and add it to the general trend of tracer equations
52      !!                 using an "Ultimate-Macho" scheme
53      !!
54      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
55      !!----------------------------------------------------------------------
56      INTEGER                     , INTENT(in   ) ::   kt         ! time step
57      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
58      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
59      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
60      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
61      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
62      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content
63      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
64      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
65      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
66      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
67      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
68      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
69      !
70      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
71      INTEGER  ::   initad                  ! number of sub-timestep for the advection
72      REAL(wp) ::   zcfl , zusnit, zdt      !   -      -
73      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box
74      !!----------------------------------------------------------------------
75      !
76      IF( kt == nit000 .AND. lwp ) THEN
77         WRITE(numout,*) 
78         WRITE(numout,*) 'ice_adv_umx : Ultimate-MACHO advection scheme'
79         WRITE(numout,*) '~~~~~~~~~~~'
80      ENDIF
81      !
82      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
83      !
84      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !       
85      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
86      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
87      IF( lk_mpp )   CALL mpp_max( zcfl )
88
89      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
90      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
91      ENDIF
92
93      zdt = rdt_ice / REAL(initad)
94
95      ! --- transport --- !
96      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
97      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
98
99      ! --- define velocity for advection: u*grad(H) --- !
100      DO jj = 2, jpjm1
101         DO ji = fs_2, fs_jpim1
102            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
103            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
104            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
105            ENDIF
106
107            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
108            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
109            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
110            ENDIF
111         END DO
112      END DO
113
114      !---------------!
115      !== advection ==!
116      !---------------!
117      DO jt = 1, initad
118         CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area
119         DO jl = 1, jpl
120            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
121            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
122            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content
123            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content
124            DO jk = 1, nlay_i
125               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
126            END DO
127            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
128            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
129            IF ( nn_pnd_scheme > 0 ) THEN
130               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
131               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
132            ENDIF
133         END DO
134      END DO
135      !
136      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
137      !
138   END SUBROUTINE ice_adv_umx
139   
140   SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
141      !!----------------------------------------------------------------------
142      !!                  ***  ROUTINE adv_umx  ***
143      !!
144      !! **  Purpose :   Compute the now trend due to total advection of
145      !!       tracers and add it to the general trend of tracer equations
146      !!
147      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
148      !!       corrected flux (monotonic correction)
149      !!       note: - this advection scheme needs a leap-frog time scheme
150      !!
151      !! ** Action : - pt  the after advective tracer
152      !!----------------------------------------------------------------------
153      INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration
154      REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step
155      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2
156      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
157      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field
158      !
159      INTEGER  ::   ji, jj           ! dummy loop indices 
160      REAL(wp) ::   ztra             ! local scalar
161      REAL(wp) ::   zfp_ui, zfp_vj   !   -      -
162      REAL(wp) ::   zfm_ui, zfm_vj   !   -      -
163      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups
164      REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd
165      !!----------------------------------------------------------------------
166      !
167      IF( nn_timing == 1 )  CALL timing_start('ice_adv_umx')
168      !
169      !  upstream advection with initial mass fluxes & intermediate update
170      ! --------------------------------------------------------------------
171      DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction
172         DO ji = 1, fs_jpim1   ! vector opt.
173            zfp_ui = puc(ji,jj) + ABS( puc(ji,jj) )
174            zfm_ui = puc(ji,jj) - ABS( puc(ji,jj) )
175            zfp_vj = pvc(ji,jj) + ABS( pvc(ji,jj) )
176            zfm_vj = pvc(ji,jj) - ABS( pvc(ji,jj) )
177            zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * ptc(ji,jj) + zfm_ui * ptc(ji+1,jj  ) )
178            zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * ptc(ji,jj) + zfm_vj * ptc(ji  ,jj+1) )
179         END DO
180      END DO
181     
182      DO jj = 2, jpjm1            ! total intermediate advective trends
183         DO ji = fs_2, fs_jpim1   ! vector opt.
184            ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
185               &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj)
186            !
187            ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ] 
188            zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme
189         END DO
190      END DO
191      CALL lbc_lnk( zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
192     
193      ! High order (_ho) fluxes
194      ! -----------------------
195      SELECT CASE( nn_limadv_ord )
196      CASE ( 20 )                          ! centered second order
197         DO jj = 2, jpjm1
198            DO ji = fs_2, fs_jpim1   ! vector opt.
199               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) )
200               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) )
201            END DO
202         END DO
203         !
204      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
205         CALL macho( kt, nn_limadv_ord, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
206         !
207         DO jj = 2, jpjm1
208            DO ji = fs_2, fs_jpim1   ! vector opt.
209               zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj)
210               zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj)
211            END DO
212         END DO
213         !
214      END SELECT
215         
216      ! antidiffusive flux : high order minus low order
217      ! --------------------------------------------------
218      DO jj = 2, jpjm1
219         DO ji = fs_2, fs_jpim1   ! vector opt.
220            zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)
221            zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj)
222         END DO
223      END DO
224      CALL lbc_lnk_multi( zfu_ho, 'U', -1., zfv_ho, 'V', -1. )         ! Lateral bondary conditions
225     
226      ! monotonicity algorithm
227      ! -------------------------
228      CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )
229     
230      ! final trend with corrected fluxes
231      ! ------------------------------------
232      DO jj = 2, jpjm1
233         DO ji = fs_2, fs_jpim1   ! vector opt. 
234            ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   &
235               &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj) 
236            ptc(ji,jj) = ptc(ji,jj) + pdt * ztra
237         END DO
238      END DO
239      CALL lbc_lnk( ptc(:,:) , 'T',  1. )
240      !
241      IF( nn_timing == 1 )  CALL timing_stop('ice_adv_umx')
242      !
243   END SUBROUTINE adv_umx
244
245
246   SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
247      !!---------------------------------------------------------------------
248      !!                    ***  ROUTINE ultimate_x  ***
249      !!     
250      !! **  Purpose :   compute 
251      !!
252      !! **  Method  :   ... ???
253      !!                 TIM = transient interpolation Modeling
254      !!
255      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
256      !!----------------------------------------------------------------------
257      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
258      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
259      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
260      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields
261      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components
262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
263      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points
264      !
265      INTEGER  ::   ji, jj    ! dummy loop indices
266      REAL(wp) ::   zc_box    !   -      -
267      REAL(wp), DIMENSION(jpi,jpj) :: zzt
268      !!----------------------------------------------------------------------
269      !
270      IF( nn_timing == 1 )  CALL timing_start('macho')
271      !
272      IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==!
273         !
274         !                                                           !--  ultimate interpolation of pt at u-point  --!
275         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )
276         !
277         !                                                           !--  advective form update in zzt  --!
278         DO jj = 2, jpjm1
279            DO ji = fs_2, fs_jpim1   ! vector opt.
280               zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
281                  &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj)
282               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
283            END DO
284         END DO
285         CALL lbc_lnk( zzt, 'T', 1. )
286         !
287         !                                                           !--  ultimate interpolation of pt at v-point  --!
288         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )
289         !
290      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==!
291         !
292         !                                                           !--  ultimate interpolation of pt at v-point  --!
293         CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v )
294         !
295         !                                                           !--  advective form update in zzt  --!
296         DO jj = 2, jpjm1
297            DO ji = fs_2, fs_jpim1
298               zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
299                  &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj)
300               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
301            END DO
302         END DO
303         CALL lbc_lnk( zzt, 'T', 1. )
304         !
305         !                                                           !--  ultimate interpolation of pt at u-point  --!
306         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )
307         !     
308      ENDIF     
309      !
310      IF( nn_timing == 1 )  CALL timing_stop('macho')
311      !
312   END SUBROUTINE macho
313
314
315   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )
316      !!---------------------------------------------------------------------
317      !!                    ***  ROUTINE ultimate_x  ***
318      !!     
319      !! **  Purpose :   compute 
320      !!
321      !! **  Method  :   ... ???
322      !!                 TIM = transient interpolation Modeling
323      !!
324      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
325      !!----------------------------------------------------------------------
326      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
327      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
328      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component
329      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
330      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point
331      !
332      INTEGER  ::   ji, jj       ! dummy loop indices
333      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
334      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4
335      !!----------------------------------------------------------------------
336      !
337      IF( nn_timing == 1 )  CALL timing_start('ultimate_x')
338      !
339      !                                                     !--  Laplacian in i-direction  --!
340      DO jj = 2, jpjm1         ! First derivative (gradient)
341         DO ji = 1, fs_jpim1
342            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
343         END DO
344         !                     ! Second derivative (Laplacian)
345         DO ji = fs_2, fs_jpim1
346            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
347         END DO
348      END DO
349      CALL lbc_lnk( ztu2, 'T', 1. )
350      !
351      !                                                     !--  BiLaplacian in i-direction  --!
352      DO jj = 2, jpjm1         ! Third derivative
353         DO ji = 1, fs_jpim1
354            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
355         END DO
356         !                     ! Fourth derivative
357         DO ji = fs_2, fs_jpim1
358            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
359         END DO
360      END DO
361      CALL lbc_lnk( ztu4, 'T', 1. )
362      !
363      !
364      SELECT CASE (k_order )
365      !
366      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
367         !       
368         DO jj = 1, jpj
369            DO ji = 1, fs_jpim1   ! vector opt.
370               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   &
371                  &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
372            END DO
373         END DO
374         !
375      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
376         !
377         DO jj = 1, jpj
378            DO ji = 1, fs_jpim1   ! vector opt.
379               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
380               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
381                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
382            END DO
383         END DO
384         CALL lbc_lnk( pt_u(:,:) , 'U',  1. )
385         ! 
386      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
387         !
388         DO jj = 1, jpj
389            DO ji = 1, fs_jpim1   ! vector opt.
390               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
391               zdx2 = e1u(ji,jj) * e1u(ji,jj)
392!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
393               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
394                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
395                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
396                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
397            END DO
398         END DO
399         !
400      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
401         !
402         DO jj = 1, jpj
403            DO ji = 1, fs_jpim1   ! vector opt.
404               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
405               zdx2 = e1u(ji,jj) * e1u(ji,jj)
406!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
407               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
408                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
409                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
410                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
411            END DO
412         END DO
413         !
414      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
415         !
416         DO jj = 1, jpj
417            DO ji = 1, fs_jpim1   ! vector opt.
418               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
419               zdx2 = e1u(ji,jj) * e1u(ji,jj)
420!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
421               zdx4 = zdx2 * zdx2
422               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
423                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
424                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
425                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
426                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
427                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
428            END DO
429         END DO
430         !
431      END SELECT
432      !
433      IF( nn_timing == 1 )  CALL timing_stop('ultimate_x')
434      !
435   END SUBROUTINE ultimate_x
436   
437 
438   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )
439      !!---------------------------------------------------------------------
440      !!                    ***  ROUTINE ultimate_y  ***
441      !!     
442      !! **  Purpose :   compute 
443      !!
444      !! **  Method  :   ... ???
445      !!                 TIM = transient interpolation Modeling
446      !!
447      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
448      !!----------------------------------------------------------------------
449      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
450      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
451      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component
452      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
453      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point
454      !
455      INTEGER  ::   ji, jj       ! dummy loop indices
456      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
457      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4
458      !!----------------------------------------------------------------------
459      !
460      IF( nn_timing == 1 )  CALL timing_start('ultimate_y')
461      !
462      !                                                     !--  Laplacian in j-direction  --!
463      DO jj = 1, jpjm1         ! First derivative (gradient)
464         DO ji = fs_2, fs_jpim1
465            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
466         END DO
467      END DO
468      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
469         DO ji = fs_2, fs_jpim1
470            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
471         END DO
472      END DO
473      CALL lbc_lnk( ztv2, 'T', 1. )
474      !
475      !                                                     !--  BiLaplacian in j-direction  --!
476      DO jj = 1, jpjm1         ! First derivative
477         DO ji = fs_2, fs_jpim1
478            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
479         END DO
480      END DO
481      DO jj = 2, jpjm1         ! Second derivative
482         DO ji = fs_2, fs_jpim1
483            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
484         END DO
485      END DO
486      CALL lbc_lnk( ztv4, 'T', 1. )
487      !
488      !
489      SELECT CASE (k_order )
490      !
491      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
492         DO jj = 1, jpjm1
493            DO ji = 1, jpi
494               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  &
495                  &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
496            END DO
497         END DO
498         !
499      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
500         DO jj = 1, jpjm1
501            DO ji = 1, jpi
502               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
503               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
504                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
505            END DO
506         END DO
507         CALL lbc_lnk( pt_v(:,:) , 'V',  1. )
508         !
509      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
510         DO jj = 1, jpjm1
511            DO ji = 1, jpi
512               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
513               zdy2 = e2v(ji,jj) * e2v(ji,jj)
514!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
515               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
516                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
517                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
518                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
519            END DO
520         END DO
521         !
522      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
523         DO jj = 1, jpjm1
524            DO ji = 1, jpi
525               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
526               zdy2 = e2v(ji,jj) * e2v(ji,jj)
527!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
528               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
529                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
530                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
531                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
532            END DO
533         END DO
534         !
535      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
536         DO jj = 1, jpjm1
537            DO ji = 1, jpi
538               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
539               zdy2 = e2v(ji,jj) * e2v(ji,jj)
540!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
541               zdy4 = zdy2 * zdy2
542               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
543                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
544                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
545                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
546                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
547                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
548            END DO
549         END DO
550         !
551      END SELECT
552      !
553      IF( nn_timing == 1 )  CALL timing_stop('ultimate_y')
554      !
555   END SUBROUTINE ultimate_y
556   
557 
558   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )
559      !!---------------------------------------------------------------------
560      !!                    ***  ROUTINE nonosc  ***
561      !!     
562      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
563      !!       scheme and the before field by a nonoscillatory algorithm
564      !!
565      !! **  Method  :   ... ???
566      !!       warning : pbef and paft must be masked, but the boundaries
567      !!       conditions on the fluxes are not necessary zalezak (1979)
568      !!       drange (1995) multi-dimensional forward-in-time and upstream-
569      !!       in-space based differencing for fluid
570      !!----------------------------------------------------------------------
571      REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step
572      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field
573      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions
574      !
575      INTEGER  ::   ji, jj    ! dummy loop indices
576      INTEGER  ::   ikm1      ! local integer
577      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars
578      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
579      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv
580      !!----------------------------------------------------------------------
581      !
582      IF( nn_timing == 1 )  CALL timing_start('nonosc_2d')
583      !
584      zbig = 1.e+40_wp
585      zsml = 1.e-15_wp
586
587      ! clem test
588      DO jj = 2, jpjm1
589         DO ji = fs_2, fs_jpim1   ! vector opt. 
590            zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   &
591               &              + pbb(ji,jj) - pbb(ji  ,jj-1) ) 
592         END DO
593      END DO
594      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
595
596      ! Determine ice masks for before and after tracers
597      WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp
598      ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1)
599      END WHERE
600
601      ! Search local extrema
602      ! --------------------
603      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
604!      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   &
605!         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  )
606!      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   &
607!         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  )
608      zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
609         &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
610      zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
611         &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
612
613      z1_dt = 1._wp / pdt
614      DO jj = 2, jpjm1
615         DO ji = fs_2, fs_jpim1   ! vector opt.
616            !
617            zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood
618               &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    )
619            zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   &
620               &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    )
621               !
622            zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux
623               & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) )
624            zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   &
625               & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) )
626               !
627            zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
628            zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt
629            zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt
630         END DO
631      END DO
632      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
633
634      ! monotonic flux in the i & j direction (paa & pbb)
635      ! -------------------------------------
636      DO jj = 2, jpjm1
637         DO ji = fs_2, fs_jpim1   ! vector opt.
638            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
639            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
640            zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) )
641            !
642            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
643            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
644            zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) )
645            !
646            paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu )
647            pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv )
648            !
649         END DO
650      END DO
651      CALL lbc_lnk_multi( paa, 'U', -1., pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
652      !
653      IF( nn_timing == 1 )  CALL timing_stop('nonosc_2d')
654      !
655   END SUBROUTINE nonosc_2d
656
657#else
658   !!----------------------------------------------------------------------
659   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
660   !!----------------------------------------------------------------------
661#endif
662
663   !!======================================================================
664END MODULE iceadv_umx
Note: See TracBrowser for help on using the repository browser.