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

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

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 18.7 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   !!   TOP 1.0 , LOCEAN-IPSL (2005)
51   !! $Header$
52   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE trc_bbl_dif( kt )
58      !!----------------------------------------------------------------------
59      !!                   ***  ROUTINE trc_bbl_dif  ***
60      !!
61      !! ** Purpose :   Compute the before tracer trend associated
62      !!      with the bottom boundary layer and add it to the general trend
63      !!      of tracer equations. The bottom boundary layer is supposed to be
64      !!      a purely diffusive bottom boundary layer.
65      !!
66      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad
67      !!      is an along bottom slope gradient) an additional lateral diffu-
68      !!      sive trend along the bottom slope is added to the general tracer
69      !!      trend, otherwise nothing is done.
70      !!      Second order operator (laplacian type) with variable coefficient
71      !!      computed as follow for temperature (idem on s):
72      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
73      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
74      !!      where ztb is a 2D array: the bottom ocean temperature and ahtb
75      !!      is a time and space varying diffusive coefficient defined by:
76      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
77      !!              = 0.       otherwise.
78      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
79      !!      is evaluated using the local density (i.e. referenced at the
80      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
81      !!      a downslope velocity of 20 cm/s if the condition for slope
82      !!      convection is satified)
83      !!      Add this before trend to the general trend tra of the
84      !!      botton ocean tracer point:
85      !!         tra = tra + difft
86      !!
87      !! ** Action  : - update tra at the bottom level with the bottom
88      !!                boundary layer trend
89      !!
90      !! References :
91      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
92      !!
93      !! History :
94      !!   8.0  !  96-06  (L. Mortier)  Original code
95      !!   8.0  !  97-11  (G. Madec)  Optimization
96      !!   8.5  !  02-08  (G. Madec)  free form + modules
97      !!   9.0  !  04-03  (C. Ethe)   Adaptation for passive tracers
98      !!----------------------------------------------------------------------
99      !! * Arguments
100      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
101
102      !! * Local declarations
103      INTEGER ::   ji, jj,jn                ! dummy loop indices
104      INTEGER ::   ik
105      INTEGER ::   ii0, ii1, ij0, ij1       ! temporary integers
106#  if defined key_partial_steps
107      INTEGER  ::   iku1, iku2, ikv1,ikv2   ! temporary intergers
108      REAL(wp) ::   ze3u, ze3v              ! temporary scalars
109#  else
110      INTEGER ::   iku, ikv
111#  endif
112      REAL(wp) ::   &
113         zsign, zt, zs, zh, zalbet,      &  ! temporary scalars
114         zgdrho, zbtr, ztra
115      REAL(wp), DIMENSION(jpi,jpj) ::    &
116        zki, zkj, zkx, zky,    &  ! temporary workspace arrays
117        ztnb, zsnb, zdep,                &
118        ztrb, zahu, zahv
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      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
207      ! --------------------------------------------
208      ! Sign of the local density gradient along the i- and j-slopes
209      ! multiplied by the slope of the ocean bottom
210
211#  if defined key_vectopt_loop   &&   ! defined key_autotasking
212      jj = 1
213      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
214#  else
215      DO jj = 1, jpjm1
216         DO ji = 1, jpim1
217#  endif
218            ! temperature, salinity anomalie and depth
219            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
220            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
221            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
222            ! masked ratio alpha/beta
223            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
224            ! local density gradient along i-bathymetric slope
225            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
226                   -          ( zsnb(ji+1,jj) - zsnb(ji,jj) )
227            ! sign of local i-gradient of density multiplied by the i-slope
228            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
229            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
230#  if ! defined key_vectopt_loop   ||   defined key_autotasking
231         END DO
232#  endif
233      END DO
234
235#  if defined key_vectopt_loop   &&   ! defined key_autotasking
236      jj = 1
237      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
238#  else
239      DO jj = 1, jpjm1
240         DO ji = 1, jpim1
241#  endif
242            ! temperature, salinity anomalie and depth
243            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
244            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
245            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
246            ! masked ratio alpha/beta
247            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
248            ! local density gradient along j-bathymetric slope
249            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
250                   -          ( zsnb(ji,jj+1) - zsnb(ji,jj) )
251            ! sign of local j-gradient of density multiplied by the j-slope
252            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
253            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
254#  if ! defined key_vectopt_loop   ||   defined key_autotasking
255         END DO
256#  endif
257      END DO
258
259
260      ! 2. Additional second order diffusive trends
261      ! -------------------------------------------
262
263      DO jn = 1, jptra
264         ! first derivative (gradient)
265         
266#  if defined key_vectopt_loop   &&   ! defined key_autotasking
267         jj = 1
268         DO ji = 1, jpij   ! vector opt. (forced unrolling)
269#  else
270         DO jj = 1, jpj
271            DO ji = 1, jpi
272#  endif
273               ik = mbkt(ji,jj) 
274               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
275#  if ! defined key_vectopt_loop   ||   defined key_autotasking
276            END DO
277#  endif
278         END DO
279
280#  if defined key_vectopt_loop   &&   ! defined key_autotasking
281         jj = 1
282         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
283#  else
284         DO jj = 1, jpjm1
285            DO ji = 1, jpim1
286#  endif
287               zkx(ji,jj) = zki(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
288               zky(ji,jj) = zkj(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
289#  if ! defined key_vectopt_loop   ||   defined key_autotasking
290            END DO
291#  endif
292         END DO
293
294         IF( cp_cfg == "orca" ) THEN
295           
296            SELECT CASE ( jp_cfg )
297               !                                           ! =======================
298            CASE ( 2 )                                  !  ORCA_R2 configuration
299               !                                        ! =======================
300               ! Gibraltar enhancement of BBL
301               ij0 = 102   ;   ij1 = 102
302               ii0 = 139   ;   ii1 = 140 
303               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
304               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
305               
306               ! Red Sea enhancement of BBL
307               ij0 =  88   ;   ij1 =  88
308               ii0 = 161   ;   ii1 = 162
309               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
310               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
311               
312               !                                        ! =======================
313            CASE ( 4 )                                  !  ORCA_R4 configuration
314               !                                        ! =======================
315               ! Gibraltar enhancement of BBL
316               ij0 =  52   ;   ij1 =  52
317               ii0 =  70   ;   ii1 =  71 
318               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
319               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
320               
321            END SELECT
322           
323         ENDIF
324
325         
326         ! second derivative (divergence) and add to the general tracer trend
327#  if defined key_vectopt_loop   &&   ! defined key_autotasking
328         jj = 1
329         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
330#  else
331         DO jj = 2, jpjm1
332            DO ji = 2, jpim1
333#  endif
334               ik = MAX( mbathy(ji,jj)-1, 1 )
335               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
336               ztra = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
337                  &    + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
338               tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
339#  if ! defined key_vectopt_loop   ||   defined key_autotasking
340            END DO
341#  endif
342         END DO
343         IF(l_ctl) THEN         ! print mean trends (used for debugging)
344            ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
345            WRITE(numout,*) ' trc/bbl  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn) 
346            tra_ctl(jn) = ztra 
347         ENDIF
348      END DO
349
350
351   END SUBROUTINE trc_bbl_dif
352
353# if defined key_trcbbl_adv
354   !!----------------------------------------------------------------------
355   !!   'key_trcbbl_adv'                    advective bottom boundary layer
356   !!----------------------------------------------------------------------
357#  include "trcbbl_adv.h90"
358# else
359   !!----------------------------------------------------------------------
360   !!   Default option :                 NO advective bottom boundary layer
361   !!----------------------------------------------------------------------
362   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
363      INTEGER, INTENT(in) :: kt
364      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
365   END SUBROUTINE trc_bbl_adv
366# endif
367
368   SUBROUTINE trc_bbl_init
369      !!----------------------------------------------------------------------
370      !!                  ***  ROUTINE trc_bbl_init  ***
371      !!
372      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
373      !!
374      !! ** Method  :   Read the namtrcbbl namelist and check the parameters
375      !!      called by tra_bbl at the first timestep (nittrc000)
376      !!
377      !! History :
378      !!    8.5  !  02-08  (G. Madec)  Original code
379      !!----------------------------------------------------------------------
380      !! * Local declarations
381      INTEGER ::   ji, jj      ! dummy loop indices
382
383      !!----------------------------------------------------------------------
384
385
386      ! Parameter control and print
387      ! ---------------------------
388      IF(lwp) THEN
389         WRITE(numout,*)
390         WRITE(numout,*) 'trc_bbl_init : * Diffusive Bottom Boundary Layer'
391         WRITE(numout,*) '~~~~~~~~~~~~'
392# if defined key_trcbbl_adv
393            WRITE(numout,*) '               * Advective Bottom Boundary Layer'
394# endif
395         WRITE(numout,*)
396      ENDIF
397 
398      DO jj = 1, jpj
399         DO ji = 1, jpi
400            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
401         END DO
402      END DO
403      DO jj = 1, jpjm1
404         DO ji = 1, jpim1
405            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
406            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
407         END DO
408      END DO
409!!bug ???
410!!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0)
411
412# if defined key_trcbbl_adv
413      w_trc_bbl(:,:,:) = 0.e0    ! initialisation of w_trc_bbl to zero
414# endif
415
416   END SUBROUTINE trc_bbl_init
417
418#else
419   !!----------------------------------------------------------------------
420   !!   Dummy module :                      No bottom boundary layer scheme
421   !!----------------------------------------------------------------------
422   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .FALSE.   !: diff bbl flag
423   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .FALSE.   !: adv  bbl flag
424CONTAINS
425   SUBROUTINE trc_bbl_dif (kt )              ! Empty routine
426      INTEGER, INTENT(in) :: kt
427      WRITE(*,*) 'trc_bbl_dif: You should not have seen this print! error?', kt
428   END SUBROUTINE trc_bbl_dif
429   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
430      INTEGER, INTENT(in) :: kt
431      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
432   END SUBROUTINE trc_bbl_adv
433#endif
434
435   !!======================================================================
436END MODULE trcbbl
Note: See TracBrowser for help on using the repository browser.