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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90 @ 10806

Last change on this file since 10806 was 10806, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.

  • Property svn:keywords set to Id
File size: 32.9 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
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
[9598]71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]72   !! $Id$
[10068]73   !! Software governed by the CeCILL license (see ./LICENSE)
[3]74   !!----------------------------------------------------------------------
75CONTAINS
76
[2715]77   INTEGER FUNCTION tra_bbl_alloc()
78      !!----------------------------------------------------------------------
79      !!                ***  FUNCTION tra_bbl_alloc  ***
80      !!----------------------------------------------------------------------
[9019]81      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
82         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
83         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
84         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
[2715]85         !
[10425]86      CALL mpp_sum ( 'trabbl', tra_bbl_alloc )
[2715]87      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
88   END FUNCTION tra_bbl_alloc
89
90
[10806]91   SUBROUTINE tra_bbl( kt, ktlev1, ktlev2, kt2lev, pts_rhs )
[3]92      !!----------------------------------------------------------------------
[2528]93      !!                  ***  ROUTINE bbl  ***
[3764]94      !!
95      !! ** Purpose :   Compute the before tracer (t & s) trend associated
[2528]96      !!              with the bottom boundary layer and add it to the general
97      !!              trend of tracer equations.
[3]98      !!
[2528]99      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
100      !!              diffusive and/or advective contribution to the tracer trend
101      !!              is added to the general tracer trend
[3764]102      !!----------------------------------------------------------------------
[10806]103      INTEGER, INTENT( in ) ::   kt              ! ocean time-step
104      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms
105      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms
106      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
[4990]107      !
[9019]108      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
[2528]109      !!----------------------------------------------------------------------
[3294]110      !
[9019]111      IF( ln_timing )   CALL timing_start( 'tra_bbl')
[3294]112      !
[9019]113      IF( l_trdtra )   THEN                         !* Save the T-S input trends
114         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
[10806]115         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)
116         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)
[2528]117      ENDIF
118
[10806]119      IF( l_bbl )   CALL bbl( kt, nit000, ktlev1, ktlev2, kt2lev, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
[3764]120
[4990]121      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
[3294]122         !
[10806]123         CALL tra_bbl_dif( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )
[2528]124         IF( ln_ctl )  &
125         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
[4990]126            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[3764]127         ! lateral boundary conditions ; just need for outputs
[10425]128         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
[3764]129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
[2528]130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
[3294]131         !
[9168]132      ENDIF
[6140]133      !
[4990]134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
[3294]135         !
[10806]136         CALL tra_bbl_adv( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )
[2528]137         IF(ln_ctl)   &
138         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
[4990]139            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
[3764]140         ! lateral boundary conditions ; just need for outputs
[10425]141         CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
[3764]142         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
[2528]143         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
[3294]144         !
[9168]145      ENDIF
[2528]146
[6140]147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
[10806]148         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)
149         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)
[4990]150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
[9019]152         DEALLOCATE( ztrdt, ztrds )
[2528]153      ENDIF
154      !
[9019]155      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
[3294]156      !
[2528]157   END SUBROUTINE tra_bbl
158
159
[10806]160   SUBROUTINE tra_bbl_dif( pt, pe3t, pt_rhs, kjpt )
[2528]161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE tra_bbl_dif  ***
[3764]163      !!
[2528]164      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
[3764]165      !!                advection terms.
[2528]166      !!
[4990]167      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
[2528]168      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
169      !!      along bottom slope gradient) an additional lateral 2nd order
170      !!      diffusion along the bottom slope is added to the general
171      !!      tracer trend, otherwise the additional trend is set to 0.
172      !!      A typical value of ahbt is 2000 m2/s (equivalent to
[3]173      !!      a downslope velocity of 20 cm/s if the condition for slope
174      !!      convection is satified)
175      !!
[10806]176      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
[3]177      !!
[503]178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
[2528]179      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]180      !!----------------------------------------------------------------------
[2715]181      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
[10806]182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! tracer fields
183      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields
184      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
[2715]185      !
[2528]186      INTEGER  ::   ji, jj, jn   ! dummy loop indices
187      INTEGER  ::   ik           ! local integers
188      REAL(wp) ::   zbtr         ! local scalars
[9019]189      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
[3]190      !!----------------------------------------------------------------------
[2528]191      !
192      DO jn = 1, kjpt                                     ! tracer loop
193         !                                                ! ===========
194         DO jj = 1, jpj
[3764]195            DO ji = 1, jpi
[5836]196               ik = mbkt(ji,jj)                             ! bottom T-level index
[10806]197               zptb(ji,jj) = pt(ji,jj,ik,jn)               ! bottom before T and S
[2528]198            END DO
199         END DO
[4990]200         !               
201         DO jj = 2, jpjm1                                    ! Compute the trend
[2528]202            DO ji = 2, jpim1
[5836]203               ik = mbkt(ji,jj)                            ! bottom T-level index
[10806]204               pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
[5836]205                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
206                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
207                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
208                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
[10806]209                  &             * r1_e1e2t(ji,jj) / pe3t(ji,jj,ik)
[2528]210            END DO
[3]211         END DO
[2528]212         !                                                  ! ===========
213      END DO                                                ! end tracer
214      !                                                     ! ===========
215   END SUBROUTINE tra_bbl_dif
[3]216
[3764]217
[10806]218   SUBROUTINE tra_bbl_adv( pt, pe3t, pt_rhs, kjpt )
[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
[10806]234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt    ! before and now tracer fields
235      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields
236      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs    ! tracer trend
[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   !   -       -
241      REAL(wp) ::   zbtr, ztra               ! local scalars
242      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
243      !!----------------------------------------------------------------------
244      !
245      !                                                          ! ===========
246      DO jn = 1, kjpt                                            ! tracer loop
[3764]247         !                                                       ! ===========
[457]248         DO jj = 1, jpjm1
[2528]249            DO ji = 1, jpim1            ! 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) )
255                  !
256                  !                                               ! up  -slope T-point (shelf bottom point)
[10806]257                  zbtr = r1_e1e2t(iis,jj) / pe3t(iis,jj,ikus)
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
[3764]260                  !
[2528]261                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
[10806]262                     zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,jk)
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
[2528]265                  END DO
[3764]266                  !
[10806]267                  zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,ikud)
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
[2528]270               ENDIF
271               !
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) )
[3764]277                  !
[2528]278                  ! up  -slope T-point (shelf bottom point)
[10806]279                  zbtr = r1_e1e2t(ji,ijs) / pe3t(ji,ijs,ikvs)
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
[3764]282                  !
[2528]283                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
[10806]284                     zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,jk)
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
[2528]287                  END DO
288                  !                                               ! down-slope T-point (deep bottom point)
[10806]289                  zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,ikvd)
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
[2528]292               ENDIF
[457]293            END DO
[2528]294            !
[457]295         END DO
[9019]296         !                                                  ! ===========
297      END DO                                                ! end tracer
298      !                                                     ! ===========
[2528]299   END SUBROUTINE tra_bbl_adv
[3]300
301
[10806]302   SUBROUTINE bbl( kt, kit000, ktlev1, ktlev2, kt2lev, cdtype )
[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
[10806]330      INTEGER         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level indices for 3-time-levelsource terms
331      INTEGER         , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms
[2528]332      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
[6140]333      !
[2528]334      INTEGER  ::   ji, jj                    ! dummy loop indices
335      INTEGER  ::   ik                        ! local integers
[4990]336      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
337      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
338      REAL(wp) ::   za, zb, zgdrho            ! local scalars
339      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
340      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
341      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
[2528]342      !!----------------------------------------------------------------------
[3294]343      !
344      IF( kt == kit000 )  THEN
[2528]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
[4990]349      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
[2528]350      DO jj = 1, jpj
351         DO ji = 1, jpi
[4990]352            ik = mbkt(ji,jj)                             ! bottom T-level index
[10806]353            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,ktlev1)    ! bottom before T and S
354            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,ktlev1)
[2528]355            !
[10806]356            zdep(ji,jj) = gdept(ji,jj,ik,kt2lev)              ! bottom T-level reference depth
357            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),ktlev2)          ! bottom velocity
358            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),ktlev2)
[2528]359         END DO
360      END DO
[4990]361      !
362      CALL eos_rab( zts, zdep, zab )
363      !
[2528]364      !                                   !-------------------!
365      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
366         !                                !-------------------!
367         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
[4990]368            DO ji = 1, fs_jpim1   ! vector opt.
369               !                                                   ! i-direction
370               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
371               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
372               !                                                         ! 2*masked bottom density gradient
373               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
374                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
[3764]375               !
[4990]376               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
377               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
[1601]378               !
[4990]379               !                                                   ! j-direction
380               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
381               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
382               !                                                         ! 2*masked bottom density gradient
383               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
384                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
[3764]385               !
[4990]386               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
[2528]387               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
[1601]388            END DO
[3]389         END DO
[1601]390         !
[2528]391      ENDIF
[6140]392      !
[2528]393      !                                   !-------------------!
394      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
395         !                                !-------------------!
396         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
[503]397         !
[2528]398         CASE( 1 )                                   != use of upper velocity
399            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
400               DO ji = 1, fs_jpim1   ! vector opt.
[4990]401                  !                                                  ! i-direction
402                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
403                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
404                  !                                                          ! 2*masked bottom density gradient
405                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
406                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
[3764]407                  !
[4990]408                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
409                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
[2528]410                  !
[4990]411                  !                                                          ! bbl velocity
[2528]412                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
413                  !
[4990]414                  !                                                  ! j-direction
415                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
416                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
417                  !                                                          ! 2*masked bottom density gradient
418                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
419                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
420                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
421                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
[2528]422                  !
[4990]423                  !                                                          ! bbl transport
[2528]424                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
425               END DO
426            END DO
[503]427            !
[2528]428         CASE( 2 )                                 != bbl velocity = F( delta rho )
429            zgbbl = grav * rn_gambbl
430            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
431               DO ji = 1, fs_jpim1   ! vector opt.
[4990]432                  !                                                  ! i-direction
[2528]433                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
[4990]434                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
435                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
[2528]436                  !
[4990]437                  ikud = mbku_d(ji,jj)
438                  ikus = mbku(ji,jj)
[2528]439                  !
[4990]440                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
441                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
442                  !                                                          !   masked bottom density gradient
443                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
444                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
445                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
446                  !
447                  !                                                          ! bbl transport (down-slope direction)
[2528]448                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
449                  !
[4990]450                  !                                                  ! j-direction
[2528]451                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
[4990]452                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
453                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
[2528]454                  !
[4990]455                  ikvd = mbkv_d(ji,jj)
456                  ikvs = mbkv(ji,jj)
[2528]457                  !
[4990]458                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
459                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
460                  !                                                          !   masked bottom density gradient
461                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
462                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
463                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
464                  !
465                  !                                                          ! bbl transport (down-slope direction)
[2528]466                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
467               END DO
468            END DO
[3]469         END SELECT
[2528]470         !
[3]471      ENDIF
[503]472      !
[2528]473   END SUBROUTINE bbl
[3]474
475
476   SUBROUTINE tra_bbl_init
477      !!----------------------------------------------------------------------
478      !!                  ***  ROUTINE tra_bbl_init  ***
479      !!
480      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
481      !!
482      !! ** Method  :   Read the nambbl namelist and check the parameters
[3294]483      !!              called by nemo_init at the first timestep (kit000)
[3]484      !!----------------------------------------------------------------------
[9019]485      INTEGER ::   ji, jj                      ! dummy loop indices
486      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
[9094]487      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
[9019]488      !!
489      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]490      !!----------------------------------------------------------------------
[3294]491      !
[4147]492      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
493      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
[6140]494901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
495      !
[4147]496      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
497      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
[9168]498902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
[4624]499      IF(lwm) WRITE ( numond, nambbl )
[2528]500      !
501      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
502      !
503      IF(lwp) THEN                   !* Parameter control and print
[3]504         WRITE(numout,*)
[2528]505         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]506         WRITE(numout,*) '~~~~~~~~~~~~'
[9019]507         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
508         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
[3]509      ENDIF
[9019]510      IF( .NOT.ln_trabbl )   RETURN
511      !
512      IF(lwp) THEN
513         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
514         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
515         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
516         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
517      ENDIF
518      !
[2715]519      !                              ! allocate trabbl arrays
520      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
[9019]521      !
[2528]522      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
523      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
[9019]524      !
[2528]525      !                             !* vertical index of  "deep" bottom u- and v-points
526      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
527         DO ji = 1, jpim1
528            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
529            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
[3]530         END DO
531      END DO
[4990]532      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[9094]533      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
[10425]534      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
[9919]535      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
[9019]536      !
[9168]537      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
[8509]538      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
[3764]539      DO jj = 1, jpjm1
[3]540         DO ji = 1, jpim1
[8509]541            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
542               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
543            ENDIF
544            !
545            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
546               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
547            ENDIF
[3]548         END DO
549      END DO
[7646]550      !
[3764]551      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
[4990]552         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
[4292]553            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
554            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
[3764]555         END DO
[2528]556      END DO
[10425]557      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1. , e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
[7646]558      !
[3764]559      !                             !* masked diffusive flux coefficients
[7753]560      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
561      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
[503]562      !
[3]563   END SUBROUTINE tra_bbl_init
564
565   !!======================================================================
566END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.