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.
icedyn_adv_umx.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90 @ 10410

Last change on this file since 10410 was 10397, checked in by smasson, 5 years ago

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