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

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 789

Last change on this file since 789 was 789, checked in by rblod, 16 years ago
Suppress jki routines and associated key_mpp_omp
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
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)
[719]	10	!! $Header$
[503]	11	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_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
	53	USE oce, ONLY : ztrds => va ! use va as 3D workspace
[3]	54	!!
[503]	55	INTEGER, INTENT( in ) :: kt ! ocean time-step
	56	!!
	57	INTEGER :: ji, jj, jk ! dummy loop indices
	58	INTEGER :: ik ! temporary integers
	59	INTEGER :: iku, iku1, iku2 ! " "
	60	INTEGER :: ikv, ikv1, ikv2 ! " "
	61	REAL(wp) :: zsign, zh, zalbet ! temporary scalars
	62	REAL(wp) :: zgdrho, zbtr ! " "
	63	REAL(wp) :: zbt, zsigna ! " "
	64	REAL(wp) :: zfui, ze3u, zt, zta ! " "
	65	REAL(wp) :: zfvj, ze3v, zs, zsa ! " "
	66	REAL(wp), DIMENSION(jpi,jpj) :: zalphax, zwu, zunb ! temporary 2D workspace
	67	REAL(wp), DIMENSION(jpi,jpj) :: zalphay, zwv, zvnb ! " "
	68	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zww, zwz ! " "
	69	REAL(wp), DIMENSION(jpi,jpj) :: ztbb, zsbb ! " "
	70	REAL(wp), DIMENSION(jpi,jpj) :: ztnb, zsnb, zdep ! " "
	71	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhdivn ! temporary 3D workspace
	72	REAL(wp) :: fsalbt, pft, pfs, pfh ! statement function
