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 trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/trabbl.F90 @ 3

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

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.5 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.   ! bottom boundary layer flag
36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
37       w_bbl,      &  ! vertical increment of velocity due to advective BBL
38       !              ! only affect tracer vertical advection
39       u_bbl, v_bbl   ! velocity involved in exhanges in the advective BBL
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#  if defined key_partial_steps
109      INTEGER  ::   iku1, iku2, ikv1,ikv2   ! temporary intergers
110      REAL(wp) ::   ze3u, ze3v              ! temporary scalars
111#  else
112      INTEGER ::   iku, ikv
113#  endif
114      REAL(wp) ::   &
115         zsign, zt, zs, zh, zalbet,      &  ! temporary scalars
116         zgdrho, zbtr, zta, zsa
117      REAL(wp), DIMENSION(jpi,jpj) ::    &
118        zki, zkj, zkw, zkx, zky, zkz,    &  ! temporary workspace arrays
119        ztnb, zsnb, zdep,                &
120        ztbb, zsbb, zahu, zahv
121      REAL(wp) ::   &
122         fsalbt, pft, pfs, pfh              ! statement function
123      !!----------------------------------------------------------------------
124      ! ratio alpha/beta
125      ! ================
126      !  fsalbt: ratio of thermal over saline expension coefficients
127      !       pft :  potential temperature in degrees celcius
128      !       pfs :  salinity anomaly (s-35) in psu
129      !       pfh :  depth in meters
130
131      fsalbt( pft, pfs, pfh ) =                                              &
132         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
133                                   - 0.203814e-03 ) * pft                    &
134                                   + 0.170907e-01 ) * pft                    &
135                                   + 0.665157e-01                            &
136         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
137         +  ( ( - 0.302285e-13 * pfh                                         &
138                - 0.251520e-11 * pfs                                         &
139                + 0.512857e-12 * pft * pft          ) * pfh                  &
140                                     - 0.164759e-06   * pfs                  &
141             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
142                                     + 0.380374e-04 ) * pfh   
143      !!----------------------------------------------------------------------
144
145      IF( kt == nit000 )   CALL tra_bbl_init
146
147      ! 0. 2D fields of bottom temperature and salinity, and bottom slope
148      ! -----------------------------------------------------------------
149      ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
150
151#  if defined key_vectopt_loop   &&   ! defined key_autotasking
152      jj = 1
153      DO ji = 1, jpij   ! vector opt. (forced unrolling)
154#  else
155      DO jj = 1, jpj
156         DO ji = 1, jpi
157#  endif
158            ik = mbkt(ji,jj)                              ! index of the bottom ocean T-level
159            ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1)   ! masked now T and S at ocean bottom
160            zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1)
161            ztbb(ji,jj) = tb(ji,jj,ik) * tmask(ji,jj,1)   ! masked before T and S at ocean bottom
162            zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1)
163            zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
164#  if ! defined key_vectopt_loop   ||   defined key_autotasking
165         END DO
166#  endif
167      END DO
168
169#  if defined key_partial_steps
170      ! partial steps correction
171#   if defined key_vectopt_loop   &&   ! defined key_autotasking
172      jj = 1
173      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
174#   else
175      DO jj = 1, jpjm1
176         DO ji = 1, jpim1
177#   endif
178            iku1 = MAX( mbathy(ji+1,jj  )-1, 1 )
179            iku2 = MAX( mbathy(ji  ,jj  )-1, 1 )
180            ikv1 = MAX( mbathy(ji  ,jj+1)-1, 1 )
181            ikv2 = MAX( mbathy(ji  ,jj  )-1, 1 )
182            ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) 
183            ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) 
184            zahu(ji,jj) = atrbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1)
185            zahv(ji,jj) = atrbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1)
186#   if ! defined key_vectopt_loop   ||   defined key_autotasking
187         END DO
188#   endif
189      END DO
190#  else
191#   if defined key_vectopt_loop   &&   ! defined key_autotasking
192      jj = 1
193      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
194#   else
195      DO jj = 1, jpjm1
196         DO ji = 1, jpim1
197#   endif
198            iku = mbku(ji,jj)
199            ikv = mbkv(ji,jj)
200            zahu(ji,jj) = atrbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1)
201            zahv(ji,jj) = atrbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1)
202#   if ! defined key_vectopt_loop   ||   defined key_autotasking
203         END DO
204#   endif
205      END DO
206#  endif
207
208      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
209      ! --------------------------------------------
210      ! Sign of the local density gradient along the i- and j-slopes
211      ! multiplied by the slope of the ocean bottom
212
213#  if defined key_vectopt_loop   &&   ! defined key_autotasking
214      jj = 1
215      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
216#  else
217      DO jj = 1, jpjm1
218         DO ji = 1, jpim1
219#  endif
220            ! temperature, salinity anomalie and depth
221            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
222            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
223            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
224            ! masked ratio alpha/beta
225            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
226            ! local density gradient along i-bathymetric slope
227            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
228                   -          ( zsnb(ji+1,jj) - zsnb(ji,jj) )
229            ! sign of local i-gradient of density multiplied by the i-slope
230            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
231            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
232#  if ! defined key_vectopt_loop   ||   defined key_autotasking
233         END DO
234#  endif
235      END DO
236
237#  if defined key_vectopt_loop   &&   ! defined key_autotasking
238      jj = 1
239      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
240#  else
241      DO jj = 1, jpjm1
242         DO ji = 1, jpim1
243#  endif
244            ! temperature, salinity anomalie and depth
245            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
246            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
247            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
248            ! masked ratio alpha/beta
249            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
250            ! local density gradient along j-bathymetric slope
251            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
252                   -          ( zsnb(ji,jj+1) - zsnb(ji,jj) )
253            ! sign of local j-gradient of density multiplied by the j-slope
254            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
255            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
256#  if ! defined key_vectopt_loop   ||   defined key_autotasking
257         END DO
258#  endif
259      END DO
260
261
262      ! 2. Additional second order diffusive trends
263      ! -------------------------------------------
264
265      ! first derivative (gradient)
266#  if defined key_vectopt_loop   &&   ! defined key_autotasking
267      jj = 1
268      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
269#  else
270      DO jj = 1, jpjm1
271         DO ji = 1, jpim1
272#  endif
273            zkx(ji,jj) = zki(ji,jj) * ( ztbb(ji+1,jj) - ztbb(ji,jj) )
274            zkz(ji,jj) = zki(ji,jj) * ( zsbb(ji+1,jj) - zsbb(ji,jj) )
275
276            zky(ji,jj) = zkj(ji,jj) * ( ztbb(ji,jj+1) - ztbb(ji,jj) )
277            zkw(ji,jj) = zkj(ji,jj) * ( zsbb(ji,jj+1) - zsbb(ji,jj) )
278#  if ! defined key_vectopt_loop   ||   defined key_autotasking
279         END DO
280#  endif
281      END DO
282
283      IF( cp_cfg == "orca" ) THEN
284
285         SELECT CASE ( jp_cfg )
286         !                                           ! =======================
287         CASE ( 2 )                                  !  ORCA_R2 configuration
288            !                                        ! =======================
289            ! Gibraltar enhancement of BBL
290            zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) )
291            zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) )
292
293            ! Red Sea enhancement of BBL
294            zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) )
295            zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) )
296
297            !                                        ! =======================
298         CASE ( 4 )                                  !  ORCA_R4 configuration
299            !                                        ! =======================
300            ! Gibraltar enhancement of BBL
301            zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) )
302            zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) )
303
304         END SELECT
305
306      ENDIF
307
308
309      ! second derivative (divergence) and add to the general tracer trend
310#  if defined key_vectopt_loop   &&   ! defined key_autotasking
311      jj = 1
312      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
313#  else
314      DO jj = 2, jpjm1
315         DO ji = 2, jpim1
316#  endif
317            ik = max( mbathy(ji,jj)-1, 1 )
318            zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
319            zta = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
320                   + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
321            zsa = (  zkz(ji,jj) - zkz(ji-1,jj  )    &
322                   + zkw(ji,jj) - zkw(ji  ,jj-1)  ) * zbtr
323            ta(ji,jj,ik) = ta(ji,jj,ik) + zta
324            sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
325#  if defined key_trdtra || defined key_trdmld
326            bbltrd(ji,jj,1) = + zta
327            bbltrd(ji,jj,2) = + zsa
328#  endif
329#  if ! defined key_vectopt_loop   ||   defined key_autotasking
330         END DO
331#  endif
332      END DO
333
334      IF( l_ctl .AND. lwp ) THEN
335         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
336         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
337         WRITE(numout,*) ' bbl  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
338         t_ctl = zta   ;   s_ctl = zsa
339      ENDIF
340
341   END SUBROUTINE tra_bbl_dif
342
343# if defined key_trabbl_adv
344   !!----------------------------------------------------------------------
345   !!   'key_trabbl_adv'                    advective bottom boundary layer
346   !!----------------------------------------------------------------------
347#  include "trabbl_adv.h90"
348# else
349   !!----------------------------------------------------------------------
350   !!   Default option :                 NO advective bottom boundary layer
351   !!----------------------------------------------------------------------
352   SUBROUTINE tra_bbl_adv (kt )              ! Empty routine
353      INTEGER, INTENT(in) :: kt
354      WRITE(*,*) kt
355   END SUBROUTINE tra_bbl_adv
356# endif
357
358   SUBROUTINE tra_bbl_init
359      !!----------------------------------------------------------------------
360      !!                  ***  ROUTINE tra_bbl_init  ***
361      !!
362      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
363      !!
364      !! ** Method  :   Read the nambbl namelist and check the parameters
365      !!      called by tra_bbl at the first timestep (nit000)
366      !!
367      !! History :
368      !!    8.5  !  02-08  (G. Madec)  Original code
369      !!----------------------------------------------------------------------
370      !! * Local declarations
371      INTEGER ::   ji, jj      ! dummy loop indices
372      NAMELIST/nambbl/ atrbbl
373      !!----------------------------------------------------------------------
374
375      ! Read Namelist nambbl : bottom boundary layer scheme
376      ! --------------------
377      REWIND ( numnam )
378      READ   ( numnam, nambbl )
379
380
381      ! Parameter control and print
382      ! ---------------------------
383      IF(lwp) THEN
384         WRITE(numout,*)
385         WRITE(numout,*) 'tra_bbl_init : * Diffusive Bottom Boundary Layer'
386         WRITE(numout,*) '~~~~~~~~~~~~'
387         IF( lk_trabbl_adv ) THEN
388            WRITE(numout,*) '               * Advective Bottom Boundary Layer'
389         ENDIF
390         WRITE(numout,*) '          Namelist nambbl : set bbl parameters'
391         WRITE(numout,*)
392         WRITE(numout,*) '          bottom boundary layer coef.    atrbbl = ', atrbbl
393         WRITE(numout,*)
394      ENDIF
395 
396      DO jj = 1, jpj
397         DO ji = 1, jpi
398            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
399         END DO
400      END DO
401      DO jj = 1, jpjm1
402         DO ji = 1, jpim1
403            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
404            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
405         END DO
406      END DO
407!!bug ???
408!!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0)
409
410# if defined key_trabbl_adv
411      w_bbl(:,:,:) = 0.e0    ! initialisation of w_bbl to zero
412# endif
413
414   END SUBROUTINE tra_bbl_init
415
416#else
417   !!----------------------------------------------------------------------
418   !!   Dummy module :                      No bottom boundary layer scheme
419   !!----------------------------------------------------------------------
420   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl_dif = .FALSE.
421   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl_adv = .FALSE.
422CONTAINS
423   SUBROUTINE tra_bbl_dif (kt )              ! Empty routine
424      INTEGER, INTENT(in) :: kt
425      WRITE(*,*) kt
426   END SUBROUTINE tra_bbl_dif
427   SUBROUTINE tra_bbl_adv (kt )              ! Empty routine
428      INTEGER, INTENT(in) :: kt
429      WRITE(*,*) kt
430   END SUBROUTINE tra_bbl_adv
431#endif
432
433   !!======================================================================
434END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.