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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2330

Last change on this file since 2330 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
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	!!==============================================================================
[2024]	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
[503]	13	!!----------------------------------------------------------------------
[2024]	14	#if defined key_trabbl \|\| defined key_esopa
[3]	15	!!----------------------------------------------------------------------
[2024]	16	!! 'key_trabbl' or bottom boundary layer
[3]	17	!!----------------------------------------------------------------------
[2024]	18	!! tra_bbl : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
	19	!! tra_bbl_dif : generic routine to compute bbl diffusive trend
	20	!! tra_bbl_adv : generic routine to compute bbl advective trend
	21	!! bbl : computation of bbl diffu. flux coef. & transport in bottom boundary layer
[2106]	22	!! tra_bbl_init : initialization, namelist read, parameters control
[503]	23	!!----------------------------------------------------------------------
[2024]	24	USE oce ! ocean dynamics and active tracers
	25	USE dom_oce ! ocean space and time domain
	26	USE phycst !
	27	USE eosbn2 ! equation of state
[2104]	28	USE trdmod_oce ! ocean space and time domain
	29	USE trdtra ! ocean active tracers trends
[2024]	30	USE iom ! IOM server
	31	USE in_out_manager ! I/O manager
	32	USE lbclnk ! ocean lateral boundary conditions
	33	USE prtctl ! Print control
[3]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
[2024]	38	PUBLIC tra_bbl ! routine called by step.F90
	39	PUBLIC tra_bbl_init ! routine called by opa.F90
	40	PUBLIC tra_bbl_dif ! routine called by tra_bbl and trc_bbl
	41	PUBLIC tra_bbl_adv ! - - - -
	42	PUBLIC bbl ! - - - -
[3]	43
[2024]	44	# if defined key_trabbl
	45	LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .TRUE. !: bottom boundary layer flag
[457]	46	# else
[2024]	47	LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .FALSE. !: bottom boundary layer flag
[457]	48	# endif
[409]	49
[2024]	50	! !!* Namelist nambbl *
	51	INTEGER , PUBLIC :: nn_bbl_ldf = 0 !: =1 : diffusive bbl or not (=0)
	52	INTEGER , PUBLIC :: nn_bbl_adv = 0 !: =1/2 : advective bbl or not (=0)
[2106]	53	! ! =1 : advective bbl using the bottom ocean velocity
	54	! ! =2 : - - using utr_bbl proportional to grad(rho)
[2024]	55	REAL(wp), PUBLIC :: rn_ahtbbl = 1.e+3 !: along slope bbl diffusive coefficient [m2/s]
	56	REAL(wp), PUBLIC :: rn_gambbl = 10.e0 !: lateral coeff. for bottom boundary layer scheme [s]
[3]	57
[2106]	58	REAL(wp), DIMENSION(jpi,jpj), PUBLIC :: utr_bbl , vtr_bbl ! u- (v-) transport in the bottom boundary layer
	59
[2024]	60	INTEGER , DIMENSION(jpi,jpj) :: mbkt ! vertical index of the bottom ocean T-level
	61	INTEGER , DIMENSION(jpi,jpj) :: mbku , mbkv ! vertical index of the (upper) bottom ocean U/V-level
	62	INTEGER , DIMENSION(jpi,jpj) :: mbku_d , mbkv_d ! vertical index of the "lower" bottom ocean U/V-level
	63	INTEGER , DIMENSION(jpi,jpj) :: mgrhu , mgrhv ! = +/-1, sign of grad(H) in u-(v-)direction
	64	REAL(wp), DIMENSION(jpi,jpj) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points
	65	REAL(wp), DIMENSION(jpi,jpj) :: ahu_bbl , ahv_bbl ! masked diffusive bbl coefficients at u and v-points
	66	REAL(wp), DIMENSION(jpi,jpj) :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points
