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_adv.h90 in branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 2208

Last change on this file since 2208 was 2208, checked in by rblod, 14 years ago
Put FCM NEMO code changes in DEV_r2191_3partymerge2010 branch
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 19.2 KB

Rev	Line
[3]	1	!!----------------------------------------------------------------------
	2	!! * trabbl_adv.h90 *
	3	!!----------------------------------------------------------------------
[503]	4	!! History : 8.5 ! 02-12 (A. de Miranda, G. Madec) Original Code
	5	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines )
	6	!!----------------------------------------------------------------------
[3]	7
	8	!!----------------------------------------------------------------------
[247]	9	!! OPA 9.0 , LOCEAN-IPSL (2005)
[1152]	10	!! $Id$
[2208]	11	!! Software governed by the CeCILL licence (NEMOGCM/License_CeCILL.txt)
[3]	12	!!----------------------------------------------------------------------
	13
	14	SUBROUTINE tra_bbl_adv( kt )
	15	!!----------------------------------------------------------------------
	16	!! * ROUTINE tra_bbl_adv *
	17	!!
	18	!! ** Purpose : Compute the before tracer (t & s) trend associated
	19	!! with the bottom boundary layer and add it to the general trend
	20	!! of tracer equations. The bottom boundary layer is supposed to be
	21	!! both an advective and diffusive bottom boundary layer.
	22	!!
	23	!! ** Method : Computes the bottom boundary horizontal and vertical
	24	!! advection terms. Add it to the general trend : ta =ta + adv_bbl.
	25	!! When the product grad( rho) * grad(h) < 0 (where grad is a
	26	!! along bottom slope gradient) an additional lateral 2nd order
	27	!! diffusion along the bottom slope is added to the general
	28	!! tracer trend, otherwise the additional trend is set to 0.
	29	!! Second order operator (laplacian type) with variable coefficient
	30	!! computed as follow for temperature (idem on s):
	31	!! difft = 1/(e1te2te3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
	32	!! + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
	33	!! where ztb is a 2D array: the bottom ocean te;perature and ahtb
	34	!! is a time and space varying diffusive coefficient defined by:
	35	!! ahbt = zahbp if grad(rho).grad(h) < 0
	36	!! = 0. otherwise.
	37	!! Note that grad(.) is the along bottom slope gradient. grad(rho)
	38	!! is evaluated using the local density (i.e. referenced at the
	39	!! local depth). Typical value of ahbt is 2000 m2/s (equivalent to
	40	!! a downslope velocity of 20 cm/s if the condition for slope
	41	!! convection is satified)
	42	!! Add this before trend to the general trend (ta,sa) of the
	43	!! botton ocean tracer point:
	44	!! ta = ta + difft
	45	!!
	46	!! ** Action : - update (ta,sa) at the bottom level with the bottom
	47	!! boundary layer trend
	48	!!
[503]	49	!! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
	50	!!----------------------------------------------------------------------
	51	USE eosbn2 ! equation of state
	52	USE oce, ONLY : ztrdt => ua ! use ua as 3D workspace
[1482]	53	USE oce, ONLY : ztrds => va ! use va as 3D workspace
	54	USE iom
[3]	55	!!
[503]	56	INTEGER, INTENT( in ) :: kt ! ocean time-step
	57	!!
	58	INTEGER :: ji, jj, jk ! dummy loop indices
	59	INTEGER :: ik ! temporary integers
	60	INTEGER :: iku, iku1, iku2 ! " "
	61	INTEGER :: ikv, ikv1, ikv2 ! " "
	62	REAL(wp) :: zsign, zh, zalbet ! temporary scalars
	63	REAL(wp) :: zgdrho, zbtr ! " "
	64	REAL(wp) :: zbt, zsigna ! " "
	65	REAL(wp) :: zfui, ze3u, zt, zta ! " "
	66	REAL(wp) :: zfvj, ze3v, zs, zsa ! " "
	67	REAL(wp), DIMENSION(jpi,jpj) :: zalphax, zwu, zunb ! temporary 2D workspace
	68	REAL(wp), DIMENSION(jpi,jpj) :: zalphay, zwv, zvnb ! " "
	69	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zww, zwz ! " "
	70	REAL(wp), DIMENSION(jpi,jpj) :: ztbb, zsbb ! " "
	71	REAL(wp), DIMENSION(jpi,jpj) :: ztnb, zsnb, zdep ! " "
	72	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhdivn ! temporary 3D workspace
	73	REAL(wp) :: fsalbt, pft, pfs, pfh ! statement function
[3]	74	!!----------------------------------------------------------------------
	75	! ratio alpha/beta
	76	! ================
	77	! fsalbt: ratio of thermal over saline expension coefficients
	78	! pft : potential temperature in degrees celcius
	79	! pfs : salinity anomaly (s-35) in psu
	80	! pfh : depth in meters
	81
	82	fsalbt( pft, pfs, pfh ) = &
	83	( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft &
	84	- 0.203814e-03 ) * pft &
	85	+ 0.170907e-01 ) * pft &
	86	+ 0.665157e-01 &
	87	+(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs &
	88	+ ( ( - 0.302285e-13 * pfh &
	89	- 0.251520e-11 * pfs &
	90	+ 0.512857e-12 * pft * pft ) * pfh &
	91	- 0.164759e-06 * pfs &
	92	+( 0.791325e-08 * pft - 0.933746e-06 ) * pft &
	93	+ 0.380374e-04 ) * pfh
	94	!!----------------------------------------------------------------------
	95
	96	IF( kt == nit000 ) CALL tra_bbl_init ! initialization at first time-step
	97
	98	! 1. 2D fields of bottom temperature and salinity, and bottom slope
	99	! -----------------------------------------------------------------
	100	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
	101
[789]	102	#if defined key_vectopt_loop
[1601]	103	DO jj = 1, 1
	104	DO ji = 1, jpij ! vector opt. (forced unrolling)
[3]	105	#else
	106	DO jj = 1, jpj
	107	DO ji = 1, jpi
	108	#endif
	109	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
[481]	110	ztnb(ji,jj) = tn(ji,jj,ik)
	111	zsnb(ji,jj) = sn(ji,jj,ik)
	112	ztbb(ji,jj) = tb(ji,jj,ik)
	113	zsbb(ji,jj) = sb(ji,jj,ik)
[3]	114	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
[481]	115
	116	zunb(ji,jj) = un(ji,jj,mbku(ji,jj))
	117	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
[3]	118	END DO
	119	END DO
	120
[481]	121
[3]	122	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
	123	! --------------------------------------------
	124	! Sign of the local density gradient along the i- and j-slopes
	125	! multiplied by the slope of the ocean bottom
	126
[1601]	127	SELECT CASE ( nn_eos )
[503]	128	!
[3]	129	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
[1601]	130	!
	131	DO jj = 1, jpjm1
	132	DO ji = 1, fs_jpim1 ! vector opt.
	133	! ... temperature, salinity anomalie and depth
	134	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
	135	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
	136	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
	137	! ... masked ratio alpha/beta
	138	zalbet = fsalbt( zt, zs, zh ) * umask(ji,jj,1)
	139	! ... local density gradient along i-bathymetric slope
	140	zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
	141	& - ( zsnb(ji+1,jj) - zsnb(ji,jj) )
	142	zgdrho = zgdrho * umask(ji,jj,1)
	143	! ... sign of local i-gradient of density multiplied by the i-slope
	144	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	145	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	146	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
	147	!
	148	! ... temperature, salinity anomalie and depth
	149	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
	150	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
	151	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
	152	! ... masked ratio alpha/beta
	153	zalbet = fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
	154	! ... local density gradient along j-bathymetric slope
	155	zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
	156	& - ( zsnb(ji,jj+1) - zsnb(ji,jj) )
	157	zgdrho = zgdrho*vmask(ji,jj,1)
	158	! ... sign of local j-gradient of density multiplied by the j-slope
	159	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	160	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	161	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
	162	END DO
[503]	163	END DO
[1601]	164	!
[3]	165	CASE ( 1 ) ! Linear formulation function of temperature only
[1601]	166	!
	167	DO jj = 1, jpjm1
	168	DO ji = 1, fs_jpim1 ! vector opt.
	169	! local 'density/temperature' gradient along i-bathymetric slope
	170	zgdrho = ( ztnb(ji+1,jj) - ztnb(ji,jj) )
	171	! sign of local i-gradient of density multiplied by the i-slope
	172	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	173	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	174	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
	175	!
	176	! local density gradient along j-bathymetric slope
	177	zgdrho = ( ztnb(ji,jj+1) - ztnb(ji,jj) )
	178	! sign of local j-gradient of density multiplied by the j-slope
	179	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	180	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	181	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
	182	END DO
[409]	183	END DO
[1601]	184	!
[481]	185	CASE ( 2 ) ! Linear formulation function of temperature and salinity
[1601]	186	!
	187	DO jj = 1, jpjm1
	188	DO ji = 1, fs_jpim1 ! vector opt.
	189	! local density gradient along i-bathymetric slope
	190	zgdrho = - ( rn_beta *( zsnb(ji+1,jj) - zsnb(ji,jj) ) &
	191	& - rn_alpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
	192	! sign of local i-gradient of density multiplied by the i-slope
	193	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	194	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	195	zalphax(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
	196	!
	197	! local density gradient along j-bathymetric slope
	198	zgdrho = - ( rn_beta *( zsnb(ji,jj+1) - zsnb(ji,jj) ) &
	199	& - rn_alpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )
	200	! sign of local j-gradient of density multiplied by the j-slope
	201	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	202	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	203	zalphay(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
	204	END DO
[409]	205	END DO
[503]	206	!
[3]	207	END SELECT
	208
	209	! lateral boundary conditions on zalphax and zalphay (unchanged sign)
	210	CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. )
	211
[21]	212
[3]	213	! 3. Velocities that are exchanged between ajacent bottom boxes.
	214	!---------------------------------------------------------------
	215
	216	! ... is equal to zero but where bbl will work.
[481]	217	u_bbl(:,:,:) = 0.e0
	218	v_bbl(:,:,:) = 0.e0
	219
	220	IF( ln_zps ) THEN ! partial steps correction
	221
[789]	222	# if defined key_vectopt_loop
[1601]	223	DO jj = 1, 1
	224	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[481]	225	# else
	226	DO jj = 1, jpjm1
	227	DO ji = 1, jpim1
[3]	228	# endif
[481]	229	iku = mbku(ji ,jj )
	230	ikv = mbkv(ji ,jj )
	231	iku1 = mbkt(ji+1,jj )
	232	iku2 = mbkt(ji ,jj )
	233	ikv1 = mbkt(ji ,jj+1)
	234	ikv2 = mbkt(ji ,jj )
	235	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
	236	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
	237
	238	IF( MAX(iku,ikv) > 1 ) THEN
	239	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * ze3u / fse3u(ji,jj,iku)
	240	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * ze3v / fse3v(ji,jj,ikv)
	241	ENDIF
	242	END DO
[3]	243	END DO
	244
[481]	245	! lateral boundary conditions on u_bbl and v_bbl (changed sign)
	246	CALL lbc_lnk( u_bbl, 'U', -1. ) ; CALL lbc_lnk( v_bbl, 'V', -1. )
	247
	248	ELSE ! if not partial step loop over the whole domain no lbc call
[3]	249
[789]	250	#if defined key_vectopt_loop
[1601]	251	DO jj = 1, 1
	252	DO ji = 1, jpij ! vector opt. (forced unrolling)
[481]	253	#else
	254	DO jj = 1, jpj
	255	DO ji = 1, jpi
	256	#endif
	257	iku = mbku(ji,jj)
	258	ikv = mbkv(ji,jj)
	259	IF( MAX(iku,ikv) > 1 ) THEN
	260	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku)
	261	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv)
	262	ENDIF
	263	END DO
	264	END DO
	265
	266	ENDIF
	267
[1601]	268
[3]	269	! 5. Along sigma advective trend
	270	! -------------------------------
	271	! ... Second order centered tracer flux at u and v-points
	272
[789]	273	# if defined key_vectopt_loop
[1601]	274	DO jj = 1, 1
	275	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	276	# else
	277	DO jj = 1, jpjm1
	278	DO ji = 1, jpim1
	279	# endif
	280	iku = mbku(ji,jj)
	281	ikv = mbkv(ji,jj)
[481]	282	zfui = e2u(ji,jj) * fse3u(ji,jj,iku) * u_bbl(ji,jj,iku)
	283	zfvj = e1v(ji,jj) * fse3v(ji,jj,ikv) * v_bbl(ji,jj,ikv)
[3]	284	! centered scheme
	285	! zwx(ji,jj) = 0.5* zfui * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
	286	! zwy(ji,jj) = 0.5* zfvj * ( ztnb(ji,jj) + ztnb(ji,jj+1) )
	287	! zww(ji,jj) = 0.5* zfui * ( zsnb(ji,jj) + zsnb(ji+1,jj) )
	288	! zwz(ji,jj) = 0.5* zfvj * ( zsnb(ji,jj) + zsnb(ji,jj+1) )
	289	! upstream scheme
[21]	290	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji ,jj ) &
	291	& +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj ) ) * 0.5
