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 NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA/trabbl.F90 @ 13409

Last change on this file since 13409 was 13409, checked in by hadcv, 4 years ago

Remaining changes prior to trunk merge

  • Property svn:keywords set to Id
File size: 32.5 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   !!            4.0  ! 2017-04  (G. Madec)  ln_trabbl namelist variable instead of a CPP key
16   !!----------------------------------------------------------------------
17
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   ! TEMP: This change not necessary after trd_tra is tiled
29   USE domain, ONLY : dom_tile
30   USE phycst         ! physical constant
31   USE eosbn2         ! equation of state
32   USE trd_oce        ! trends: ocean variables
33   USE trdtra         ! trends: active tracers
34   !
35   USE iom            ! IOM library               
36   USE in_out_manager ! I/O manager
37   USE lbclnk         ! ocean lateral boundary conditions
38   USE prtctl         ! Print control
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 nemogcm.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   !                                !!* Namelist nambbl *
52   LOGICAL , PUBLIC ::   ln_trabbl   !: bottom boundary layer flag
53   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
54   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
55   !                                            !  =1 : advective bbl using the bottom ocean velocity
56   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
57   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
58   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
59
60   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
61
62   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
64
65   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
66   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
68   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)
69
70   !! * Substitutions
71#  include "do_loop_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
74   !! $Id$
75   !! Software governed by the CeCILL license (see ./LICENSE)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   INTEGER FUNCTION tra_bbl_alloc()
80      !!----------------------------------------------------------------------
81      !!                ***  FUNCTION tra_bbl_alloc  ***
82      !!----------------------------------------------------------------------
83      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
84         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
85         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
86         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
87         !
88      CALL mpp_sum ( 'trabbl', tra_bbl_alloc )
89      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
90   END FUNCTION tra_bbl_alloc
91
92
93   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
94      !!----------------------------------------------------------------------
95      !!                  ***  ROUTINE bbl  ***
96      !!
97      !! ** Purpose :   Compute the before tracer (t & s) trend associated
98      !!              with the bottom boundary layer and add it to the general
99      !!              trend of tracer equations.
100      !!
101      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
102      !!              diffusive and/or advective contribution to the tracer trend
103      !!              is added to the general tracer trend
104      !!----------------------------------------------------------------------
105      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step
106      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
107      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
108      !
109      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
110      REAL(wp), SAVE :: zsum1, zsum2, zsum3, zsum4
111      ! TEMP: This change not necessary after trd_tra is tiled
112      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ztrdt, ztrds
113      !!----------------------------------------------------------------------
114      !
115      IF( ln_timing )   CALL timing_start( 'tra_bbl')
116      !
117      IF( l_trdtra )   THEN                         !* Save the T-S input trends
118         IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
119            ! TEMP: This can be A2D after trd_tra is tiled
120            ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
121         ENDIF
122
123         DO_3D_00_00( 1, jpkm1 )
124            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs)
125            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs)
126         END_3D
127      ENDIF
128
129      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
130
131      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
132         !
133         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
134         IF( sn_cfctl%l_prtctl )  &
135         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, psum1=zsum1, &
136            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, psum2=zsum2, &
137            &          clinfo3='tra' )
138         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
139            ! lateral boundary conditions ; just need for outputs
140            CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
141            CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
142            CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
143         ENDIF
144         !
145      ENDIF
146      !
147      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
148         !
149         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
150         IF(sn_cfctl%l_prtctl)   &
151         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, psum1=zsum3, &
152            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, psum2=zsum4, &
153            &          clinfo3='tra' )
154         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
155            ! lateral boundary conditions ; just need for outputs
156            CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
157            CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
158            CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
159         ENDIF
160         !
161      ENDIF
162
163      ! TEMP: These changes not necessary after trd_tra is tiled
164      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
165         DO_3D_00_00( 1, jpkm1 )
166            ztrdt(ji,jj,jk) = pts(ji,jj,jk,jp_tem,Krhs) - ztrdt(ji,jj,jk)
167            ztrds(ji,jj,jk) = pts(ji,jj,jk,jp_sal,Krhs) - ztrds(ji,jj,jk)
168         END_3D
169
170         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
171            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
172
173            ! TODO: TO BE TILED- trd_tra
174            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
175            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
176            DEALLOCATE( ztrdt, ztrds )
177
178            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
179         ENDIF
180      ENDIF
181      !
182      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
183      !
184   END SUBROUTINE tra_bbl
185
186
187   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
188      !!----------------------------------------------------------------------
189      !!                  ***  ROUTINE tra_bbl_dif  ***
190      !!
191      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
192      !!                advection terms.
193      !!
194      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
195      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
196      !!      along bottom slope gradient) an additional lateral 2nd order
197      !!      diffusion along the bottom slope is added to the general
198      !!      tracer trend, otherwise the additional trend is set to 0.
199      !!      A typical value of ahbt is 2000 m2/s (equivalent to
200      !!      a downslope velocity of 20 cm/s if the condition for slope
201      !!      convection is satified)
202      !!
203      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
204      !!
205      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
206      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
207      !!----------------------------------------------------------------------
208      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
209      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
210      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
211      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
212      !
213      INTEGER  ::   ji, jj, jn   ! dummy loop indices
214      INTEGER  ::   ik           ! local integers
215      REAL(wp) ::   zbtr         ! local scalars
216      REAL(wp), DIMENSION(A2D) ::   zptb   ! workspace
217      !!----------------------------------------------------------------------
218      !
219      DO jn = 1, kjpt                                     ! tracer loop
220         !                                                ! ===========
221         DO_2D_11_11
222            ik = mbkt(ji,jj)                             ! bottom T-level index
223            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
224         END_2D
225         !               
226         DO_2D_00_00
227            ik = mbkt(ji,jj)                            ! bottom T-level index
228            pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
229               &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
230               &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
231               &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
232               &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
233               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
234         END_2D
235         !                                                  ! ===========
236      END DO                                                ! end tracer
237      !                                                     ! ===========
238   END SUBROUTINE tra_bbl_dif
239
240
241   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
242      !!----------------------------------------------------------------------
243      !!                  ***  ROUTINE trc_bbl  ***
244      !!
245      !! ** Purpose :   Compute the before passive tracer trend associated
246      !!     with the bottom boundary layer and add it to the general trend
247      !!     of tracer equations.
248      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
249      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
250      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
251      !!                       transport proportional to the along-slope density gradient
252      !!
253      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
254      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
255      !!----------------------------------------------------------------------
256      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
257      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
258      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
259      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
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      !                                                          ! ===========
269      DO jn = 1, kjpt                                            ! tracer loop
270         !                                                       ! ===========
271         DO_2D_10_10            ! CAUTION start from i=1 to update i=2 when cyclic east-west
272            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
273               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
274               iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
275               ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
276               zu_bbl = ABS( utr_bbl(ji,jj) )
277               !
278               !                                               ! up  -slope T-point (shelf bottom point)
279               zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
280               ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
281               pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
282               !
283               DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
284                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
285                  ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
286                  pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
287               END DO
288               !
289               zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
290               ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
291               pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
292            ENDIF
293            !
294            IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
295               ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
296               ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
297               ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
298               zv_bbl = ABS( vtr_bbl(ji,jj) )
299               !
300               ! up  -slope T-point (shelf bottom point)
301               zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
302               ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
303               pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
304               !
305               DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
306                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
307                  ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
308                  pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
309               END DO
310               !                                               ! down-slope T-point (deep bottom point)
311               zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
312               ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
313               pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
314            ENDIF
315         END_2D
316         !                                                       ! ===========
317      END DO                                                     ! end tracer
318      !                                                          ! ===========
319   END SUBROUTINE tra_bbl_adv
320
321
322   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
323      !!----------------------------------------------------------------------
324      !!                  ***  ROUTINE bbl  ***
325      !!
326      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
327      !!                advection terms.
328      !!
329      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
330      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
331      !!      along bottom slope gradient) an additional lateral 2nd order
332      !!      diffusion along the bottom slope is added to the general
333      !!      tracer trend, otherwise the additional trend is set to 0.
334      !!      A typical value of ahbt is 2000 m2/s (equivalent to
335      !!      a downslope velocity of 20 cm/s if the condition for slope
336      !!      convection is satified)
337      !!              * advective bbl (nn_bbl_adv=1 or 2) :
338      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
339      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
340      !!        i.e. transport proportional to the along-slope density gradient
341      !!
342      !!      NB: the along slope density gradient is evaluated using the
343      !!      local density (i.e. referenced at a common local depth).
344      !!
345      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
346      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
347      !!----------------------------------------------------------------------
348      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
349      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
350      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
351      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
352      !
353      INTEGER  ::   ji, jj                    ! dummy loop indices
354      INTEGER  ::   ik                        ! local integers
355      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
356      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
357      REAL(wp) ::   za, zb, zgdrho            ! local scalars
358      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
359      REAL(wp), DIMENSION(A2D,jpts)   :: zts, zab         ! 3D workspace
360      REAL(wp), DIMENSION(A2D)        :: zub, zvb, zdep   ! 2D workspace
361      !!----------------------------------------------------------------------
362      !
363      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
364         IF( kt == kit000 )  THEN
365            IF(lwp)  WRITE(numout,*)
366            IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
367            IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
368         ENDIF
369      ENDIF
370      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
371      DO_2D_11_11
372         ik = mbkt(ji,jj)                             ! bottom T-level index
373         zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
374         zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
375         !
376         zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
377         zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
378         zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
379      END_2D
380      !
381      CALL eos_rab( zts, zdep, zab, Kmm )
382      !
383      !                                   !-------------------!
384      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
385         !                                !-------------------!
386         DO_2D_10_10
387            !                                                   ! i-direction
388            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
389            zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
390            !                                                         ! 2*masked bottom density gradient
391            zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
392               &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
393            !
394            zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
395            ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
396            !
397            !                                                   ! j-direction
398            za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
399            zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
400            !                                                         ! 2*masked bottom density gradient
401            zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
402               &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
403            !
404            zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
405            ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
406         END_2D
407         !
408      ENDIF
409      !
410      !                                   !-------------------!
411      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
412         !                                !-------------------!
413         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
414         !
415         CASE( 1 )                                   != use of upper velocity
416            DO_2D_10_10
417               !                                                  ! i-direction
418               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
419               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
420               !                                                          ! 2*masked bottom density gradient
421               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
422                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
423               !
424               zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
425               zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
426               !
427               !                                                          ! bbl velocity
428               utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
429               !
430               !                                                  ! j-direction
431               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
432               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
433               !                                                          ! 2*masked bottom density gradient
434               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
435                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
436               zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
437               zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
438               !
439               !                                                          ! bbl transport
440               vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
441            END_2D
442            !
443         CASE( 2 )                                 != bbl velocity = F( delta rho )
444            zgbbl = grav * rn_gambbl
445            DO_2D_10_10
446               !                                                  ! i-direction
447               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
448               iid  = ji + MAX( 0, mgrhu(ji,jj) )
449               iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
450               !
451               ikud = mbku_d(ji,jj)
452               ikus = mbku(ji,jj)
453               !
454               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
455               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
456               !                                                          !   masked bottom density gradient
457               zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
458                  &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
459               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
460               !
461               !                                                          ! bbl transport (down-slope direction)
462               utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
463               !
464               !                                                  ! j-direction
465               !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
466               ijd  = jj + MAX( 0, mgrhv(ji,jj) )
467               ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
468               !
469               ikvd = mbkv_d(ji,jj)
470               ikvs = mbkv(ji,jj)
471               !
472               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
473               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
474               !                                                          !   masked bottom density gradient
475               zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
476                  &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
477               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
478               !
479               !                                                          ! bbl transport (down-slope direction)
480               vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
481            END_2D
482         END SELECT
483         !
484      ENDIF
485      !
486   END SUBROUTINE bbl
487
488
489   SUBROUTINE tra_bbl_init
490      !!----------------------------------------------------------------------
491      !!                  ***  ROUTINE tra_bbl_init  ***
492      !!
493      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
494      !!
495      !! ** Method  :   Read the nambbl namelist and check the parameters
496      !!              called by nemo_init at the first timestep (kit000)
497      !!----------------------------------------------------------------------
498      INTEGER ::   ji, jj                      ! dummy loop indices
499      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
500      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
501      !!
502      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
503      !!----------------------------------------------------------------------
504      !
505      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
506901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
507      !
508      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
509902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
510      IF(lwm) WRITE ( numond, nambbl )
511      !
512      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
513      !
514      IF(lwp) THEN                   !* Parameter control and print
515         WRITE(numout,*)
516         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
517         WRITE(numout,*) '~~~~~~~~~~~~'
518         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
519         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
520      ENDIF
521      IF( .NOT.ln_trabbl )   RETURN
522      !
523      IF(lwp) THEN
524         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
525         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
526         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
527         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
528      ENDIF
529      !
530      !                              ! allocate trabbl arrays
531      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
532      !
533      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
534      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
535      !
536      !                             !* vertical index of  "deep" bottom u- and v-points
537      DO_2D_10_10
538         mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
539         mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
540      END_2D
541      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
542      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
543      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
544      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
545      !
546      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
547      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
548      DO_2D_10_10
549         IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
550            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
551         ENDIF
552         !
553         IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
554            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
555         ENDIF
556      END_2D
557      !
558      DO_2D_10_10
559         e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
560         e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
561      END_2D
562      CALL lbc_lnk_multi( 'trabbl', e3u_bbl_0, 'U', 1. , e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
563      !
564      !                             !* masked diffusive flux coefficients
565      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
566      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
567      !
568   END SUBROUTINE tra_bbl_init
569
570   !!======================================================================
571END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.