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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90 @ 9421

Last change on this file since 9421 was 9421, checked in by clem, 6 years ago

cosmetics

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