source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 4726

Last change on this file since 4726 was 4726, checked in by mathiot, 6 years ago

ISF branch: change name of 2 variables (icedep ⇒ risfdep and lmask ⇒ ssmask), cosmetic changes and add ldfslp key

  • Property svn:keywords set to Id
File size: 37.7 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 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
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  !: =1   : diffusive bbl or not (=0)
53   INTEGER , PUBLIC ::   nn_bbl_adv  !: =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   !: along slope bbl diffusive coefficient [m2/s]
57   REAL(wp), PUBLIC ::   rn_gambbl   !: 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
69   !! * Substitutions
70#  include "domzgr_substitute.h90"
71#  include "vectopt_loop_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
74   !! $Id$
75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   INTEGER FUNCTION tra_bbl_alloc()
80      !!----------------------------------------------------------------------
81      !!                ***  FUNCTION tra_bbl_alloc  ***
82      !!----------------------------------------------------------------------
83      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
84         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
85         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
86         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc                )
87         !
88      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
89      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
90   END FUNCTION tra_bbl_alloc
91
92
93   SUBROUTINE tra_bbl( kt )
94      !!----------------------------------------------------------------------
95      !!                  ***  ROUTINE bbl  ***
96      !!
97      !! ** Purpose :   Compute the before tracer (t & s) trend associated
98      !!              with the bottom boundary layer and add it to the general
99      !!              trend of tracer equations.
100      !!
101      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
102      !!              diffusive and/or advective contribution to the tracer trend
103      !!              is added to the general tracer trend
104      !!----------------------------------------------------------------------
105      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
106      !!
107      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
108      !!----------------------------------------------------------------------
109      !
110      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
111      !
112      IF( l_trdtra )   THEN                        !* Save ta and sa trends
113         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
114         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
115         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
116      ENDIF
117
118      IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
119
120      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
121         !
122         CALL tra_bbl_dif( tsb, tsa, jpts )
123         IF( ln_ctl )  &
124         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
125         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
126         ! lateral boundary conditions ; just need for outputs
127         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
128         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
129         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
130         !
131      END IF
132
133      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
134         !
135         CALL tra_bbl_adv( tsb, tsa, jpts )
136         IF(ln_ctl)   &
137         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
138         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
139         ! lateral boundary conditions ; just need for outputs
140         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
141         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
142         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
143         !
144      END IF
145
146      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
147         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
148         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
149         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
150         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
151         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
152      ENDIF
153      !
154      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
155      !
156   END SUBROUTINE tra_bbl
157
158
159   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
160      !!----------------------------------------------------------------------
161      !!                  ***  ROUTINE tra_bbl_dif  ***
162      !!
163      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
164      !!                advection terms.
165      !!
166      !! ** Method  :
167      !!        * diffusive bbl (nn_bbl_ldf=1) :
168      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
169      !!      along bottom slope gradient) an additional lateral 2nd order
170      !!      diffusion along the bottom slope is added to the general
171      !!      tracer trend, otherwise the additional trend is set to 0.
172      !!      A typical value of ahbt is 2000 m2/s (equivalent to
173      !!      a downslope velocity of 20 cm/s if the condition for slope
174      !!      convection is satified)
175      !!
176      !! ** Action  :   pta   increased by the bbl diffusive trend
177      !!
178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
179      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
180      !!----------------------------------------------------------------------
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_e12t(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_e12t(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_e12t(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_e12t(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_e12t(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_e12t(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_e12t(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      !
402      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
403      !
404      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
405      !
406
407      IF( kt == kit000 )  THEN
408         IF(lwp)  WRITE(numout,*)
409         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
410         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
411      ENDIF
412
413      !                                        !* bottom temperature, salinity, velocity and depth
414#if defined key_vectopt_loop
415      DO jj = 1, 1   ! vector opt. (forced unrolling)
416         DO ji = 1, jpij
417#else
418      DO jj = 1, jpj
419         DO ji = 1, jpi
420#endif
421            ik = mbkt(ji,jj)                        ! bottom T-level index
422            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * ssmask(ji,jj)      ! bottom before T and S
423            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * ssmask(ji,jj)
424            zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth
425            !
426            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
427            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
428         END DO
429      END DO
430
431      !                                   !-------------------!
432      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
433         !                                !-------------------!
434         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
435            DO ji = 1, jpim1
436               !                                                ! i-direction
437               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
438               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
439               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
440               !                                                         ! masked bbl i-gradient of density
441               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
442                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
443               !
444               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
445               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
446               !
447               !                                                ! j-direction
448               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
449               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
450               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
451               !                                                         ! masked bbl j-gradient of density
452               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
453                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
454               !
455               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
456               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
457               !
458            END DO
459         END DO
460         !
461      ENDIF
462
463      !                                   !-------------------!
464      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
465         !                                !-------------------!
466         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
467         !
468         CASE( 1 )                                   != use of upper velocity
469            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
470               DO ji = 1, fs_jpim1   ! vector opt.
471                  !                                               ! i-direction
472                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
473                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
474                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
475                  !                                                           ! masked bbl i-gradient of density
476                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
477                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
478                  !
479                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
480                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
481                  !
482                  !                                                           ! bbl velocity
483                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
484                  !
485                  !                                               ! j-direction
486                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
487                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
488                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
489                  !                                                           ! masked bbl j-gradient of density
490                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
491                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
492                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
493                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
494                  !
495                  !                                                           ! bbl velocity
496                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
497               END DO
498            END DO
499            !
500         CASE( 2 )                                 != bbl velocity = F( delta rho )
501            zgbbl = grav * rn_gambbl
502            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
503               DO ji = 1, fs_jpim1   ! vector opt.
504                  !                                         ! i-direction
505                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
506                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
507                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
508                  !
509                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
510                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
511                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
512                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
513                  zgdrho =    fsbeta( zt, zs, zh )                                    &
514                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
515                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
516                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
517                  !
518                  !                                             ! bbl transport (down-slope direction)
519                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
520                  !
521                  !                                         ! j-direction
522                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
523                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
524                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
525                  !
526                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
527                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
528                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
529                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
530                  zgdrho =    fsbeta( zt, zs, zh )                                    &
531                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
532                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
533                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
534                  !
535                  !                                             ! bbl transport (down-slope direction)
536                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
537               END DO
538            END DO
539         END SELECT
540         !
541      ENDIF
542      !
543      CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
544      !
545      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
546      !
547   END SUBROUTINE bbl
548
549
550   SUBROUTINE tra_bbl_init
551      !!----------------------------------------------------------------------
552      !!                  ***  ROUTINE tra_bbl_init  ***
553      !!
554      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
555      !!
556      !! ** Method  :   Read the nambbl namelist and check the parameters
557      !!              called by nemo_init at the first timestep (kit000)
558      !!----------------------------------------------------------------------
559      INTEGER ::   ji, jj               ! dummy loop indices
560      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
561      INTEGER  ::   ios                 ! Local integer output status for namelist read
562      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
563      !!
564      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
565      !!----------------------------------------------------------------------
566      !
567      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
568      !
569      CALL wrk_alloc( jpi, jpj, zmbk )
570      !
571
572      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
573      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
574901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
575
576      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
577      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
578902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
579      IF(lwm) WRITE ( numond, nambbl )
580      !
581      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
582      !
583      IF(lwp) THEN                   !* Parameter control and print
584         WRITE(numout,*)
585         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
586         WRITE(numout,*) '~~~~~~~~~~~~'
587         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
588         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
589         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
590         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
591         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
592      ENDIF
593
594      !                              ! allocate trabbl arrays
595      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
596
597      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
598      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
599
600      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
601
602      !                             !* vertical index of  "deep" bottom u- and v-points
603      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
604         DO ji = 1, jpim1
605            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
606            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
607         END DO
608      END DO
609      ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
610      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
611      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
612
613                                     !* sign of grad(H) at u- and v-points
614      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
615      DO jj = 1, jpjm1
616         DO ji = 1, jpim1
617            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
618            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
619         END DO
620      END DO
621
622      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
623         DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0)
624            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
625            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
626         END DO
627      END DO
628      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
629
630      !                             !* masked diffusive flux coefficients
631      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
632      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
633
634
635      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
636         !
637         SELECT CASE ( jp_cfg )
638         CASE ( 2 )                          ! ORCA_R2
639            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
640            ii0 = 139   ;   ii1 = 140
641            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
642            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
643            !
644            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
645            ii0 = 161   ;   ii1 = 162
646            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
647            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
648            !
649         CASE ( 4 )                          ! ORCA_R4
650            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
651            ii0 =  70   ;   ii1 =  71
652            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
653            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
654         END SELECT
655         !
656      ENDIF
657      !
658      CALL wrk_dealloc( jpi, jpj, zmbk )
659      !
660      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
661      !
662   END SUBROUTINE tra_bbl_init
663
664#else
665   !!----------------------------------------------------------------------
666   !!   Dummy module :                      No bottom boundary layer scheme
667   !!----------------------------------------------------------------------
668   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
669CONTAINS
670   SUBROUTINE tra_bbl_init               ! Dummy routine
671   END SUBROUTINE tra_bbl_init
672   SUBROUTINE tra_bbl( kt )              ! Dummy routine
673      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
674   END SUBROUTINE tra_bbl
675#endif
676
677   !!======================================================================
678END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.