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

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90 @ 13551

Last change on this file since 13551 was 13551, checked in by hadcv, 4 years ago

#2365: Replace trd_tra workarounds with ctl_warn if using tiling

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