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 @ 10954

Last change on this file since 10954 was 10954, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TRA modules and all knock on effects of these conversions. SETTE tested

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