source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

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