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 @ 340

Last change on this file since 340 was 334, checked in by opalod, 19 years ago

nemo_v1_update_022 : CE + RB + CT : add print control possibility
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.8 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   !!   OPA 9.0 , LODYC-IPSL (2003)
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE trc_bbl_dif( kt )
57      !!----------------------------------------------------------------------
58      !!                   ***  ROUTINE trc_bbl_dif  ***
59      !!
60      !! ** Purpose :   Compute the before tracer trend associated
61      !!      with the bottom boundary layer and add it to the general trend
62      !!      of tracer equations. The bottom boundary layer is supposed to be
63      !!      a purely diffusive bottom boundary layer.
64      !!
65      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad
66      !!      is an along bottom slope gradient) an additional lateral diffu-
67      !!      sive trend along the bottom slope is added to the general tracer
68      !!      trend, otherwise nothing is done.
69      !!      Second order operator (laplacian type) with variable coefficient
70      !!      computed as follow for temperature (idem on s):
71      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
72      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
73      !!      where ztb is a 2D array: the bottom ocean temperature and ahtb
74      !!      is a time and space varying diffusive coefficient defined by:
75      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
76      !!              = 0.       otherwise.
77      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
78      !!      is evaluated using the local density (i.e. referenced at the
79      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
80      !!      a downslope velocity of 20 cm/s if the condition for slope
81      !!      convection is satified)
82      !!      Add this before trend to the general trend tra of the
83      !!      botton ocean tracer point:
84      !!         tra = tra + difft
85      !!
86      !! ** Action  : - update tra at the bottom level with the bottom
87      !!                boundary layer trend
88      !!
89      !! References :
90      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
91      !!
92      !! History :
93      !!   8.0  !  96-06  (L. Mortier)  Original code
94      !!   8.0  !  97-11  (G. Madec)  Optimization
95      !!   8.5  !  02-08  (G. Madec)  free form + modules
96      !!   9.0  !  04-03  (C. Ethe)   Adaptation for passive tracers
97      !!----------------------------------------------------------------------
98      !! * Arguments
99      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
100
101      !! * Local declarations
102      INTEGER ::   ji, jj,jn                ! dummy loop indices
103      INTEGER ::   ik
104      INTEGER ::   ii0, ii1, ij0, ij1       ! temporary integers
105#  if defined key_partial_steps
106      INTEGER  ::   iku1, iku2, ikv1,ikv2   ! temporary intergers
107      REAL(wp) ::   ze3u, ze3v              ! temporary scalars
108#  else
109      INTEGER ::   iku, ikv
110#  endif
111      REAL(wp) ::   &
112         zsign, zt, zs, zh, zalbet,      &  ! temporary scalars
113         zgdrho, zbtr, ztra
114      REAL(wp), DIMENSION(jpi,jpj) ::    &
115        zki, zkj, zkx, zky,    &  ! temporary workspace arrays
116        ztnb, zsnb, zdep,                &
117        ztrb, zahu, zahv
118      CHARACTER (len=22) :: charout
119      REAL(wp) ::   &
120         fsalbt, pft, pfs, pfh              ! statement function
121      !!----------------------------------------------------------------------
122      ! ratio alpha/beta
123      ! ================
124      !  fsalbt: ratio of thermal over saline expension coefficients
125      !       pft :  potential temperature in degrees celcius
126      !       pfs :  salinity anomaly (s-35) in psu
127      !       pfh :  depth in meters
128
129      fsalbt( pft, pfs, pfh ) =                                              &
130         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
131                                   - 0.203814e-03 ) * pft                    &
132                                   + 0.170907e-01 ) * pft                    &
133                                   + 0.665157e-01                            &
134         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
135         +  ( ( - 0.302285e-13 * pfh                                         &
136                - 0.251520e-11 * pfs                                         &
137                + 0.512857e-12 * pft * pft          ) * pfh                  &
138                                     - 0.164759e-06   * pfs                  &
139             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
140                                     + 0.380374e-04 ) * pfh   
141      !!----------------------------------------------------------------------
142
143
144      IF( kt == nittrc000 )   CALL trc_bbl_init
145
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            zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
162#  if ! defined key_vectopt_loop   ||   defined key_autotasking
163         END DO
164#  endif
165      END DO
166
167#  if defined key_partial_steps
168      ! partial steps correction
169#   if defined key_vectopt_loop   &&   ! defined key_autotasking
170      jj = 1
171      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
172#   else
173      DO jj = 1, jpjm1
174         DO ji = 1, jpim1
175#   endif
176            iku1 = MAX( mbathy(ji+1,jj  )-1, 1 )
177            iku2 = MAX( mbathy(ji  ,jj  )-1, 1 )
178            ikv1 = MAX( mbathy(ji  ,jj+1)-1, 1 )
179            ikv2 = MAX( mbathy(ji  ,jj  )-1, 1 )
180            ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) 
181            ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) 
182            zahu(ji,jj) = atrbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1)
183            zahv(ji,jj) = atrbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1)
184#   if ! defined key_vectopt_loop   ||   defined key_autotasking
185         END DO
186#   endif
187      END DO
188#  else
189#   if defined key_vectopt_loop   &&   ! defined key_autotasking
190      jj = 1
191      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
192#   else
193      DO jj = 1, jpjm1
194         DO ji = 1, jpim1
195#   endif
196            iku = mbku(ji,jj)
197            ikv = mbkv(ji,jj)
198            zahu(ji,jj) = atrbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1)
199            zahv(ji,jj) = atrbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1)
200#   if ! defined key_vectopt_loop   ||   defined key_autotasking
201         END DO
202#   endif
203      END DO
204#  endif
205
206!!
207!!     OFFLINE VERSION OF DIFFUSIVE BBL
208!!
209#if defined key_off_tra
210
211      ! 2. Additional second order diffusive trends
212      ! -------------------------------------------
213
214      DO jn = 1, jptra
215         ! first derivative (gradient)
216         
217#  if defined key_vectopt_loop   &&   ! defined key_autotasking
218         jj = 1
219         DO ji = 1, jpij   ! vector opt. (forced unrolling)
220#  else
221         DO jj = 1, jpj
222            DO ji = 1, jpi
223#  endif
224               ik = mbkt(ji,jj) 
225               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
226#  if ! defined key_vectopt_loop   ||   defined key_autotasking
227            END DO
228#  endif
229         END DO
230
231#  if defined key_vectopt_loop   &&   ! defined key_autotasking
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               zkx(ji,jj) = bblx(ji,jj) * zahu(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
239               zky(ji,jj) = bbly(ji,jj) * zahv(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
240#  if ! defined key_vectopt_loop   ||   defined key_autotasking
241            END DO
242#  endif
243         END DO
244!!
245!!  ONLINE VERSION OF DIFFUSIVE BBL
246!!
247#else
248      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
249      ! --------------------------------------------
250      ! Sign of the local density gradient along the i- and j-slopes
251      ! multiplied by the slope of the ocean bottom
252
253#  if defined key_vectopt_loop   &&   ! defined key_autotasking
254      jj = 1
255      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
256#  else
257      DO jj = 1, jpjm1
258         DO ji = 1, jpim1
259#  endif
260            ! temperature, salinity anomalie and depth
261            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
262            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
263            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
264            ! masked ratio alpha/beta
265            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
266            ! local density gradient along i-bathymetric slope
267            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
268                   -          ( zsnb(ji+1,jj) - zsnb(ji,jj) )
269            ! sign of local i-gradient of density multiplied by the i-slope
270            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
271            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
272#  if ! defined key_vectopt_loop   ||   defined key_autotasking
273         END DO
274#  endif
275      END DO
276
277#  if defined key_vectopt_loop   &&   ! defined key_autotasking
278      jj = 1
279      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
280#  else
281      DO jj = 1, jpjm1
282         DO ji = 1, jpim1
283#  endif
284            ! temperature, salinity anomalie and depth
285            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
286            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
287            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
288            ! masked ratio alpha/beta
289            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
290            ! local density gradient along j-bathymetric slope
291            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
292                   -          ( zsnb(ji,jj+1) - zsnb(ji,jj) )
293            ! sign of local j-gradient of density multiplied by the j-slope
294            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
295            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
296#  if ! defined key_vectopt_loop   ||   defined key_autotasking
297         END DO
298#  endif
299      END DO
300
301      ! 2. Additional second order diffusive trends
302      ! -------------------------------------------
303
304      DO jn = 1, jptra
305         ! first derivative (gradient)
306
307#  if defined key_vectopt_loop   &&   ! defined key_autotasking
308         jj = 1
309         DO ji = 1, jpij   ! vector opt. (forced unrolling)
310#  else
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313#  endif
314               ik = mbkt(ji,jj)
315               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
316#  if ! defined key_vectopt_loop   ||   defined key_autotasking
317            END DO
318#  endif
319         END DO
320#  if defined key_vectopt_loop   &&   ! defined key_autotasking
321         jj = 1
322         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
323#  else
324         DO jj = 1, jpjm1
325            DO ji = 1, jpim1
326#  endif
327               zkx(ji,jj) = zki(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
328               zky(ji,jj) = zkj(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
329#  if ! defined key_vectopt_loop   ||   defined key_autotasking
330            END DO
331#  endif
332         END DO
333
334#endif
335
336         IF( cp_cfg == "orca" ) THEN
337           
338            SELECT CASE ( jp_cfg )
339               !                                           ! =======================
340            CASE ( 2 )                                  !  ORCA_R2 configuration
341               !                                        ! =======================
342               ! Gibraltar enhancement of BBL
343               ij0 = 102   ;   ij1 = 102
344               ii0 = 139   ;   ii1 = 140 
345               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
346               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
347               
348               ! Red Sea enhancement of BBL
349               ij0 =  88   ;   ij1 =  88
350               ii0 = 161   ;   ii1 = 162
351               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
352               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
353               
354               !                                        ! =======================
355            CASE ( 4 )                                  !  ORCA_R4 configuration
356               !                                        ! =======================
357               ! Gibraltar enhancement of BBL
358               ij0 =  52   ;   ij1 =  52
359               ii0 =  70   ;   ii1 =  71 
360               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
361               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
362               
363            END SELECT
364           
365         ENDIF
366
367         
368         ! second derivative (divergence) and add to the general tracer trend
369#  if defined key_vectopt_loop   &&   ! defined key_autotasking
370         jj = 1
371         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
372#  else
373         DO jj = 2, jpjm1
374            DO ji = 2, jpim1
375#  endif
376               ik = MAX( mbathy(ji,jj)-1, 1 )
377               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
378               ztra = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
379                  &    + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
380               tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
381#  if ! defined key_vectopt_loop   ||   defined key_autotasking
382            END DO
383#  endif
384         END DO
385
386      END DO
387
388      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
389         WRITE(charout, FMT="('bbl - dif')")
390         CALL prt_ctl_trc_info(charout)
391         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
392      ENDIF
393
394   END SUBROUTINE trc_bbl_dif
395
396# if defined key_trcbbl_adv
397   !!----------------------------------------------------------------------
398   !!   'key_trcbbl_adv'                    advective bottom boundary layer
399   !!----------------------------------------------------------------------
400#  include "trcbbl_adv.h90"
401# else
402   !!----------------------------------------------------------------------
403   !!   Default option :                 NO advective bottom boundary layer
404   !!----------------------------------------------------------------------
405   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
406      INTEGER, INTENT(in) :: kt
407      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
408   END SUBROUTINE trc_bbl_adv
409# endif
410
411   SUBROUTINE trc_bbl_init
412      !!----------------------------------------------------------------------
413      !!                  ***  ROUTINE trc_bbl_init  ***
414      !!
415      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
416      !!
417      !! ** Method  :   Read the namtrcbbl namelist and check the parameters
418      !!      called by tra_bbl at the first timestep (nittrc000)
419      !!
420      !! History :
421      !!    8.5  !  02-08  (G. Madec)  Original code
422      !!----------------------------------------------------------------------
423      !! * Local declarations
424      INTEGER ::   ji, jj      ! dummy loop indices
425
426      !!----------------------------------------------------------------------
427
428
429      ! Parameter control and print
430      ! ---------------------------
431      IF(lwp) THEN
432         WRITE(numout,*)
433         WRITE(numout,*) 'trc_bbl_init : * Diffusive Bottom Boundary Layer'
434         WRITE(numout,*) '~~~~~~~~~~~~'
435# if defined key_trcbbl_adv
436            WRITE(numout,*) '               * Advective Bottom Boundary Layer'
437# endif
438         WRITE(numout,*)
439      ENDIF
440 
441      DO jj = 1, jpj
442         DO ji = 1, jpi
443            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
444         END DO
445      END DO
446      DO jj = 1, jpjm1
447         DO ji = 1, jpim1
448            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
449            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
450         END DO
451      END DO
452!!bug ???
453!!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0)
454
455# if defined key_trcbbl_adv
456      w_trc_bbl(:,:,:) = 0.e0    ! initialisation of w_trc_bbl to zero
457# endif
458
459   END SUBROUTINE trc_bbl_init
460
461#else
462   !!----------------------------------------------------------------------
463   !!   Dummy module :                      No bottom boundary layer scheme
464   !!----------------------------------------------------------------------
465   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .FALSE.   !: diff bbl flag
466   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .FALSE.   !: adv  bbl flag
467CONTAINS
468   SUBROUTINE trc_bbl_dif (kt )              ! Empty routine
469      INTEGER, INTENT(in) :: kt
470      WRITE(*,*) 'trc_bbl_dif: You should not have seen this print! error?', kt
471   END SUBROUTINE trc_bbl_dif
472   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
473      INTEGER, INTENT(in) :: kt
474      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
475   END SUBROUTINE trc_bbl_adv
476#endif
477
478   !!======================================================================
479END MODULE trcbbl
Note: See TracBrowser for help on using the repository browser.