source: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS_SAL/MY_SRC/trabbl_crs.F90 @ 4315

Last change on this file since 4315 was 4315, checked in by cetlod, 7 years ago

dev_IO_CRS : update trabbl_crs module

  • Property svn:executable set to *
File size: 38.9 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         ! 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                  ! up  -slope T-point (shelf bottom point)
309                  zbtr = r1_bt_crs(ji,ijs,ikvs)
310                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
311                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
312                  !                   
313                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
314                     zbtr = r1_bt_crs(ji,ijd,jk)
315                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
316                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
317                  END DO
318                  !                                               ! down-slope T-point (deep bottom point)
319                  zbtr = r1_bt_crs(ji,ijd,ikvd)
320                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
321                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
322               ENDIF
323            END DO
324            !
325         END DO
326         !                                                  ! ===========
327      END DO                                                ! end tracer
328      !                                                     ! ===========
329      !
330      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
331      !
332   END SUBROUTINE tra_bbl_adv_crs
333
334
335   SUBROUTINE bbl_crs( kt, kit000, cdtype )
336      !!----------------------------------------------------------------------
337      !!                  ***  ROUTINE bbl  ***
338      !!                   
339      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
340      !!                advection terms.
341      !!
342      !! ** Method  :   
343      !!        * diffusive bbl (nn_bbl_ldf=1) :
344      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
345      !!      along bottom slope gradient) an additional lateral 2nd order
346      !!      diffusion along the bottom slope is added to the general
347      !!      tracer trend, otherwise the additional trend is set to 0.
348      !!      A typical value of ahbt is 2000 m2/s (equivalent to
349      !!      a downslope velocity of 20 cm/s if the condition for slope
350      !!      convection is satified)
351      !!        * advective bbl (nn_bbl_adv=1 or 2) :
352      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
353      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
354      !!        i.e. transport proportional to the along-slope density gradient
355      !!
356      !!      NB: the along slope density gradient is evaluated using the
357      !!      local density (i.e. referenced at a common local depth).
358      !!
359      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
360      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
361      !!---------------------------------------------------------------------- 
362      !
363      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
364      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
365      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
366      !!
367      INTEGER  ::   ji, jj                    ! dummy loop indices
368      INTEGER  ::   ik                        ! local integers
369      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
370      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
371      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
372      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
373      !!
374      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
375      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep
376      !!----------------------- zv_bbl-----------------------------------------------
377      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
378      ! ================            pft :  potential temperature in degrees celcius
379      !                             pfs :  salinity anomaly (s-35) in psu
380      !                             pfh :  depth in meters
381      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
382      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
383         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
384                                   - 0.203814e-03 ) * pft                    &
385                                   + 0.170907e-01 ) * pft                    &
386                                   + 0.665157e-01                            &
387         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
388         +  ( ( - 0.302285e-13 * pfh                                         &
389                - 0.251520e-11 * pfs                                         &
390                + 0.512857e-12 * pft * pft          ) * pfh                  &
391                                     - 0.164759e-06   * pfs                  &
392             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
393                                     + 0.380374e-04 ) * pfh
394      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
395         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
396                                 - 0.301985e-05 ) * pft                      &
397                                 + 0.785567e-03                              &
398         + (     0.515032e-08 * pfs                                          &
399               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
400               +(  (   0.121551e-17 * pfh                                    &
401                     - 0.602281e-15 * pfs                                    &
402                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
403                                          + 0.408195e-10   * pfs             &
404                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
405                                          - 0.121555e-07 ) * pfh
406      !!----------------------------------------------------------------------
407     
408      !
409      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
410      !
411      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
412      !
413     
414      IF( kt == kit000 )  THEN
415         IF(lwp)  WRITE(numout,*)
416         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
417         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
418      ENDIF
419     
420      !                                        !* bottom temperature, salinity, velocity and depth
421#if defined key_vectopt_loop
422      DO jj = 1, 1   ! vector opt. (forced unrolling)
423         DO ji = 1, jpij
424#else
425      DO jj = 1, jpj
426         DO ji = 1, jpi
427#endif
428            ik = mbkt_crs(ji,jj)                        ! bottom T-level index
429            ztb (ji,jj) = tsb_crs(ji,jj,ik,jp_tem) * tmask_crs(ji,jj,1)      ! bottom before T and S
430            zsb (ji,jj) = tsb_crs(ji,jj,ik,jp_sal) * tmask_crs(ji,jj,1)
431            zdep(ji,jj) = gdept_crs(ji,jj,ik)        ! bottom T-level reference depth !!
432           
433            zub(ji,jj) = un_crs(ji,jj,mbku_crs(ji,jj))      ! bottom velocity
434            zvb(ji,jj) = vn_crs(ji,jj,mbkv_crs(ji,jj)) 
435         END DO
436      END DO
437     
438      !                                   !-------------------!
439      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
440         !                                !-------------------!
441         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
442            DO ji = 1, jpim1             
443               !                                                ! i-direction
444               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
445               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
446               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
447               !                                                         ! masked bbl i-gradient of density
448               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
449                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask_crs(ji,jj,1)
450               !                                                     
451               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
452               ahu_bbl_crs(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
453               !
454               !                                                ! j-direction
455               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
456               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
457               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
458               !                                                         ! masked bbl j-gradient of density
459               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
460                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask_crs(ji,jj,1)
461               !                                                   
462               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
463               ahv_bbl_crs(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
464               !
465            END DO
466         END DO
467         !
468      ENDIF
469
470      !                                   !-------------------!
471      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
472         !                                !-------------------!
473         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
474         !
475         CASE( 1 )                                   != use of upper velocity
476            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
477               DO ji = 1, fs_jpim1   ! vector opt.
478                  !                                               ! i-direction
479                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
480                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
481                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
482                  !                                                           ! masked bbl i-gradient of density
483                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
484                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask_crs(ji,jj,1)
485                  !                                                         
486                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
487                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
488                  !
489                  !                                                           ! bbl velocity
490                  utr_bbl_crs(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u_crs(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
491                  !
492                  !                                               ! j-direction
493                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
494                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
495                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
496                  !                                                           ! masked bbl j-gradient of density
497                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
498                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask_crs(ji,jj,1)
499                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
500                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
501                  !
502                  !                                                           ! bbl velocity
503                  vtr_bbl_crs(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v_crs(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
504               END DO
505            END DO
506            !
507         CASE( 2 )                                 != bbl velocity = F( delta rho )
508            zgbbl = grav * rn_gambbl
509            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
510               DO ji = 1, fs_jpim1   ! vector opt.
511                  !                                         ! i-direction
512                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
513                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
514                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku_crs(ji,jj)
515                  !
516                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
517                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
518                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
519                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
520                  zgdrho =    fsbeta( zt, zs, zh )                                    &
521                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
522                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask_crs(ji,jj,1)
523                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
524                  !
525                  !                                             ! bbl transport (down-slope direction)
526                  utr_bbl_crs(ji,jj) = e2u_crs(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
527                  !
528                  !                                         ! j-direction
529                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
530                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
531                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv_crs(ji,jj)
532                  !
533                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
534                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
535                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
536                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
537                  zgdrho =    fsbeta( zt, zs, zh )                                    &
538                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
539                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask_crs(ji,jj,1)
540                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
541                  !
542                  !                                             ! bbl transport (down-slope direction)
543                  vtr_bbl_crs(ji,jj) = e1v_crs(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
544               END DO
545            END DO
546         END SELECT
547         !
548      ENDIF
549      !
550      CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
551      !
552      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
553      !
554   END SUBROUTINE bbl_crs
555
556
557   SUBROUTINE tra_bbl_init_crs
558      !!----------------------------------------------------------------------
559      !!                  ***  ROUTINE tra_bbl_init  ***
560      !!
561      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
562      !!
563      !! ** Method  :   Read the nambbl namelist and check the parameters
564      !!              called by nemo_init at the first timestep (kit000)
565      !!----------------------------------------------------------------------
566      INTEGER ::   ji, jj               ! dummy loop indices
567      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
568      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
569      !!
570      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
571      !!----------------------------------------------------------------------
572      !
573      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
574      !
575      CALL wrk_alloc( jpi, jpj, zmbk ) 
576      !
577
578      REWIND ( numnam )              !* Read Namelist nambbl : bottom boundary layer scheme
579      READ   ( numnam, nambbl )
580      !
581      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
582      !
583      IF(lwp) THEN                   !* Parameter control and print
584         WRITE(numout,*)
585         WRITE(numout,*) 'tra_bbl_init_crs : bottom boundary layer initialisation'
586         WRITE(numout,*) '~~~~~~~~~~~~'
587         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
588         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
589         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
590         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
591         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
592      ENDIF
593
594      !                              ! allocate trabbl arrays
595      IF( tra_bbl_alloc_crs() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
596     
597      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
598      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
599
600      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
601
602
603      !                             !* inverse of surface of T-cells
604      r1_e1e2t(:,:) = 1._wp / e1e2w_msk(:,:,1)
605     
606      !                             !* vertical index of  "deep" bottom u- and v-points
607      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
608         DO ji = 1, jpim1
609            mbku_d(ji,jj) = MAX(  mbkt_crs(ji+1,jj  ) , mbkt_crs(ji,jj)  )   ! >= 1 as mbkt=1 over land
610            mbkv_d(ji,jj) = MAX(  mbkt_crs(ji  ,jj+1) , mbkt_crs(ji,jj)  )
611         END DO
612      END DO
613      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
614      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL crs_lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
615      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL crs_lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
616
617                                        !* sign of grad(H) at u- and v-points
618      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
619      DO jj = 1, jpjm1               
620         DO ji = 1, jpim1
621            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
622            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)) )  )
623         END DO
624      END DO
625
626      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
627         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
628            e3u_bbl_0(ji,jj) = MIN( e3u_crs(ji,jj,mbkt_crs(ji+1,jj  )), e3u_crs(ji,jj,mbkt_crs(ji,jj)) ) 
629            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
630         END DO
631      END DO
632      CALL crs_lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL crs_lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
633
634      !                             !* masked diffusive flux coefficients
635      ahu_bbl_0(:,:) = rn_ahtbbl_crs * e2u_crs(:,:) * e3u_bbl_0(:,:) / e1u_crs(:,:)  * umask_crs(:,:,1)
636      ahv_bbl_0(:,:) = rn_ahtbbl_crs * e1v_crs(:,:) * e3v_bbl_0(:,:) / e2v_crs(:,:)  * vmask_crs(:,:,1)
637
638!! cc reflexion for the coarsening
639
640!      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
641         !
642!         SELECT CASE ( jp_cfg )
643!         CASE ( 2 )                          ! ORCA_R2
644!            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
645!            ii0 = 139   ;   ii1 = 140 
646!            ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) &
647!          & =  4.e0*ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1))
648!            ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) &
649!          & =  4.e0*ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1))
650!            !
651!            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
652!            ii0 = 161   ;   ii1 = 162
653!            ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1(ij1)) &
654!          & = 10.e0*ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1))
655!            ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1(ij1)) &
656!          & = 10.e0*ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1))
657!            !
658!         CASE ( 4 )                          ! ORCA_R4
659!            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL  !!cc
660!            ii0 =  70   ;   ii1 =  71 
661!            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
662!            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
663!         END SELECT
664!         !
665!      ENDIF
666      !
667      CALL wrk_dealloc( jpi, jpj, zmbk ) 
668      !
669      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
670      !
671   END SUBROUTINE tra_bbl_init_crs
672
673#else
674   !!----------------------------------------------------------------------
675   !!   Dummy module :                      No bottom boundary layer scheme
676   !!----------------------------------------------------------------------
677   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl_crs = .FALSE.   !: bbl flag
678CONTAINS
679   SUBROUTINE tra_bbl_init_crs               ! Dummy routine
680   END SUBROUTINE tra_bbl_init_crs
681   SUBROUTINE tra_bbl_crs( kt )              ! Dummy routine
682      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
683   END SUBROUTINE tra_bbl_crs
684#endif
685
686   !!======================================================================
687END MODULE trabbl_crs
Note: See TracBrowser for help on using the repository browser.