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 NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90 @ 14537

Last change on this file since 14537 was 14537, checked in by hadcv, 3 years ago

#2600: Reorganise dom_tile code

  • Property svn:keywords set to Id
File size: 31.9 KB
RevLine 
[3]1MODULE trabbl
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
[2528]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
[3764]10   !!            3.3  ! 2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   ! 2010-04  (G. Madec)  Campin & Goosse advective bbl
[2528]12   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
[4990]14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
[9019]15   !!            4.0  ! 2017-04  (G. Madec)  ln_trabbl namelist variable instead of a CPP key
[503]16   !!----------------------------------------------------------------------
[9019]17
[3]18   !!----------------------------------------------------------------------
[2715]19   !!   tra_bbl_alloc : allocate trabbl arrays
[2528]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
[503]25   !!----------------------------------------------------------------------
[2528]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
[5836]30   USE trd_oce        ! trends: ocean variables
[2528]31   USE trdtra         ! trends: active tracers
[4990]32   !
[14072]33   USE iom            ! IOM library
[2528]34   USE in_out_manager ! I/O manager
35   USE lbclnk         ! ocean lateral boundary conditions
36   USE prtctl         ! Print control
[3294]37   USE timing         ! Timing
[14072]38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3]39
40   IMPLICIT NONE
41   PRIVATE
42
[2528]43   PUBLIC   tra_bbl       !  routine called by step.F90
[9124]44   PUBLIC   tra_bbl_init  !  routine called by nemogcm.F90
[2528]45   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
[9124]46   PUBLIC   tra_bbl_adv   !     -      -          -
[2528]47   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
[3]48
[4147]49   !                                !!* Namelist nambbl *
[9019]50   LOGICAL , PUBLIC ::   ln_trabbl   !: bottom boundary layer flag
[4147]51   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
52   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
[2528]53   !                                            !  =1 : advective bbl using the bottom ocean velocity
54   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
[4147]55   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
56   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
[409]57
[4990]58   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
[3764]59
[2715]60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
[3]62
[3764]63   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
64   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
[3]67
68   !! * Substitutions
[12377]69#  include "do_loop_substitute.h90"
[13237]70#  include "domzgr_substitute.h90"
[3]71   !!----------------------------------------------------------------------
[9598]72   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]73   !! $Id$
[10068]74   !! Software governed by the CeCILL license (see ./LICENSE)
[3]75   !!----------------------------------------------------------------------
76CONTAINS
77
[2715]78   INTEGER FUNCTION tra_bbl_alloc()
79      !!----------------------------------------------------------------------
80      !!                ***  FUNCTION tra_bbl_alloc  ***
81      !!----------------------------------------------------------------------
[9019]82      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
83         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(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) ,                                    STAT=tra_bbl_alloc )
[2715]86         !
[10425]87      CALL mpp_sum ( 'trabbl', tra_bbl_alloc )
[2715]88      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
89   END FUNCTION tra_bbl_alloc
90
91
[12377]92   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
[3]93      !!----------------------------------------------------------------------
[2528]94      !!                  ***  ROUTINE bbl  ***
[3764]95      !!
96      !! ** Purpose :   Compute the before tracer (t & s) trend associated
[2528]97      !!              with the bottom boundary layer and add it to the general
98      !!              trend of tracer equations.
[3]99      !!
[2528]100      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
101      !!              diffusive and/or advective contribution to the tracer trend
102      !!              is added to the general tracer trend
[3764]103      !!----------------------------------------------------------------------
[12377]104      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step
105      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
106      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
[4990]107      !
[13982]108      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
[9019]109      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
[2528]110      !!----------------------------------------------------------------------
[3294]111      !
[9019]112      IF( ln_timing )   CALL timing_start( 'tra_bbl')
[3294]113      !
[9019]114      IF( l_trdtra )   THEN                         !* Save the T-S input trends
[13982]115         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
[12377]116         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
117         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
[2528]118      ENDIF
119
[12377]120      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
[3764]121
[4990]122      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
[3294]123         !
[12377]124         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
125         IF( sn_cfctl%l_prtctl )  &
126         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
127            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[13982]128         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
[14537]129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
[13982]131         ENDIF
[3294]132         !
[9168]133      ENDIF
[6140]134      !
[4990]135      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
[3294]136         !
[12377]137         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
138         IF(sn_cfctl%l_prtctl)   &
[13982]139         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, &
[12377]140            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[14537]141         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
[13982]142            ! lateral boundary conditions ; just need for outputs
143            CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
[14537]144         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
145         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
[13982]146         ENDIF
[3294]147         !
[9168]148      ENDIF
[2528]149
[6140]150      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
[12377]151         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
152         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
153         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
154         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
[9019]155         DEALLOCATE( ztrdt, ztrds )
[2528]156      ENDIF
157      !
[9019]158      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
[3294]159      !
[2528]160   END SUBROUTINE tra_bbl
161
162
[12377]163   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
[2528]164      !!----------------------------------------------------------------------
165      !!                  ***  ROUTINE tra_bbl_dif  ***
[3764]166      !!
[2528]167      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
[3764]168      !!                advection terms.
[2528]169      !!
[4990]170      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
[2528]171      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
172      !!      along bottom slope gradient) an additional lateral 2nd order
173      !!      diffusion along the bottom slope is added to the general
174      !!      tracer trend, otherwise the additional trend is set to 0.
175      !!      A typical value of ahbt is 2000 m2/s (equivalent to
[3]176      !!      a downslope velocity of 20 cm/s if the condition for slope
177      !!      convection is satified)
178      !!
[12377]179      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
[3]180      !!
[503]181      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
[2528]182      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]183      !!----------------------------------------------------------------------
[2715]184      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
[12377]185      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
186      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
187      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
[2715]188      !
[2528]189      INTEGER  ::   ji, jj, jn   ! dummy loop indices
190      INTEGER  ::   ik           ! local integers
191      REAL(wp) ::   zbtr         ! local scalars
[13982]192      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zptb   ! workspace
[3]193      !!----------------------------------------------------------------------
[2528]194      !
195      DO jn = 1, kjpt                                     ! tracer loop
196         !                                                ! ===========
[13295]197         DO_2D( 1, 1, 1, 1 )
[12377]198            ik = mbkt(ji,jj)                             ! bottom T-level index
199            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
200         END_2D
[14072]201         !
[13497]202         DO_2D( 0, 0, 0, 0 )                               ! Compute the trend
[12377]203            ik = mbkt(ji,jj)                            ! bottom T-level index
204            pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
205               &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
206               &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
207               &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
208               &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
209               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
210         END_2D
[2528]211         !                                                  ! ===========
212      END DO                                                ! end tracer
213      !                                                     ! ===========
214   END SUBROUTINE tra_bbl_dif
[3]215
[3764]216
[12377]217   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
[2528]218      !!----------------------------------------------------------------------
219      !!                  ***  ROUTINE trc_bbl  ***
220      !!
[3764]221      !! ** Purpose :   Compute the before passive tracer trend associated
[2528]222      !!     with the bottom boundary layer and add it to the general trend
223      !!     of tracer equations.
224      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
225      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
[3764]226      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
227      !!                       transport proportional to the along-slope density gradient
[2528]228      !!
229      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
230      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]231      !!----------------------------------------------------------------------
[2715]232      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
[12377]233      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
235      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
[2715]236      !
[2528]237      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
238      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
239      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
240      REAL(wp) ::   zbtr, ztra               ! local scalars
241      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
242      !!----------------------------------------------------------------------
243      !                                                          ! ===========
244      DO jn = 1, kjpt                                            ! tracer loop
[3764]245         !                                                       ! ===========
[14537]246         DO_2D_OVR( 1, 0, 1, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
[13982]247            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
248               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
249               iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
250               ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
251               zu_bbl = ABS( utr_bbl(ji,jj) )
[2528]252               !
[13982]253               !                                               ! up  -slope T-point (shelf bottom point)
254               zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
255               ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
256               pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
257               !
258               DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
259                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
260                  ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
261                  pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
262               END DO
263               !
264               zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
265               ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
266               pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
267            ENDIF
[2528]268            !
[13982]269            IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
270               ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
271               ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
272               ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
273               zv_bbl = ABS( vtr_bbl(ji,jj) )
274               !
275               ! up  -slope T-point (shelf bottom point)
276               zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
277               ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
278               pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
279               !
280               DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
281                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
282                  ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
283                  pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
284               END DO
285               !                                               ! down-slope T-point (deep bottom point)
286               zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
287               ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
288               pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
289            ENDIF
290         END_2D
291         !                                                       ! ===========
292      END DO                                                     ! end tracer
293      !                                                          ! ===========
[2528]294   END SUBROUTINE tra_bbl_adv
[3]295
296
[12377]297   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
[2528]298      !!----------------------------------------------------------------------
299      !!                  ***  ROUTINE bbl  ***
[3764]300      !!
[2528]301      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
[3764]302      !!                advection terms.
[2528]303      !!
[4990]304      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
[2528]305      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
306      !!      along bottom slope gradient) an additional lateral 2nd order
307      !!      diffusion along the bottom slope is added to the general
308      !!      tracer trend, otherwise the additional trend is set to 0.
309      !!      A typical value of ahbt is 2000 m2/s (equivalent to
310      !!      a downslope velocity of 20 cm/s if the condition for slope
311      !!      convection is satified)
[4990]312      !!              * advective bbl (nn_bbl_adv=1 or 2) :
[2528]313      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
314      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
315      !!        i.e. transport proportional to the along-slope density gradient
316      !!
317      !!      NB: the along slope density gradient is evaluated using the
318      !!      local density (i.e. referenced at a common local depth).
319      !!
320      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
321      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]322      !!----------------------------------------------------------------------
[2528]323      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
[4990]324      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
[2528]325      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
[12377]326      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
[6140]327      !
[2528]328      INTEGER  ::   ji, jj                    ! dummy loop indices
329      INTEGER  ::   ik                        ! local integers
[4990]330      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
331      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
332      REAL(wp) ::   za, zb, zgdrho            ! local scalars
333      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
[13982]334      REAL(wp), DIMENSION(A2D(nn_hls),jpts)   :: zts, zab         ! 3D workspace
335      REAL(wp), DIMENSION(A2D(nn_hls))        :: zub, zvb, zdep   ! 2D workspace
[2528]336      !!----------------------------------------------------------------------
[3294]337      !
[14537]338      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
[13982]339         IF( kt == kit000 )  THEN
340            IF(lwp)  WRITE(numout,*)
341            IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
342            IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
343         ENDIF
[2528]344      ENDIF
[4990]345      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
[13295]346      DO_2D( 1, 1, 1, 1 )
[12377]347         ik = mbkt(ji,jj)                             ! bottom T-level index
348         zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
349         zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
350         !
351         zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
352         zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
353         zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
354      END_2D
[4990]355      !
[12377]356      CALL eos_rab( zts, zdep, zab, Kmm )
[4990]357      !
[2528]358      !                                   !-------------------!
359      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
360         !                                !-------------------!
[13497]361         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
[12377]362            !                                                   ! i-direction
363            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
364            zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
365            !                                                         ! 2*masked bottom density gradient
366            zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
367               &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
368            !
[13226]369            zsign  = SIGN(  0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
[12377]370            ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
371            !
372            !                                                   ! j-direction
373            za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
374            zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
375            !                                                         ! 2*masked bottom density gradient
376            zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
377               &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
378            !
[13226]379            zsign = SIGN(  0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
[12377]380            ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
381         END_2D
[1601]382         !
[2528]383      ENDIF
[6140]384      !
[2528]385      !                                   !-------------------!
386      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
387         !                                !-------------------!
388         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
[503]389         !
[2528]390         CASE( 1 )                                   != use of upper velocity
[13497]391            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
[12377]392               !                                                  ! i-direction
393               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
394               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
[14072]395               !                                                          ! 2*masked bottom density gradient
[12377]396               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
397                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
398               !
[13226]399               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
400               zsigna= SIGN(  0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
[12377]401               !
402               !                                                          ! bbl velocity
403               utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
404               !
405               !                                                  ! j-direction
406               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
407               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
408               !                                                          ! 2*masked bottom density gradient
409               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
410                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
[13226]411               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
412               zsigna= SIGN(  0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
[12377]413               !
414               !                                                          ! bbl transport
415               vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
416            END_2D
[503]417            !
[2528]418         CASE( 2 )                                 != bbl velocity = F( delta rho )
419            zgbbl = grav * rn_gambbl
[13497]420            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
[12377]421               !                                                  ! i-direction
422               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
423               iid  = ji + MAX( 0, mgrhu(ji,jj) )
424               iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
425               !
426               ikud = mbku_d(ji,jj)
427               ikus = mbku(ji,jj)
428               !
429               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
430               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
431               !                                                          !   masked bottom density gradient
432               zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
433                  &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
434               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
435               !
436               !                                                          ! bbl transport (down-slope direction)
437               utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
438               !
439               !                                                  ! j-direction
440               !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
441               ijd  = jj + MAX( 0, mgrhv(ji,jj) )
442               ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
443               !
444               ikvd = mbkv_d(ji,jj)
445               ikvs = mbkv(ji,jj)
446               !
447               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
448               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
449               !                                                          !   masked bottom density gradient
450               zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
451                  &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
452               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
453               !
454               !                                                          ! bbl transport (down-slope direction)
455               vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
456            END_2D
[3]457         END SELECT
[2528]458         !
[3]459      ENDIF
[503]460      !
[2528]461   END SUBROUTINE bbl
[3]462
463
464   SUBROUTINE tra_bbl_init
465      !!----------------------------------------------------------------------
466      !!                  ***  ROUTINE tra_bbl_init  ***
467      !!
468      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
469      !!
470      !! ** Method  :   Read the nambbl namelist and check the parameters
[3294]471      !!              called by nemo_init at the first timestep (kit000)
[3]472      !!----------------------------------------------------------------------
[9019]473      INTEGER ::   ji, jj                      ! dummy loop indices
474      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
[9094]475      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
[9019]476      !!
477      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]478      !!----------------------------------------------------------------------
[3294]479      !
[4147]480      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
[11536]481901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
[6140]482      !
[4147]483      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
[11536]484902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
[4624]485      IF(lwm) WRITE ( numond, nambbl )
[2528]486      !
487      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
488      !
489      IF(lwp) THEN                   !* Parameter control and print
[3]490         WRITE(numout,*)
[2528]491         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]492         WRITE(numout,*) '~~~~~~~~~~~~'
[9019]493         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
494         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
[3]495      ENDIF
[9019]496      IF( .NOT.ln_trabbl )   RETURN
497      !
498      IF(lwp) THEN
499         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
500         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
501         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
502         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
503      ENDIF
504      !
[2715]505      !                              ! allocate trabbl arrays
506      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
[9019]507      !
[13532]508      IF(lwp) THEN
509         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
510         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
511      ENDIF
[9019]512      !
[2528]513      !                             !* vertical index of  "deep" bottom u- and v-points
[13497]514      DO_2D( 1, 0, 1, 0 )                 ! (the "shelf" bottom k-indices are mbku and mbkv)
[12377]515         mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
516         mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
517      END_2D
[4990]518      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[14072]519      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )
520      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp)
[9919]521      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
[9019]522      !
[9168]523      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
[8509]524      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
[13295]525      DO_2D( 1, 0, 1, 0 )
[12377]526         IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
[13226]527            mgrhu(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
[12377]528         ENDIF
529         !
530         IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
[13226]531            mgrhv(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
[12377]532         ENDIF
533      END_2D
[7646]534      !
[13497]535      DO_2D( 1, 0, 1, 0 )           !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0)
[12377]536         e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
537         e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
538      END_2D
[13226]539      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp )      ! lateral boundary conditions
[7646]540      !
[3764]541      !                             !* masked diffusive flux coefficients
[7753]542      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
543      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
[503]544      !
[3]545   END SUBROUTINE tra_bbl_init
546
547   !!======================================================================
548END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.