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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8214

Last change on this file since 8214 was 7954, checked in by gm, 7 years ago

#1880 (HPC-09): OPA & TOP replace key_trabbl by ln_trabbl + in CONFIG, remove reference to key_zdfXXX and key_trabbl

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