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 branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

Updates to average physics variables for TOP substepping.

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