[3]	73	!!----------------------------------------------------------------------
	74	! ratio alpha/beta
	75	! ================
	76	! fsalbt: ratio of thermal over saline expension coefficients
	77	! pft : potential temperature in degrees celcius
	78	! pfs : salinity anomaly (s-35) in psu
	79	! pfh : depth in meters
	80
	81	fsalbt( pft, pfs, pfh ) = &
	82	( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft &
	83	- 0.203814e-03 ) * pft &
	84	+ 0.170907e-01 ) * pft &
	85	+ 0.665157e-01 &
	86	+(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs &
	87	+ ( ( - 0.302285e-13 * pfh &
	88	- 0.251520e-11 * pfs &
	89	+ 0.512857e-12 * pft * pft ) * pfh &
	90	- 0.164759e-06 * pfs &
	91	+( 0.791325e-08 * pft - 0.933746e-06 ) * pft &
	92	+ 0.380374e-04 ) * pfh
	93	!!----------------------------------------------------------------------
	94
	95	IF( kt == nit000 ) CALL tra_bbl_init ! initialization at first time-step
	96
	97	! 1. 2D fields of bottom temperature and salinity, and bottom slope
	98	! -----------------------------------------------------------------
	99	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
	100
[789]	101	#if defined key_vectopt_loop
[3]	102	jj = 1
	103	DO ji = 1, jpij ! vector opt. (forced unrolling)
	104	#else
	105	DO jj = 1, jpj
	106	DO ji = 1, jpi
	107	#endif
	108	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
[481]	109	ztnb(ji,jj) = tn(ji,jj,ik)
	110	zsnb(ji,jj) = sn(ji,jj,ik)
	111	ztbb(ji,jj) = tb(ji,jj,ik)
	112	zsbb(ji,jj) = sb(ji,jj,ik)
[3]	113	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
[481]	114
	115	zunb(ji,jj) = un(ji,jj,mbku(ji,jj))
	116	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
[789]	117	#if ! defined key_vectopt_loop
[3]	118	END DO
	119	#endif
	120	END DO
	121
[481]	122
[3]	123	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
	124	! --------------------------------------------
	125	! Sign of the local density gradient along the i- and j-slopes
	126	! multiplied by the slope of the ocean bottom
	127
	128	SELECT CASE ( neos )
[503]	129	!
[3]	130	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
[503]	131	!
[3]	132	DO jj = 1, jpjm1
[503]	133	DO ji = 1, fs_jpim1 ! vector opt.
	134	! ... temperature, salinity anomalie and depth
	135	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
	136	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
	137	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
	138	! ... masked ratio alpha/beta
	139	zalbet = fsalbt( zt, zs, zh ) * umask(ji,jj,1)
	140	! ... local density gradient along i-bathymetric slope
	141	zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
	142	& - ( zsnb(ji+1,jj) - zsnb(ji,jj) )
	143	zgdrho = zgdrho * umask(ji,jj,1)
	144	! ... sign of local i-gradient of density multiplied by the i-slope
	145	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	146	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	147	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
	148	END DO
[3]	149	END DO
[503]	150	!
[3]	151	DO jj = 1, jpjm1
[503]	152	DO ji = 1, fs_jpim1 ! vector opt.
	153	! ... temperature, salinity anomalie and depth
	154	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
	155	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
	156	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
	157	! ... masked ratio alpha/beta
	158	zalbet = fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
	159	! ... local density gradient along j-bathymetric slope
	160	zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
	161	& - ( zsnb(ji,jj+1) - zsnb(ji,jj) )
	162	zgdrho = zgdrho*vmask(ji,jj,1)
	163	! ... sign of local j-gradient of density multiplied by the j-slope
	164	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	165	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	166	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
	167	END DO
[3]	168	END DO
[503]	169	!
[3]	170	CASE ( 1 ) ! Linear formulation function of temperature only
[503]	171	!
[409]	172	DO jj = 1, jpjm1
[481]	173	DO ji = 1, fs_jpim1 ! vector opt.
[409]	174	! local 'density/temperature' gradient along i-bathymetric slope
[481]	175	zgdrho = ( ztnb(ji+1,jj) - ztnb(ji,jj) )
[409]	176	! sign of local i-gradient of density multiplied by the i-slope
[481]	177	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	178	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	179	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
[3]	180
[481]	181	! local density gradient along j-bathymetric slope
	182	zgdrho = ( ztnb(ji,jj+1) - ztnb(ji,jj) )
	183	! sign of local j-gradient of density multiplied by the j-slope
	184	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	185	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	186	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
[409]	187	END DO
	188	END DO
[503]	189	!
[481]	190	CASE ( 2 ) ! Linear formulation function of temperature and salinity
[503]	191	!
[409]	192	DO jj = 1, jpjm1
[481]	193	DO ji = 1, fs_jpim1 ! vector opt.
	194	! local density gradient along i-bathymetric slope
	195	zgdrho = - ( rbeta*( zsnb(ji+1,jj) - zsnb(ji,jj) ) &
	196	& - ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
	197	! sign of local i-gradient of density multiplied by the i-slope
	198	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	199	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
	200	zalphax(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
	201
[409]	202	! local density gradient along j-bathymetric slope
[481]	203	zgdrho = - ( rbeta*( zsnb(ji,jj+1) - zsnb(ji,jj) ) &
	204	- ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )
[409]	205	! sign of local j-gradient of density multiplied by the j-slope
[481]	206	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	207	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
	208	zalphay(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
[409]	209	END DO
	210	END DO
[503]	211	!
[3]	212	CASE DEFAULT
[474]	213	WRITE(ctmp1,*) ' bad flag value for neos = ', neos
	214	CALL ctl_stop( ctmp1 )
[503]	215	!
[3]	216	END SELECT
	217
	218	! lateral boundary conditions on zalphax and zalphay (unchanged sign)
	219	CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. )
	220
[21]	221
[3]	222	! 3. Velocities that are exchanged between ajacent bottom boxes.
	223	!---------------------------------------------------------------
	224
	225	! ... is equal to zero but where bbl will work.
[481]	226	u_bbl(:,:,:) = 0.e0
	227	v_bbl(:,:,:) = 0.e0
	228
	229	IF( ln_zps ) THEN ! partial steps correction
	230
[789]	231	# if defined key_vectopt_loop
[481]	232	jj = 1
	233	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
	234	# else
	235	DO jj = 1, jpjm1
	236	DO ji = 1, jpim1
[3]	237	# endif
[481]	238	iku = mbku(ji ,jj )
	239	ikv = mbkv(ji ,jj )
	240	iku1 = mbkt(ji+1,jj )
	241	iku2 = mbkt(ji ,jj )
	242	ikv1 = mbkt(ji ,jj+1)
	243	ikv2 = mbkt(ji ,jj )
	244	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
	245	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
	246
	247	IF( MAX(iku,ikv) > 1 ) THEN
	248	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * ze3u / fse3u(ji,jj,iku)
	249	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * ze3v / fse3v(ji,jj,ikv)
	250	ENDIF
[789]	251	# if ! defined key_vectopt_loop
[481]	252	END DO
	253	# endif
[3]	254	END DO
	255
[481]	256	! lateral boundary conditions on u_bbl and v_bbl (changed sign)
	257	CALL lbc_lnk( u_bbl, 'U', -1. ) ; CALL lbc_lnk( v_bbl, 'V', -1. )
	258
	259	ELSE ! if not partial step loop over the whole domain no lbc call
[3]	260
[789]	261	#if defined key_vectopt_loop
[481]	262	jj = 1
	263	DO ji = 1, jpij ! vector opt. (forced unrolling)
	264	#else
	265	DO jj = 1, jpj
	266	DO ji = 1, jpi
	267	#endif
	268	iku = mbku(ji,jj)
	269	ikv = mbkv(ji,jj)
	270	IF( MAX(iku,ikv) > 1 ) THEN
	271	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku)
	272	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv)
	273	ENDIF
[789]	274	#if ! defined key_vectopt_loop
[481]	275	END DO
	276	#endif
	277	END DO
	278
	279	ENDIF
	280
	281
[3]	282	! 5. Along sigma advective trend
	283	! -------------------------------
	284	! ... Second order centered tracer flux at u and v-points
	285
[789]	286	# if defined key_vectopt_loop
[3]	287	jj = 1
	288	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
	289	# else
	290	DO jj = 1, jpjm1
	291	DO ji = 1, jpim1
	292	# endif
	293	iku = mbku(ji,jj)
	294	ikv = mbkv(ji,jj)
[481]	295	zfui = e2u(ji,jj) * fse3u(ji,jj,iku) * u_bbl(ji,jj,iku)
	296	zfvj = e1v(ji,jj) * fse3v(ji,jj,ikv) * v_bbl(ji,jj,ikv)
[3]	297	! centered scheme
	298	! zwx(ji,jj) = 0.5* zfui * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
	299	! zwy(ji,jj) = 0.5* zfvj * ( ztnb(ji,jj) + ztnb(ji,jj+1) )
	300	! zww(ji,jj) = 0.5* zfui * ( zsnb(ji,jj) + zsnb(ji+1,jj) )
	301	! zwz(ji,jj) = 0.5* zfvj * ( zsnb(ji,jj) + zsnb(ji,jj+1) )
	302	! upstream scheme
[21]	303	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji ,jj ) &
	304	& +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj ) ) * 0.5
