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 @ 2971

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

Latest changes for substepping - add key_trabbl.

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