Home

Context Navigation

source: palm/trunk/SOURCE/prandtl_fluxes.f90 @ 190

Last change on this file since 190 was 187, checked in by letzel, 16 years ago
new: descriptions of plant canopy model and prandtl layer (trunk/DOC/misc) changed: more consistent flux definitions; modification of the integrated version of the profile function for momentum for unstable stratification (wall_fluxes, production_e) bugfix: change definition of us_wall from 1D to 2D (prandtl_fluxes, wall_fluxes)
Property svn:keywords set to `Id`
File size: 12.9 KB

Rev	Line
[1]	1	SUBROUTINE prandtl_fluxes
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[187]	6	! Bugfix: change definition of us_wall from 1D to 2D:
	7	! Modification of the evaluation of the vertical turbulent momentum
	8	! fluxes u'w' and v'w'; the first usws that is computed corresponds
	9	! to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first
	10	! vsws that is computed corresponds to -v'w'/u* and not as priorily
	11	! assumed to (-v'w')**0.5. Therefore, the intermediate result for
	12	! usws has to be multiplied by -u* instead by itself in order to
	13	! get u'w'. Accordingly, the intermediate result for vsws has to be
	14	! multiplied by -u* instead by itself in order to get v'w'. As u*
	15	! is calculated for the position of a scalar an additional
	16	! interpolation of u* to the position of u and v, respectively,
	17	! is necessary. As u* is not determined for the ghost points on
	18	! each PE, an additional exchange of information from neighbouring
	19	! PEs is necessary.!
	20	! Change: Modification of the integrated version of the profile function for
	21	! momentum for unstable stratification
[1]	22	!
	23	! Former revisions:
	24	! -----------------
[3]	25	! $Id: prandtl_fluxes.f90 187 2008-08-06 16:25:09Z letzel $
[77]	26	!
[110]	27	! 108 2007-08-24 15:10:38Z letzel
	28	! assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean
	29	!
[77]	30	! 75 2007-03-22 09:54:05Z raasch
	31	! moisture renamed humidity
	32	!
[3]	33	! RCS Log replace by Id keyword, revision history cleaned up
	34	!
[1]	35	! Revision 1.19 2006/04/26 12:24:35 raasch
	36	! +OpenMP directives and optimization (array assignments replaced by DO loops)
	37	!
	38	! Revision 1.1 1998/01/23 10:06:06 raasch
	39	! Initial revision
	40	!
	41	!
	42	! Description:
	43	! ------------
	44	! Diagnostic computation of vertical fluxes in the Prandtl layer from the
	45	! values of the variables at grid point k=1
	46	!------------------------------------------------------------------------------!
	47
	48	USE arrays_3d
	49	USE control_parameters
	50	USE grid_variables
	51	USE indices
	52
	53	IMPLICIT NONE
	54
	55	INTEGER :: i, j, k
[108]	56	REAL :: a, b, e_q, rifm, uv_total, z_p
[1]	57
	58	!
	59	!-- Compute theta*
	60	IF ( constant_heatflux ) THEN
	61	!
	62	!-- For a given heat flux in the Prandtl layer:
	63	!-- for u* use the value from the previous time step
	64	!$OMP PARALLEL DO
	65	DO i = nxl-1, nxr+1
	66	DO j = nys-1, nyn+1
	67	ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
	68	!
	69	!-- ts must be limited, because otherwise overflow may occur in case of
	70	!-- us=0 when computing rif further below
	71	IF ( ts(j,i) < -1.05E5 ) ts = -1.0E5
	72	IF ( ts(j,i) > 1.0E5 ) ts = 1.0E5
	73	ENDDO
	74	ENDDO
	75
	76	ELSE
	77	!
	78	!-- For a given surface temperature:
	79	!-- (the Richardson number is still the one from the previous time step)
	80	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
	81	DO i = nxl-1, nxr+1
	82	DO j = nys-1, nyn+1
	83
	84	k = nzb_s_inner(j,i)
	85	z_p = zu(k+1) - zw(k)
	86
	87	IF ( rif(j,i) >= 0.0 ) THEN
	88	!
	89	!-- Stable stratification
	90	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
	91	LOG( z_p / z0(j,i) ) + &
	92	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
	93	)
	94	ELSE
	95	!
	96	!-- Unstable stratification
	97	a = SQRT( 1.0 - 16.0 * rif(j,i) )
[187]	98	b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
	99
	100	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
	101	LOG( z_p / z0(j,i) ) - &
	102	2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
