source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

  • Property svn:keywords set to Id
File size: 22.3 KB
Line 
1MODULE traadv_ubs
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_ubs  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  1.0  !  2006-08  (L. Debreu, R. Benshila)  Original code
7   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_adv_ubs : update the tracer trend with the horizontal
12   !!                 advection trends using a third order biaised scheme 
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and active tracers
15   USE dom_oce         ! ocean space and time domain
16   USE trdmod_oce         ! ocean space and time domain
17   USE trdtra
18   USE lib_mpp
19   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
20   USE in_out_manager  ! I/O manager
21   USE diaptr          ! poleward transport diagnostics
22   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
23   USE trc_oce         ! share passive tracers/Ocean variables
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   tra_adv_ubs   ! routine called by traadv module
29
30   LOGICAL :: l_trd  ! flag to compute trends or not
31
32   !! * Control permutation of array indices
33#  include "oce_ftrans.h90"
34#  include "dom_oce_ftrans.h90"
35#  include "trc_oce_ftrans.h90"
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn,      &
48      &                                       ptb, ptn, pta, kjpt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE tra_adv_ubs  ***
51      !!                 
52      !! ** Purpose :   Compute the now trend due to the advection of tracers
53      !!      and add it to the general trend of passive tracer equations.
54      !!
55      !! ** Method  :   The upstream biased third (UBS) is order scheme based
56      !!      on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
57      !!      It is only used in the horizontal direction.
58      !!      For example the i-component of the advective fluxes are given by :
59      !!                !  e1u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
60      !!          zwx = !  or
61      !!                !  e1u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
62      !!      where zltu is the second derivative of the before temperature field:
63      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
64      !!      This results in a dissipatively dominant (i.e. hyper-diffusive)
65      !!      truncation error. The overall performance of the advection scheme
66      !!      is similar to that reported in (Farrow and Stevens, 1995).
67      !!      For stability reasons, the first term of the fluxes which corresponds
68      !!      to a second order centered scheme is evaluated using the now velocity
69      !!      (centered in time) while the second term which is the diffusive part
70      !!      of the scheme, is evaluated using the before velocity (forward in time).
71      !!      Note that UBS is not positive. Do not use it on passive tracers.
72      !!      On the vertical, the advection is evaluated using a TVD scheme, as
73      !!      the UBS have been found to be too diffusive.
74      !!
75      !! ** Action : - update (pta) with the now advective tracer trends
76      !!
77      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
78      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
79      !!----------------------------------------------------------------------
80      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
81      USE oce     , ONLY:   zwx  => ua       , zwy  => va         ! (ua,va) used as workspace
82      USE wrk_nemo, ONLY:   ztu  => wrk_3d_1 , ztv  => wrk_3d_2   ! 3D workspace
83      USE wrk_nemo, ONLY:   zltu => wrk_3d_3 , zltv => wrk_3d_4   !  -      -
84      USE wrk_nemo, ONLY:   zti  => wrk_3d_5 , ztw  => wrk_3d_6   !  -      -
85
86      !! DCSE_NEMO: need additional directives for renamed module variables
87!FTRANS zwx zwy ztu ztv zltu zltv zti ztw :I :I :z
88
89      !
90      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
91      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
92      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
93      REAL(wp), DIMENSION(       jpkorig  ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
94
95!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
96!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
97!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
98
99!FTRANS pun pvn pwn :I :I :z
100!FTRANS ptb ptn pta :I :I :z :
101      REAL(wp), INTENT(in   ) ::   pun(jpi,jpj,jpkorig)        ! ocean velocity component (u)
102      REAL(wp), INTENT(in   ) ::   pvn(jpi,jpj,jpkorig)        ! ocean velocity component (v)
103      REAL(wp), INTENT(in   ) ::   pwn(jpi,jpj,jpkorig)        ! ocean velocity component (w)
104!! DCSE_NEMO: Next two arguments made inout to silence the cray compile,
105!! which rightly complains about the call to nonosc_v (which also has them
106!! as inout)
107      REAL(wp), INTENT(inout) ::   ptb(jpi,jpj,jpkorig,kjpt)   ! tracer fields (before)
108      REAL(wp), INTENT(inout) ::   ptn(jpi,jpj,jpkorig,kjpt)   ! tracer fields (now)
109      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpkorig,kjpt)   ! tracer trend
110
111      !
112      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
113      REAL(wp) ::   ztra, zbtr, zcoef, z2dtt                       ! local scalars
114      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
115      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
116      !!----------------------------------------------------------------------
117
118      IF( wrk_in_use(3, 1,2,3,4,5,6) )THEN
119         CALL ctl_stop('tra_adv_ubs: requested workspace arrays unavailable')   ;   RETURN
120      ENDIF
121
122      IF( kt == nit000 )  THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
125         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
126         !
127         l_trd = .FALSE.
128         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
129      ENDIF
130      !
131      !                                                          ! ===========
132      DO jn = 1, kjpt                                            ! tracer loop
133         !                                                       ! ===========
134         ! 1. Bottom value : flux set to zero
135         ! ----------------------------------
136         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
137         !                                             
138#if defined key_z_first
139         DO jj = 1, jpjm1
140            DO ji = 1, jpim1
141               DO jk = 1, jpkm1
142#else
143         DO jk = 1, jpkm1                                 ! Horizontal slab
144            !                                   
145            !  Laplacian
146            DO jj = 1, jpjm1            ! First derivative (gradient)
147               DO ji = 1, fs_jpim1   ! vector opt.
148#endif
149                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
150                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
151                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
152                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
153               END DO
154            END DO
155#if defined key_z_first
156         END DO
157         DO jj = 2, jpjm1               ! Second derivative (divergence)
158            DO ji = 2, jpim1
159               DO jk = 1, jpkm1
160#else
161            DO jj = 2, jpjm1            ! Second derivative (divergence)
162               DO ji = fs_2, fs_jpim1   ! vector opt.
163#endif
164                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
165                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
166                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
167               END DO
168            END DO
169            !                                   
170         END DO                                           ! End of slab         
171         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
172
173         !   
174         !  Horizontal advective fluxes               
175#if defined key_z_first
176         DO jj = 1, jpjm1
177            DO ji = 1, jpim1
178               DO jk = 1, jpkm1
179#else
180         DO jk = 1, jpkm1                                 ! Horizontal slab
181            DO jj = 1, jpjm1
182               DO ji = 1, fs_jpim1   ! vector opt.
183#endif
184                  ! upstream transport
185                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
186                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
187                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
188                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
189                  ! centered scheme
190                  zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
191                  zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
192                  ! UBS scheme
193                  zwx(ji,jj,jk) =  zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 
194                  zwy(ji,jj,jk) =  zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) 
195               END DO
196            END DO
197         END DO                                           ! End of slab         
198
199         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
200
201         ! Horizontal advective trends
202#if defined key_z_first
203         DO jj = 2, jpjm1
204            DO ji = 2, jpim1
205               DO jk = 1, jpkm1
206#else
207         DO jk = 1, jpkm1
208            !  Tracer flux divergence at t-point added to the general trend
209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1   ! vector opt.
211#endif
212                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
213                  ! horizontal advective
214                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
215                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
216                  ! add it to the general tracer trends
217                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
218               END DO
219            END DO
220            !                                             
221         END DO                                           !   End of slab
222
223         ! Horizontal trend used in tra_adv_ztvd subroutine
224         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
225
226         ! 3. Save the horizontal advective trends for diagnostic
227         ! ------------------------------------------------------
228         !                                 ! trend diagnostics (contribution of upstream fluxes)
229         IF( l_trd ) THEN
230             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
231             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
232         END IF
233         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
234         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
235            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
236            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
237         ENDIF
238         
239         ! TVD scheme for the vertical direction 
240         ! ----------------------
241         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
242
243         !  Bottom value : flux set to zero
244         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
245
246         ! Surface value
247         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
248         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
249         ENDIF
250         !  upstream advection with initial mass fluxes & intermediate update
251         ! -------------------------------------------------------------------
252         ! Interior value
253#if defined key_z_first
254         DO jj = 1, jpj
255            DO ji = 1, jpi
256               DO jk = 2, jpkm1
257#else
258         DO jk = 2, jpkm1
259            DO jj = 1, jpj
260               DO ji = 1, jpi
261#endif
262                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
263                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
264                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
265               END DO
266            END DO
267         END DO 
268         ! update and guess with monotonic sheme
269#if defined key_z_first
270         DO jj = 2, jpjm1
271            DO ji = 2, jpim1
272               DO jk = 1, jpkm1
273                  z2dtt = p2dt(jk)
274#else
275         DO jk = 1, jpkm1
276            z2dtt = p2dt(jk)
277            DO jj = 2, jpjm1
278               DO ji = fs_2, fs_jpim1   ! vector opt.
279#endif
280                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
281                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
282                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
283                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
284               END DO
285            END DO
286         END DO
287         !
288         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
289
290         !  antidiffusive flux : high order minus low order
291#if defined key_z_first
292         DO jj = 1, jpj
293            DO ji = 1, jpi
294               ztw(ji,jj,1) = 0.e0   ! Surface value
295               DO jk = 2, jpkm1      ! Interior value
296#else
297         ztw(:,:,1) = 0.e0       ! Surface value
298         DO jk = 2, jpkm1        ! Interior value
299            DO jj = 1, jpj
300               DO ji = 1, jpi
301#endif
302                  ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk)
303               END DO
304            END DO
305         END DO
306         !
307         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
308
309         !  final trend with corrected fluxes
310#if defined key_z_first
311         DO jj = 2, jpjm1 
312            DO ji = 2, jpim1
313               DO jk = 1, jpkm1
314#else
315         DO jk = 1, jpkm1
316            DO jj = 2, jpjm1 
317               DO ji = fs_2, fs_jpim1   ! vector opt.   
318#endif
319                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
320                  ! k- vertical advective trends 
321                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
322                  ! added to the general tracer trends
323                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
324               END DO
325            END DO
326         END DO
327
328         !  Save the final vertical advective trends
329         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
330            ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
331#if defined key_z_first
332            DO jj = 2, jpjm1
333               DO ji = 2, jpim1
334                  DO jk = 1, jpkm1
335#else
336            DO jk = 1, jpkm1
337               DO jj = 2, jpjm1
338                  DO ji = fs_2, fs_jpim1   ! vector opt.
339#endif
340                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
341                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
342                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
343                  END DO
344               END DO
345            END DO
346            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv )
347         ENDIF
348         !
349      ENDDO
350      !
351      IF( wrk_not_released(3, 1,2,3,4,5,6) )   CALL ctl_stop('tra_adv_ubs: failed to release workspace arrays')
352      !
353
354      !! * Reset control of array index permutation
355!FTRANS CLEAR
356#  include "oce_ftrans.h90"
357#  include "dom_oce_ftrans.h90"
358#  include "trc_oce_ftrans.h90"
359
360   END SUBROUTINE tra_adv_ubs
361
362
363   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
364      !!---------------------------------------------------------------------
365      !!                    ***  ROUTINE nonosc_z  ***
366      !!     
367      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
368      !!       scheme and the before field by a nonoscillatory algorithm
369      !!
370      !! **  Method  :   ... ???
371      !!       warning : pbef and paft must be masked, but the boundaries
372      !!       conditions on the fluxes are not necessary zalezak (1979)
373      !!       drange (1995) multi-dimensional forward-in-time and upstream-
374      !!       in-space based differencing for fluid
375      !!----------------------------------------------------------------------
376      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
377      USE wrk_nemo, ONLY:   zbetup => wrk_3d_1, zbetdo => wrk_3d_2   ! 3D workspace
378
379      !! DCSE_NEMO: need additional directives for renamed module variables
380!FTRANS zbetup zbetdo :I :I :z
381
382      !
383      REAL(wp), INTENT(in   ), DIMENSION(jpkorig)          ::   p2dt   ! vertical profile of tracer time-step
384
385      !! DCSE_NEMO: This style defeats ftrans
386!     REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
387!     REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
388!     REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
389
390!FTRANS pbef paft pcc :I :I :z
391      REAL(wp), INTENT(inout) ::   pbef(jpi,jpj,jpkorig)   ! before field
392      REAL(wp), INTENT(inout) ::   paft(jpi,jpj,jpkorig)   ! after field
393      REAL(wp), INTENT(inout) ::   pcc(jpi,jpj,jpkorig)    ! monotonic flux in the k direction
394      !
395      INTEGER  ::   ji, jj, jk   ! dummy loop indices
396      INTEGER  ::   ikm1         ! local integer
397      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
398      !!----------------------------------------------------------------------
399
400      IF( wrk_in_use(3, 1,2) ) THEN
401         CALL ctl_stop('nonosc_z: requested workspace arrays unavailable')   ;   RETURN
402      ENDIF
403
404      zbig  = 1.e+40_wp
405      zrtrn = 1.e-15_wp
406      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
407
408      ! Search local extrema
409      ! --------------------
410      ! large negative value (-zbig) inside land
411      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
412      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
413      ! search maximum in neighbourhood
414#if defined key_z_first
415      DO jj = 2, jpjm1
416         DO ji = 2, jpim1
417            DO jk = 1, jpkm1
418               ikm1 = MAX(jk-1,1)
419#else
420      DO jk = 1, jpkm1
421         ikm1 = MAX(jk-1,1)
422         DO jj = 2, jpjm1
423            DO ji = fs_2, fs_jpim1   ! vector opt.
424#endif
425               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
426                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
427                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
428            END DO
429         END DO
430      END DO
431      ! large positive value (+zbig) inside land
432      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
433      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
434      ! search minimum in neighbourhood
435#if defined key_z_first
436      DO jj = 2, jpjm1
437         DO ji = 2, jpim1
438            DO jk = 1, jpkm1
439               ikm1 = MAX(jk-1,1)
440#else
441      DO jk = 1, jpkm1
442         ikm1 = MAX(jk-1,1)
443         DO jj = 2, jpjm1
444            DO ji = fs_2, fs_jpim1   ! vector opt.
445#endif
446               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
447                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
448                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
449            END DO
450         END DO
451      END DO
452
453      ! restore masked values to zero
454      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
455      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
456
457
458      ! 2. Positive and negative part of fluxes and beta terms
459      ! ------------------------------------------------------
460
461#if defined key_z_first
462      DO jj = 2, jpjm1
463         DO ji = 2, jpim1
464            DO jk = 1, jpkm1
465               z2dtt = p2dt(jk)
466#else
467      DO jk = 1, jpkm1
468         z2dtt = p2dt(jk)
469         DO jj = 2, jpjm1
470            DO ji = fs_2, fs_jpim1   ! vector opt.
471#endif
472               ! positive & negative part of the flux
473               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
474               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
475               ! up & down beta terms
476               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
477               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
478               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
479            END DO
480         END DO
481      END DO
482      ! monotonic flux in the k direction, i.e. pcc
483      ! -------------------------------------------
484#if defined key_z_first
485      DO jj = 2, jpjm1
486         DO ji = 2, jpim1
487            DO jk = 2, jpkm1
488#else
489      DO jk = 2, jpkm1
490         DO jj = 2, jpjm1
491            DO ji = fs_2, fs_jpim1   ! vector opt.
492#endif
493               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
494               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
495               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
496               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
497            END DO
498         END DO
499      END DO
500      !
501      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('nonosc_z: failed to release workspace arrays')
502      !
503   END SUBROUTINE nonosc_z
504
505   !!======================================================================
506END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.