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_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/TRA/trabbl.F90 @ 14002

Last change on this file since 14002 was 14002, checked in by emanuelaclementi, 4 years ago

phasing with trunk rev14001 - ticket #2152 #2339

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