[2106]	67	REAL(wp), DIMENSION(jpi,jpj) :: e1e2t_r ! thichness of the bbl (e3) at u and v-points
	68	LOGICAL, PUBLIC :: l_bbl !: flag to compute bbl diffu. flux coef and transport
[3]	69
	70	!! * Substitutions
	71	# include "domzgr_substitute.h90"
	72	# include "vectopt_loop_substitute.h90"
	73	!!----------------------------------------------------------------------
[2287]	74	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[2104]	75	!! $Id$
[2287]	76	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	77	!!----------------------------------------------------------------------
	78
	79	CONTAINS
	80
[2106]	81
[2024]	82	SUBROUTINE tra_bbl( kt )
[3]	83	!!----------------------------------------------------------------------
[2024]	84	!! * ROUTINE bbl *
	85	!!
[3]	86	!! ** Purpose : Compute the before tracer (t & s) trend associated
[2024]	87	!! with the bottom boundary layer and add it to the general trend
	88	!! of tracer equations.
[3]	89	!!
[2106]	90	!! ** Method : Depending on namtra_bbl namelist parameters the bbl
	91	!! diffusive and/or advective contribution to the tracer trend
	92	!! is added to the general tracer trend
[2024]	93	!!----------------------------------------------------------------------
	94	INTEGER, INTENT( in ) :: kt ! ocean time-step
[2104]	95	!!
[2024]	96	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
	97	!!----------------------------------------------------------------------
	98
	99	IF( l_trdtra ) THEN !* Save ta and sa trends
	100	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
	101	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
	102	ENDIF
	103
[2123]	104	IF( l_bbl ) CALL bbl( kt, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl)
[2024]	105
[2106]	106
	107	IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl
[2024]	108	CALL tra_bbl_dif( tsb, tsa, jpts )
	109	IF( ln_ctl ) &
	110	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, &
	111	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	112	! lateral boundary conditions ; just need for outputs
	113	CALL lbc_lnk( ahu_bbl, 'U', 1. ) ; CALL lbc_lnk( ahv_bbl, 'V', 1. )
	114	CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef
	115	CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef
	116	END IF
	117
[2106]	118	IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl
	119	CALL tra_bbl_adv( tsb, tsa, jpts )
[2024]	120	IF(ln_ctl) &
	121	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv - Ta: ', mask1=tmask, &
	122	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	123	! lateral boundary conditions ; just need for outputs
[2106]	124	CALL lbc_lnk( utr_bbl, 'U', 1. ) ; CALL lbc_lnk( vtr_bbl, 'V', 1. )
[2024]	125	CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport
	126	CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport
	127	END IF
	128
	129	IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics
	130	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
	131	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
	132	CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
	133	CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
	134	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
	135	ENDIF
	136	!
	137	END SUBROUTINE tra_bbl
	138
	139
[2148]	140	SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
[2024]	141	!!----------------------------------------------------------------------
	142	!! * ROUTINE tra_bbl_dif *
	143	!!
	144	!! ** Purpose : Computes the bottom boundary horizontal and vertical
	145	!! advection terms.
	146	!!
	147	!! ** Method :
	148	!! * diffusive bbl (nn_bbl_ldf=1) :
	149	!! When the product grad( rho) * grad(h) < 0 (where grad is an
	150	!! along bottom slope gradient) an additional lateral 2nd order
	151	!! diffusion along the bottom slope is added to the general
	152	!! tracer trend, otherwise the additional trend is set to 0.
	153	!! A typical value of ahbt is 2000 m2/s (equivalent to
[3]	154	!! a downslope velocity of 20 cm/s if the condition for slope
	155	!! convection is satified)
	156	!!
[2148]	157	!! ** Action : pta increased by the bbl diffusive trend
[2106]	158	!!
[503]	159	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
[2024]	160	!! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
	161	!!----------------------------------------------------------------------
[2104]	162	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[2148]	163	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	164	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2104]	165	!!
	166	INTEGER :: ji, jj, jn ! dummy loop indices
