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_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago

Merge branch 'dynamic_memory' into master-svn-dyn

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