[409]	292	zwy(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * ztbb(ji ,jj ) &
	293	& +( zfvj - ABS( zfvj ) ) * ztbb(ji ,jj+1) ) * 0.5
[21]	294	zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji ,jj ) &
	295	& +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj ) ) * 0.5
[409]	296	zwz(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * zsbb(ji ,jj ) &
	297	& +( zfvj - ABS( zfvj ) ) * zsbb(ji ,jj+1) ) * 0.5
[3]	298	END DO
[1601]	299	END DO
[789]	300	# if defined key_vectopt_loop
[1601]	301	DO jj = 1, 1
	302	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	303	# else
	304	DO jj = 2, jpjm1
	305	DO ji = 2, jpim1
	306	# endif
	307	ik = mbkt(ji,jj)
	308	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
[216]	309	! horizontal advective trends
[3]	310	zta = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
	311	& + zwy(ji,jj) - zwy(ji ,jj-1) )
	312	zsa = - zbtr * ( zww(ji,jj) - zww(ji-1,jj ) &
	313	& + zwz(ji,jj) - zwz(ji ,jj-1) )
	314
[216]	315	! add it to the general tracer trends
[3]	316	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
	317	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
	318	END DO
	319	END DO
	320
[503]	321	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' bbl - Ta: ', mask1=tmask, &
	322	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[3]	323
	324	! 6. Vertical advection velocities
	325	! --------------------------------
	326	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
	327	DO jk= 1, jpkm1
	328	DO jj=1, jpjm1
