New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

trabbl.F90 in branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 5 years ago
Introduction of stochastic physics in NEMO, based on Andrea Storto's code. For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.
File size: 36.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: