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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization

  • Property svn:keywords set to Id
File size: 34.2 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   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
15   !!----------------------------------------------------------------------
16#if   defined key_trabbl
17   !!----------------------------------------------------------------------
18   !!   'key_trabbl'   or                             bottom boundary layer
19   !!----------------------------------------------------------------------
20   !!   tra_bbl_alloc : allocate trabbl arrays
21   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
22   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
23   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
24   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
25   !!   tra_bbl_init  : initialization, namelist read, parameters control
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and active tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constant
30   USE eosbn2         ! equation of state
31   USE trd_oce        ! trends: ocean variables
32   USE trdtra         ! trends: active tracers
33   !
34   USE iom            ! IOM library               
35   USE in_out_manager ! I/O manager
36   USE lbclnk         ! ocean lateral boundary conditions
37   USE prtctl         ! Print control
38   USE wrk_nemo       ! Memory Allocation
39   USE timing         ! Timing
40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   tra_bbl       !  routine called by step.F90
46   PUBLIC   tra_bbl_init  !  routine called by opa.F90
47   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
48   PUBLIC   tra_bbl_adv   !  -          -          -              -
49   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
50
51   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
52
53   !                                !!* Namelist nambbl *
54   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
55   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
56   !                                            !  =1 : advective bbl using the bottom ocean velocity
57   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
58   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
59   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
60
61   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
62
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
65
66   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
67   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
70
71   !! * Substitutions
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) ,                                      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      INTEGER  ::   ji, jj, jk        ! dummy loop indices
108      !
109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
113      !
114      IF( l_trdtra )   THEN                         !* Save the input trends
115         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
116!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
117         DO jk = 1, jpk
118            DO jj = 1, jpj
119               DO ji = 1, jpi
120                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
121                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
122               END DO
123            END DO
124         END DO
125      ENDIF
126
127      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
128
129      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
130         !
131         CALL tra_bbl_dif( tsb, tsa, jpts )
132         IF( ln_ctl )  &
133         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - 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( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
137         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
138         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
139         !
140      END IF
141      !
142      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
143         !
144         CALL tra_bbl_adv( tsb, tsa, jpts )
145         IF(ln_ctl)   &
146         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
147            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
148         ! lateral boundary conditions ; just need for outputs
149         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
150         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
151         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
152         !
153      END IF
154
155      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
156!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
157         DO jk = 1, jpk
158            DO jj = 1, jpj
159               DO ji = 1, jpi
160                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
161                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk)
162               END DO
163            END DO
164         END DO
165         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
166         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
167         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
168      ENDIF
169      !
170      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
171      !
172   END SUBROUTINE tra_bbl
173
174
175   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
176      !!----------------------------------------------------------------------
177      !!                  ***  ROUTINE tra_bbl_dif  ***
178      !!
179      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
180      !!                advection terms.
181      !!
182      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
183      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
184      !!      along bottom slope gradient) an additional lateral 2nd order
185      !!      diffusion along the bottom slope is added to the general
186      !!      tracer trend, otherwise the additional trend is set to 0.
187      !!      A typical value of ahbt is 2000 m2/s (equivalent to
188      !!      a downslope velocity of 20 cm/s if the condition for slope
189      !!      convection is satified)
190      !!
191      !! ** Action  :   pta   increased by the bbl diffusive trend
192      !!
193      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
194      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
195      !!----------------------------------------------------------------------
196      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
197      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
198      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
199      !
200      INTEGER  ::   ji, jj, jn   ! dummy loop indices
201      INTEGER  ::   ik           ! local integers
202      REAL(wp) ::   zbtr         ! local scalars
203      REAL(wp), POINTER, DIMENSION(:,:) :: zptb
204      !!----------------------------------------------------------------------
205      !
206      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
207      !
208      CALL wrk_alloc( jpi, jpj, zptb )
209      !
210      DO jn = 1, kjpt                                     ! tracer loop
211         !                                                ! ===========
212!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
213         DO jj = 1, jpj
214            DO ji = 1, jpi
215               ik = mbkt(ji,jj)                             ! bottom T-level index
216               zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S
217            END DO
218         END DO
219         !               
220!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
221         DO jj = 2, jpjm1                                    ! Compute the trend
222            DO ji = 2, jpim1
223               ik = mbkt(ji,jj)                            ! bottom T-level index
224               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
225                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
226                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
227                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
228                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
229                  &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik)
230            END DO
231         END DO
232         !                                                  ! ===========
233      END DO                                                ! end tracer
234      !                                                     ! ===========
235      CALL wrk_dealloc( jpi, jpj, zptb )
236      !
237      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
238      !
239   END SUBROUTINE tra_bbl_dif
240
241
242   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
243      !!----------------------------------------------------------------------
244      !!                  ***  ROUTINE trc_bbl  ***
245      !!
246      !! ** Purpose :   Compute the before passive tracer trend associated
247      !!     with the bottom boundary layer and add it to the general trend
248      !!     of tracer equations.
249      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
250      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
251      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
252      !!                       transport proportional to the along-slope density gradient
253      !!
254      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
255      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
256      !!----------------------------------------------------------------------
257      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
258      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
259      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
260      !
261      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
262      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
263      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
264      REAL(wp) ::   zbtr, ztra               ! local scalars
265      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
266      !!----------------------------------------------------------------------
267      !
268      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
269      !                                                          ! ===========
270      DO jn = 1, kjpt                                            ! tracer loop
271         !                                                       ! ===========
272         DO jj = 1, jpjm1
273            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
274               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
275                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
276                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
277                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
278                  zu_bbl = ABS( utr_bbl(ji,jj) )
279                  !
280                  !                                               ! up  -slope T-point (shelf bottom point)
281                  zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus)
282                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
283                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
284                  !
285                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
286                     zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk)
287                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
288                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
289                  END DO
290                  !
291                  zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud)
292                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
293                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
294               ENDIF
295               !
296               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
297                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
298                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
299                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
300                  zv_bbl = ABS( vtr_bbl(ji,jj) )
301                  !
302                  ! up  -slope T-point (shelf bottom point)
303                  zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs)
304                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
305                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
306                  !
307                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
308                     zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk)
309                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
310                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
311                  END DO
312                  !                                               ! down-slope T-point (deep bottom point)
313                  zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd)
314                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
315                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
316               ENDIF
317            END DO
318            !
319         END DO
320         !                                                       ! ===========
321      END DO                                                     ! end tracer
322      !                                                          ! ===========
323      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
324      !
325   END SUBROUTINE tra_bbl_adv
326
327
328   SUBROUTINE bbl( kt, kit000, cdtype )
329      !!----------------------------------------------------------------------
330      !!                  ***  ROUTINE bbl  ***
331      !!
332      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
333      !!                advection terms.
334      !!
335      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
336      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
337      !!      along bottom slope gradient) an additional lateral 2nd order
338      !!      diffusion along the bottom slope is added to the general
339      !!      tracer trend, otherwise the additional trend is set to 0.
340      !!      A typical value of ahbt is 2000 m2/s (equivalent to
341      !!      a downslope velocity of 20 cm/s if the condition for slope
342      !!      convection is satified)
343      !!              * advective bbl (nn_bbl_adv=1 or 2) :
344      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
345      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
346      !!        i.e. transport proportional to the along-slope density gradient
347      !!
348      !!      NB: the along slope density gradient is evaluated using the
349      !!      local density (i.e. referenced at a common local depth).
350      !!
351      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
352      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
353      !!----------------------------------------------------------------------
354      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
355      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
356      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
357      !
358      INTEGER  ::   ji, jj                    ! dummy loop indices
359      INTEGER  ::   ik                        ! local integers
360      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
361      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
362      REAL(wp) ::   za, zb, zgdrho            ! local scalars
363      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
364      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
365      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
366      !!----------------------------------------------------------------------
367      !
368      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
369      !
370      IF( kt == kit000 )  THEN
371         IF(lwp)  WRITE(numout,*)
372         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
373         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
374      ENDIF
375      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
376!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
377      DO jj = 1, jpj
378         DO ji = 1, jpi
379            ik = mbkt(ji,jj)                             ! bottom T-level index
380            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
381            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
382            !
383            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth
384            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
385            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
386         END DO
387      END DO
388      !
389      CALL eos_rab( zts, zdep, zab )
390      !
391      !                                   !-------------------!
392      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
393         !                                !-------------------!
394!$OMP PARALLEL DO schedule(static) private(jj,ji,za,zb,zgdrho,zsign)
395         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
396            DO ji = 1, fs_jpim1   ! vector opt.
397               !                                                   ! i-direction
398               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
399               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
400               !                                                         ! 2*masked bottom density gradient
401               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
402                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
403               !
404               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
405               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
406               !
407               !                                                   ! j-direction
408               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
409               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
410               !                                                         ! 2*masked bottom density gradient
411               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
412                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
413               !
414               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
415               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
416            END DO
417         END DO
418         !
419      ENDIF
420      !
421      !                                   !-------------------!
422      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
423         !                                !-------------------!
424         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
425         !
426         CASE( 1 )                                   != use of upper velocity
427!$OMP PARALLEL DO schedule(static) private(jj,ji,za,zb,zgdrho,zsign,zsigna)
428            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
429               DO ji = 1, fs_jpim1   ! vector opt.
430                  !                                                  ! i-direction
431                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
432                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
433                  !                                                          ! 2*masked bottom density gradient
434                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
435                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
436                  !
437                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
438                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
439                  !
440                  !                                                          ! bbl velocity
441                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
442                  !
443                  !                                                  ! j-direction
444                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
445                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
446                  !                                                          ! 2*masked bottom density gradient
447                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
448                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
449                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
450                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
451                  !
452                  !                                                          ! bbl transport
453                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
454               END DO
455            END DO
456            !
457         CASE( 2 )                                 != bbl velocity = F( delta rho )
458            zgbbl = grav * rn_gambbl
459!$OMP PARALLEL DO schedule(static) private(jj,ji,iid,iis,ikud,ikus,za,zb,zgdrho,ijd,ijs,ikvd,ikvs)
460            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
461               DO ji = 1, fs_jpim1   ! vector opt.
462                  !                                                  ! i-direction
463                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
464                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
465                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
466                  !
467                  ikud = mbku_d(ji,jj)
468                  ikus = mbku(ji,jj)
469                  !
470                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
471                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
472                  !                                                          !   masked bottom density gradient
473                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
474                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
475                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
476                  !
477                  !                                                          ! bbl transport (down-slope direction)
478                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
479                  !
480                  !                                                  ! j-direction
481                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
482                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
483                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
484                  !
485                  ikvd = mbkv_d(ji,jj)
486                  ikvs = mbkv(ji,jj)
487                  !
488                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
489                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
490                  !                                                          !   masked bottom density gradient
491                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
492                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
493                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
494                  !
495                  !                                                          ! bbl transport (down-slope direction)
496                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
497               END DO
498            END DO
499         END SELECT
500         !
501      ENDIF
502      !
503      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
504      !
505   END SUBROUTINE bbl
506
507
508   SUBROUTINE tra_bbl_init
509      !!----------------------------------------------------------------------
510      !!                  ***  ROUTINE tra_bbl_init  ***
511      !!
512      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
513      !!
514      !! ** Method  :   Read the nambbl namelist and check the parameters
515      !!              called by nemo_init at the first timestep (kit000)
516      !!----------------------------------------------------------------------
517      INTEGER ::   ji, jj               ! dummy loop indices
518      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
519      INTEGER ::   ios                  !   -      -
520      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
521      !
522      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
523      !!----------------------------------------------------------------------
524      !
525      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
526      !
527      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
528      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
529901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
530      !
531      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
532      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
533902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
534      IF(lwm) WRITE ( numond, nambbl )
535      !
536      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
537      !
538      IF(lwp) THEN                   !* Parameter control and print
539         WRITE(numout,*)
540         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
541         WRITE(numout,*) '~~~~~~~~~~~~'
542         WRITE(numout,*) '   Namelist nambbl : set bbl parameters'
543         WRITE(numout,*) '      diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
544         WRITE(numout,*) '      advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
545         WRITE(numout,*) '      diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
546         WRITE(numout,*) '      advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
547      ENDIF
548
549      !                              ! allocate trabbl arrays
550      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
551
552      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
553      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
554
555      !                             !* vertical index of  "deep" bottom u- and v-points
556!$OMP PARALLEL DO schedule(static) private(jj,ji)
557      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
558         DO ji = 1, jpim1
559            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
560            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
561         END DO
562      END DO
563      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
564      CALL wrk_alloc( jpi, jpj, zmbk )
565      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
566      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
567      CALL wrk_dealloc( jpi, jpj, zmbk )
568
569      !                                 !* sign of grad(H) at u- and v-points
570      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
571!$OMP PARALLEL DO schedule(static) private(jj,ji)
572      DO jj = 1, jpjm1
573         DO ji = 1, jpim1
574            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
575            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
576         END DO
577      END DO
578      !
579!$OMP PARALLEL DO schedule(static) private(jj,ji)
580      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
581         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
582            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
583            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
584         END DO
585      END DO
586      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
587      !
588      !                             !* masked diffusive flux coefficients
589!$OMP PARALLEL DO schedule(static) private(jj,ji)
590      DO jj = 1, jpj
591         DO ji = 1, jpi
592            ahu_bbl_0(ji,jj) = rn_ahtbbl * e2_e1u(ji,jj) * e3u_bbl_0(ji,jj) * umask(ji,jj,1)
593            ahv_bbl_0(ji,jj) = rn_ahtbbl * e1_e2v(ji,jj) * e3v_bbl_0(ji,jj) * vmask(ji,jj,1)
594         END DO
595      END DO
596
597      !
598      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
599      !
600   END SUBROUTINE tra_bbl_init
601
602#else
603   !!----------------------------------------------------------------------
604   !!   Dummy module :                      No bottom boundary layer scheme
605   !!----------------------------------------------------------------------
606   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
607CONTAINS
608   SUBROUTINE tra_bbl_init               ! Dummy routine
609   END SUBROUTINE tra_bbl_init
610   SUBROUTINE tra_bbl( kt )              ! Dummy routine
611      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
612   END SUBROUTINE tra_bbl
613#endif
614
615   !!======================================================================
616END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.