New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_umx.F90 in NEMO/branches/UKMO/dev_r10037_GPU/src/ICE – NEMO

source: NEMO/branches/UKMO/dev_r10037_GPU/src/ICE/icedyn_adv_umx.F90 @ 11467

Last change on this file since 11467 was 11467, checked in by andmirek, 5 years ago

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