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

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

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