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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA/trabbl.F90 @ 14644

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

Merge trunk -r14642:HEAD

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