[409]	329	DO ji = 1, fs_jpim1 ! vector opt.
	330	zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk) * fse3u(ji,jj,jk)
	331	zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk) * fse3v(ji,jj,jk)
[3]	332	END DO
	333	END DO
	334
	335	! ... horizontal divergence
	336	DO jj = 2, jpjm1
	337	DO ji = fs_2, fs_jpim1 ! vector opt.
[409]	338	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
[3]	339	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
	340	+ zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
	341	END DO
	342	END DO
	343	END DO
	344
	345
	346	! ... horizontal bottom divergence
[481]	347
	348	IF( ln_zps ) THEN
[1601]	349
[789]	350	# if defined key_vectopt_loop
[1601]	351	DO jj = 1, 1
	352	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	353	# else
[481]	354	DO jj = 1, jpjm1
	355	DO ji = 1, jpim1
[3]	356	# endif
[481]	357	iku = mbku(ji ,jj )
	358	ikv = mbkv(ji ,jj )
	359	iku1 = mbkt(ji+1,jj )
	360	iku2 = mbkt(ji ,jj )
	361	ikv1 = mbkt(ji ,jj+1)
	362	ikv2 = mbkt(ji ,jj )
	363	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
	364	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
[1601]	365
[481]	366	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * ze3u
	367	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * ze3v
	368	END DO
	369	END DO