[409]	305	zwy(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * ztbb(ji ,jj ) &
	306	& +( zfvj - ABS( zfvj ) ) * ztbb(ji ,jj+1) ) * 0.5
[21]	307	zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji ,jj ) &
	308	& +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj ) ) * 0.5
[409]	309	zwz(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * zsbb(ji ,jj ) &
	310	& +( zfvj - ABS( zfvj ) ) * zsbb(ji ,jj+1) ) * 0.5
[789]	311	#if ! defined key_vectopt_loop
[3]	312	END DO
	313	#endif
	314	END DO
[789]	315	# if defined key_vectopt_loop
[3]	316	jj = 1
	317	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
	318	# else
	319	DO jj = 2, jpjm1
	320	DO ji = 2, jpim1
	321	# endif
	322	ik = mbkt(ji,jj)
	323	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
[216]	324	! horizontal advective trends
[3]	325	zta = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
	326	& + zwy(ji,jj) - zwy(ji ,jj-1) )
	327	zsa = - zbtr * ( zww(ji,jj) - zww(ji-1,jj ) &
	328	& + zwz(ji,jj) - zwz(ji ,jj-1) )
	329
[216]	330	! add it to the general tracer trends
[3]	331	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
	332	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
[789]	333	#if ! defined key_vectopt_loop
[3]	334	END DO
	335	#endif
	336	END DO
	337
[503]	338	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' bbl - Ta: ', mask1=tmask, &
	339	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[3]	340
	341	! 6. Vertical advection velocities
	342	! --------------------------------
	343	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
	344	DO jk= 1, jpkm1
	345	DO jj=1, jpjm1
[409]	346	DO ji = 1, fs_jpim1 ! vector opt.
	347	zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk) * fse3u(ji,jj,jk)
	348	zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk) * fse3v(ji,jj,jk)
[3]	349	END DO
	350	END DO
	351
	352	! ... horizontal divergence
	353	DO jj = 2, jpjm1
	354	DO ji = fs_2, fs_jpim1 ! vector opt.
[409]	355	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
[3]	356	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
	357	+ zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
	358	END DO
	359	END DO
	360	END DO
	361
	362
	363	! ... horizontal bottom divergence
[481]	364
	365	IF( ln_zps ) THEN
	366
[789]	367	# if defined key_vectopt_loop
[481]	368	jj = 1
	369	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[3]	370	# else
[481]	371	DO jj = 1, jpjm1
	372	DO ji = 1, jpim1
[3]	373	# endif
[481]	374	iku = mbku(ji ,jj )
	375	ikv = mbkv(ji ,jj )
	376	iku1 = mbkt(ji+1,jj )
	377	iku2 = mbkt(ji ,jj )
	378	ikv1 = mbkt(ji ,jj+1)
	379	ikv2 = mbkt(ji ,jj )
	380	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
	381	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
	382
	383	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * ze3u
	384	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * ze3v
[789]	385	#if ! defined key_vectopt_loop
[481]	386	END DO
[3]	387	#endif
[481]	388	END DO
	389
	390	ELSE
[3]	391
[789]	392	# if defined key_vectopt_loop
[481]	393	jj = 1
	394	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
	395	# else
	396	DO jj = 1, jpjm1
	397	DO ji = 1, jpim1
	398	# endif
	399	iku = mbku(ji,jj)
	400	ikv = mbkv(ji,jj)
	401	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
	402	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
[789]	403	#if ! defined key_vectopt_loop
[481]	404	END DO
	405	#endif
	406	END DO
	407
	408	ENDIF
	409
	410
[789]	411	# if defined key_vectopt_loop
[3]	412	jj = 1
	413	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
	414	# else
	415	DO jj = 2, jpjm1
	416	DO ji = 2, jpim1
	417	# endif
	418	ik = mbkt(ji,jj)
	419	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
	420	zhdivn(ji,jj,ik) = &
[481]	421	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) ) &
	422	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) ) &
	423	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) ) &
	424	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) ) &
[3]	425	& ) / zbt
	426
[789]	427	# if ! defined key_vectopt_loop
[3]	428	END DO
	429	# endif
	430	END DO
	431
	432	! 7. compute additional vertical velocity to be used in t boxes
	433	! -------------------------------------------------------------
	434	! ... Computation from the bottom
	435	! Note that w_bbl(:,:,jpk) has been set to 0 in tra_bbl_init
	436	DO jk = jpkm1, 1, -1
	437	DO jj= 2, jpjm1
	438	DO ji = fs_2, fs_jpim1 ! vector opt.
	439	w_bbl(ji,jj,jk) = w_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
	440	END DO
	441	END DO
	442	END DO
	443
	444	! Boundary condition on w_bbl (unchanged sign)
	445	CALL lbc_lnk( w_bbl, 'W', 1. )
[503]	446	!
[3]	447	END SUBROUTINE tra_bbl_adv

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

Download in other formats:

Original Format