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

Last change on this file since 10180 was 10180, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 4b: reduce communications in si3, see #2133

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