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

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation

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