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.
trabbl_adv.h90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 789

Last change on this file since 789 was 789, checked in by rblod, 16 years ago

Suppress jki routines and associated key_mpp_omp

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