source: trunk/NEMO/TOP_SRC/TRP/trcbbl_adv.h90 @ 941

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

  • Property svn:executable set to *
File size: 20.3 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                     ***  trcbbl_adv.h90  ***
3   !!----------------------------------------------------------------------
4
5   !!----------------------------------------------------------------------
6   !!   TOP 1.0 , LOCEAN-IPSL (2005)
7   !! $Header$
8   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
9   !!----------------------------------------------------------------------
10
11   SUBROUTINE trc_bbl_adv( kt )
12      !!----------------------------------------------------------------------
13      !!                  ***  ROUTINE trc_bbl_adv  ***
14      !!                   
15      !! ** Purpose :   Compute the before tracer trend associated
16      !!     with the bottom boundary layer and add it to the general trend
17      !!     of tracer equations. The bottom boundary layer is supposed to be
18      !!     both an advective and diffusive bottom boundary layer.
19      !!
20      !! ** Method  :   Computes the bottom boundary horizontal and vertical
21      !!      advection terms. Add it to the general trend : tra =tra + adv_bbl.
22      !!        When the product grad( rho) * grad(h) < 0 (where grad is a
23      !!      along bottom slope gradient) an additional lateral 2nd order
24      !!      diffusion along the bottom slope is added to the general
25      !!      tracer trend, otherwise the additional trend is set to 0.
26      !!      Second order operator (laplacian type) with variable coefficient
27      !!      computed as follow for temperature (idem on s):
28      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
29      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
30      !!      where ztb is a 2D array: the bottom ocean te;perature and ahtb
31      !!      is a time and space varying diffusive coefficient defined by:
32      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
33      !!              = 0.       otherwise.
34      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
35      !!      is evaluated using the local density (i.e. referenced at the
36      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
37      !!      a downslope velocity of 20 cm/s if the condition for slope
38      !!      convection is satified)
39      !!        Add this before trend to the general trend tra of the
40      !!      botton ocean tracer point:
41      !!              tra = tra + difft
42      !!
43      !! ** Action  : - update tra at the bottom level with the bottom
44      !!                boundary layer trend
45      !!
46      !! References :
47      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
48      !!
49      !! History :
50      !!   8.5  !  02-12  (A. de Miranda, G. Madec)  Original Code
51      !!   9.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines )
52      !!   9.0  !  04-03  (C. Ethe) Adaptation for Passive tracers
53      !!----------------------------------------------------------------------     
54      !gh
55     
56      !! * Arguments
57      INTEGER, INTENT( in ) ::   kt        ! ocean time-step
58     
59      !! * Local declarations
60      INTEGER :: ji, jj, jk, jn            ! dummy loop indices
61      INTEGER :: ik, iku, ikv              ! temporary integers
62
63      REAL(wp) ::   &
64         zsign, zt, zs, zh, zalbet,     &  ! temporary scalars
65         zgdrho, zbtr, ztra                !    "         "
66      REAL(wp), DIMENSION(jpi,jpj) ::   & 
67          ztnb, zsnb, zdep, ztrb       ! temporary workspace arrays
68      REAL(wp), DIMENSION(jpi,jpj) ::   &  ! temporary workspace arrays
69         zalphax, zwu, zunb,            &  !    "                  "
70         zalphay, zwv, zvnb,            &  !    "                  "
71         zwx, zwy, zww, zwz,            &  !    "                  "
72         zti, zsi ,ztmin,ztmax, zsmin,zsmax!    "                  "
73                              !    "                  "
74      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
75         zhdivn                            ! temporary workspace arrays
76      REAL(wp) ::   &
77         zfui, zfvj, zbt,  zsigna,     &  ! temporary scalars
78    iku1,iku2,ikv1,ikv2,      &  ! temporary scalars
79    ze3u,ze3v,          &  ! temporary scalars
80    z2,z2dtt                  ! temporary scalars
81      REAL(wp) ::   &
82         fsalbt, pft, pfs, pfh             ! statement function
83      CHARACTER (len=22) :: charout
84      !!----------------------------------------------------------------------
85      ! ratio alpha/beta
86      ! ================
87      !  fsalbt: ratio of thermal over saline expension coefficients
88      !       pft :  potential temperature in degrees celcius
89      !       pfs :  salinity anomaly (s-35) in psu
90      !       pfh :  depth in meters
91
92      fsalbt( pft, pfs, pfh ) =                                              &
93         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
94                                   - 0.203814e-03 ) * pft                    &
95                                   + 0.170907e-01 ) * pft                    &
96                                   + 0.665157e-01                            &
97         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
98         +  ( ( - 0.302285e-13 * pfh                                         &
99                - 0.251520e-11 * pfs                                         &
100                + 0.512857e-12 * pft * pft          ) * pfh                  &
101                                     - 0.164759e-06   * pfs                  &
102             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
103                                     + 0.380374e-04 ) * pfh
104      !!----------------------------------------------------------------------
105
106   IF( kt == nit000 )   CALL trc_bbl_init    ! initialization at first time-step
107   
108      ! 1. 2D fields of bottom temperature and salinity, and bottom slope
109      ! -----------------------------------------------------------------
110      ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
111
112#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
113      jj = 1
114      DO ji = 1, jpij   ! vector opt. (forced unrolling)
115#else
116      DO jj = 1, jpj
117         DO ji = 1, jpi
118#endif
119            ik = mbkt(ji,jj)                    ! index of the bottom ocean T-level
120            ztnb(ji,jj) = tn(ji,jj,ik)    ! masked now T at the ocean bottom
121            zsnb(ji,jj) = sn(ji,jj,ik)    ! masked now S at the ocean bottom
122            zdep(ji,jj) = fsdept(ji,jj,ik)                 ! depth of the ocean bottom T-level
123!gh
124            zunb(ji,jj) = un(ji,jj,mbku(ji,jj))
125            zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
126#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
127         END DO
128#endif
129      END DO
130
131      ! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
132      ! --------------------------------------------
133      ! Sign of the local density gradient along the i- and j-slopes
134      ! multiplied by the slope of the ocean bottom
135
136      SELECT CASE ( neos )
137
138      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
139
140      DO jj = 1, jpjm1
141        DO ji = 1, fs_jpim1   ! vector opt.
142      !   ... temperature, salinity anomalie and depth
143          zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
144          zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
145          zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
146      !   ... masked ratio alpha/beta
147          zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
148      !   ... local density gradient along i-bathymetric slope
149          zgdrho = zalbet*( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
150                     -    ( zsnb(ji+1,jj) - zsnb(ji,jj) )
151          zgdrho = zgdrho * umask(ji,jj,1)
152      !   ... sign of local i-gradient of density multiplied by the i-slope
153          zsign = SIGN( 0.5, -zgdrho     * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
154     zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
155          zalphax(ji,jj) = ( 0.5 + zsigna ) * ( 0.5-zsign ) * umask(ji,jj,1)
156        END DO
157      END DO
158
159      DO jj = 1, jpjm1
160        DO ji = 1, fs_jpim1   ! vector opt.
161      !   ... temperature, salinity anomalie and depth
162          zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
163          zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
164          zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
165      !   ... masked ratio alpha/beta
166          zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
167      !   ... local density gradient along j-bathymetric slope
168          zgdrho = zalbet*( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
169                     -    ( zsnb(ji,jj+1) - zsnb(ji,jj) )
170          zgdrho = zgdrho * vmask(ji,jj,1)
171      !   ... sign of local j-gradient of density multiplied by the j-slope
172          zsign = SIGN( 0.5, -zgdrho     * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
173     zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
174          zalphay(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
175        END DO
176      END DO
177
178
179      CASE ( 1 )               ! Linear formulation function of temperature only
180
181      DO jj = 1, jpjm1
182        DO ji = 1, fs_jpim1   ! vector opt.
183            ! temperature, salinity anomalie and depth
184            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
185            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
186            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
187!gh            ! masked ratio alpha/beta
188            ! local density gradient along i-bathymetric slope
189            zgdrho =  ( ztnb(ji+1,jj) - ztnb(ji,jj) )   
190            ! sign of local i-gradient of density multiplied by the i-slope
191            zsign = SIGN( 0.5, - zgdrho    * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
192       zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
193            zalphax(ji,jj) = ( 0.5 - zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
194         END DO
195      END DO
196
197      DO jj = 1, jpjm1
198        DO ji = 1, fs_jpim1   ! vector opt.
199            ! temperature, salinity anomalie and depth
200            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
201            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
202            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
203!gh             ! masked ratio alpha/beta
204            ! local density gradient along j-bathymetric slope
205            zgdrho =  ( ztnb(ji,jj+1) - ztnb(ji,jj) )   
206            ! sign of local j-gradient of density multiplied by the j-slope
207            zsign = SIGN( 0.5, -zgdrho     * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
208       zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
209            zalphay(ji,jj) = ( 0.5 - zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
210         END DO
211      END DO
212       
213      CASE ( 2 )               ! Linear formulation function of temperature and salinity
214 DO jj = 1, jpjm1
215        DO ji = 1, fs_jpim1   ! vector opt.
216       ! local density gradient along i-bathymetric slope
217            zgdrho = - ( rbeta*( zsnb(ji+1,jj) - zsnb(ji,jj) )   &
218                     -  ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
219            ! sign of local i-gradient of density multiplied by the i-slope
220            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
221       zsigna= SIGN( 0.5, zunb(ji,jj)*(  zdep(ji+1,jj) - zdep(ji,jj) ))
222            zalphax(ji,jj)=(0.5-zsigna)*(0.5-zsign)*umask(ji,jj,1)
223        END DO
224      END DO
225
226      DO jj = 1, jpjm1
227        DO ji = 1, fs_jpim1   ! vector opt.
228      ! local density gradient along j-bathymetric slope
229            zgdrho = - ( rbeta*( zsnb(ji,jj+1) - zsnb(ji,jj) )   &
230                   -    ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )   
231            ! sign of local j-gradient of density multiplied by the j-slope
232            zsign = SIGN( 0.5, -zgdrho     * ( zdep(ji,jj+1) - zdep(ji,jj) ) ) 
233       zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
234            zalphay(ji,jj) = ( 0.5 - zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
235         END DO
236      END DO
237
238       
239      CASE DEFAULT
240
241         IF(lwp) WRITE(numout,cform_err)
242         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
243         nstop = nstop + 1
244
245      END SELECT
246
247      ! lateral boundary conditions on zalphax and zalphay   (unchanged sign)
248       CALL lbc_lnk( zalphax, 'U', 1. )   ;   CALL lbc_lnk( zalphay, 'V', 1. )
249
250
251      ! 3. Velocities that are exchanged between ajacent bottom boxes.
252      !---------------------------------------------------------------
253      ! ... is equal to zero but where bbl will work.
254       u_trc_bbl(:,:,:) = 0.e0
255       v_trc_bbl(:,:,:) = 0.e0
256       
257       
258!gh     
259       IF( ln_zps ) THEN
260      ! partial steps correction   
261     
262#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
263       jj = 1
264       DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
265#else
266       DO jj = 1, jpjm1
267          DO ji = 1, jpim1
268#endif
269             iku  = mbku(ji  ,jj  )
270             ikv  = mbkv(ji  ,jj  ) 
271             iku1 = mbkt(ji+1,jj  )
272             iku2 = mbkt(ji  ,jj  )
273             ikv1 = mbkt(ji  ,jj+1)
274             ikv2 = mbkt(ji  ,jj  )
275             ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
276             ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
277             
278             IF( MAX(iku,ikv) >  1 ) THEN
279                u_trc_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * ze3u / fse3u(ji,jj,iku)
280                v_trc_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * ze3v / fse3v(ji,jj,ikv)       
281             ENDIF
282#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
283          END DO
284#endif
285       END DO
286
287      ! lateral boundary conditions on u_trc_bbl and v_trc_bbl   (changed sign)
288       CALL lbc_lnk( u_trc_bbl, 'U', -1. )   ;   CALL lbc_lnk( v_trc_bbl, 'V', -1. )
289       
290    ELSE    ! z-coordinate - full steps or s-coordinate   
291       ! if not partial step loop over the whole domain no lbc call
292
293#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
294      jj = 1
295      DO ji = 1, jpij   ! vector opt. (forced unrolling)
296#else
297      DO jj = 1, jpj
298         DO ji = 1, jpi
299#endif
300            iku = mbku(ji,jj)
301            ikv = mbkv(ji,jj)
302            IF( MAX(iku,ikv) >  1 ) THEN
303               u_trc_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku)
304               v_trc_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv)       
305            ENDIF
306#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
307         END DO
308#endif
309          END DO
310       
311       ENDIF
312       
313       DO jn = 1, jptra
314
315#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
316          jj = 1
317          DO ji = 1, jpij   ! vector opt. (forced unrolling)
318#else
319          DO jj = 1, jpj
320             DO ji = 1, jpi
321#endif
322                ik = mbkt(ji,jj)                               ! index of the bottom ocean T-level
323                ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)    ! masked now T at the ocean bottom
324#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
325             END DO
326#endif
327          END DO
328
329
330
331          ! 5. Along sigma advective trend
332          ! -------------------------------
333          ! ... Second order centered tracer flux at u and v-points
334       
335# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
336          jj = 1
337          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
338# else
339          DO jj = 1, jpjm1
340             DO ji = 1, jpim1
341# endif
342                iku = mbku(ji,jj)
343                ikv = mbkv(ji,jj)
344                zfui = e2u(ji,jj) * fse3u(ji,jj,iku) * u_trc_bbl(ji,jj,iku)
345                zfvj = e1v(ji,jj) * fse3v(ji,jj,ikv) * v_trc_bbl(ji,jj,ikv)
346                ! upstream scheme
347                zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztrb(ji  ,jj  )   &
348                   &          +( zfui - ABS( zfui ) ) * ztrb(ji+1,jj  ) ) * 0.5
349                zwy(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * ztrb(ji  ,jj  )   &
350                   &          +( zfvj - ABS( zfvj ) ) * ztrb(ji  ,jj+1) ) * 0.5
351#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
352             END DO
353#endif
354          END DO
355     
356# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
357          jj = 1
358          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
359# else
360          DO jj = 2, jpjm1
361             DO ji = 2, jpim1
362# endif
363                ik = mbkt(ji,jj)
364                zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,ik) )
365                ! horizontal advective trends
366                ztra = - zbtr * (  zwx(ji,jj) - zwx(ji-1,jj  )   &
367                   &             + zwy(ji,jj) - zwy(ji  ,jj-1)  )
368 
369                ! add it to the general tracer trends
370                tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
371#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
372             END DO
373#endif
374          END DO
375     
376        END DO
377   
378        IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
379         WRITE(charout, FMT="('bbl - adv')")
380         CALL prt_ctl_trc_info(charout)
381         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
382      ENDIF         
383       
384       ! 6. Vertical advection velocities
385       ! --------------------------------
386       ! ... computes divergence perturbation (velocties to be removed from upper t boxes :
387       DO jk= 1, jpkm1
388
389          DO jj=1, jpjm1
390             DO ji = 1, fs_jpim1   ! vertor opt.
391                zwu(ji,jj) = -e2u(ji,jj) * u_trc_bbl(ji,jj,jk) * fse3u(ji,jj,jk)
392                zwv(ji,jj) = -e1v(ji,jj) * v_trc_bbl(ji,jj,jk) * fse3v(ji,jj,jk)
393             END DO
394          END DO
395
396          ! ... horizontal divergence
397          DO jj = 2, jpjm1
398             DO ji = fs_2, fs_jpim1   ! vector opt.
399                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
400                zhdivn(ji,jj,jk) = (  zwu(ji,jj) - zwu(ji-1,jj  )   &
401                   &                + zwv(ji,jj) - zwv(ji  ,jj-1)  ) / zbt
402             END DO
403          END DO
404       END DO
405
406
407       ! ... horizontal bottom divergence
408!gh
409       IF( ln_zps ) THEN
410     
411# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
412          jj = 1
413          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
414# else
415          DO jj = 1, jpjm1
416             DO ji = 1, jpim1
417# endif
418                iku  = mbku(ji  ,jj  )
419                ikv  = mbkv(ji  ,jj  ) 
420                iku1 = mbkt(ji+1,jj  )
421                iku2 = mbkt(ji  ,jj  )
422                ikv1 = mbkt(ji  ,jj+1)
423                ikv2 = mbkt(ji  ,jj  )
424                ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
425                ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
426               
427                zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * ze3u 
428                zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * ze3v
429#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
430             END DO
431#endif
432          END DO
433   
434       ELSE
435
436# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
437          jj = 1
438          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
439# else
440          DO jj = 1, jpjm1
441             DO ji = 1, jpim1
442# endif
443                iku = mbku(ji,jj)
444                ikv = mbkv(ji,jj)
445                zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
446                zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
447#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
448             END DO
449#endif
450          END DO
451       ENDIF
452
453# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
454       jj = 1
455       DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
456# else
457       DO jj = 2, jpjm1
458          DO ji = 2, jpim1
459# endif
460             ik = mbkt(ji,jj)
461             zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
462             zhdivn(ji,jj,ik) =   &
463                &   (  zwu(ji  ,jj  ) * ( zunb(ji  ,jj  ) - un(ji  ,jj  ,ik) )   &
464                &    - zwu(ji-1,jj  ) * ( zunb(ji-1,jj  ) - un(ji-1,jj  ,ik) )   &
465                &    + zwv(ji  ,jj  ) * ( zvnb(ji  ,jj  ) - vn(ji  ,jj  ,ik) )   &
466                &    - zwv(ji  ,jj-1) * ( zvnb(ji  ,jj-1) - vn(ji  ,jj-1,ik) )   &
467                &   ) / zbt
468
469# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
470          END DO
471# endif
472       END DO
473
474       ! 7. compute additional vertical velocity to be used in t boxes
475       ! -------------------------------------------------------------
476       
477       ! ... Computation from the bottom
478       ! Note that w_trc_bbl(:,:,jpk) has been set to 0 in trc_bbl_init
479       DO jk = jpkm1, 1, -1
480          DO jj= 2, jpjm1
481             DO ji = fs_2, fs_jpim1   ! vector opt.
482                w_trc_bbl(ji,jj,jk) = w_trc_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
483             END DO
484          END DO
485       END DO
486
487      ! Boundary condition on w_bbl   (unchanged sign)
488      CALL lbc_lnk( w_trc_bbl, 'W', 1. )
489
490   END SUBROUTINE trc_bbl_adv
Note: See TracBrowser for help on using the repository browser.