[2106]	167	INTEGER :: ik ! local integers
[2123]	168	REAL(wp) :: zbtr ! local scalars
[2148]	169	REAL(wp), DIMENSION(jpi,jpj) :: zptb ! tracer trend
[2104]	170	!!----------------------------------------------------------------------
[2024]	171	!
	172	DO jn = 1, kjpt ! tracer loop
	173	! ! ===========
	174	# if defined key_vectopt_loop
	175	DO jj = 1, 1 ! vector opt. (forced unrolling)
[2148]	176	DO ji = 1, jpij
	177	#else
	178	DO jj = 1, jpj
	179	DO ji = 1, jpi
	180	#endif
	181	ik = mbkt(ji,jj) ! bottom T-level index
	182	zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S
	183	END DO
	184	END DO
	185	! ! Compute the trend
	186	# if defined key_vectopt_loop
	187	DO jj = 1, 1 ! vector opt. (forced unrolling)
[2024]	188	DO ji = jpi+1, jpij-jpi-1
	189	# else
	190	DO jj = 2, jpjm1
	191	DO ji = 2, jpim1
	192	# endif
[2106]	193	ik = mbkt(ji,jj) ! bottom T-level index
[2024]	194	zbtr = e1e2t_r(ji,jj) / fse3t(ji,jj,ik)
[2148]	195	pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) &
	196	& + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) &
	197	& - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) &
	198	& + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) &
	199	& - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) * zbtr
[2024]	200	END DO
	201	END DO
[2106]	202	! ! ===========
	203	END DO ! end tracer
	204	! ! ===========
[2024]	205	END SUBROUTINE tra_bbl_dif
[2104]	206
[2024]	207
[2148]	208	SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
[2024]	209	!!----------------------------------------------------------------------
	210	!! * ROUTINE trc_bbl *
[2106]	211	!!
[2024]	212	!! ** Purpose : Compute the before passive tracer trend associated
	213	!! with the bottom boundary layer and add it to the general trend
	214	!! of tracer equations.
[2106]	215	!! ** Method : advective bbl (nn_bbl_adv = 1 or 2) :
	216	!! nn_bbl_adv = 1 use of the ocean near bottom velocity as bbl velocity
	217	!! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e.
	218	!! transport proportional to the along-slope density gradient
[3]	219	!!
[503]	220	!!
[2024]	221	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
	222	!! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
	223	!!
	224	!!----------------------------------------------------------------------
[2104]	225	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[2148]	226	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	227	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2104]	228	!!
[2024]	229	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2123]	230	INTEGER :: iis , iid , ijs , ijd ! local integers
[2104]	231	INTEGER :: ikus, ikud, ikvs, ikvd ! - -
	232	REAL(wp) :: zbtr, ztra ! local scalars
	233	REAL(wp) :: zu_bbl, zv_bbl ! - -
[3]	234	!!----------------------------------------------------------------------
	235
[2024]	236	! ! ===========
	237	DO jn = 1, kjpt ! tracer loop
	238	! ! ===========
	239	# if defined key_vectopt_loop
	240	DO jj = 1, 1
[2123]	241	DO ji = 1, jpij-jpi-1 ! vector opt. (forced unrolling)
[2024]	242	# else
[2123]	243	DO jj = 1, jpjm1
[2024]	244	DO ji = 1, jpim1 ! CAUTION start from i=1 to update i=2 when cyclic east-west
	245	# endif
	246	IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection
	247	! down-slope i/k-indices (deep) & up-slope i/k indices (shelf)
	248	iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) )
	249	ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj)
	250	zu_bbl = ABS( utr_bbl(ji,jj) )
	251	!
	252	! ! up -slope T-point (shelf bottom point)
	253	zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus)
[2148]	254	ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
	255	pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
[2024]	256	!
	257	DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column)
	258	zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk)
[2148]	259	ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
	260	pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
[2024]	261	END DO
	262	!
	263	zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud)
[2148]	264	ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
	265	pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
[2024]	266	ENDIF
[2104]	267	!
[2024]	268	IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection
	269	! down-slope j/k-indices (deep) & up-slope j/k indices (shelf)
	270	ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) )
	271	ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj)
	272	zv_bbl = ABS( vtr_bbl(ji,jj) )
	273	!
	274	! up -slope T-point (shelf bottom point)
	275	zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs)
[2148]	276	ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
	277	pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
[2024]	278	!
	279	DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column)
	280	zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk)
[2148]	281	ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
	282	pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra
[2024]	283	END DO
	284	! ! down-slope T-point (deep bottom point)
	285	zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd)
[2148]	286	ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
	287	pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
[2024]	288	ENDIF
	289	END DO
	290	!
	291	END DO
[2106]	292	! ! ===========
	293	END DO ! end tracer
	294	! ! ===========
[2024]	295	END SUBROUTINE tra_bbl_adv
	296
[2104]	297
[2024]	298	SUBROUTINE bbl( kt, cdtype )
	299	!!----------------------------------------------------------------------
	300	!! * ROUTINE bbl *
	301	!!
	302	!! ** Purpose : Computes the bottom boundary horizontal and vertical
	303	!! advection terms.
	304	!!
	305	!! ** Method :
	306	!! * diffusive bbl (nn_bbl_ldf=1) :
	307	!! When the product grad( rho) * grad(h) < 0 (where grad is an
	308	!! along bottom slope gradient) an additional lateral 2nd order
	309	!! diffusion along the bottom slope is added to the general
	310	!! tracer trend, otherwise the additional trend is set to 0.
	311	!! A typical value of ahbt is 2000 m2/s (equivalent to
	312	!! a downslope velocity of 20 cm/s if the condition for slope
	313	!! convection is satified)
	314	!! * advective bbl (nn_bbl_adv=1 or 2) :
	315	!! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity
	316	!! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation
	317	!! i.e. transport proportional to the along-slope density gradient
	318	!!
	319	!! NB: the along slope density gradient is evaluated using the
	320	!! local density (i.e. referenced at a common local depth).
	321	!!
	322	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
	323	!! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
	324	!!----------------------------------------------------------------------
[2104]	325	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	326	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	327	!!
[2024]	328	INTEGER :: ji, jj ! dummy loop indices
[2104]	329	INTEGER :: ik ! local integers
	330	INTEGER :: iis , iid , ijs , ijd ! - -
	331	INTEGER :: ikus, ikud, ikvs, ikvd ! - -
	332	REAL(wp) :: zsign, zsigna, zgbbl ! local scalars
	333	REAL(wp) :: zgdrho, zt, zs, zh ! - -
	334	REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb, ztb, zsb, zdep ! 2D workspace
[2024]	335	!!
	336	REAL(wp) :: fsalbt, fsbeta, pft, pfs, pfh ! statement function
	337	!!----------------------- zv_bbl-----------------------------------------------
	338	! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
	339	! ================ pft : potential temperature in degrees celcius
	340	! pfs : salinity anomaly (s-35) in psu
	341	! pfh : depth in meters
	342	! nn_eos = 0 (Jackett and McDougall 1994 formulation)
	343	fsalbt( pft, pfs, pfh ) = & ! alpha/beta
[3]	344	( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft &
[2024]	345	- 0.203814e-03 ) * pft &
	346	+ 0.170907e-01 ) * pft &
	347	+ 0.665157e-01 &
[3]	348	+(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs &
	349	+ ( ( - 0.302285e-13 * pfh &
[2024]	350	- 0.251520e-11 * pfs &
	351	+ 0.512857e-12 * pft * pft ) * pfh &
	352	- 0.164759e-06 * pfs &
	353	+( 0.791325e-08 * pft - 0.933746e-06 ) * pft &
	354	+ 0.380374e-04 ) * pfh
	355	fsbeta( pft, pfs, pfh ) = & ! beta
	356	( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft &
	357	- 0.301985e-05 ) * pft &
	358	+ 0.785567e-03 &
	359	+ ( 0.515032e-08 * pfs &
	360	+ 0.788212e-08 * pft - 0.356603e-06 ) * pfs &
	361	+( ( 0.121551e-17 * pfh &
	362	- 0.602281e-15 * pfs &
	363	- 0.175379e-14 * pft + 0.176621e-12 ) * pfh &
	364	+ 0.408195e-10 * pfs &
	365	+ ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft &
	366	- 0.121555e-07 ) * pfh
[3]	367	!!----------------------------------------------------------------------
[2024]	368
[2104]	369	IF( kt == nit000 ) THEN
	370	IF(lwp) WRITE(numout,*)
	371	IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
	372	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	373	ENDIF
	374
[2024]	375	! !* bottom temperature, salinity, velocity and depth
	376	#if defined key_vectopt_loop
	377	DO jj = 1, 1 ! vector opt. (forced unrolling)
	378	DO ji = 1, jpij
	379	#else
[3]	380	DO jj = 1, jpj
	381	DO ji = 1, jpi
[2024]	382	#endif
	383	ik = mbkt(ji,jj) ! bottom T-level index
[2148]	384	ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1) ! bottom before T and S
	385	zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
[2024]	386	zdep(ji,jj) = fsdept_0(ji,jj,ik) ! bottom T-level reference depth
	387	!
	388	zub(ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity
	389	zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
[3]	390	END DO
	391	END DO
[2024]	392
	393	! !-------------------!
	394	IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl !
	395	! !-------------------!
[2106]	396	DO jj = 1, jpjm1 ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
	397	DO ji = 1, jpim1
[2024]	398	! ! i-direction
	399	zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth
	400	zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
[1601]	401	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
[2024]	402	! ! masked bbl i-gradient of density
	403	zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) &
	404	& - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1)
	405	!
	406	zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope )
	407	ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff.
[1601]	408	!
[2024]	409	! ! j-direction
	410	zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth
	411	zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
[1601]	412	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
[2024]	413	! ! masked bbl j-gradient of density
	414	zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) &
	415	& - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1)
	416	!
	417	zsign = SIGN( 0.5, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope )
	418	ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
[1601]	419	!
	420	END DO
[3]	421	END DO
[1601]	422	!
[2024]	423	ENDIF
[409]	424
[2024]	425	! !-------------------!
	426	IF( nn_bbl_adv /= 0 ) THEN ! advective bbl !
	427	! !-------------------!
	428	SELECT CASE ( nn_bbl_adv ) !* bbl transport type
[503]	429	!
[2024]	430	CASE( 1 ) != use of upper velocity
[2123]	431	DO jj = 1, jpjm1 ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0
[2024]	432	DO ji = 1, fs_jpim1 ! vector opt.
	433	! ! i-direction
	434	zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth
	435	zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
	436	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
	437	! ! masked bbl i-gradient of density
	438	zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) &
	439	& - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1)
	440	!
	441	zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope
	442	zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope
	443	!
	444	! ! bbl velocity
	445	utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
	446	!
	447	! ! j-direction
	448	zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth
	449	zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
	450	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
	451	! ! masked bbl j-gradient of density
	452	zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) &
	453	& - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1)
	454	zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope
	455	zsigna= SIGN( 0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope
	456	!
	457	! ! bbl velocity
	458	vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
	459	END DO
	460	END DO
[2104]	461	!
[2024]	462	CASE( 2 ) != bbl velocity = F( delta rho )
	463	zgbbl = grav * rn_gambbl
[2123]	464	DO jj = 1, jpjm1 ! criteria: rho_up > rho_down
[2024]	465	DO ji = 1, fs_jpim1 ! vector opt.
	466	! ! i-direction
	467	! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf)
	468	iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) )
	469	ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj)
	470	!
	471	! ! mid-depth density anomalie (up-slope minus down-slope)
	472	zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! mid slope depth of T, S, and depth
	473	zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
	474	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
	475	zgdrho = fsbeta( zt, zs, zh ) &
	476	& * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) ) &
	477	& - ( zsb(iid,jj) - zsb(iis,jj) ) ) * umask(ji,jj,1)
	478	zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep
	479	!
	480	! ! bbl transport (down-slope direction)
	481	utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
	482	!
	483	! ! j-direction
	484	! down-slope T-point j/k-index (deep) & of the up -slope T-point j/k-index (shelf)
	485	ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) )
	486	ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj)
	487	!
	488	! ! mid-depth density anomalie (up-slope minus down-slope)
	489	zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) ) ! mid slope depth of T, S, and depth
	490	zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
	491	zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
	492	zgdrho = fsbeta( zt, zs, zh ) &
	493	& * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) ) &
	494	& - ( zsb(ji,ijd) - zsb(ji,ijs) ) ) * vmask(ji,jj,1)
	495	zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep
	496	!
	497	! ! bbl transport (down-slope direction)
	498	vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
	499	END DO
	500	END DO
[3]	501	END SELECT
[2024]	502	!
[3]	503	ENDIF
[503]	504	!
[2024]	505	END SUBROUTINE bbl
[3]	506
	507
	508	SUBROUTINE tra_bbl_init
	509	!!----------------------------------------------------------------------
	510	!! * ROUTINE tra_bbl_init *
	511	!!
	512	!! ** Purpose : Initialization for the bottom boundary layer scheme.
	513	!!
	514	!! ** Method : Read the nambbl namelist and check the parameters
[2104]	515	!! called by tra_bbl at the first timestep (nit000)
[3]	516	!!----------------------------------------------------------------------
[2024]	517	INTEGER :: ji, jj ! dummy loop indices
	518	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integer
	519	REAL(wp), DIMENSION(jpi,jpj) :: zmbk ! 2D workspace
[2104]	520	!!
[2024]	521	NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]	522	!!----------------------------------------------------------------------
	523
[2024]	524	REWIND ( numnam ) !* Read Namelist nambbl : bottom boundary layer scheme
[3]	525	READ ( numnam, nambbl )
[2106]	526	!
	527	l_bbl = .TRUE. !* flag to compute bbl coef and transport
	528	!
[2024]	529	IF(lwp) THEN !* Parameter control and print
[3]	530	WRITE(numout,*)
[2024]	531	WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]	532	WRITE(numout,*) '~~~~~~~~~~~~'
[503]	533	WRITE(numout,*) ' Namelist nambbl : set bbl parameters'
[2024]	534	WRITE(numout,*) ' diffusive bbl (=1) or not (=0) nn_bbl_ldf = ', nn_bbl_ldf
	535	WRITE(numout,*) ' advective bbl (=1/2) or not (=0) nn_bbl_adv = ', nn_bbl_adv
	536	WRITE(numout,*) ' diffusive bbl coefficient rn_ahtbbl = ', rn_ahtbbl, ' m2/s'
	537	WRITE(numout,*) ' advective bbl coefficient rn_gambbl = ', rn_gambbl, ' s'
[3]	538	ENDIF
[2024]	539
	540	IF( nn_bbl_adv == 1 ) WRITE(numout,) ' Advective BBL using upper velocity'
	541	IF( nn_bbl_adv == 2 ) WRITE(numout,) ' Advective BBL using velocity = F( delta rho)'
	542
	543	IF( nn_eos /= 0 ) CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
	544
	545
	546	! !* inverse of surface of T-cells
	547	e1e2t_r(:,:) = 1.0 / ( e1t(:,:) * e2t(:,:) )
	548
	549	! !* vertical index of bottom t-, u- and v-points
	550	DO jj = 1, jpj ! bottom k-index of T-level
