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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

  • Property svn:keywords set to Id
File size: 38.0 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
10   !!            3.3  ! 2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   ! 2010-04  (G. Madec)  Campin & Goosse advective bbl
12   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
[503]14   !!----------------------------------------------------------------------
[2528]15#if   defined key_trabbl   ||   defined key_esopa
[3]16   !!----------------------------------------------------------------------
[2528]17   !!   'key_trabbl'   or                             bottom boundary layer
[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
30   USE trdmod_oce     ! trends: ocean variables
31   USE trdtra         ! trends: active tracers
32   USE iom            ! IOM server               
33   USE in_out_manager ! I/O manager
34   USE lbclnk         ! ocean lateral boundary conditions
35   USE prtctl         ! Print control
[3294]36   USE wrk_nemo       ! Memory Allocation
37   USE timing         ! Timing
[3]38
[3294]39
[3]40   IMPLICIT NONE
41   PRIVATE
42
[2528]43   PUBLIC   tra_bbl       !  routine called by step.F90
44   PUBLIC   tra_bbl_init  !  routine called by opa.F90
45   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
46   PUBLIC   tra_bbl_adv   !  -          -          -              -
47   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
[3]48
[2528]49   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
[503]50
[2528]51   !                                           !!* Namelist nambbl *
52   INTEGER , PUBLIC ::   nn_bbl_ldf = 0         !: =1   : diffusive bbl or not (=0)
53   INTEGER , PUBLIC ::   nn_bbl_adv = 0         !: =1/2 : advective bbl or not (=0)
54   !                                            !  =1 : advective bbl using the bottom ocean velocity
55   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
56   REAL(wp), PUBLIC ::   rn_ahtbbl  = 1.e3_wp   !: along slope bbl diffusive coefficient [m2/s]
57   REAL(wp), PUBLIC ::   rn_gambbl  = 10.0_wp   !: lateral coeff. for bottom boundary layer scheme [s]
[409]58
[2715]59   LOGICAL , PUBLIC ::   l_bbl                  !: flag to compute bbl diffu. flux coef and transport
60   
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
62   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
[3]63
[2715]64   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level
65   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_e1e2t               ! inverse of the cell surface at t-point      [1/m2]
[3]69
70   !! * Substitutions
71#  include "domzgr_substitute.h90"
72#  include "vectopt_loop_substitute.h90"
73   !!----------------------------------------------------------------------
[2528]74   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
75   !! $Id$
76   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]77   !!----------------------------------------------------------------------
78CONTAINS
79
[2715]80   INTEGER FUNCTION tra_bbl_alloc()
81      !!----------------------------------------------------------------------
82      !!                ***  FUNCTION tra_bbl_alloc  ***
83      !!----------------------------------------------------------------------
84      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
85         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
86         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
87         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , r1_e1e2t(jpi,jpj)                  , STAT= tra_bbl_alloc )
88         !
89      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
90      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
91   END FUNCTION tra_bbl_alloc
92
93
[2528]94   SUBROUTINE tra_bbl( kt )
[3]95      !!----------------------------------------------------------------------
[2528]96      !!                  ***  ROUTINE bbl  ***
97      !!                   
[3]98      !! ** Purpose :   Compute the before tracer (t & s) trend associated
[2528]99      !!              with the bottom boundary layer and add it to the general
100      !!              trend of tracer equations.
[3]101      !!
[2528]102      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
103      !!              diffusive and/or advective contribution to the tracer trend
104      !!              is added to the general tracer trend
105      !!---------------------------------------------------------------------- 
106      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
107      !!
[3294]108      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
[2528]109      !!----------------------------------------------------------------------
[3294]110      !
111      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
112      !
113      IF( l_trdtra )   THEN                        !* Save ta and sa trends
114         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
115         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
116         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[2528]117      ENDIF
118
[3294]119      IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
120 
121      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
122         !
[2528]123         CALL tra_bbl_dif( tsb, tsa, jpts )
124         IF( ln_ctl )  &
125         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
126         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
127         ! lateral boundary conditions ; just need for outputs                         
128         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef     
130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
[3294]131         !
[2528]132      END IF
133
134      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
[3294]135         !
[2528]136         CALL tra_bbl_adv( tsb, tsa, jpts )
137         IF(ln_ctl)   &
138         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
139         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
140         ! lateral boundary conditions ; just need for outputs                         
141         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
142         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport     
143         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
[3294]144         !
[2528]145      END IF
146
147      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
148         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
149         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
[3294]152         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
[2528]153      ENDIF
154      !
[3294]155      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
156      !
[2528]157   END SUBROUTINE tra_bbl
158
159
160   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE tra_bbl_dif  ***
163      !!                   
164      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
165      !!                advection terms.
166      !!
167      !! ** Method  :   
168      !!        * diffusive bbl (nn_bbl_ldf=1) :
169      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
170      !!      along bottom slope gradient) an additional lateral 2nd order
171      !!      diffusion along the bottom slope is added to the general
172      !!      tracer trend, otherwise the additional trend is set to 0.
173      !!      A typical value of ahbt is 2000 m2/s (equivalent to
[3]174      !!      a downslope velocity of 20 cm/s if the condition for slope
175      !!      convection is satified)
176      !!
[2528]177      !! ** Action  :   pta   increased by the bbl diffusive trend
[3]178      !!
[503]179      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
[2528]180      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
181      !!---------------------------------------------------------------------- 
[2715]182      !
183      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
184      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
185      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
186      !
[2528]187      INTEGER  ::   ji, jj, jn   ! dummy loop indices
188      INTEGER  ::   ik           ! local integers
189      REAL(wp) ::   zbtr         ! local scalars
[3294]190      REAL(wp), POINTER, DIMENSION(:,:) :: zptb
[3]191      !!----------------------------------------------------------------------
[2528]192      !
[3294]193      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
[2715]194      !
[3294]195      CALL wrk_alloc( jpi, jpj, zptb )
196      !
[2528]197      DO jn = 1, kjpt                                     ! tracer loop
198         !                                                ! ===========
[789]199#  if defined key_vectopt_loop
[2528]200         DO jj = 1, 1   ! vector opt. (forced unrolling)
201            DO ji = 1, jpij
202#else
203         DO jj = 1, jpj
204            DO ji = 1, jpi 
205#endif
206               ik = mbkt(ji,jj)                        ! bottom T-level index
207               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S
208            END DO
209         END DO
210         !                                                ! Compute the trend
211#  if defined key_vectopt_loop
212         DO jj = 1, 1   ! vector opt. (forced unrolling)
213            DO ji = jpi+1, jpij-jpi-1
[3]214#  else
[2528]215         DO jj = 2, jpjm1
216            DO ji = 2, jpim1
[3]217#  endif
[2528]218               ik = mbkt(ji,jj)                            ! bottom T-level index
[2715]219               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik)
[2528]220               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
221                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
222                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
223                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
224                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
225            END DO
[3]226         END DO
[2528]227         !                                                  ! ===========
228      END DO                                                ! end tracer
229      !                                                     ! ===========
[3294]230      CALL wrk_dealloc( jpi, jpj, zptb )
[2715]231      !
[3294]232      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
233      !
[2528]234   END SUBROUTINE tra_bbl_dif
235   
[3]236
[2528]237   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE trc_bbl  ***
240      !!
241      !! ** Purpose :   Compute the before passive tracer trend associated
242      !!     with the bottom boundary layer and add it to the general trend
243      !!     of tracer equations.
244      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
245      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
246      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
247      !!                       transport proportional to the along-slope density gradient                   
248      !!
249      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
250      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
251      !!---------------------------------------------------------------------- 
[2715]252      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
253      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
254      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
255      !
[2528]256      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
257      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
258      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
259      REAL(wp) ::   zbtr, ztra               ! local scalars
260      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
261      !!----------------------------------------------------------------------
262      !
[3294]263      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
[2528]264      !                                                          ! ===========
265      DO jn = 1, kjpt                                            ! tracer loop
266         !                                                       ! ===========         
[789]267# if defined key_vectopt_loop
[1601]268         DO jj = 1, 1
[2528]269            DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling)
[457]270# else
271         DO jj = 1, jpjm1
[2528]272            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
[457]273# endif
[2528]274               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
275                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
276                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
277                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
278                  zu_bbl = ABS( utr_bbl(ji,jj) )
279                  !
280                  !                                               ! up  -slope T-point (shelf bottom point)
[2715]281                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
[2528]282                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
283                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
284                  !                   
285                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
[2715]286                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
[2528]287                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
288                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
289                  END DO
290                  !
[2715]291                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
[2528]292                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
293                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
294               ENDIF
295               !
296               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
297                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
298                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
299                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
300                  zv_bbl = ABS( vtr_bbl(ji,jj) )
301                  !
302                  ! up  -slope T-point (shelf bottom point)
[2715]303                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
[2528]304                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
305                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
306                  !                   
307                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
[2715]308                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
[2528]309                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
310                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
311                  END DO
312                  !                                               ! down-slope T-point (deep bottom point)
[2715]313                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
[2528]314                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
315                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
316               ENDIF
[457]317            END DO
[2528]318            !
[457]319         END DO
[2528]320         !                                                  ! ===========
321      END DO                                                ! end tracer
322      !                                                     ! ===========
[3294]323      !
324      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
325      !
[2528]326   END SUBROUTINE tra_bbl_adv
[3]327
328
[3294]329   SUBROUTINE bbl( kt, kit000, cdtype )
[2528]330      !!----------------------------------------------------------------------
331      !!                  ***  ROUTINE bbl  ***
332      !!                   
333      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
334      !!                advection terms.
335      !!
336      !! ** Method  :   
337      !!        * diffusive bbl (nn_bbl_ldf=1) :
338      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
339      !!      along bottom slope gradient) an additional lateral 2nd order
340      !!      diffusion along the bottom slope is added to the general
341      !!      tracer trend, otherwise the additional trend is set to 0.
342      !!      A typical value of ahbt is 2000 m2/s (equivalent to
343      !!      a downslope velocity of 20 cm/s if the condition for slope
344      !!      convection is satified)
345      !!        * advective bbl (nn_bbl_adv=1 or 2) :
346      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
347      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
348      !!        i.e. transport proportional to the along-slope density gradient
349      !!
350      !!      NB: the along slope density gradient is evaluated using the
351      !!      local density (i.e. referenced at a common local depth).
352      !!
353      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
354      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
355      !!---------------------------------------------------------------------- 
[2715]356      !
[2528]357      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
[3294]358      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
[2528]359      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
360      !!
361      INTEGER  ::   ji, jj                    ! dummy loop indices
362      INTEGER  ::   ik                        ! local integers
363      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
364      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
365      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
366      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
367      !!
368      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
[3294]369      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep
[2528]370      !!----------------------- zv_bbl-----------------------------------------------
371      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
372      ! ================            pft :  potential temperature in degrees celcius
373      !                             pfs :  salinity anomaly (s-35) in psu
374      !                             pfh :  depth in meters
375      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
376      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
377         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
378                                   - 0.203814e-03 ) * pft                    &
379                                   + 0.170907e-01 ) * pft                    &
380                                   + 0.665157e-01                            &
381         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
382         +  ( ( - 0.302285e-13 * pfh                                         &
383                - 0.251520e-11 * pfs                                         &
384                + 0.512857e-12 * pft * pft          ) * pfh                  &
385                                     - 0.164759e-06   * pfs                  &
386             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
387                                     + 0.380374e-04 ) * pfh
388      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
389         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
390                                 - 0.301985e-05 ) * pft                      &
391                                 + 0.785567e-03                              &
392         + (     0.515032e-08 * pfs                                          &
393               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
394               +(  (   0.121551e-17 * pfh                                    &
395                     - 0.602281e-15 * pfs                                    &
396                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
397                                          + 0.408195e-10   * pfs             &
398                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
399                                          - 0.121555e-07 ) * pfh
400      !!----------------------------------------------------------------------
[3294]401     
402      !
403      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
404      !
405      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
406      !
[2715]407     
[3294]408      IF( kt == kit000 )  THEN
[2528]409         IF(lwp)  WRITE(numout,*)
410         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
411         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
412      ENDIF
413     
414      !                                        !* bottom temperature, salinity, velocity and depth
415#if defined key_vectopt_loop
416      DO jj = 1, 1   ! vector opt. (forced unrolling)
417         DO ji = 1, jpij
418#else
419      DO jj = 1, jpj
420         DO ji = 1, jpi
421#endif
422            ik = mbkt(ji,jj)                        ! bottom T-level index
423            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
424            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
425            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
426            !
427            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
428            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
429         END DO
430      END DO
431     
432      !                                   !-------------------!
433      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
434         !                                !-------------------!
435         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
436            DO ji = 1, jpim1             
437               !                                                ! i-direction
438               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
439               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
[1601]440               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
[2528]441               !                                                         ! masked bbl i-gradient of density
442               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
443                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
444               !                                                     
445               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
446               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
[1601]447               !
[2528]448               !                                                ! j-direction
449               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
450               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
[1601]451               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
[2528]452               !                                                         ! masked bbl j-gradient of density
453               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
454                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
455               !                                                   
456               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
457               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
[1601]458               !
459            END DO
[3]460         END DO
[1601]461         !
[2528]462      ENDIF
[409]463
[2528]464      !                                   !-------------------!
465      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
466         !                                !-------------------!
467         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
[503]468         !
[2528]469         CASE( 1 )                                   != use of upper velocity
470            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
471               DO ji = 1, fs_jpim1   ! vector opt.
472                  !                                               ! i-direction
473                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
474                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
475                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
476                  !                                                           ! masked bbl i-gradient of density
477                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
478                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
479                  !                                                         
480                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
481                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
482                  !
483                  !                                                           ! bbl velocity
484                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
485                  !
486                  !                                               ! j-direction
487                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
488                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
489                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
490                  !                                                           ! masked bbl j-gradient of density
491                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
492                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
493                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
494                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
495                  !
496                  !                                                           ! bbl velocity
497                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
498               END DO
499            END DO
[503]500            !
[2528]501         CASE( 2 )                                 != bbl velocity = F( delta rho )
502            zgbbl = grav * rn_gambbl
503            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
504               DO ji = 1, fs_jpim1   ! vector opt.
505                  !                                         ! i-direction
506                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
507                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
508                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
509                  !
510                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
511                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
512                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
513                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
514                  zgdrho =    fsbeta( zt, zs, zh )                                    &
515                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
516                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
517                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
518                  !
519                  !                                             ! bbl transport (down-slope direction)
520                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
521                  !
522                  !                                         ! j-direction
523                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
524                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
525                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
526                  !
527                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
528                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
529                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
530                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
531                  zgdrho =    fsbeta( zt, zs, zh )                                    &
532                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
533                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
534                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
535                  !
536                  !                                             ! bbl transport (down-slope direction)
537                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
538               END DO
539            END DO
[3]540         END SELECT
[2528]541         !
[3]542      ENDIF
[503]543      !
[3294]544      CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
[2715]545      !
[3294]546      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
547      !
[2528]548   END SUBROUTINE bbl
[3]549
550
551   SUBROUTINE tra_bbl_init
552      !!----------------------------------------------------------------------
553      !!                  ***  ROUTINE tra_bbl_init  ***
554      !!
555      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
556      !!
557      !! ** Method  :   Read the nambbl namelist and check the parameters
[3294]558      !!              called by nemo_init at the first timestep (kit000)
[3]559      !!----------------------------------------------------------------------
[2528]560      INTEGER ::   ji, jj               ! dummy loop indices
561      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
[3294]562      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
[2528]563      !!
564      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]565      !!----------------------------------------------------------------------
[3294]566      !
567      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
568      !
569      CALL wrk_alloc( jpi, jpj, zmbk ) 
570      !
[3]571
[2528]572      REWIND ( numnam )              !* Read Namelist nambbl : bottom boundary layer scheme
[3]573      READ   ( numnam, nambbl )
[2528]574      !
575      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
576      !
577      IF(lwp) THEN                   !* Parameter control and print
[3]578         WRITE(numout,*)
[2528]579         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]580         WRITE(numout,*) '~~~~~~~~~~~~'
[503]581         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
[2528]582         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
583         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
584         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
585         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
[3]586      ENDIF
[2715]587
588      !                              ! allocate trabbl arrays
589      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
590     
[2528]591      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
592      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
593
594      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
595
596
597      !                             !* inverse of surface of T-cells
[2715]598      r1_e1e2t(:,:) = 1._wp / ( e1t(:,:) * e2t(:,:) )
[2528]599     
600      !                             !* vertical index of  "deep" bottom u- and v-points
601      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
602         DO ji = 1, jpim1
603            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
604            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
[3]605         END DO
606      END DO
[2528]607      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
608      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
609      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
610
611                                        !* sign of grad(H) at u- and v-points
612      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
613      DO jj = 1, jpjm1               
[3]614         DO ji = 1, jpim1
[2528]615            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
616            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
[3]617         END DO
618      END DO
619
[2528]620      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
621         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
622            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 
623            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 
624         END DO
625      END DO
626      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
[481]627
[2528]628      !                             !* masked diffusive flux coefficients
629      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
630      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
[481]631
[2528]632
633      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
634         !
635         SELECT CASE ( jp_cfg )
636         CASE ( 2 )                          ! ORCA_R2
637            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
638            ii0 = 139   ;   ii1 = 140 
639            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
640            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
641            !
642            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
643            ii0 = 161   ;   ii1 = 162
644            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
645            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
646            !
647         CASE ( 4 )                          ! ORCA_R4
648            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
649            ii0 =  70   ;   ii1 =  71 
650            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
651            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
652         END SELECT
653         !
654      ENDIF
[503]655      !
[3294]656      CALL wrk_dealloc( jpi, jpj, zmbk ) 
[2715]657      !
[3294]658      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
659      !
[3]660   END SUBROUTINE tra_bbl_init
661
662#else
663   !!----------------------------------------------------------------------
664   !!   Dummy module :                      No bottom boundary layer scheme
665   !!----------------------------------------------------------------------
[2528]666   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
[3]667CONTAINS
[2528]668   SUBROUTINE tra_bbl_init               ! Dummy routine
669   END SUBROUTINE tra_bbl_init
670   SUBROUTINE tra_bbl( kt )              ! Dummy routine
671      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
672   END SUBROUTINE tra_bbl
[3]673#endif
674
675   !!======================================================================
676END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.