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.
trcbbl.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcbbl.F90 @ 349

Last change on this file since 349 was 349, checked in by opalod, 18 years ago

nemo_v1_update_031 : CT : change header names
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 KB
Line 
1MODULE trcbbl
2   !!==============================================================================
3   !!                       ***  MODULE  trcbbl  ***
4   !! Ocean passive tracers physics :  advective and/or diffusive bottom boundary
5   !!                                  layer scheme
6   !!==============================================================================
7#if  defined key_passivetrc && ( defined key_trcbbl_dif   ||   defined key_trcbbl_adv )
8   !!----------------------------------------------------------------------
9   !!   'key_trcbbl_dif'   or            diffusive bottom boundary layer
10   !!   'key_trcbbl_adv'                 advective bottom boundary layer
11   !!----------------------------------------------------------------------
12   !!   trc_bbl_dif  : update the passive tracer trends due to the bottom
13   !!                  boundary layer (diffusive only)
14   !!   trc_bbl_adv  : update the passive tracer trends due to the bottom
15   !!                  boundary layer (advective and/or diffusive)
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce_trc             ! ocean dynamics and active tracers variables
19   USE trc                 ! ocean passive tracers variables
20   USE prtctl_trc          ! Print control for debbuging
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC trc_bbl_dif    ! routine called by step.F90
27   PUBLIC trc_bbl_adv    ! routine called by step.F90
28
29   !! * Shared module variables
30   LOGICAL, PUBLIC, PARAMETER ::   &  !:
31      lk_trcbbl_dif = .TRUE.   !: diffusive bottom boundary layer flag
32# if defined key_trcbbl_adv
33   LOGICAL, PUBLIC, PARAMETER ::    &  !:
34      lk_trcbbl_adv = .TRUE.   !: advective bottom boundary layer flag
35   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
36       u_trc_bbl, v_trc_bbl, &  !: velocity involved in exhanges in the advective BBL
37       w_trc_bbl                !: vertical increment of velocity due to advective BBL
38       !                        !  only affect tracer vertical advection
39# else
40   LOGICAL, PUBLIC, PARAMETER ::    &  !:
41      lk_trcbbl_adv = .FALSE.  !: advective bottom boundary layer flag
42# endif
43
44   !! * Module variables
45   INTEGER, DIMENSION(jpi,jpj) ::   &  !:
46      mbkt, mbku, mbkv                 ! ???
47
48   !! * Substitutions
49#  include "passivetrc_substitute.h90"
50   !!----------------------------------------------------------------------
51   !!   TOP 1.0 , LOCEAN-IPSL (2005)
52   !! $Header$
53   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE trc_bbl_dif( kt )
59      !!----------------------------------------------------------------------
60      !!                   ***  ROUTINE trc_bbl_dif  ***
61      !!
62      !! ** Purpose :   Compute the before tracer trend associated
63      !!      with the bottom boundary layer and add it to the general trend
64      !!      of tracer equations. The bottom boundary layer is supposed to be
65      !!      a purely diffusive bottom boundary layer.
66      !!
67      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad
68      !!      is an along bottom slope gradient) an additional lateral diffu-
69      !!      sive trend along the bottom slope is added to the general tracer
70      !!      trend, otherwise nothing is done.
71      !!      Second order operator (laplacian type) with variable coefficient
72      !!      computed as follow for temperature (idem on s):
73      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
74      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
75      !!      where ztb is a 2D array: the bottom ocean temperature and ahtb
76      !!      is a time and space varying diffusive coefficient defined by:
77      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
78      !!              = 0.       otherwise.
79      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
80      !!      is evaluated using the local density (i.e. referenced at the
81      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
82      !!      a downslope velocity of 20 cm/s if the condition for slope
83      !!      convection is satified)
84      !!      Add this before trend to the general trend tra of the
85      !!      botton ocean tracer point:
86      !!         tra = tra + difft
87      !!
88      !! ** Action  : - update tra at the bottom level with the bottom
89      !!                boundary layer trend
90      !!
91      !! References :
92      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
93      !!
94      !! History :
95      !!   8.0  !  96-06  (L. Mortier)  Original code
96      !!   8.0  !  97-11  (G. Madec)  Optimization
97      !!   8.5  !  02-08  (G. Madec)  free form + modules
98      !!   9.0  !  04-03  (C. Ethe)   Adaptation for passive tracers
99      !!----------------------------------------------------------------------
100      !! * Arguments
101      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
102
103      !! * Local declarations
104      INTEGER ::   ji, jj,jn                ! dummy loop indices
105      INTEGER ::   ik
106      INTEGER ::   ii0, ii1, ij0, ij1       ! temporary integers
107#  if defined key_partial_steps
108      INTEGER  ::   iku1, iku2, ikv1,ikv2   ! temporary intergers
109      REAL(wp) ::   ze3u, ze3v              ! temporary scalars
110#  else
111      INTEGER ::   iku, ikv
112#  endif
113      REAL(wp) ::   &
114         zsign, zt, zs, zh, zalbet,      &  ! temporary scalars
115         zgdrho, zbtr, ztra
116      REAL(wp), DIMENSION(jpi,jpj) ::    &
117        zki, zkj, zkx, zky,    &  ! temporary workspace arrays
118        ztnb, zsnb, zdep,                &
119        ztrb, zahu, zahv
120      CHARACTER (len=22) :: charout
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
146      IF( kt == nittrc000 )   CALL trc_bbl_init
147
148
149      ! 0. 2D fields of bottom temperature and salinity, and bottom slope
150      ! -----------------------------------------------------------------
151      ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
152
153#  if defined key_vectopt_loop   &&   ! defined key_autotasking
154      jj = 1
155      DO ji = 1, jpij   ! vector opt. (forced unrolling)
156#  else
157      DO jj = 1, jpj
158         DO ji = 1, jpi
159#  endif
160            ik = mbkt(ji,jj)                              ! index of the bottom ocean T-level
161            ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1)   ! masked now T and S at ocean bottom
162            zsnb(ji,jj) = sn(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!!
209!!     OFFLINE VERSION OF DIFFUSIVE BBL
210!!
211#if defined key_off_tra
212
213      ! 2. Additional second order diffusive trends
214      ! -------------------------------------------
215
216      DO jn = 1, jptra
217         ! first derivative (gradient)
218         
219#  if defined key_vectopt_loop   &&   ! defined key_autotasking
220         jj = 1
221         DO ji = 1, jpij   ! vector opt. (forced unrolling)
222#  else
223         DO jj = 1, jpj
224            DO ji = 1, jpi
225#  endif
226               ik = mbkt(ji,jj) 
227               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
228#  if ! defined key_vectopt_loop   ||   defined key_autotasking
229            END DO
230#  endif
231         END DO
232
233#  if defined key_vectopt_loop   &&   ! defined key_autotasking
234         jj = 1
235         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
236#  else
237         DO jj = 1, jpjm1
238            DO ji = 1, jpim1
239#  endif
240               zkx(ji,jj) = bblx(ji,jj) * zahu(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
241               zky(ji,jj) = bbly(ji,jj) * zahv(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
242#  if ! defined key_vectopt_loop   ||   defined key_autotasking
243            END DO
244#  endif
245         END DO
246!!
247!!  ONLINE VERSION OF DIFFUSIVE BBL
248!!
249#else
250      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
251      ! --------------------------------------------
252      ! Sign of the local density gradient along the i- and j-slopes
253      ! multiplied by the slope of the ocean bottom
254
255#  if defined key_vectopt_loop   &&   ! defined key_autotasking
256      jj = 1
257      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
258#  else
259      DO jj = 1, jpjm1
260         DO ji = 1, jpim1
261#  endif
262            ! temperature, salinity anomalie and depth
263            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
264            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
265            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
266            ! masked ratio alpha/beta
267            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
268            ! local density gradient along i-bathymetric slope
269            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
270                   -          ( zsnb(ji+1,jj) - zsnb(ji,jj) )
271            ! sign of local i-gradient of density multiplied by the i-slope
272            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
273            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
274#  if ! defined key_vectopt_loop   ||   defined key_autotasking
275         END DO
276#  endif
277      END DO
278
279#  if defined key_vectopt_loop   &&   ! defined key_autotasking
280      jj = 1
281      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
282#  else
283      DO jj = 1, jpjm1
284         DO ji = 1, jpim1
285#  endif
286            ! temperature, salinity anomalie and depth
287            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
288            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
289            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
290            ! masked ratio alpha/beta
291            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
292            ! local density gradient along j-bathymetric slope
293            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
294                   -          ( zsnb(ji,jj+1) - zsnb(ji,jj) )
295            ! sign of local j-gradient of density multiplied by the j-slope
296            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
297            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
298#  if ! defined key_vectopt_loop   ||   defined key_autotasking
299         END DO
300#  endif
301      END DO
302
303      ! 2. Additional second order diffusive trends
304      ! -------------------------------------------
305
306      DO jn = 1, jptra
307         ! first derivative (gradient)
308
309#  if defined key_vectopt_loop   &&   ! defined key_autotasking
310         jj = 1
311         DO ji = 1, jpij   ! vector opt. (forced unrolling)
312#  else
313         DO jj = 1, jpj
314            DO ji = 1, jpi
315#  endif
316               ik = mbkt(ji,jj)
317               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
318#  if ! defined key_vectopt_loop   ||   defined key_autotasking
319            END DO
320#  endif
321         END DO
322#  if defined key_vectopt_loop   &&   ! defined key_autotasking
323         jj = 1
324         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
325#  else
326         DO jj = 1, jpjm1
327            DO ji = 1, jpim1
328#  endif
329               zkx(ji,jj) = zki(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
330               zky(ji,jj) = zkj(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
331#  if ! defined key_vectopt_loop   ||   defined key_autotasking
332            END DO
333#  endif
334         END DO
335
336#endif
337
338         IF( cp_cfg == "orca" ) THEN
339           
340            SELECT CASE ( jp_cfg )
341               !                                           ! =======================
342            CASE ( 2 )                                  !  ORCA_R2 configuration
343               !                                        ! =======================
344               ! Gibraltar enhancement of BBL
345               ij0 = 102   ;   ij1 = 102
346               ii0 = 139   ;   ii1 = 140 
347               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
348               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
349               
350               ! Red Sea enhancement of BBL
351               ij0 =  88   ;   ij1 =  88
352               ii0 = 161   ;   ii1 = 162
353               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
354               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
355               
356               !                                        ! =======================
357            CASE ( 4 )                                  !  ORCA_R4 configuration
358               !                                        ! =======================
359               ! Gibraltar enhancement of BBL
360               ij0 =  52   ;   ij1 =  52
361               ii0 =  70   ;   ii1 =  71 
362               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
363               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
364               
365            END SELECT
366           
367         ENDIF
368
369         
370         ! second derivative (divergence) and add to the general tracer trend
371#  if defined key_vectopt_loop   &&   ! defined key_autotasking
372         jj = 1
373         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
374#  else
375         DO jj = 2, jpjm1
376            DO ji = 2, jpim1
377#  endif
378               ik = MAX( mbathy(ji,jj)-1, 1 )
379               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
380               ztra = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
381                  &    + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
382               tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
383#  if ! defined key_vectopt_loop   ||   defined key_autotasking
384            END DO
385#  endif
386         END DO
387
388      END DO
389
390      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
391         WRITE(charout, FMT="('bbl - dif')")
392         CALL prt_ctl_trc_info(charout)
393         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
394      ENDIF
395
396   END SUBROUTINE trc_bbl_dif
397
398# if defined key_trcbbl_adv
399   !!----------------------------------------------------------------------
400   !!   'key_trcbbl_adv'                    advective bottom boundary layer
401   !!----------------------------------------------------------------------
402#  include "trcbbl_adv.h90"
403# else
404   !!----------------------------------------------------------------------
405   !!   Default option :                 NO advective bottom boundary layer
406   !!----------------------------------------------------------------------
407   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
408      INTEGER, INTENT(in) :: kt
409      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
410   END SUBROUTINE trc_bbl_adv
411# endif
412
413   SUBROUTINE trc_bbl_init
414      !!----------------------------------------------------------------------
415      !!                  ***  ROUTINE trc_bbl_init  ***
416      !!
417      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
418      !!
419      !! ** Method  :   Read the namtrcbbl namelist and check the parameters
420      !!      called by tra_bbl at the first timestep (nittrc000)
421      !!
422      !! History :
423      !!    8.5  !  02-08  (G. Madec)  Original code
424      !!----------------------------------------------------------------------
425      !! * Local declarations
426      INTEGER ::   ji, jj      ! dummy loop indices
427
428      !!----------------------------------------------------------------------
429
430
431      ! Parameter control and print
432      ! ---------------------------
433      IF(lwp) THEN
434         WRITE(numout,*)
435         WRITE(numout,*) 'trc_bbl_init : * Diffusive Bottom Boundary Layer'
436         WRITE(numout,*) '~~~~~~~~~~~~'
437# if defined key_trcbbl_adv
438            WRITE(numout,*) '               * Advective Bottom Boundary Layer'
439# endif
440         WRITE(numout,*)
441      ENDIF
442 
443      DO jj = 1, jpj
444         DO ji = 1, jpi
445            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
446         END DO
447      END DO
448      DO jj = 1, jpjm1
449         DO ji = 1, jpim1
450            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
451            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
452         END DO
453      END DO
454!!bug ???
455!!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0)
456
457# if defined key_trcbbl_adv
458      w_trc_bbl(:,:,:) = 0.e0    ! initialisation of w_trc_bbl to zero
459# endif
460
461   END SUBROUTINE trc_bbl_init
462
463#else
464   !!----------------------------------------------------------------------
465   !!   Dummy module :                      No bottom boundary layer scheme
466   !!----------------------------------------------------------------------
467   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .FALSE.   !: diff bbl flag
468   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .FALSE.   !: adv  bbl flag
469CONTAINS
470   SUBROUTINE trc_bbl_dif (kt )              ! Empty routine
471      INTEGER, INTENT(in) :: kt
472      WRITE(*,*) 'trc_bbl_dif: You should not have seen this print! error?', kt
473   END SUBROUTINE trc_bbl_dif
474   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
475      INTEGER, INTENT(in) :: kt
476      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
477   END SUBROUTINE trc_bbl_adv
478#endif
479
480   !!======================================================================
481END MODULE trcbbl
Note: See TracBrowser for help on using the repository browser.