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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2104

Last change on this file since 2104 was 2104, checked in by cetlod, 14 years ago

update DEV_r2006_merge_TRA_TRC according to review

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 36.3 KB
Line 
1MODULE trabbl
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
6   !! History :  OPA  !  1996-06  (L. Mortier)  Original code
7   !!            8.0  !  1997-11  (G. Madec)    Optimization
8   !!   NEMO     1.0  !  2002-08  (G. Madec)  free form + modules
9   !!             -   !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
10   !!            3.3  !  2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   !  2010-04  (G. Madec)  Campin & Goosse advective bbl
12   !!             -   !  2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!----------------------------------------------------------------------
14#if   defined key_trabbl   ||   defined key_esopa
15   !!----------------------------------------------------------------------
16   !!   'key_trabbl'   or                             bottom boundary layer
17   !!----------------------------------------------------------------------
18   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
19   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
20   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
21   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
22   !!   tra_bbl_init  : initialization, namlist read, parameters control
23   !!----------------------------------------------------------------------
24   USE oce            ! ocean dynamics and active tracers
25   USE dom_oce        ! ocean space and time domain
26   USE phycst         !
27   USE eosbn2         ! equation of state
28   USE trdmod_oce     ! ocean space and time domain
29   USE trdtra         ! ocean active tracers trends
30   USE iom            ! IOM server               
31   USE in_out_manager ! I/O manager
32   USE lbclnk         ! ocean lateral boundary conditions
33   USE prtctl         ! Print control
34   USE trc_oce        ! share passive tracers/Ocean variables
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_bbl       !  routine called by step.F90
40   PUBLIC   tra_bbl_init  !  routine called by opa.F90
41   PUBLIC   tra_bbl_dif   !  routine called by tra_bbl and trc_bbl
42   PUBLIC   tra_bbl_adv   !  -          -          -              -
43   PUBLIC   bbl           !  -          -          -              -
44
45# if defined key_trabbl
46   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
47# else
48   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bottom boundary layer flag
49# endif
50
51   LOGICAL, PUBLIC              ::   l_bbl               !: flag to compute bbl diffu. flux coef and transport
52   
53   !                                         !!* Namelist nambbl *
54   INTEGER , PUBLIC ::   nn_bbl_ldf = 0       !: =1   : diffusive bbl or not (=0)
55   INTEGER , PUBLIC ::   nn_bbl_adv = 0       !: =1/2 : advective bbl or not (=0)
56   !                                          !  =1 : advective bbl using the model velocity
57   !                                          !  =2 :   -        -  using utr_bbl proportional to grad(rho)
58   REAL(wp), PUBLIC ::   rn_ahtbbl  = 1.e+3   !: along slope bbl diffusive coefficient [m2/s]
59   REAL(wp), PUBLIC ::   rn_gambbl  = 10.e0   !: lateral coeff. for bottom boundary layer scheme [s]
60
61   REAL(wp), DIMENSION(jpi,jpj), PUBLIC ::   utr_bbl, vtr_bbl   ! u- (v-) transport in the bottom boundary layer
62   
63   INTEGER , DIMENSION(jpi,jpj) ::   mbkt                   ! vertical index of the bottom ocean T-level
64   INTEGER , DIMENSION(jpi,jpj) ::   mbku     , mbkv        ! vertical index of the (upper) bottom ocean U/V-level
65   INTEGER , DIMENSION(jpi,jpj) ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level
66   INTEGER , DIMENSION(jpi,jpj) ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction
67   REAL(wp), DIMENSION(jpi,jpj) ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
68   REAL(wp), DIMENSION(jpi,jpj) ::   ahu_bbl  , ahv_bbl     ! masked diffusive bbl coefficients at u and v-points
69   REAL(wp), DIMENSION(jpi,jpj) ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points
70   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2t_r                ! thichness of the bbl (e3) at u and v-points
71
72   !! * Substitutions
73#  include "domzgr_substitute.h90"
74#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
77   !! $Id$
78   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80
81CONTAINS
82
83   SUBROUTINE tra_bbl( kt )
84      !!----------------------------------------------------------------------
85      !!                  ***  ROUTINE bbl  ***
86      !!                   
87      !! ** Purpose :   Compute the before tracer (t & s) trend associated
88      !!     with the bottom boundary layer and add it to the general trend
89      !!     of tracer equations.
90      !!
91      !!---------------------------------------------------------------------- 
92      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
93      !!
94      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
95      !!----------------------------------------------------------------------
96
97      IF( l_trdtra )   THEN                    !* Save ta and sa trends
98         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
99         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
100      ENDIF
101
102      !* bbl coef and transport are computed only if not already done in passive tracers routine
103      IF( l_bbl )      CALL bbl( kt, 'TRA' ) 
104
105      !* Diffusive bbl :
106      IF( nn_bbl_ldf == 1 ) THEN
107         CALL tra_bbl_dif( tsb, tsa, jpts )
108         IF( ln_ctl )  &
109         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
110         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
111         ! lateral boundary conditions ; just need for outputs                         
112         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
113         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef     
114         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
115      END IF
116
117      !* Advective bbl : bbl upstream advective trends added to the tracer trends
118      IF( nn_bbl_adv /= 0 ) THEN
119         CALL tra_bbl_adv( tsb, tsa, jpts ) 
120         IF(ln_ctl)   &
121         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
122         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
123         ! lateral boundary conditions ; just need for outputs                         
124         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )   
125         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport     
126         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
127      END IF
128
129      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
130         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
131         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
132         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
133         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
134         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
135      ENDIF
136      !
137   END SUBROUTINE tra_bbl
138
139
140   SUBROUTINE tra_bbl_dif( ptrab, ptraa, kjpt )
141      !!----------------------------------------------------------------------
142      !!                  ***  ROUTINE tra_bbl_dif  ***
143      !!                   
144      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
145      !!                advection terms.
146      !!
147      !! ** Method  :   
148      !!        * diffusive bbl (nn_bbl_ldf=1) :
149      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
150      !!      along bottom slope gradient) an additional lateral 2nd order
151      !!      diffusion along the bottom slope is added to the general
152      !!      tracer trend, otherwise the additional trend is set to 0.
153      !!      A typical value of ahbt is 2000 m2/s (equivalent to
154      !!      a downslope velocity of 20 cm/s if the condition for slope
155      !!      convection is satified)
156      !!
157      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
158      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
159      !!---------------------------------------------------------------------- 
160      INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers
161      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields
162      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend
163      !!
164      INTEGER  ::   ji, jj, jn   ! dummy loop indices
165      INTEGER  ::   ik           ! local integer
166      REAL(wp) ::   zbtr, ztra   ! local scalars
167      REAL(wp), DIMENSION(jpi,jpj) :: ztrb, zkx, zky   ! 2D workspace
168      !!----------------------------------------------------------------------
169      !
170      !                                                   ! ===========
171      DO jn = 1, kjpt                                     ! tracer loop
172         !                                                ! ===========
173#if defined key_vectopt_loop
174         DO jj = 1, 1   ! vector opt. (forced unrolling)
175            DO ji = 1, jpij
176#else
177         DO jj = 1, jpj
178            DO ji = 1, jpi
179#endif
180               ik = mbkt(ji,jj)                        ! bottom T-level index
181               ztrb(ji,jj) = ptrab(ji,jj,ik,jn)              ! bottom before T and S
182            END DO
183         END DO
184         !
185!!gm  forced unrolling should be uuseless in the loop below (no indirect adressing)
186#  if defined key_vectopt_loop
187         DO jj = 1, 1   ! vector opt. (forced unrolling)
188            DO ji = 1, jpij-jpi
189#  else
190         DO jj = 1, jpjm1
191            DO ji = 1, jpim1
192#  endif
193               zkx(ji,jj) = ahu_bbl(ji,jj) * ( ztrb(ji+1,jj  ) - ztrb(ji,jj) )  ! diffusive i-flux
194               zky(ji,jj) = ahv_bbl(ji,jj) * ( ztrb(ji  ,jj+1) - ztrb(ji,jj) )  ! diffusive j-flux
195            END DO
196         END DO
197            !                                        ! Add the diffusive trends
198#  if defined key_vectopt_loop
199         DO jj = 1, 1   ! vector opt. (forced unrolling)
200            DO ji = jpi+1, jpij-jpi-1
201#  else
202         DO jj = 2, jpjm1
203            DO ji = 2, jpim1
204#  endif
205               ik = mbkt(ji,jj)
206               zbtr = e1e2t_r(ji,jj)  / fse3t(ji,jj,ik)
207               ztra = ( zkx(ji,jj) - zkx(ji-1,jj) + zky(ji,jj) - zky(ji,jj-1) ) * zbtr
208               ptraa(ji,jj,ik,jn) = ptraa(ji,jj,ik,jn) + ztra
209            END DO
210         END DO
211         !
212      END DO
213      !
214   END SUBROUTINE tra_bbl_dif
215   
216
217   SUBROUTINE tra_bbl_adv( ptrab, ptraa, kjpt )
218      !!----------------------------------------------------------------------
219      !!                  ***  ROUTINE trc_bbl  ***
220      !!                   
221      !! ** Purpose :   Compute the before passive tracer trend associated
222      !!     with the bottom boundary layer and add it to the general trend
223      !!     of tracer equations.
224      !!        * advective bbl (nn_bbl_adv=1 or 2) :
225      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
226      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
227      !!        i.e. transport proportional to the along-slope density gradient
228      !!
229      !!      NB: the along slope density gradient is evaluated using the
230      !!      local density (i.e. referenced at a common local depth).
231      !!
232      !!
233      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
234      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
235      !!
236      !!---------------------------------------------------------------------- 
237      INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers
238      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields
239      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend
240      !!
241      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
242      INTEGER  ::   ik                       ! local integers
243      INTEGER  ::   iis , iid , ijs , ijd    !   -       -
244      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
245      REAL(wp) ::   zbtr, ztra               ! local scalars
246      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
247      !!----------------------------------------------------------------------
248
249      !                                                          ! ===========
250      DO jn = 1, kjpt                                            ! tracer loop
251         !                                                       ! ===========         
252# if defined key_vectopt_loop
253         DO jj = 1, 1
254            DO ji = jpi+1, jpij-jpi-1   ! vector opt. (forced unrolling)
255# else
256         DO jj = 2, jpjm1
257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
258# endif
259               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
260                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
261                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
262                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
263                  zu_bbl = ABS( utr_bbl(ji,jj) )
264                  !
265                  !                                               ! up  -slope T-point (shelf bottom point)
266                  zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus)
267                  ztra = zu_bbl * ( ptrab(iid,jj,ikus,jn) - ptrab(iis,jj,ikus,jn) ) * zbtr
268                  ptraa(iis,jj,ikus,jn) = ptraa(iis,jj,ikus,jn) + ztra
269                  !                   
270                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
271                     zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk)
272                     ztra = zu_bbl * ( ptrab(iid,jj,jk+1,jn) - ptrab(iid,jj,jk,jn) ) * zbtr
273                     ptraa(iid,jj,jk,jn) = ptraa(iid,jj,jk,jn) + ztra
274                  END DO
275                  !
276                  zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud)
277                  ztra = zu_bbl * ( ptrab(iis,jj,ikus,jn) - ptrab(iid,jj,ikud,jn) ) * zbtr
278                  ptraa(iid,jj,ikud,jn) = ptraa(iid,jj,ikud,jn) + ztra
279               ENDIF
280               !
281               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
282                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
283                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
284                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
285                  zv_bbl = ABS( vtr_bbl(ji,jj) )
286                  !
287                  ! up  -slope T-point (shelf bottom point)
288                  zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs)
289                  ztra = zv_bbl * ( ptrab(ji,ijd,ikvs,jn) - ptrab(ji,ijs,ikvs,jn) ) * zbtr
290                  ptraa(ji,ijs,ikvs,jn) = ptraa(ji,ijs,ikvs,jn) + ztra
291                  !                   
292                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
293                     zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk)
294                     ztra = zv_bbl * ( ptrab(ji,ijd,jk+1,jn) - ptrab(ji,ijd,jk,jn) ) * zbtr
295                     ptraa(ji,ijd,jk,jn) = ptraa(ji,ijd,jk,jn)  + ztra
296                  END DO
297                  !                                               ! down-slope T-point (deep bottom point)
298                  zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd)
299                  ztra = zv_bbl * ( ptrab(ji,ijs,ikvs,jn) - ptrab(ji,ijd,ikvd,jn) ) * zbtr
300                  ptraa(ji,ijd,ikvd,jn) = ptraa(ji,ijd,ikvd,jn) + ztra
301               ENDIF
302            END DO
303            !
304         END DO
305            !
306      END DO
307      !
308   END SUBROUTINE tra_bbl_adv
309
310
311   SUBROUTINE bbl( kt, cdtype )
312      !!----------------------------------------------------------------------
313      !!                  ***  ROUTINE bbl  ***
314      !!                   
315      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
316      !!                advection terms.
317      !!
318      !! ** Method  :   
319      !!        * diffusive bbl (nn_bbl_ldf=1) :
320      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
321      !!      along bottom slope gradient) an additional lateral 2nd order
322      !!      diffusion along the bottom slope is added to the general
323      !!      tracer trend, otherwise the additional trend is set to 0.
324      !!      A typical value of ahbt is 2000 m2/s (equivalent to
325      !!      a downslope velocity of 20 cm/s if the condition for slope
326      !!      convection is satified)
327      !!        * advective bbl (nn_bbl_adv=1 or 2) :
328      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
329      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
330      !!        i.e. transport proportional to the along-slope density gradient
331      !!
332      !!      NB: the along slope density gradient is evaluated using the
333      !!      local density (i.e. referenced at a common local depth).
334      !!
335      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
336      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
337      !!---------------------------------------------------------------------- 
338      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
339      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
340      !!
341      INTEGER  ::   ji, jj                    ! dummy loop indices
342      INTEGER  ::   ik                        ! local integers
343      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
344      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
345      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
346      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
347      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb, ztb, zsb, zdep  !  2D workspace
348      !!
349      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
350      !!----------------------- zv_bbl-----------------------------------------------
351      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
352      ! ================            pft :  potential temperature in degrees celcius
353      !                             pfs :  salinity anomaly (s-35) in psu
354      !                             pfh :  depth in meters
355      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
356      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
357         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
358                                   - 0.203814e-03 ) * pft                    &
359                                   + 0.170907e-01 ) * pft                    &
360                                   + 0.665157e-01                            &
361         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
362         +  ( ( - 0.302285e-13 * pfh                                         &
363                - 0.251520e-11 * pfs                                         &
364                + 0.512857e-12 * pft * pft          ) * pfh                  &
365                                     - 0.164759e-06   * pfs                  &
366             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
367                                     + 0.380374e-04 ) * pfh
368      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
369         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
370                                 - 0.301985e-05 ) * pft                      &
371                                 + 0.785567e-03                              &
372         + (     0.515032e-08 * pfs                                          &
373               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
374               +(  (   0.121551e-17 * pfh                                    &
375                     - 0.602281e-15 * pfs                                    &
376                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
377                                          + 0.408195e-10   * pfs             &
378                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
379                                          - 0.121555e-07 ) * pfh
380      !!----------------------------------------------------------------------
381     
382      IF( kt == nit000 )  THEN
383         IF(lwp)  WRITE(numout,*)
384         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
385         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
386      ENDIF
387     
388      !                                        !* bottom temperature, salinity, velocity and depth
389#if defined key_vectopt_loop
390      DO jj = 1, 1   ! vector opt. (forced unrolling)
391         DO ji = 1, jpij
392#else
393      DO jj = 1, jpj
394         DO ji = 1, jpi
395#endif
396            ik = mbkt(ji,jj)                        ! bottom T-level index
397            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)      ! bottom before T and S
398            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal)
399            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
400            !
401            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
402            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
403         END DO
404      END DO
405     
406      !                                   !-------------------!
407      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
408         !                                !-------------------!
409         !                                        ! bbl diffusive fluxes
410         !                                             ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
411#  if defined key_vectopt_loop
412         DO jj = 1, 1   ! vector opt. (forced unrolling)
413            DO ji = 1, jpij-jpi
414#  else
415         DO jj = 1, jpjm1
416            DO ji = 1, jpim1
417#  endif
418               !                                                ! i-direction
419               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
420               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
421               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
422               !                                                         ! masked bbl i-gradient of density
423               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
424                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
425               !                                                     
426               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
427               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
428               !
429               !                                                ! j-direction
430               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
431               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
432               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
433               !                                                         ! masked bbl j-gradient of density
434               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
435                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
436               !                                                   
437               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
438               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
439               !
440            END DO
441         END DO
442         !
443      ENDIF
444
445      !                                   !-------------------!
446      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
447         !                                !-------------------!
448         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
449         !
450         CASE( 1 )                                   != use of upper velocity
451            DO jj = 2, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
452               DO ji = 1, fs_jpim1   ! vector opt.
453                  !                                               ! i-direction
454                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
455                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
456                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
457                  !                                                           ! masked bbl i-gradient of density
458                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
459                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
460                  !                                                         
461                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
462                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
463                  !
464                  !                                                           ! bbl velocity
465                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
466                  !
467                  !                                               ! j-direction
468                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
469                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
470                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
471                  !                                                           ! masked bbl j-gradient of density
472                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
473                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
474                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
475                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
476                  !
477                  !                                                           ! bbl velocity
478                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
479               END DO
480            END DO
481            !
482         CASE( 2 )                                 != bbl velocity = F( delta rho )
483            zgbbl = grav * rn_gambbl
484            DO jj = 2, jpjm1                            ! criteria: rho_up > rho_down
485               DO ji = 1, fs_jpim1   ! vector opt.
486                  !                                         ! i-direction
487                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
488                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
489                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
490                  !
491                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
492                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
493                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
494                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
495                  zgdrho =    fsbeta( zt, zs, zh )                                    &
496                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
497                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
498                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
499                  !
500                  !                                             ! bbl transport (down-slope direction)
501                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
502                  !
503                  !                                         ! j-direction
504                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
505                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
506                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
507                  !
508                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
509                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
510                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
511                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
512                  zgdrho =    fsbeta( zt, zs, zh )                                    &
513                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
514                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
515                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
516                  !
517                  !                                             ! bbl transport (down-slope direction)
518                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
519               END DO
520            END DO
521         END SELECT
522         !
523      ENDIF
524      !
525   END SUBROUTINE bbl
526
527
528   SUBROUTINE tra_bbl_init
529      !!----------------------------------------------------------------------
530      !!                  ***  ROUTINE tra_bbl_init  ***
531      !!
532      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
533      !!
534      !! ** Method  :   Read the nambbl namelist and check the parameters
535      !!              called by tra_bbl at the first timestep (nit000)
536      !!----------------------------------------------------------------------
537      INTEGER ::   ji, jj               ! dummy loop indices
538      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
539      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
540      !!
541      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
542      !!----------------------------------------------------------------------
543
544      REWIND ( numnam )              !* Read Namelist nambbl : bottom boundary layer scheme
545      READ   ( numnam, nambbl )
546
547      IF(lwp) THEN                   !* Parameter control and print
548         WRITE(numout,*)
549         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
550         WRITE(numout,*) '~~~~~~~~~~~~'
551         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
552         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
553         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
554         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
555         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
556      ENDIF
557     
558      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
559      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
560
561      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
562
563
564      !                             !* inverse of surface of T-cells
565      e1e2t_r(:,:) = 1.0 / ( e1t(:,:) * e2t(:,:) )
566     
567      !                             !* vertical index of bottom t-, u- and v-points
568      DO jj = 1, jpj                      ! bottom k-index of T-level
569         DO ji = 1, jpi
570            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )
571         END DO
572      END DO
573      DO jj = 1, jpjm1                    ! bottom k-index of u- (v-) level (shelf and deep)
574         DO ji = 1, jpim1
575            mbku  (ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )   ! "shelf"
576            mbkv  (ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
577            mbku_d(ji,jj) = MAX( MAX( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )   ! "deep"
578            mbkv_d(ji,jj) = MAX( MAX( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
579         END DO
580      END DO
581      zmbk(:,:) = FLOAT( mbku  (:,:) )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
582      zmbk(:,:) = FLOAT( mbkv  (:,:) )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
583      zmbk(:,:) = FLOAT( mbku_d(:,:) )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
584      zmbk(:,:) = FLOAT( mbkv_d(:,:) )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
585
586      DO jj = 1, jpj                !* sign of grad(H) at u- and v-points
587         DO ji = 1, jpi
588            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
589            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
590         END DO
591      END DO
592
593      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
594         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
595            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 
596            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 
597         END DO
598      END DO
599      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
600
601      !                             !* masked diffusive flux coefficients
602      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
603      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
604
605      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
606         !
607         SELECT CASE ( jp_cfg )
608         CASE ( 2 )                          ! ORCA_R2
609            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
610            ii0 = 139   ;   ii1 = 140 
611            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
612            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
613            !
614            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
615            ii0 = 161   ;   ii1 = 162
616            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
617            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
618            !
619         CASE ( 4 )                          ! ORCA_R4
620            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
621            ii0 =  70   ;   ii1 =  71 
622            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
623            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
624         END SELECT
625         !
626      ENDIF
627      !
628      l_bbl = .TRUE.     !: flag to compute bbl coef and transport
629      !
630   END SUBROUTINE tra_bbl_init
631
632#else
633   !!----------------------------------------------------------------------
634   !!   Dummy module :                      No bottom boundary layer scheme
635   !!----------------------------------------------------------------------
636   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
637CONTAINS
638   SUBROUTINE tra_bbl_init               ! Dummy routine
639   END SUBROUTINE tra_bbl_init
640   SUBROUTINE tra_bbl( kt )              ! Dummy routine
641      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
642   END SUBROUTINE tra_bbl
643#endif
644
645   !!======================================================================
646END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.