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.F90 in tags/nemo_dev_x8/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_dev_x8/NEMO/OPA_SRC/TRA/trabbl.F90 @ 280

Last change on this file since 280 was 106, checked in by opalod, 20 years ago

CT : UPDATE067 : Add control indices nictl, njctl used in SUM function output to compare mono versus multi procs runs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.1 KB
Line 
1MODULE trabbl
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
6#if   defined key_trabbl_dif   ||   defined key_trabbl_adv   || defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_trabbl_dif'   or            diffusive bottom boundary layer
9   !!   'key_trabbl_adv'                 advective bottom boundary layer
10   !!----------------------------------------------------------------------
11   !!   tra_bbl_dif  : update the active tracer trends due to the bottom
12   !!                  boundary layer (diffusive only)
13   !!   tra_bbl_adv  : update the active tracer trends due to the bottom
14   !!                  boundary layer (advective and/or diffusive)
15   !!   tra_bbl_init : initialization, namlist read, parameters control
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and active tracers
19   USE dom_oce         ! ocean space and time domain
20   USE trdtra_oce     ! ocean active tracer trends
21   USE in_out_manager  ! I/O manager
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC tra_bbl_dif    ! routine called by step.F90
28   PUBLIC tra_bbl_adv    ! routine called by step.F90
29
30   !! * Shared module variables
31   LOGICAL, PUBLIC, PARAMETER ::   &  !:
32      lk_trabbl_dif = .TRUE.   !: diffusive bottom boundary layer flag
33# if defined key_trabbl_adv
34   LOGICAL, PUBLIC, PARAMETER ::    &  !:
35      lk_trabbl_adv = .TRUE.   !: advective bottom boundary layer flag
36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
37       u_bbl, v_bbl,  &  !: velocity involved in exhanges in the advective BBL
38       w_bbl             !: vertical increment of velocity due to advective BBL
39       !                 !  only affect tracer vertical advection
40# else
41   LOGICAL, PUBLIC, PARAMETER ::    &  !:
42      lk_trabbl_adv = .FALSE.  !: advective bottom boundary layer flag
43# endif
44
45   !! * Module variables
46   INTEGER, DIMENSION(jpi,jpj) ::   &  !:
47      mbkt, mbku, mbkv                 ! ???
48   REAL(wp) ::        &  !!! * bbl namelist *
49      atrbbl = 1.e+3      ! lateral coeff. for bottom boundary layer scheme (m2/s)
50
51   !! * Substitutions
52#  include "domzgr_substitute.h90"
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !!   OPA 9.0 , LODYC-IPSL (2003)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE tra_bbl_dif( kt )
61      !!----------------------------------------------------------------------
62      !!                   ***  ROUTINE tra_bbl_dif  ***
63      !!
64      !! ** Purpose :   Compute the before tracer (t & s) trend associated
65      !!      with the bottom boundary layer and add it to the general trend
66      !!      of tracer equations. The bottom boundary layer is supposed to be
67      !!      a purely diffusive bottom boundary layer.
68      !!
69      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad
70      !!      is an along bottom slope gradient) an additional lateral diffu-
71      !!      sive trend along the bottom slope is added to the general tracer
72      !!      trend, otherwise nothing is done.
73      !!      Second order operator (laplacian type) with variable coefficient
74      !!      computed as follow for temperature (idem on s):
75      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
76      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
77      !!      where ztb is a 2D array: the bottom ocean temperature and ahtb
78      !!      is a time and space varying diffusive coefficient defined by:
79      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
80      !!              = 0.       otherwise.
81      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
82      !!      is evaluated using the local density (i.e. referenced at the
83      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
84      !!      a downslope velocity of 20 cm/s if the condition for slope
85      !!      convection is satified)
86      !!      Add this before trend to the general trend (ta,sa) of the
87      !!      botton ocean tracer point:
88      !!         ta = ta + difft
89      !!
90      !! ** Action  : - update (ta,sa) at the bottom level with the bottom
91      !!                boundary layer trend
92      !!              - save the trends in bbltrd ('key_diatrends')
93      !!
94      !! References :
95      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
96      !!
97      !! History :
98      !!   8.0  !  96-06  (L. Mortier)  Original code
99      !!   8.0  !  97-11  (G. Madec)  Optimization
100      !!   8.5  !  02-08  (G. Madec)  free form + modules
101      !!----------------------------------------------------------------------
102      !! * Arguments
103      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
104
105      !! * Local declarations
106      INTEGER ::   ji, jj                   ! dummy loop indices
107      INTEGER ::   ik
108      INTEGER ::   ii0, ii1, ij0, ij1       ! temporary integers
109#  if defined key_partial_steps
110      INTEGER  ::   iku1, iku2, ikv1,ikv2   ! temporary intergers
111      REAL(wp) ::   ze3u, ze3v              ! temporary scalars
112#  else
113      INTEGER ::   iku, ikv
114#  endif
115      REAL(wp) ::   &
116         zsign, zt, zs, zh, zalbet,      &  ! temporary scalars
117         zgdrho, zbtr, zta, zsa
118      REAL(wp), DIMENSION(jpi,jpj) ::    &
119        zki, zkj, zkw, zkx, zky, zkz,    &  ! temporary workspace arrays
120        ztnb, zsnb, zdep,                &
121        ztbb, zsbb, zahu, zahv
122      REAL(wp) ::   &
123         fsalbt, pft, pfs, pfh              ! statement function
124      !!----------------------------------------------------------------------
125      ! ratio alpha/beta
126      ! ================
127      !  fsalbt: ratio of thermal over saline expension coefficients
128      !       pft :  potential temperature in degrees celcius
129      !       pfs :  salinity anomaly (s-35) in psu
130      !       pfh :  depth in meters
131
132      fsalbt( pft, pfs, pfh ) =                                              &
133         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
134                                   - 0.203814e-03 ) * pft                    &
135                                   + 0.170907e-01 ) * pft                    &
136                                   + 0.665157e-01                            &
137         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
138         +  ( ( - 0.302285e-13 * pfh                                         &
139                - 0.251520e-11 * pfs                                         &
140                + 0.512857e-12 * pft * pft          ) * pfh                  &
141                                     - 0.164759e-06   * pfs                  &
142             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
143                                     + 0.380374e-04 ) * pfh   
144      !!----------------------------------------------------------------------
145
146      IF( kt == nit000 )   CALL tra_bbl_init
147
148      ! 0. 2D fields of bottom temperature and salinity, and bottom slope
149      ! -----------------------------------------------------------------
150      ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
151
152#  if defined key_vectopt_loop   &&   ! defined key_autotasking
153      jj = 1
154      DO ji = 1, jpij   ! vector opt. (forced unrolling)
155#  else
156      DO jj = 1, jpj
157         DO ji = 1, jpi
158#  endif
159            ik = mbkt(ji,jj)                              ! index of the bottom ocean T-level
160            ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1)   ! masked now T and S at ocean bottom
161            zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1)
162            ztbb(ji,jj) = tb(ji,jj,ik) * tmask(ji,jj,1)   ! masked before T and S at ocean bottom
163            zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1)
164            zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
165#  if ! defined key_vectopt_loop   ||   defined key_autotasking
166         END DO
167#  endif
168      END DO
169
170#  if defined key_partial_steps
171      ! partial steps correction
172#   if defined key_vectopt_loop   &&   ! defined key_autotasking
173      jj = 1
174      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
175#   else
176      DO jj = 1, jpjm1
177         DO ji = 1, jpim1
178#   endif
179            iku1 = MAX( mbathy(ji+1,jj  )-1, 1 )
180            iku2 = MAX( mbathy(ji  ,jj  )-1, 1 )
181            ikv1 = MAX( mbathy(ji  ,jj+1)-1, 1 )
182            ikv2 = MAX( mbathy(ji  ,jj  )-1, 1 )
183            ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) 
184            ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) 
185            zahu(ji,jj) = atrbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1)
186            zahv(ji,jj) = atrbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1)
187#   if ! defined key_vectopt_loop   ||   defined key_autotasking
188         END DO
189#   endif
190      END DO
191#  else
192#   if defined key_vectopt_loop   &&   ! defined key_autotasking
193      jj = 1
194      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
195#   else
196      DO jj = 1, jpjm1
197         DO ji = 1, jpim1
198#   endif
199            iku = mbku(ji,jj)
200            ikv = mbkv(ji,jj)
201            zahu(ji,jj) = atrbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1)
202            zahv(ji,jj) = atrbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1)
203#   if ! defined key_vectopt_loop   ||   defined key_autotasking
204         END DO
205#   endif
206      END DO
207#  endif
208
209      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
210      ! --------------------------------------------
211      ! Sign of the local density gradient along the i- and j-slopes
212      ! multiplied by the slope of the ocean bottom
213
214#  if defined key_vectopt_loop   &&   ! defined key_autotasking
215      jj = 1
216      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
217#  else
218      DO jj = 1, jpjm1
219         DO ji = 1, jpim1
220#  endif
221            ! temperature, salinity anomalie and depth
222            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
223            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
224            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
225            ! masked ratio alpha/beta
226            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
227            ! local density gradient along i-bathymetric slope
228            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
229                   -          ( zsnb(ji+1,jj) - zsnb(ji,jj) )
230            ! sign of local i-gradient of density multiplied by the i-slope
231            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
232            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
233#  if ! defined key_vectopt_loop   ||   defined key_autotasking
234         END DO
235#  endif
236      END DO
237
238#  if defined key_vectopt_loop   &&   ! defined key_autotasking
239      jj = 1
240      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
241#  else
242      DO jj = 1, jpjm1
243         DO ji = 1, jpim1
244#  endif
245            ! temperature, salinity anomalie and depth
246            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
247            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
248            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
249            ! masked ratio alpha/beta
250            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
251            ! local density gradient along j-bathymetric slope
252            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
253                   -          ( zsnb(ji,jj+1) - zsnb(ji,jj) )
254            ! sign of local j-gradient of density multiplied by the j-slope
255            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
256            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
257#  if ! defined key_vectopt_loop   ||   defined key_autotasking
258         END DO
259#  endif
260      END DO
261
262
263      ! 2. Additional second order diffusive trends
264      ! -------------------------------------------
265
266      ! first derivative (gradient)
267#  if defined key_vectopt_loop   &&   ! defined key_autotasking
268      jj = 1
269      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
270#  else
271      DO jj = 1, jpjm1
272         DO ji = 1, jpim1
273#  endif
274            zkx(ji,jj) = zki(ji,jj) * ( ztbb(ji+1,jj) - ztbb(ji,jj) )
275            zkz(ji,jj) = zki(ji,jj) * ( zsbb(ji+1,jj) - zsbb(ji,jj) )
276
277            zky(ji,jj) = zkj(ji,jj) * ( ztbb(ji,jj+1) - ztbb(ji,jj) )
278            zkw(ji,jj) = zkj(ji,jj) * ( zsbb(ji,jj+1) - zsbb(ji,jj) )
279#  if ! defined key_vectopt_loop   ||   defined key_autotasking
280         END DO
281#  endif
282      END DO
283
284      IF( cp_cfg == "orca" ) THEN
285
286         SELECT CASE ( jp_cfg )
287         !                                           ! =======================
288         CASE ( 2 )                                  !  ORCA_R2 configuration
289            !                                        ! =======================
290            ! Gibraltar enhancement of BBL
291            ij0 = 102   ;   ij1 = 102
292            ii0 = 139   ;   ii1 = 140 
293            zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
294            zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
295
296            ! Red Sea enhancement of BBL
297            ij0 =  88   ;   ij1 =  88
298            ii0 = 161   ;   ii1 = 162
299            zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
300            zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
301
302            !                                        ! =======================
303         CASE ( 4 )                                  !  ORCA_R4 configuration
304            !                                        ! =======================
305            ! Gibraltar enhancement of BBL
306            ij0 =  52   ;   ij1 =  52
307            ii0 =  70   ;   ii1 =  71 
308            zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
309            zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
310
311         END SELECT
312
313      ENDIF
314
315
316      ! second derivative (divergence) and add to the general tracer trend
317#  if defined key_vectopt_loop   &&   ! defined key_autotasking
318      jj = 1
319      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
320#  else
321      DO jj = 2, jpjm1
322         DO ji = 2, jpim1
323#  endif
324            ik = max( mbathy(ji,jj)-1, 1 )
325            zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
326            zta = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
327                   + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
328            zsa = (  zkz(ji,jj) - zkz(ji-1,jj  )    &
329                   + zkw(ji,jj) - zkw(ji  ,jj-1)  ) * zbtr
330            ta(ji,jj,ik) = ta(ji,jj,ik) + zta
331            sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
332#  if defined key_trdtra || defined key_trdmld
333            bbltrd(ji,jj,1) = + zta
334            bbltrd(ji,jj,2) = + zsa
335#  endif
336#  if ! defined key_vectopt_loop   ||   defined key_autotasking
337         END DO
338#  endif
339      END DO
340
341      IF(l_ctl) THEN
342         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
343         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
344         WRITE(numout,*) ' bbl  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
345         t_ctl = zta   ;   s_ctl = zsa
346      ENDIF
347
348   END SUBROUTINE tra_bbl_dif
349
350# if defined key_trabbl_adv
351   !!----------------------------------------------------------------------
352   !!   'key_trabbl_adv'                    advective bottom boundary layer
353   !!----------------------------------------------------------------------
354#  include "trabbl_adv.h90"
355# else
356   !!----------------------------------------------------------------------
357   !!   Default option :                 NO advective bottom boundary layer
358   !!----------------------------------------------------------------------
359   SUBROUTINE tra_bbl_adv (kt )              ! Empty routine
360      INTEGER, INTENT(in) :: kt
361      WRITE(*,*) 'tra_bbl_adv: You should not have seen this print! error?', kt
362   END SUBROUTINE tra_bbl_adv
363# endif
364
365   SUBROUTINE tra_bbl_init
366      !!----------------------------------------------------------------------
367      !!                  ***  ROUTINE tra_bbl_init  ***
368      !!
369      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
370      !!
371      !! ** Method  :   Read the nambbl namelist and check the parameters
372      !!      called by tra_bbl at the first timestep (nit000)
373      !!
374      !! History :
375      !!    8.5  !  02-08  (G. Madec)  Original code
376      !!----------------------------------------------------------------------
377      !! * Local declarations
378      INTEGER ::   ji, jj      ! dummy loop indices
379      NAMELIST/nambbl/ atrbbl
380      !!----------------------------------------------------------------------
381
382      ! Read Namelist nambbl : bottom boundary layer scheme
383      ! --------------------
384      REWIND ( numnam )
385      READ   ( numnam, nambbl )
386
387
388      ! Parameter control and print
389      ! ---------------------------
390      IF(lwp) THEN
391         WRITE(numout,*)
392         WRITE(numout,*) 'tra_bbl_init : * Diffusive Bottom Boundary Layer'
393         WRITE(numout,*) '~~~~~~~~~~~~'
394         IF( lk_trabbl_adv ) THEN
395            WRITE(numout,*) '               * Advective Bottom Boundary Layer'
396         ENDIF
397         WRITE(numout,*) '          Namelist nambbl : set bbl parameters'
398         WRITE(numout,*)
399         WRITE(numout,*) '          bottom boundary layer coef.    atrbbl = ', atrbbl
400         WRITE(numout,*)
401      ENDIF
402 
403      DO jj = 1, jpj
404         DO ji = 1, jpi
405            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
406         END DO
407      END DO
408      DO jj = 1, jpjm1
409         DO ji = 1, jpim1
410            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
411            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
412         END DO
413      END DO
414!!bug ???
415!!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0)
416
417# if defined key_trabbl_adv
418      w_bbl(:,:,:) = 0.e0    ! initialisation of w_bbl to zero
419# endif
420
421   END SUBROUTINE tra_bbl_init
422
423#else
424   !!----------------------------------------------------------------------
425   !!   Dummy module :                      No bottom boundary layer scheme
426   !!----------------------------------------------------------------------
427   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl_dif = .FALSE.   !: diff bbl flag
428   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl_adv = .FALSE.   !: adv  bbl flag
429CONTAINS
430   SUBROUTINE tra_bbl_dif (kt )              ! Empty routine
431      INTEGER, INTENT(in) :: kt
432      WRITE(*,*) 'tra_bbl_dif: You should not have seen this print! error?', kt
433   END SUBROUTINE tra_bbl_dif
434   SUBROUTINE tra_bbl_adv (kt )              ! Empty routine
435      INTEGER, INTENT(in) :: kt
436      WRITE(*,*) 'tra_bbl_adv: You should not have seen this print! error?', kt
437   END SUBROUTINE tra_bbl_adv
438#endif
439
440   !!======================================================================
441END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.