[1601]	370	!
[481]	371	ELSE
[1601]	372	!
[789]	373	# if defined key_vectopt_loop
[1601]	374	DO jj = 1, 1
	375	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[481]	376	# else
	377	DO jj = 1, jpjm1
	378	DO ji = 1, jpim1
	379	# endif
	380	iku = mbku(ji,jj)
	381	ikv = mbkv(ji,jj)
	382	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
	383	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
	384	END DO
	385	END DO
[1601]	386	!
[481]	387	ENDIF
	388
	389
[789]	390	# if defined key_vectopt_loop
[1601]	391	DO jj = 1, 1
	392	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3]	393	# else
	394	DO jj = 2, jpjm1
	395	DO ji = 2, jpim1
	396	# endif
	397	ik = mbkt(ji,jj)
	398	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
	399	zhdivn(ji,jj,ik) = &
[481]	400	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) ) &
	401	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) ) &
	402	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) ) &
	403	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) ) &
[3]	404	& ) / zbt
	405
	406	END DO
[1601]	407	END DO
[3]	408
	409	! 7. compute additional vertical velocity to be used in t boxes
	410	! -------------------------------------------------------------
	411	! ... Computation from the bottom
	412	! Note that w_bbl(:,:,jpk) has been set to 0 in tra_bbl_init
	413	DO jk = jpkm1, 1, -1
	414	DO jj= 2, jpjm1
	415	DO ji = fs_2, fs_jpim1 ! vector opt.
	416	w_bbl(ji,jj,jk) = w_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
	417	END DO
	418	END DO
	419	END DO
[1601]	420	CALL lbc_lnk( w_bbl, 'W', 1. ) ! Boundary condition on w_bbl (unchanged sign)
[3]	421
[1482]	422	CALL iom_put( "uoce_bbl", u_bbl ) ! bbl i-current
	423	CALL iom_put( "voce_bbl", v_bbl ) ! bbl j-current
	424	CALL iom_put( "woce_bbl", w_bbl ) ! bbl vert. current
[503]	425	!
[3]	426	END SUBROUTINE tra_bbl_adv

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

Download in other formats:

Original Format