[1]	103	ENDIF
	104
	105	ENDDO
	106	ENDDO
	107	ENDIF
	108
	109	!
	110	!-- Compute z_p/L (corresponds to the Richardson-flux number)
[75]	111	IF ( .NOT. humidity ) THEN
[1]	112	!$OMP PARALLEL DO PRIVATE( k, z_p )
	113	DO i = nxl-1, nxr+1
	114	DO j = nys-1, nyn+1
	115	k = nzb_s_inner(j,i)
	116	z_p = zu(k+1) - zw(k)
	117	rif(j,i) = z_p * kappa * g * ts(j,i) / &
	118	( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
	119	!
	120	!-- Limit the value range of the Richardson numbers.
	121	!-- This is necessary for very small velocities (u,v --> 0), because
	122	!-- the absolute value of rif can then become very large, which in
	123	!-- consequence would result in very large shear stresses and very
	124	!-- small momentum fluxes (both are generally unrealistic).
	125	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
	126	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
	127	ENDDO
	128	ENDDO
	129	ELSE
	130	!$OMP PARALLEL DO PRIVATE( k, z_p )
	131	DO i = nxl-1, nxr+1
	132	DO j = nys-1, nyn+1
	133	k = nzb_s_inner(j,i)
	134	z_p = zu(k+1) - zw(k)
	135	rif(j,i) = z_p * kappa * g * &
	136	( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
	137	( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
	138	!
	139	!-- Limit the value range of the Richardson numbers.
	140	!-- This is necessary for very small velocities (u,v --> 0), because
	141	!-- the absolute value of rif can then become very large, which in
	142	!-- consequence would result in very large shear stresses and very
	143	!-- small momentum fluxes (both are generally unrealistic).
	144	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
	145	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
	146	ENDDO
	147	ENDDO
	148	ENDIF
	149
	150	!
	151	!-- Compute u* at the scalars' grid points
	152	!$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
	153	DO i = nxl, nxr
	154	DO j = nys, nyn
	155
	156	k = nzb_s_inner(j,i)
	157	z_p = zu(k+1) - zw(k)
	158
	159	!
	160	!-- Compute the absolute value of the horizontal velocity
	161	uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
	162	( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2 &
	163	)
	164
	165	IF ( rif(j,i) >= 0.0 ) THEN
	166	!
	167	!-- Stable stratification
	168	us(j,i) = kappa * uv_total / ( &
	169	LOG( z_p / z0(j,i) ) + &
	170	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
	171	)
	172	ELSE
	173	!
	174	!-- Unstable stratification
[187]	175	a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	176	b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
	177
	178	us(j,i) = kappa * uv_total / ( &
	179	LOG( z_p / z0(j,i) ) - &
	180	LOG( ( 1.0 + a )*2 ( 1.0 + a**2 ) / ( &
	181	( 1.0 + b )*2 ( 1.0 + b**2 ) ) ) + &
	182	2.0 * ( ATAN( a ) - ATAN( b ) ) &
	183	)
[1]	184	ENDIF
	185	ENDDO
	186	ENDDO
	187
	188	!
[187]	189	!-- Values of us at ghost point locations are needed for the evaluation of usws
	190	!-- and vsws.
	191	CALL exchange_horiz_2d( us )
	192	!
[1]	193	!-- Compute u'w' for the total model domain.
	194	!-- First compute the corresponding component of u* and square it.
	195	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
	196	DO i = nxl, nxr
	197	DO j = nys, nyn
	198
	199	k = nzb_u_inner(j,i)
	200	z_p = zu(k+1) - zw(k)
	201
	202	!
	203	!-- Compute Richardson-flux number for this point
	204	rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
	205	IF ( rifm >= 0.0 ) THEN
	206	!
	207	!-- Stable stratification
	208	usws(j,i) = kappa * u(k+1,j,i) / ( &
	209	LOG( z_p / z0(j,i) ) + &
	210	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
	211	)
	212	ELSE
	213	!
	214	!-- Unstable stratification
[187]	215	a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
	216	b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
	217
	218	usws(j,i) = kappa * u(k+1,j,i) / ( &
	219	LOG( z_p / z0(j,i) ) - &
	220	LOG( (1.0 + a )*2 ( 1.0 + a**2 ) / ( &
	221	(1.0 + b )*2 ( 1.0 + b**2 ) ) ) + &
	222	2.0 * ( ATAN( a ) - ATAN( b ) ) &
[1]	223	)
	224	ENDIF
[187]	225	usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
[1]	226	ENDDO
	227	ENDDO
	228
	229	!
	230	!-- Compute v'w' for the total model domain.
	231	!-- First compute the corresponding component of u* and square it.
	232	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
	233	DO i = nxl, nxr
	234	DO j = nys, nyn
	235
	236	k = nzb_v_inner(j,i)
	237	z_p = zu(k+1) - zw(k)
	238
	239	!
	240	!-- Compute Richardson-flux number for this point
	241	rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
	242	IF ( rifm >= 0.0 ) THEN
	243	!
	244	!-- Stable stratification
	245	vsws(j,i) = kappa * v(k+1,j,i) / ( &
	246	LOG( z_p / z0(j,i) ) + &
	247	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
	248	)
	249	ELSE
	250	!
	251	!-- Unstable stratification
[187]	252	a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
	253	b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
	254
	255	vsws(j,i) = kappa * v(k+1,j,i) / ( &
	256	LOG( z_p / z0(j,i) ) - &
	257	LOG( (1.0 + a )*2 ( 1.0 + a**2 ) / ( &
	258	(1.0 + b )*2 ( 1.0 + b**2 ) ) ) + &
	259	2.0 * ( ATAN( a ) - ATAN( b ) ) &
[1]	260	)
	261	ENDIF
[187]	262	vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
[1]	263	ENDDO
	264	ENDDO
	265
	266	!
	267	!-- If required compute q*
[75]	268	IF ( humidity .OR. passive_scalar ) THEN
[1]	269	IF ( constant_waterflux ) THEN
	270	!
	271	!-- For a given water flux in the Prandtl layer:
	272	!$OMP PARALLEL DO
	273	DO i = nxl-1, nxr+1
	274	DO j = nys-1, nyn+1
	275	qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
	276	ENDDO
	277	ENDDO
	278
	279	ELSE
	280	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
	281	DO i = nxl-1, nxr+1
	282	DO j = nys-1, nyn+1
	283
	284	k = nzb_s_inner(j,i)
	285	z_p = zu(k+1) - zw(k)
	286
[108]	287	!
	288	!-- assume saturation for atmosphere coupled to ocean
	289	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	290	e_q = 6.1 * &
	291	EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
	292	q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
	293	ENDIF
[1]	294	IF ( rif(j,i) >= 0.0 ) THEN
	295	!
	296	!-- Stable stratification
	297	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
	298	LOG( z_p / z0(j,i) ) + &
	299	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
	300	)
	301	ELSE
	302	!
	303	!-- Unstable stratification
[187]	304	a = SQRT( 1.0 - 16.0 * rif(j,i) )
	305	b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
	306
	307	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
	308	LOG( z_p / z0(j,i) ) - &
	309	2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
[1]	310	ENDIF
	311
	312	ENDDO
	313	ENDDO
	314	ENDIF
	315	ENDIF
	316
	317	!
[187]	318	!-- Exchange the boundaries for the momentum fluxes (only for sake of
	319	!-- completeness)
[1]	320	CALL exchange_horiz_2d( usws )
	321	CALL exchange_horiz_2d( vsws )
[75]	322	IF ( humidity .OR. passive_scalar ) CALL exchange_horiz_2d( qsws )
[1]	323
	324	!
	325	!-- Compute the vertical kinematic heat flux
	326	IF ( .NOT. constant_heatflux ) THEN
	327	!$OMP PARALLEL DO
	328	DO i = nxl-1, nxr+1
	329	DO j = nys-1, nyn+1
	330	shf(j,i) = -ts(j,i) * us(j,i)
	331	ENDDO
	332	ENDDO
	333	ENDIF
	334
	335	!
	336	!-- Compute the vertical water/scalar flux
[75]	337	IF ( .NOT. constant_heatflux .AND. ( humidity .OR. passive_scalar ) ) THEN
[1]	338	!$OMP PARALLEL DO
	339	DO i = nxl-1, nxr+1
	340	DO j = nys-1, nyn+1
	341	qsws(j,i) = -qs(j,i) * us(j,i)
	342	ENDDO
	343	ENDDO
	344	ENDIF
	345
	346	!
	347	!-- Bottom boundary condition for the TKE
	348	IF ( ibc_e_b == 2 ) THEN
	349	!$OMP PARALLEL DO
	350	DO i = nxl-1, nxr+1
	351	DO j = nys-1, nyn+1
	352	e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
	353	!
	354	!-- As a test: cm = 0.4
	355	! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
	356	e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
	357	ENDDO
	358	ENDDO
	359	ENDIF
	360
	361
	362	END SUBROUTINE prandtl_fluxes

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |