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

Context Navigation

source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 13320

Last change on this file since 13320 was 13191, checked in by jwhile, 4 years ago
Updates for 1d runnig
File size: 36.7 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	!
[13191]	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
[13191]	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(:,:)
[13191]	201	#if defined key_traldf_c2d \|\| key_traldf_c3d
[11442]	202	IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
	203	CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
	204	ENDIF
[13191]	205	#else
	206	IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) &
	207	& CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
	208	'key_traldf_c2d or key_traldf_c3d')
	209	#endif
	210
	211
[11442]	212	ahv_bbl_1(:,:) = ahv_bbl(:,:)
[13191]	213	#if defined key_traldf_c2d \|\| key_traldf_c3d
[11442]	214	IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
	215	CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
	216	ENDIF
[13191]	217	#else
	218	IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) &
	219	& CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
	220	'key_traldf_c2d or key_traldf_c3d')
	221	#endif
	222
[11442]	223	!
[2528]	224	DO jn = 1, kjpt ! tracer loop
	225	! ! ===========
	226	DO jj = 1, jpj
[3764]	227	DO ji = 1, jpi
[4990]	228	ik = mbkt(ji,jj) ! bottom T-level index
	229	zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S
[2528]	230	END DO
	231	END DO
[13191]	232	!
[4990]	233	DO jj = 2, jpjm1 ! Compute the trend
[2528]	234	DO ji = 2, jpim1
[4990]	235	ik = mbkt(ji,jj) ! bottom T-level index
[4292]	236	zbtr = r1_e12t(ji,jj) / fse3t(ji,jj,ik)
[2528]	237	pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) &
[11442]	238	& + ( ahu_bbl_1(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) &
	239	& - ahu_bbl_1(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) &
	240	& + ahv_bbl_1(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) &
	241	& - ahv_bbl_1(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) * zbtr
[2528]	242	END DO
[3]	243	END DO
[2528]	244	! ! ===========
	245	END DO ! end tracer
	246	! ! ===========
[7771]	247	DEALLOCATE( zptb )
[2715]	248	!
[3294]	249	IF( nn_timing == 1 ) CALL timing_stop('tra_bbl_dif')
	250	!
[2528]	251	END SUBROUTINE tra_bbl_dif
[3]	252
[3764]	253
[2528]	254	SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
	255	!!----------------------------------------------------------------------
	256	!! * ROUTINE trc_bbl *
	257	!!
[3764]	258	!! ** Purpose : Compute the before passive tracer trend associated
[2528]	259	!! with the bottom boundary layer and add it to the general trend
	260	!! of tracer equations.
	261	!! ** Method : advective bbl (nn_bbl_adv = 1 or 2) :
	262	!! nn_bbl_adv = 1 use of the ocean near bottom velocity as bbl velocity
[3764]	263	!! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e.
	264	!! transport proportional to the along-slope density gradient
[2528]	265	!!
	266	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
	267	!! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]	268	!!----------------------------------------------------------------------
[2715]	269	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	270	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
[3764]	271	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	272	!
[2528]	273	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	274	INTEGER :: iis , iid , ijs , ijd ! local integers
	275	INTEGER :: ikus, ikud, ikvs, ikvd ! - -
	276	REAL(wp) :: zbtr, ztra ! local scalars
	277	REAL(wp) :: zu_bbl, zv_bbl ! - -
	278	!!----------------------------------------------------------------------
	279	!
[3294]	280	IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_adv')
[2528]	281	! ! ===========
	282	DO jn = 1, kjpt ! tracer loop
[3764]	283	! ! ===========
[457]	284	DO jj = 1, jpjm1
[2528]	285	DO ji = 1, jpim1 ! CAUTION start from i=1 to update i=2 when cyclic east-west
	286	IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection
	287	! down-slope i/k-indices (deep) & up-slope i/k indices (shelf)
	288	iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) )
	289	ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj)
	290	zu_bbl = ABS( utr_bbl(ji,jj) )
	291	!
	292	! ! up -slope T-point (shelf bottom point)
[4292]	293	zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
[2528]	294	ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
	295	pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
[3764]	296	!
[2528]	297	DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column)
[4292]	298	zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
[2528]	299	ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
	300	pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
	301	END DO
[3764]	302	!
[4292]	303	zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
[2528]	304	ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
	305	pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
	306	ENDIF
	307	!
	308	IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection
	309	! down-slope j/k-indices (deep) & up-slope j/k indices (shelf)
	310	ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) )
	311	ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj)
	312	zv_bbl = ABS( vtr_bbl(ji,jj) )
[3764]	313	!
[2528]	314	! up -slope T-point (shelf bottom point)
[4292]	315	zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
[2528]	316	ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
	317	pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
[3764]	318	!
[2528]	319	DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column)
[4292]	320	zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
[2528]	321	ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
	322	pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra
	323	END DO
	324	! ! down-slope T-point (deep bottom point)
[4292]	325	zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
[2528]	326	ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
	327	pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
	328	ENDIF
[457]	329	END DO
[2528]	330	!
[457]	331	END DO
[2528]	332	! ! ===========
	333	END DO ! end tracer
	334	! ! ===========
[3294]	335	!
	336	IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_adv')
	337	!
[2528]	338	END SUBROUTINE tra_bbl_adv
[3]	339
	340
[3294]	341	SUBROUTINE bbl( kt, kit000, cdtype )
[2528]	342	!!----------------------------------------------------------------------
	343	!! * ROUTINE bbl *
[3764]	344	!!
[2528]	345	!! ** Purpose : Computes the bottom boundary horizontal and vertical
[3764]	346	!! advection terms.
[2528]	347	!!
[4990]	348	!! ** Method : * diffusive bbl (nn_bbl_ldf=1) :
[2528]	349	!! When the product grad( rho) * grad(h) < 0 (where grad is an
	350	!! along bottom slope gradient) an additional lateral 2nd order
	351	!! diffusion along the bottom slope is added to the general
	352	!! tracer trend, otherwise the additional trend is set to 0.
	353	!! A typical value of ahbt is 2000 m2/s (equivalent to
	354	!! a downslope velocity of 20 cm/s if the condition for slope
	355	!! convection is satified)
[4990]	356	!! * advective bbl (nn_bbl_adv=1 or 2) :
[2528]	357	!! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity
	358	!! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation
	359	!! i.e. transport proportional to the along-slope density gradient
	360	!!
	361	!! NB: the along slope density gradient is evaluated using the
	362	!! local density (i.e. referenced at a common local depth).
	363	!!
	364	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
	365	!! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
[3764]	366	!!----------------------------------------------------------------------
[2528]	367	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[4990]	368	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	369	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	370	!!
	371	INTEGER :: ji, jj ! dummy loop indices
	372	INTEGER :: ik ! local integers
[4990]	373	INTEGER :: iis, iid, ikus, ikud ! - -
	374	INTEGER :: ijs, ijd, ikvs, ikvd ! - -
	375	REAL(wp) :: za, zb, zgdrho ! local scalars
	376	REAL(wp) :: zsign, zsigna, zgbbl ! - -
	377	REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts, zab ! 3D workspace
	378	REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb, zdep ! 2D workspace
[2528]	379	!!----------------------------------------------------------------------
[3294]	380	!
	381	IF( nn_timing == 1 ) CALL timing_start( 'bbl')
	382	!
	383	IF( kt == kit000 ) THEN
[2528]	384	IF(lwp) WRITE(numout,*)
	385	IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
	386	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	387	ENDIF
[4990]	388	! !* bottom variables (T, S, alpha, beta, depth, velocity)
[2528]	389	DO jj = 1, jpj
	390	DO ji = 1, jpi
[4990]	391	ik = mbkt(ji,jj) ! bottom T-level index
	392	zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem) ! bottom before T and S
	393	zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
[2528]	394	!
[4990]	395	zdep(ji,jj) = fsdept(ji,jj,ik) ! bottom T-level reference depth
	396	zub (ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity
	397	zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
[2528]	398	END DO
	399	END DO
[4990]	400	!
	401	CALL eos_rab( zts, zdep, zab )
	402	!
[2528]	403	! !-------------------!
	404	IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl !
	405	! !-------------------!
	406	DO jj = 1, jpjm1 ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
[4990]	407	DO ji = 1, fs_jpim1 ! vector opt.
	408	! ! i-direction
	409	za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point
	410	zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
	411	! ! 2*masked bottom density gradient
	412	zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) &
	413	& - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1)
[3764]	414	!
[4990]	415	zsign = SIGN( 0.5, -zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope )
	416	ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff.
[1601]	417	!
[4990]	418	! ! j-direction
	419	za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point
	420	zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
	421	! ! 2*masked bottom density gradient
	422	zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) &
	423	& - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1)
[3764]	424	!
[4990]	425	zsign = SIGN( 0.5, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope )
[2528]	426	ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
[1601]	427	END DO
[3]	428	END DO
[1601]	429	!
[2528]	430	ENDIF
[409]	431
[2528]	432	! !-------------------!
	433	IF( nn_bbl_adv /= 0 ) THEN ! advective bbl !
	434	! !-------------------!
	435	SELECT CASE ( nn_bbl_adv ) !* bbl transport type
[503]	436	!
[2528]	437	CASE( 1 ) != use of upper velocity
	438	DO jj = 1, jpjm1 ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0
	439	DO ji = 1, fs_jpim1 ! vector opt.
[4990]	440	! ! i-direction
	441	za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point
	442	zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
[13191]	443	! ! 2*masked bottom density gradient
[4990]	444	zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) &
	445	- zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1)
[3764]	446	!
[4990]	447	zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope
	448	zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope
[2528]	449	!
[4990]	450	! ! bbl velocity
[2528]	451	utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
	452	!
[4990]	453	! ! j-direction
	454	za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point
	455	zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
	456	! ! 2*masked bottom density gradient
	457	zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) &
	458	& - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1)
	459	zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope
	460	zsigna= SIGN( 0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope
[2528]	461	!
[4990]	462	! ! bbl transport
[2528]	463	vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
	464	END DO
	465	END DO
[503]	466	!
[2528]	467	CASE( 2 ) != bbl velocity = F( delta rho )
	468	zgbbl = grav * rn_gambbl
	469	DO jj = 1, jpjm1 ! criteria: rho_up > rho_down
	470	DO ji = 1, fs_jpim1 ! vector opt.
[4990]	471	! ! i-direction
[2528]	472	! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf)
[4990]	473	iid = ji + MAX( 0, mgrhu(ji,jj) )
	474	iis = ji + 1 - MAX( 0, mgrhu(ji,jj) )
[2528]	475	!
[4990]	476	ikud = mbku_d(ji,jj)
	477	ikus = mbku(ji,jj)
[2528]	478	!
[4990]	479	za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point
	480	zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
	481	! ! masked bottom density gradient
	482	zgdrho = 0.5 * ( za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) ) &
	483	& - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) ) ) * umask(ji,jj,1)
	484	zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep
	485	!
	486	! ! bbl transport (down-slope direction)
[2528]	487	utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
	488	!
[4990]	489	! ! j-direction
[2528]	490	! down-slope T-point j/k-index (deep) & of the up -slope T-point j/k-index (shelf)
[4990]	491	ijd = jj + MAX( 0, mgrhv(ji,jj) )
	492	ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) )
[2528]	493	!
[4990]	494	ikvd = mbkv_d(ji,jj)
	495	ikvs = mbkv(ji,jj)
[2528]	496	!
[4990]	497	za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point
	498	zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
	499	! ! masked bottom density gradient
	500	zgdrho = 0.5 * ( za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) ) &
	501	& - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) ) ) * vmask(ji,jj,1)
	502	zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep
	503	!
	504	! ! bbl transport (down-slope direction)
[2528]	505	vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
	506	END DO
	507	END DO
[3]	508	END SELECT
[2528]	509	!
[3]	510	ENDIF
[503]	511	!
[3294]	512	IF( nn_timing == 1 ) CALL timing_stop( 'bbl')
	513	!
[2528]	514	END SUBROUTINE bbl
[3]	515
	516
	517	SUBROUTINE tra_bbl_init
	518	!!----------------------------------------------------------------------
	519	!! * ROUTINE tra_bbl_init *
	520	!!
	521	!! ** Purpose : Initialization for the bottom boundary layer scheme.
	522	!!
	523	!! ** Method : Read the nambbl namelist and check the parameters
[3294]	524	!! called by nemo_init at the first timestep (kit000)
[3]	525	!!----------------------------------------------------------------------
[2528]	526	INTEGER :: ji, jj ! dummy loop indices
[4990]	527	INTEGER :: ii0, ii1, ij0, ij1 ! local integer
	528	INTEGER :: ios ! - -
[3294]	529	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
[2528]	530	!!
	531	NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
[3]	532	!!----------------------------------------------------------------------
[3294]	533	!
	534	IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_init')
	535	!
[3764]	536	CALL wrk_alloc( jpi, jpj, zmbk )
[3294]	537	!
[3]	538
[4147]	539	REWIND( numnam_ref ) ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
	540	READ ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
	541	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
	542
	543	REWIND( numnam_cfg ) ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
	544	READ ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
	545	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
[4624]	546	IF(lwm) WRITE ( numond, nambbl )
[2528]	547	!
	548	l_bbl = .TRUE. !* flag to compute bbl coef and transport
	549	!
	550	IF(lwp) THEN !* Parameter control and print
[3]	551	WRITE(numout,*)
[2528]	552	WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
[3]	553	WRITE(numout,*) '~~~~~~~~~~~~'
[503]	554	WRITE(numout,*) ' Namelist nambbl : set bbl parameters'
[2528]	555	WRITE(numout,*) ' diffusive bbl (=1) or not (=0) nn_bbl_ldf = ', nn_bbl_ldf
	556	WRITE(numout,*) ' advective bbl (=1/2) or not (=0) nn_bbl_adv = ', nn_bbl_adv
	557	WRITE(numout,*) ' diffusive bbl coefficient rn_ahtbbl = ', rn_ahtbbl, ' m2/s'
	558	WRITE(numout,*) ' advective bbl coefficient rn_gambbl = ', rn_gambbl, ' s'
