New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trabbl.F90 in NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/trabbl.F90 @ 14986

Last change on this file since 14986 was 14986, checked in by sparonuz, 3 years ago

Merge trunk -r14984:HEAD

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