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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

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