[3]	559	ENDIF
[2715]	560
	561	! ! allocate trabbl arrays
	562	IF( tra_bbl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
[3764]	563
[2528]	564	IF( nn_bbl_adv == 1 ) WRITE(numout,) ' Advective BBL using upper velocity'
	565	IF( nn_bbl_adv == 2 ) WRITE(numout,) ' Advective BBL using velocity = F( delta rho)'
	566
	567	! !* vertical index of "deep" bottom u- and v-points
	568	DO jj = 1, jpjm1 ! (the "shelf" bottom k-indices are mbku and mbkv)
	569	DO ji = 1, jpim1
	570	mbku_d(ji,jj) = MAX( mbkt(ji+1,jj ) , mbkt(ji,jj) ) ! >= 1 as mbkt=1 over land
	571	mbkv_d(ji,jj) = MAX( mbkt(ji ,jj+1) , mbkt(ji,jj) )
[3]	572	END DO
	573	END DO
[4990]	574	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[2528]	575	zmbk(:,:) = REAL( mbku_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
	576	zmbk(:,:) = REAL( mbkv_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
	577
[10302]	578	!! AXY (16/08/17): remove the following per George and Andrew bug-hunt
	579	!! !* sign of grad(H) at u- and v-points
	580	!! mgrhu(jpi,:) = 0 ; mgrhu(:,jpj) = 0 ; mgrhv(jpi,:) = 0 ; mgrhv(:,jpj) = 0
	581	!! DO jj = 1, jpjm1
	582	!! DO ji = 1, jpim1
	583	!! mgrhu(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	584	!! mgrhv(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	585	!! END DO
	586	!! END DO
	587
	588	!! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
	589	!* sign of grad(H) at u- and v-points; zero if grad(H) = 0
	590	mgrhu(:,:) = 0 ; mgrhv(:,:) = 0
[3764]	591	DO jj = 1, jpjm1
[3]	592	DO ji = 1, jpim1
[10302]	593	#if defined key_bbl_old_nonconserve
	594	! This key allows old (non conservative version) to be used for continuity of results
	595	mgrhu(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	596	mgrhv(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	597	#else
	598	IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
	599	mgrhu(ji,jj) = INT( SIGN( 1.e0, &
	600	gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	601	ENDIF
[13191]	602	!
[10302]	603	IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
	604	mgrhv(ji,jj) = INT( SIGN( 1.e0, &
	605	gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) )
	606	ENDIF
	607	#endif
[3]	608	END DO
	609	END DO
	610
[3764]	611	DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point
[4990]	612	DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0)
[4292]	613	e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj )), e3u_0(ji,jj,mbkt(ji,jj)) )
	614	e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
[3764]	615	END DO
[2528]	616	END DO
	617	CALL lbc_lnk( e3u_bbl_0, 'U', 1. ) ; CALL lbc_lnk( e3v_bbl_0, 'V', 1. ) ! lateral boundary conditions
[481]	618
[3764]	619	! !* masked diffusive flux coefficients
[2528]	620	ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:) * umask(:,:,1)
	621	ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:) * vmask(:,:,1)
[481]	622
[2528]	623	IF( cp_cfg == "orca" ) THEN !* ORCA configuration : regional enhancement of ah_bbl
	624	!
	625	SELECT CASE ( jp_cfg )
	626	CASE ( 2 ) ! ORCA_R2
	627	ij0 = 102 ; ij1 = 102 ! Gibraltar enhancement of BBL
[3764]	628	ii0 = 139 ; ii1 = 140
[2528]	629	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	630	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	631	!
	632	ij0 = 88 ; ij1 = 88 ! Red Sea enhancement of BBL
	633	ii0 = 161 ; ii1 = 162
	634	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	635	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	636	!
	637	CASE ( 4 ) ! ORCA_R4
	638	ij0 = 52 ; ij1 = 52 ! Gibraltar enhancement of BBL
[3764]	639	ii0 = 70 ; ii1 = 71
[2528]	640	ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	641	ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
	642	END SELECT
	643	!
	644	ENDIF
[503]	645	!
[3764]	646	CALL wrk_dealloc( jpi, jpj, zmbk )
[2715]	647	!
[3294]	648	IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_init')
	649	!
[3]	650	END SUBROUTINE tra_bbl_init
	651
	652	#else
	653	!!----------------------------------------------------------------------
	654	!! Dummy module : No bottom boundary layer scheme
	655	!!----------------------------------------------------------------------
[2528]	656	LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .FALSE. !: bbl flag
[3]	657	CONTAINS
[2528]	658	SUBROUTINE tra_bbl_init ! Dummy routine
	659	END SUBROUTINE tra_bbl_init
	660	SUBROUTINE tra_bbl( kt ) ! Dummy routine
	661	WRITE(,) 'tra_bbl: You should not have seen this print! error?', kt
	662	END SUBROUTINE tra_bbl
[3]	663	#endif
	664
	665	!!======================================================================
	666	END MODULE trabbl

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

Download in other formats: