source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90 @ 9604

Last change on this file since 9604 was 9604, checked in by clem, 3 years ago

change history of the ice routines

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