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_crs.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl_crs.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by cbricaud, 9 years ago

bug correction

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