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/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/TRA/trabbl.F90 @ 13248

Last change on this file since 13248 was 13248, checked in by francesca, 4 years ago

dev_r12558_HPC-08_epico_Extra_Halo: merge with trunk@13237, see #2366

  • Property svn:keywords set to Id
File size: 31.3 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   !
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
[4990]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"
[13248]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      !
[9019]108      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
[2528]109      !!----------------------------------------------------------------------
[3294]110      !
[9019]111      IF( ln_timing )   CALL timing_start( 'tra_bbl')
[3294]112      !
[9019]113      IF( l_trdtra )   THEN                         !* Save the T-S input trends
114         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
[12377]115         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
116         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
[2528]117      ENDIF
118
[12377]119      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
[3764]120
[4990]121      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
[3294]122         !
[12377]123         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
124         IF( sn_cfctl%l_prtctl )  &
125         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
126            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[3764]127         ! lateral boundary conditions ; just need for outputs
[13247]128         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1.0_wp , ahv_bbl, 'V', 1.0_wp )
[3764]129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
[2528]130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
[3294]131         !
[9168]132      ENDIF
[6140]133      !
[4990]134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
[3294]135         !
[12377]136         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
137         IF(sn_cfctl%l_prtctl)   &
138         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
139            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[3764]140         ! lateral boundary conditions ; just need for outputs
[13247]141         CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
[3764]142         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
[2528]143         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
[3294]144         !
[9168]145      ENDIF
[2528]146
[6140]147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
[12377]148         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
149         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
150         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
151         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
[9019]152         DEALLOCATE( ztrdt, ztrds )
[2528]153      ENDIF
154      !
[9019]155      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
[3294]156      !
[2528]157   END SUBROUTINE tra_bbl
158
159
[12377]160   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
[2528]161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE tra_bbl_dif  ***
[3764]163      !!
[2528]164      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
[3764]165      !!                advection terms.
[2528]166      !!
[4990]167      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
[2528]168      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
169      !!      along bottom slope gradient) an additional lateral 2nd order
170      !!      diffusion along the bottom slope is added to the general
171      !!      tracer trend, otherwise the additional trend is set to 0.
172      !!      A typical value of ahbt is 2000 m2/s (equivalent to
[3]173      !!      a downslope velocity of 20 cm/s if the condition for slope
174      !!      convection is satified)
175      !!
[12377]176      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
[3]177      !!
[503]178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
[2528]179      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]180      !!----------------------------------------------------------------------
[2715]181      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
[12377]182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
184      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
[2715]185      !
[2528]186      INTEGER  ::   ji, jj, jn   ! dummy loop indices
187      INTEGER  ::   ik           ! local integers
188      REAL(wp) ::   zbtr         ! local scalars
[9019]189      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
[3]190      !!----------------------------------------------------------------------
[2528]191      !
192      DO jn = 1, kjpt                                     ! tracer loop
193         !                                                ! ===========
[12377]194         DO_2D_11_11
195            ik = mbkt(ji,jj)                             ! bottom T-level index
196            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
197         END_2D
[4990]198         !               
[12377]199         DO_2D_00_00
200            ik = mbkt(ji,jj)                            ! bottom T-level index
201            pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
202               &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
203               &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
204               &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
205               &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
206               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
207         END_2D
[2528]208         !                                                  ! ===========
209      END DO                                                ! end tracer
210      !                                                     ! ===========
211   END SUBROUTINE tra_bbl_dif
[3]212
[3764]213
[12377]214   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
[2528]215      !!----------------------------------------------------------------------
216      !!                  ***  ROUTINE trc_bbl  ***
217      !!
[3764]218      !! ** Purpose :   Compute the before passive tracer trend associated
[2528]219      !!     with the bottom boundary layer and add it to the general trend
220      !!     of tracer equations.
221      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
222      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
[3764]223      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
224      !!                       transport proportional to the along-slope density gradient
[2528]225      !!
226      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
227      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]228      !!----------------------------------------------------------------------
[2715]229      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
[12377]230      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
231      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
232      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
[2715]233      !
[2528]234      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
235      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
236      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
237      REAL(wp) ::   zbtr, ztra               ! local scalars
238      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
239      !!----------------------------------------------------------------------
240      !
241      !                                                          ! ===========
242      DO jn = 1, kjpt                                            ! tracer loop
[3764]243         !                                                       ! ===========
[457]244         DO jj = 1, jpjm1
[2528]245            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
246               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
247                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
248                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
249                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
250                  zu_bbl = ABS( utr_bbl(ji,jj) )
251                  !
252                  !                                               ! up  -slope T-point (shelf bottom point)
[12377]253                  zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
254                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
255                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
[3764]256                  !
[2528]257                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
[12377]258                     zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
259                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
260                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
[2528]261                  END DO
[3764]262                  !
[12377]263                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
264                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
265                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
[2528]266               ENDIF
267               !
268               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
269                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
270                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
271                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
272                  zv_bbl = ABS( vtr_bbl(ji,jj) )
[3764]273                  !
[2528]274                  ! up  -slope T-point (shelf bottom point)
[12377]275                  zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
276                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
277                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
[3764]278                  !
[2528]279                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
[12377]280                     zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
281                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
282                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
[2528]283                  END DO
284                  !                                               ! down-slope T-point (deep bottom point)
[12377]285                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
286                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
287                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
[2528]288               ENDIF
[457]289            END DO
[2528]290            !
[457]291         END DO
[9019]292         !                                                  ! ===========
293      END DO                                                ! end tracer
294      !                                                     ! ===========
[2528]295   END SUBROUTINE tra_bbl_adv
[3]296
297
[12377]298   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
[2528]299      !!----------------------------------------------------------------------
300      !!                  ***  ROUTINE bbl  ***
[3764]301      !!
[2528]302      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
[3764]303      !!                advection terms.
[2528]304      !!
[4990]305      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
[2528]306      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
307      !!      along bottom slope gradient) an additional lateral 2nd order
308      !!      diffusion along the bottom slope is added to the general
309      !!      tracer trend, otherwise the additional trend is set to 0.
310      !!      A typical value of ahbt is 2000 m2/s (equivalent to
311      !!      a downslope velocity of 20 cm/s if the condition for slope
312      !!      convection is satified)
[4990]313      !!              * advective bbl (nn_bbl_adv=1 or 2) :
[2528]314      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
315      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
316      !!        i.e. transport proportional to the along-slope density gradient
317      !!
318      !!      NB: the along slope density gradient is evaluated using the
319      !!      local density (i.e. referenced at a common local depth).
320      !!
321      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
322      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]323      !!----------------------------------------------------------------------
[2528]324      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
[4990]325      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
[2528]326      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
[12377]327      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
[6140]328      !
[2528]329      INTEGER  ::   ji, jj                    ! dummy loop indices
330      INTEGER  ::   ik                        ! local integers
[4990]331      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
332      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
333      REAL(wp) ::   za, zb, zgdrho            ! local scalars
334      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
335      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
336      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
[2528]337      !!----------------------------------------------------------------------
[3294]338      !
339      IF( kt == kit000 )  THEN
[2528]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
[4990]344      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
[12377]345      DO_2D_11_11
346         ik = mbkt(ji,jj)                             ! bottom T-level index
347         zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
348         zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
349         !
350         zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
351         zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
352         zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
353      END_2D
[4990]354      !
[12377]355      CALL eos_rab( zts, zdep, zab, Kmm )
[4990]356      !
[2528]357      !                                   !-------------------!
358      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
359         !                                !-------------------!
[12377]360         DO_2D_10_10
361            !                                                   ! i-direction
362            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
363            zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
364            !                                                         ! 2*masked bottom density gradient
365            zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
366               &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
367            !
[13247]368            zsign  = SIGN(  0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
[12377]369            ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
370            !
371            !                                                   ! j-direction
372            za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
373            zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
374            !                                                         ! 2*masked bottom density gradient
375            zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
376               &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
377            !
[13247]378            zsign = SIGN(  0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
[12377]379            ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
380         END_2D
[1601]381         !
[2528]382      ENDIF
[6140]383      !
[2528]384      !                                   !-------------------!
385      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
386         !                                !-------------------!
387         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
[503]388         !
[2528]389         CASE( 1 )                                   != use of upper velocity
[12377]390            DO_2D_10_10
391               !                                                  ! i-direction
392               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
393               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
394               !                                                          ! 2*masked bottom density gradient
395               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
396                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
397               !
[13247]398               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
399               zsigna= SIGN(  0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
[12377]400               !
401               !                                                          ! bbl velocity
402               utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
403               !
404               !                                                  ! j-direction
405               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
406               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
407               !                                                          ! 2*masked bottom density gradient
408               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
409                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
[13247]410               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
411               zsigna= SIGN(  0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
[12377]412               !
413               !                                                          ! bbl transport
414               vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
415            END_2D
[503]416            !
[2528]417         CASE( 2 )                                 != bbl velocity = F( delta rho )
418            zgbbl = grav * rn_gambbl
[12377]419            DO_2D_10_10
420               !                                                  ! i-direction
421               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
422               iid  = ji + MAX( 0, mgrhu(ji,jj) )
423               iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
424               !
425               ikud = mbku_d(ji,jj)
426               ikus = mbku(ji,jj)
427               !
428               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
429               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
430               !                                                          !   masked bottom density gradient
431               zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
432                  &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
433               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
434               !
435               !                                                          ! bbl transport (down-slope direction)
436               utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
437               !
438               !                                                  ! j-direction
439               !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
440               ijd  = jj + MAX( 0, mgrhv(ji,jj) )
441               ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
442               !
443               ikvd = mbkv_d(ji,jj)
444               ikvs = mbkv(ji,jj)
445               !
446               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
447               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
448               !                                                          !   masked bottom density gradient
449               zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
450                  &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
451               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
452               !
453               !                                                          ! bbl transport (down-slope direction)
454               vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
455            END_2D
[3]456         END SELECT
[2528]457         !
[3]458      ENDIF
[503]459      !
[2528]460   END SUBROUTINE bbl
[3]461
462
463   SUBROUTINE tra_bbl_init
464      !!----------------------------------------------------------------------
465      !!                  ***  ROUTINE tra_bbl_init  ***
466      !!
467      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
468      !!
469      !! ** Method  :   Read the nambbl namelist and check the parameters
[3294]470      !!              called by nemo_init at the first timestep (kit000)
[3]471      !!----------------------------------------------------------------------
[9019]472      INTEGER ::   ji, jj                      ! dummy loop indices
473      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
[9094]474      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
[9019]475      !!
476      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]477      !!----------------------------------------------------------------------
[3294]478      !
[4147]479      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
[11536]480901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
[6140]481      !
[4147]482      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
[11536]483902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
[4624]484      IF(lwm) WRITE ( numond, nambbl )
[2528]485      !
486      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
487      !
488      IF(lwp) THEN                   !* Parameter control and print
[3]489         WRITE(numout,*)
[2528]490         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]491         WRITE(numout,*) '~~~~~~~~~~~~'
[9019]492         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
493         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
[3]494      ENDIF
[9019]495      IF( .NOT.ln_trabbl )   RETURN
496      !
497      IF(lwp) THEN
498         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
499         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
500         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
501         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
502      ENDIF
503      !
[2715]504      !                              ! allocate trabbl arrays
505      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
[9019]506      !
[2528]507      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
508      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
[9019]509      !
[2528]510      !                             !* vertical index of  "deep" bottom u- and v-points
[12377]511      DO_2D_10_10
512         mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
513         mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
514      END_2D
[4990]515      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[9094]516      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
[13247]517      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp) 
[9919]518      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
[9019]519      !
[9168]520      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
[8509]521      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
[12377]522      DO_2D_10_10
523         IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
[13247]524            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]525         ENDIF
526         !
527         IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
[13247]528            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]529         ENDIF
530      END_2D
[7646]531      !
[12377]532      DO_2D_10_10
533         e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
534         e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
535      END_2D
[13247]536      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp )      ! lateral boundary conditions
[7646]537      !
[3764]538      !                             !* masked diffusive flux coefficients
[7753]539      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
540      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
[503]541      !
[3]542   END SUBROUTINE tra_bbl_init
543
544   !!======================================================================
545END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.