[3]	551	DO ji = 1, jpi
[2024]	552	mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )
[3]	553	END DO
	554	END DO
[2024]	555	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level (shelf and deep)
[3]	556	DO ji = 1, jpim1
[2024]	557	mbku (ji,jj) = MAX( MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) - 1, 1 ) ! "shelf"
	558	mbkv (ji,jj) = MAX( MIN( mbathy(ji ,jj+1), mbathy(ji,jj) ) - 1, 1 )
	559	mbku_d(ji,jj) = MAX( MAX( mbathy(ji+1,jj ), mbathy(ji,jj) ) - 1, 1 ) ! "deep"
	560	mbkv_d(ji,jj) = MAX( MAX( mbathy(ji ,jj+1), mbathy(ji,jj) ) - 1, 1 )
[3]	561	END DO
	562	END DO
[2024]	563	zmbk(:,:) = FLOAT( mbku (:,:) ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	564	zmbk(:,:) = FLOAT( mbkv (:,:) ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	565	zmbk(:,:) = FLOAT( mbku_d(:,:) ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
	566	zmbk(:,:) = FLOAT( mbkv_d(:,:) ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
[3]	567
[2123]	568	!* sign of grad(H) at u- and v-points
	569	mgrhu(jpi,:) = 0. ; mgrhu(:,jpj) = 0. ; mgrhv(jpi,:) = 0. ; mgrhv(:,jpj) = 0.
	570	DO jj = 1, jpjm1
	571	DO ji = 1, jpim1
[2024]	572	mgrhu(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )
	573	mgrhv(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )
	574	END DO
	575	END DO
[481]	576
[2024]	577	DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point
	578	DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0)
	579	e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj )), fse3u_0(ji,jj,mbkt(ji,jj)) )
	580	e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )
	581	END DO
	582	END DO
	583	CALL lbc_lnk( e3u_bbl_0, 'U', 1. ) ; CALL lbc_lnk( e3v_bbl_0, 'V', 1. ) ! lateral boundary conditions
[481]	584
[2024]	585	! !* masked diffusive flux coefficients
	586	ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:) * umask(:,:,1)
	587	ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:) * vmask(:,:,1)
	588
[2148]	589
[2024]	590	IF( cp_cfg == "orca" ) THEN !* ORCA configuration : regional enhancement of ah_bbl
	591	!
	592	SELECT CASE ( jp_cfg )
	593	CASE ( 2 ) ! ORCA_R2
	594	ij0 = 102 ; ij1 = 102 ! Gibraltar enhancement of BBL
	595	ii0 = 139 ; ii1 = 140
	596	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	597	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	598	!
	599	ij0 = 88 ; ij1 = 88 ! Red Sea enhancement of BBL
	600	ii0 = 161 ; ii1 = 162
	601	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	602	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	603	!
	604	CASE ( 4 ) ! ORCA_R4
	605	ij0 = 52 ; ij1 = 52 ! Gibraltar enhancement of BBL
	606	ii0 = 70 ; ii1 = 71
	607	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	608	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	609	END SELECT
	610	!
	611	ENDIF
[503]	612	!
[3]	613	END SUBROUTINE tra_bbl_init
	614
	615	#else
	616	!!----------------------------------------------------------------------
	617	!! Dummy module : No bottom boundary layer scheme
	618	!!----------------------------------------------------------------------
[2024]	619	LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .FALSE. !: bbl flag
[3]	620	CONTAINS
[2104]	621	SUBROUTINE tra_bbl_init ! Dummy routine
	622	END SUBROUTINE tra_bbl_init
	623	SUBROUTINE tra_bbl( kt ) ! Dummy routine
[2024]	624	WRITE(,) 'tra_bbl: You should not have seen this print! error?', kt
	625	END SUBROUTINE tra_bbl
[3]	626	#endif
	627
	628	!!======================================================================
	629	END MODULE trabbl

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

Download in other formats: