Home

Context Navigation

source: palm/trunk/SOURCE/model_1d_mod.f90 @ 2636

Last change on this file since 2636 was 2339, checked in by gronemeier, 7 years ago
corrected timestamp in header
Property svn:keywords set to `Id`
File size: 37.1 KB

Rev	Line
[2338]	1	!> @file model_1d_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[1961]	22	!
[2299]	23	!
[1961]	24	! Former revisions:
	25	! -----------------
	26	! $Id: model_1d_mod.f90 2339 2017-08-07 13:55:26Z raasch $
[2339]	27	! corrected timestamp in header
	28	!
	29	! 2338 2017-08-07 12:15:38Z gronemeier
[2338]	30	! renamed init_1d_model to model_1d_mod and and formatted it as a module;
	31	! reformatted output of profiles
	32	!
[2339]	33	! 2337 2017-08-07 08:59:53Z gronemeier
[2337]	34	! revised calculation of mixing length
	35	! removed rounding of time step
	36	! corrected calculation of virtual potential temperature
	37	!
	38	! 2334 2017-08-04 11:57:04Z gronemeier
[2334]	39	! set c_m = 0.4 according to Detering and Etling (1985)
	40	!
	41	! 2299 2017-06-29 10:14:38Z maronga
[2299]	42	! Removed german text
[1961]	43	!
[2299]	44	! 2101 2017-01-05 16:42:31Z suehring
	45	!
[2060]	46	! 2059 2016-11-10 14:20:40Z maronga
	47	! Corrected min/max values of Rif.
	48	!
[2001]	49	! 2000 2016-08-20 18:09:15Z knoop
	50	! Forced header and separation lines into 80 columns
	51	!
[1961]	52	! 1960 2016-07-12 16:34:24Z suehring
[1960]	53	! Remove passive_scalar from IF-statements, as 1D-scalar profile is effectively
	54	! not used.
	55	! Formatting adjustment
[1809]	56	!
	57	! 1808 2016-04-05 19:44:00Z raasch
	58	! routine local_flush replaced by FORTRAN statement
	59	!
[1710]	60	! 1709 2015-11-04 14:47:01Z maronga
	61	! Set initial time step to 10 s to avoid instability of the 1d model for small
	62	! grid spacings
	63	!
[1698]	64	! 1697 2015-10-28 17:14:10Z raasch
	65	! small E- and F-FORMAT changes to avoid informative compiler messages about
	66	! insufficient field width
	67	!
[1692]	68	! 1691 2015-10-26 16:17:44Z maronga
	69	! Renamed prandtl_layer to constant_flux_layer. rif is replaced by ol and zeta.
	70	!
[1683]	71	! 1682 2015-10-07 23:56:08Z knoop
	72	! Code annotations made doxygen readable
	73	!
[1354]	74	! 1353 2014-04-08 15:21:23Z heinze
	75	! REAL constants provided with KIND-attribute
	76	!
[1347]	77	! 1346 2014-03-27 13:18:20Z heinze
	78	! Bugfix: REAL constants provided with KIND-attribute especially in call of
	79	! intrinsic function like MAX, MIN, SIGN
	80	!
[1323]	81	! 1322 2014-03-20 16:38:49Z raasch
	82	! REAL functions provided with KIND-attribute
	83	!
[1321]	84	! 1320 2014-03-20 08:40:49Z raasch
[1320]	85	! ONLY-attribute added to USE-statements,
	86	! kind-parameters added to all INTEGER and REAL declaration statements,
	87	! kinds are defined in new module kinds,
	88	! revision history before 2012 removed,
	89	! comment fields (!:) to be used for variable explanations added to
	90	! all variable declaration statements
[1321]	91	!
[1037]	92	! 1036 2012-10-22 13:43:42Z raasch
	93	! code put under GPL (PALM 3.9)
	94	!
[1017]	95	! 1015 2012-09-27 09:23:24Z raasch
	96	! adjustment of mixing length to the Prandtl mixing length at first grid point
	97	! above ground removed
	98	!
[1002]	99	! 1001 2012-09-13 14:08:46Z raasch
	100	! all actions concerning leapfrog scheme removed
	101	!
[997]	102	! 996 2012-09-07 10:41:47Z raasch
	103	! little reformatting
	104	!
[979]	105	! 978 2012-08-09 08:28:32Z fricke
	106	! roughness length for scalar quantities z0h1d added
	107	!
[1]	108	! Revision 1.1 1998/03/09 16:22:10 raasch
	109	! Initial revision
	110	!
	111	!
	112	! Description:
	113	! ------------
[1682]	114	!> 1D-model to initialize the 3D-arrays.
	115	!> The temperature profile is set as steady and a corresponding steady solution
	116	!> of the wind profile is being computed.
	117	!> All subroutines required can be found within this file.
[1691]	118	!>
	119	!> @todo harmonize code with new surface_layer_fluxes module
[1709]	120	!> @bug 1D model crashes when using small grid spacings in the order of 1 m
[1]	121	!------------------------------------------------------------------------------!
[2338]	122	MODULE model_1d_mod
[1]	123
[1320]	124	USE arrays_3d, &
[2338]	125	ONLY: dd2zu, ddzu, ddzw, dzu, l_grid, pt_init, q_init, ug, u_init, &
	126	vg, v_init, zu
[1320]	127
[2338]	128	USE control_parameters, &
	129	ONLY: constant_diffusion, constant_flux_layer, dissipation_1d, f, g, &
	130	humidity, ibc_e_b, intermediate_timestep_count, &
	131	intermediate_timestep_count_max, kappa, km_constant, &
	132	message_string, mixing_length_1d, prandtl_number, &
	133	roughness_length, run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
	134
[1320]	135	USE indices, &
[2338]	136	ONLY: nzb, nzb_diff, nzt
[1320]	137
	138	USE kinds
[1]	139
[2338]	140	USE pegrid
	141
	142
[1]	143	IMPLICIT NONE
	144
[2338]	145
	146	INTEGER(iwp) :: current_timestep_number_1d = 0 !< current timestep number (1d-model)
	147	INTEGER(iwp) :: damp_level_ind_1d !< lower grid index of damping layer (1d-model) !#
	148
	149	LOGICAL :: run_control_header_1d = .FALSE. !< flag for output of run control header (1d-model)
	150	LOGICAL :: stop_dt_1d = .FALSE. !< termination flag, used in case of too small timestep (1d-model)
	151
	152	REAL(wp) :: c_m = 0.4_wp !< model constant, 0.4 according to Detering and Etling (1985)
	153	REAL(wp) :: damp_level_1d = -1.0_wp !< namelist parameter !#
	154	REAL(wp) :: dt_1d = 60.0_wp !< dynamic timestep (1d-model)
	155	REAL(wp) :: dt_max_1d = 300.0_wp !< timestep limit (1d-model)
	156	REAL(wp) :: dt_pr_1d = 9999999.9_wp !< namelist parameter !#
	157	REAL(wp) :: dt_run_control_1d = 60.0_wp !< namelist parameter !#
	158	REAL(wp) :: end_time_1d = 864000.0_wp !< namelist parameter !#
	159	REAL(wp) :: qs1d !< characteristic humidity scale (1d-model)
	160	REAL(wp) :: simulated_time_1d = 0.0_wp !< updated simulated time (1d-model)
	161	REAL(wp) :: time_pr_1d = 0.0_wp !< updated simulated time for profile output (1d-model)
	162	REAL(wp) :: time_run_control_1d = 0.0_wp !< updated simulated time for run-control output (1d-model)
	163	REAL(wp) :: ts1d !< characteristic temperature scale (1d-model)
	164	REAL(wp) :: us1d !< friction velocity (1d-model) !#
	165	REAL(wp) :: usws1d !< u-component of the momentum flux (1d-model) !#
	166	REAL(wp) :: vsws1d !< v-component of the momentum flux (1d-model) !#
	167	REAL(wp) :: z01d !< roughness length for momentum (1d-model)
	168	REAL(wp) :: z0h1d !< roughness length for scalars (1d-model)
	169
	170
	171	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d !< tke (1d-model) !#
	172	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d_p !< prognostic value of tke (1d-model)
	173	REAL(wp), DIMENSION(:), ALLOCATABLE :: kh1d !< turbulent diffusion coefficient for heat (1d-model) !#
	174	REAL(wp), DIMENSION(:), ALLOCATABLE :: km1d !< turbulent diffusion coefficient for momentum (1d-model)!#
	175	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length Blackadar (1d-model)
	176	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d !< mixing length for turbulent diffusion coefficients (1d-model) !#
	177	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_diss !< mixing length for dissipation (1d-model)
	178	REAL(wp), DIMENSION(:), ALLOCATABLE :: rif1d !< Richardson flux number (1d-model) !#
	179	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_e !< tendency of e (1d-model)
	180	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_em !< weighted tendency of e for previous sub-timestep (1d-model)
	181	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_u !< tendency of u (1d-model)
	182	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_um !< weighted tendency of u for previous sub-timestep (1d-model)
	183	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_v !< tendency of v (1d-model)
	184	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_vm !< weighted tendency of v for previous sub-timestep (1d-model)
	185	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d !< u-velocity component (1d-model) !#
	186	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d_p !< prognostic value of u-velocity component (1d-model)
	187	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d !< v-velocity component (1d-model) !#
	188	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d_p !< prognostic value of v-velocity component (1d-model)
	189
	190	!
	191	!-- Initialize 1D model
	192	INTERFACE init_1d_model
	193	MODULE PROCEDURE init_1d_model
	194	END INTERFACE init_1d_model
	195
	196	!
	197	!-- Print profiles
	198	INTERFACE print_1d_model
	199	MODULE PROCEDURE print_1d_model
	200	END INTERFACE print_1d_model
	201
	202	!
	203	!-- Print run control information
	204	INTERFACE run_control_1d
	205	MODULE PROCEDURE run_control_1d
	206	END INTERFACE run_control_1d
	207
	208	!
	209	!-- Main procedure
	210	INTERFACE time_integration_1d
	211	MODULE PROCEDURE time_integration_1d
	212	END INTERFACE time_integration_1d
	213
	214	!
	215	!-- Calculate time step
	216	INTERFACE timestep_1d
	217	MODULE PROCEDURE timestep_1d
	218	END INTERFACE timestep_1d
	219
	220	SAVE
	221
	222	PRIVATE
	223	!
	224	!-- Public interfaces
	225	PUBLIC init_1d_model
	226
	227	!
	228	!-- Public variables
	229	PUBLIC damp_level_1d, damp_level_ind_1d, dt_pr_1d, dt_run_control_1d, &
	230	e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d, &
	231	vsws1d
	232
	233
	234	CONTAINS
	235
	236	SUBROUTINE init_1d_model
	237
	238	IMPLICIT NONE
	239
[1682]	240	CHARACTER (LEN=9) :: time_to_string !<
[1320]	241
[1682]	242	INTEGER(iwp) :: k !<
[1320]	243
[1682]	244	REAL(wp) :: lambda !<
[1]	245
	246	!
	247	!-- Allocate required 1D-arrays
[2338]	248	ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), &
	249	km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1), &
	250	l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1), &
	251	te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &
	252	te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &
	253	u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
[1]	254
	255	!
	256	!-- Initialize arrays
	257	IF ( constant_diffusion ) THEN
[1001]	258	km1d = km_constant
	259	kh1d = km_constant / prandtl_number
[1]	260	ELSE
[1353]	261	e1d = 0.0_wp; e1d_p = 0.0_wp
	262	kh1d = 0.0_wp; km1d = 0.0_wp
	263	rif1d = 0.0_wp
[1]	264	!
	265	!-- Compute the mixing length
[1353]	266	l_black(nzb) = 0.0_wp
[1]	267
	268	IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN
	269	!
	270	!-- Blackadar mixing length
[1353]	271	IF ( f /= 0.0_wp ) THEN
	272	lambda = 2.7E-4_wp * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
	273	ABS( f ) + 1E-10_wp
[1]	274	ELSE
[1353]	275	lambda = 30.0_wp
[1]	276	ENDIF
	277
	278	DO k = nzb+1, nzt+1
[1353]	279	l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
[1]	280	ENDDO
	281
	282	ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN
	283	!
	284	!-- Use the same mixing length as in 3D model
	285	l_black(1:nzt) = l_grid
	286	l_black(nzt+1) = l_black(nzt)
	287
	288	ENDIF
	289	ENDIF
[2337]	290	l1d = l_black
	291	l1d_diss = l_black
	292	u1d = u_init
	293	u1d_p = u_init
	294	v1d = v_init
	295	v1d_p = v_init
[1]	296
	297	!
	298	!-- Set initial horizontal velocities at the lowest grid levels to a very small
	299	!-- value in order to avoid too small time steps caused by the diffusion limit
	300	!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
[1353]	301	u1d(0:1) = 0.1_wp
	302	u1d_p(0:1) = 0.1_wp
	303	v1d(0:1) = 0.1_wp
	304	v1d_p(0:1) = 0.1_wp
[1]	305
	306	!
	307	!-- For u, theta and the momentum fluxes plausible values are set
[1691]	308	IF ( constant_flux_layer ) THEN
[1353]	309	us1d = 0.1_wp ! without initial friction the flow would not change
[1]	310	ELSE
[1353]	311	e1d(nzb+1) = 1.0_wp
	312	km1d(nzb+1) = 1.0_wp
	313	us1d = 0.0_wp
[1]	314	ENDIF
[1353]	315	ts1d = 0.0_wp
	316	usws1d = 0.0_wp
	317	vsws1d = 0.0_wp
[996]	318	z01d = roughness_length
[978]	319	z0h1d = z0h_factor * z01d
[1960]	320	IF ( humidity ) qs1d = 0.0_wp
[1]	321
	322	!
[46]	323	!-- Tendencies must be preset in order to avoid runtime errors within the
	324	!-- first Runge-Kutta step
[1353]	325	te_em = 0.0_wp
	326	te_um = 0.0_wp
	327	te_vm = 0.0_wp
[46]	328
	329	!
[2338]	330	!-- Set model constant
	331	IF ( dissipation_1d == 'as_in_3d_model' ) c_m = 0.1_wp
	332
	333	!
[1]	334	!-- Set start time in hh:mm:ss - format
	335	simulated_time_chr = time_to_string( simulated_time_1d )
	336
	337	!
[2337]	338	!-- Integrate the 1D-model equations using the Runge-Kutta scheme
[1]	339	CALL time_integration_1d
	340
	341
	342	END SUBROUTINE init_1d_model
	343
	344
	345
	346	!------------------------------------------------------------------------------!
	347	! Description:
	348	! ------------
[2338]	349	!> Runge-Kutta time differencing scheme for the 1D-model.
[1]	350	!------------------------------------------------------------------------------!
[1682]	351
	352	SUBROUTINE time_integration_1d
[1]	353
	354	IMPLICIT NONE
	355
[1682]	356	CHARACTER (LEN=9) :: time_to_string !<
[1320]	357
[2338]	358	INTEGER(iwp) :: k !< loop index
[1320]	359
[2338]	360	REAL(wp) :: a !< auxiliary variable
	361	REAL(wp) :: b !< auxiliary variable
	362	REAL(wp) :: dissipation !< dissipation of TKE
	363	REAL(wp) :: dpt_dz !< vertical temperature gradient
	364	REAL(wp) :: flux !< vertical temperature gradient
	365	REAL(wp) :: kmzm !< Km(z-dz/2)
	366	REAL(wp) :: kmzp !< Km(z+dz/2)
	367	REAL(wp) :: l_stable !< mixing length for stable case
	368	REAL(wp) :: pt_0 !< reference temperature
	369	REAL(wp) :: uv_total !< horizontal wind speed
[1]	370
	371	!
	372	!-- Determine the time step at the start of a 1D-simulation and
	373	!-- determine and printout quantities used for run control
[1709]	374	dt_1d = 10.0_wp
[1]	375	CALL run_control_1d
	376
	377	!
	378	!-- Start of time loop
	379	DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d )
	380
	381	!
	382	!-- Depending on the timestep scheme, carry out one or more intermediate
	383	!-- timesteps
	384
	385	intermediate_timestep_count = 0
	386	DO WHILE ( intermediate_timestep_count < &
	387	intermediate_timestep_count_max )
	388
	389	intermediate_timestep_count = intermediate_timestep_count + 1
	390
	391	CALL timestep_scheme_steering
	392
	393	!
	394	!-- Compute all tendency terms. If a Prandtl-layer is simulated, k starts
	395	!-- at nzb+2.
	396	DO k = nzb_diff, nzt
	397
[1353]	398	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	399	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[1]	400	!
	401	!-- u-component
	402	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
[1001]	403	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
	404	- kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) &
	405	) * ddzw(k)
[1]	406	!
	407	!-- v-component
[1001]	408	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	409	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
	410	- kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) &
	411	) * ddzw(k)
[1]	412	ENDDO
	413	IF ( .NOT. constant_diffusion ) THEN
	414	DO k = nzb_diff, nzt
	415	!
	416	!-- TKE
[1353]	417	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	418	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[75]	419	IF ( .NOT. humidity ) THEN
[1]	420	pt_0 = pt_init(k)
	421	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	422	ELSE
[1353]	423	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
	424	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
[2337]	425	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
	426	pt_init(k-1) * q_init(k-1) ) &
	427	) * dd2zu(k)
[1]	428	ENDIF
	429
	430	IF ( dissipation_1d == 'detering' ) THEN
	431	!
[2337]	432	!-- According to Detering, c_e=c_m^3
	433	dissipation = c_m*3 e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	434	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
[2337]	435	dissipation = ( 0.19_wp &
	436	+ 0.74_wp * l1d_diss(k) / l_grid(k) &
	437	) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	438	ENDIF
	439
	440	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
	441	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
	442	) &
	443	- g / pt_0 * kh1d(k) * flux &
	444	+ ( &
[1001]	445	kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &
	446	- kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &
[1]	447	) * ddzw(k) &
[1001]	448	- dissipation
[1]	449	ENDDO
	450	ENDIF
	451
	452	!
	453	!-- Tendency terms at the top of the Prandtl-layer.
	454	!-- Finite differences of the momentum fluxes are computed using half the
	455	!-- normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
[1691]	456	IF ( constant_flux_layer ) THEN
[1]	457
	458	k = nzb+1
[1353]	459	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	460	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[75]	461	IF ( .NOT. humidity ) THEN
[1]	462	pt_0 = pt_init(k)
	463	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	464	ELSE
[1353]	465	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
	466	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
[2337]	467	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
	468	pt_init(k-1) * q_init(k-1) ) &
[1]	469	) * dd2zu(k)
	470	ENDIF
	471
	472	IF ( dissipation_1d == 'detering' ) THEN
	473	!
	474	!-- According to Detering, c_e=0.064
[2337]	475	dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	476	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
[2337]	477	dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) ) &
	478	* e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	479	ENDIF
	480
	481	!
	482	!-- u-component
[1001]	483	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
	484	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d &
[1353]	485	) * 2.0_wp * ddzw(k)
[1]	486	!
	487	!-- v-component
[1001]	488	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	489	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d &
[1353]	490	) * 2.0_wp * ddzw(k)
[1]	491	!
	492	!-- TKE
	493	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	494	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	495	) &
	496	- g / pt_0 * kh1d(k) * flux &
	497	+ ( &
[1001]	498	kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &
	499	- kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &
[1]	500	) * ddzw(k) &
[1001]	501	- dissipation
[1]	502	ENDIF
	503
	504	!
	505	!-- Prognostic equations for all 1D variables
	506	DO k = nzb+1, nzt
	507
[1001]	508	u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
	509	tsc(3) * te_um(k) )
	510	v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
	511	tsc(3) * te_vm(k) )
[1]	512
	513	ENDDO
	514	IF ( .NOT. constant_diffusion ) THEN
	515	DO k = nzb+1, nzt
	516
[1001]	517	e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
	518	tsc(3) * te_em(k) )
[1]	519
	520	ENDDO
	521	!
	522	!-- Eliminate negative TKE values, which can result from the
	523	!-- integration due to numerical inaccuracies. In such cases the TKE
	524	!-- value is reduced to 10 percent of its old value.
[1353]	525	WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d
[1]	526	ENDIF
	527
	528	!
	529	!-- Calculate tendencies for the next Runge-Kutta step
	530	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	531	IF ( intermediate_timestep_count == 1 ) THEN
	532
	533	DO k = nzb+1, nzt
	534	te_um(k) = te_u(k)
	535	te_vm(k) = te_v(k)
	536	ENDDO
	537
	538	IF ( .NOT. constant_diffusion ) THEN
	539	DO k = nzb+1, nzt
	540	te_em(k) = te_e(k)
	541	ENDDO
	542	ENDIF
	543
	544	ELSEIF ( intermediate_timestep_count < &
	545	intermediate_timestep_count_max ) THEN
	546
	547	DO k = nzb+1, nzt
[1353]	548	te_um(k) = -9.5625_wp * te_u(k) + 5.3125_wp * te_um(k)
	549	te_vm(k) = -9.5625_wp * te_v(k) + 5.3125_wp * te_vm(k)
[1]	550	ENDDO
	551
	552	IF ( .NOT. constant_diffusion ) THEN
	553	DO k = nzb+1, nzt
[1353]	554	te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
[1]	555	ENDDO
	556	ENDIF
	557
	558	ENDIF
	559	ENDIF
	560
	561
	562	!
	563	!-- Boundary conditions for the prognostic variables.
	564	!-- At the top boundary (nzt+1) u,v and e keep their initial values
[2334]	565	!-- (ug(nzt+1), vg(nzt+1), 0).
	566	!-- At the bottom boundary, Dirichlet condition is used for u and v (0)
	567	!-- and Neumann condition for e (e(nzb)=e(nzb+1)).
[1353]	568	u1d_p(nzb) = 0.0_wp
	569	v1d_p(nzb) = 0.0_wp
[667]	570
[1]	571	!
	572	!-- Swap the time levels in preparation for the next time step.
	573	u1d = u1d_p
	574	v1d = v1d_p
	575	IF ( .NOT. constant_diffusion ) THEN
	576	e1d = e1d_p
	577	ENDIF
	578
	579	!
	580	!-- Compute diffusion quantities
	581	IF ( .NOT. constant_diffusion ) THEN
	582
	583	!
	584	!-- First compute the vertical fluxes in the Prandtl-layer
[1691]	585	IF ( constant_flux_layer ) THEN
[1]	586	!
	587	!-- Compute theta* using Rif numbers of the previous time step
[2334]	588	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	589	!
	590	!-- Stable stratification
[1353]	591	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	592	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
	593	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
[1]	594	)
	595	ELSE
	596	!
	597	!-- Unstable stratification
[1353]	598	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
	599	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	600	zu(nzb+1) * z0h1d )
[2337]	601
	602	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	603	LOG( (a-1.0_wp) / (a+1.0_wp) * &
	604	(b+1.0_wp) / (b-1.0_wp) )
[1]	605	ENDIF
	606
[1691]	607	ENDIF ! constant_flux_layer
[1]	608
	609	!
	610	!-- Compute the Richardson-flux numbers,
	611	!-- first at the top of the Prandtl-layer using u* of the previous
	612	!-- time step (+1E-30, if u* = 0), then in the remaining area. There
	613	!-- the rif-numbers of the previous time step are used.
	614
[1691]	615	IF ( constant_flux_layer ) THEN
[75]	616	IF ( .NOT. humidity ) THEN
[1]	617	pt_0 = pt_init(nzb+1)
	618	flux = ts1d
	619	ELSE
[1353]	620	pt_0 = pt_init(nzb+1) * ( 1.0_wp + 0.61_wp * q_init(nzb+1) )
	621	flux = ts1d + 0.61_wp * pt_init(k) * qs1d
[1]	622	ENDIF
	623	rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
[1353]	624	( pt_0 * ( us1d**2 + 1E-30_wp ) )
[1]	625	ENDIF
	626
	627	DO k = nzb_diff, nzt
[75]	628	IF ( .NOT. humidity ) THEN
[1]	629	pt_0 = pt_init(k)
	630	flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
	631	ELSE
[1353]	632	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
[1]	633	flux = ( ( pt_init(k+1) - pt_init(k-1) ) &
[2337]	634	+ 0.61_wp &
	635	* ( pt_init(k+1) * q_init(k+1) &
	636	- pt_init(k-1) * q_init(k-1) ) &
[1]	637	) * dd2zu(k)
	638	ENDIF
[1353]	639	IF ( rif1d(k) >= 0.0_wp ) THEN
	640	rif1d(k) = g / pt_0 * flux / &
	641	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	642	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	643	+ 1E-30_wp &
[1]	644	)
	645	ELSE
[1353]	646	rif1d(k) = g / pt_0 * flux / &
	647	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	648	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	649	+ 1E-30_wp &
	650	) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
[1]	651	ENDIF
	652	ENDDO
	653	!
	654	!-- Richardson-numbers must remain restricted to a realistic value
	655	!-- range. It is exceeded excessively for very small velocities
	656	!-- (u,v --> 0).
[2059]	657	WHERE ( rif1d < -5.0_wp ) rif1d = -5.0_wp
	658	WHERE ( rif1d > 1.0_wp ) rif1d = 1.0_wp
[1]	659
	660	!
	661	!-- Compute u* from the absolute velocity value
[1691]	662	IF ( constant_flux_layer ) THEN
[1]	663	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
	664
[1353]	665	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	666	!
	667	!-- Stable stratification
	668	us1d = kappa * uv_total / ( &
[1353]	669	LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * &
[1]	670	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	671	)
	672	ELSE
	673	!
	674	!-- Unstable stratification
[1353]	675	a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) )
	676	b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	677	zu(nzb+1) * z01d ) )
[2337]	678	us1d = kappa * uv_total / ( &
	679	LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / &
	680	(1.0_wp+a) ) + &
	681	2.0_wp * ( ATAN( b ) - ATAN( a ) ) &
	682	)
[1]	683	ENDIF
	684
	685	!
	686	!-- Compute the momentum fluxes for the diffusion terms
	687	usws1d = - u1d(nzb+1) / uv_total * us1d**2
	688	vsws1d = - v1d(nzb+1) / uv_total * us1d**2
	689
	690	!
	691	!-- Boundary condition for the turbulent kinetic energy at the top
	692	!-- of the Prandtl-layer. c_m = 0.4 according to Detering.
	693	!-- Additional Neumann condition de/dz = 0 at nzb is set to ensure
	694	!-- compatibility with the 3D model.
	695	IF ( ibc_e_b == 2 ) THEN
[2334]	696	e1d(nzb+1) = ( us1d / c_m )**2
[1]	697	ENDIF
	698	e1d(nzb) = e1d(nzb+1)
	699
[1960]	700	IF ( humidity ) THEN
[1]	701	!
	702	!-- Compute q*
[2334]	703	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	704	!
[1960]	705	!-- Stable stratification
	706	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
[1353]	707	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
	708	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
[1]	709	)
[1960]	710	ELSE
[1]	711	!
[1960]	712	!-- Unstable stratification
	713	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
	714	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	715	zu(nzb+1) * z0h1d )
[2337]	716	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
	717	LOG( (a-1.0_wp) / (a+1.0_wp) * &
	718	(b+1.0_wp) / (b-1.0_wp) )
	719	ENDIF
[1]	720	ELSE
[1353]	721	qs1d = 0.0_wp
[2337]	722	ENDIF
[1]	723
[1691]	724	ENDIF ! constant_flux_layer
[1]	725
	726	!
[2337]	727	!-- Compute the diabatic mixing length. The unstable stratification
	728	!-- must not be considered for l1d (km1d) as it is already considered
	729	!-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
	730	!-- too large.
[1]	731	IF ( mixing_length_1d == 'blackadar' ) THEN
	732	DO k = nzb+1, nzt
[1353]	733	IF ( rif1d(k) >= 0.0_wp ) THEN
	734	l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
[2337]	735	l1d_diss(k) = l1d(k)
[1]	736	ELSE
[2337]	737	l1d(k) = l_black(k)
	738	l1d_diss(k) = l_black(k) * &
	739	SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
[1]	740	ENDIF
	741	ENDDO
	742	ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN
	743	DO k = nzb+1, nzt
	744	dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
[1353]	745	IF ( dpt_dz > 0.0_wp ) THEN
	746	l_stable = 0.76_wp * SQRT( e1d(k) ) / &
	747	SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
[1]	748	ELSE
	749	l_stable = l_grid(k)
	750	ENDIF
	751	l1d(k) = MIN( l_grid(k), l_stable )
[2337]	752	l1d_diss(k) = l1d(k)
[1]	753	ENDDO
	754	ENDIF
	755
	756	!
	757	!-- Compute the diffusion coefficients for momentum via the
	758	!-- corresponding Prandtl-layer relationship and according to
[2337]	759	!-- Prandtl-Kolmogorov, respectively
[1691]	760	IF ( constant_flux_layer ) THEN
[1353]	761	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
	762	km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &
	763	( 1.0_wp + 5.0_wp * rif1d(nzb+1) )
[1]	764	ELSE
[1353]	765	km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &
	766	( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp
[1]	767	ENDIF
	768	ENDIF
	769	DO k = nzb_diff, nzt
[2337]	770	km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
[1]	771	ENDDO
	772
	773	!
	774	!-- Add damping layer
	775	DO k = damp_level_ind_1d+1, nzt+1
[1353]	776	km1d(k) = 1.1_wp * km1d(k-1)
[1346]	777	km1d(k) = MIN( km1d(k), 10.0_wp )
[1]	778	ENDDO
	779
	780	!
	781	!-- Compute the diffusion coefficient for heat via the relationship
	782	!-- kh = phim / phih * km
	783	DO k = nzb+1, nzt
[1353]	784	IF ( rif1d(k) >= 0.0_wp ) THEN
[1]	785	kh1d(k) = km1d(k)
	786	ELSE
[1353]	787	kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
[1]	788	ENDIF
	789	ENDDO
	790
	791	ENDIF ! .NOT. constant_diffusion
	792
	793	ENDDO ! intermediate step loop
	794
	795	!
	796	!-- Increment simulated time and output times
	797	current_timestep_number_1d = current_timestep_number_1d + 1
	798	simulated_time_1d = simulated_time_1d + dt_1d
	799	simulated_time_chr = time_to_string( simulated_time_1d )
	800	time_pr_1d = time_pr_1d + dt_1d
	801	time_run_control_1d = time_run_control_1d + dt_1d
	802
	803	!
	804	!-- Determine and print out quantities for run control
	805	IF ( time_run_control_1d >= dt_run_control_1d ) THEN
	806	CALL run_control_1d
	807	time_run_control_1d = time_run_control_1d - dt_run_control_1d
	808	ENDIF
	809
	810	!
	811	!-- Profile output on file
	812	IF ( time_pr_1d >= dt_pr_1d ) THEN
	813	CALL print_1d_model
	814	time_pr_1d = time_pr_1d - dt_pr_1d
	815	ENDIF
	816
	817	!
	818	!-- Determine size of next time step
	819	CALL timestep_1d
	820
	821	ENDDO ! time loop
	822
	823
	824	END SUBROUTINE time_integration_1d
	825
	826
	827	!------------------------------------------------------------------------------!
	828	! Description:
	829	! ------------
[1682]	830	!> Compute and print out quantities for run control of the 1D model.
[1]	831	!------------------------------------------------------------------------------!
[1682]	832
	833	SUBROUTINE run_control_1d
[1]	834
[1682]	835
[1320]	836	USE constants, &
	837	ONLY: pi
[1]	838
	839	IMPLICIT NONE
	840
[2338]	841	INTEGER(iwp) :: k !< loop index
[1320]	842
[2338]	843	REAL(wp) :: alpha !< angle of wind vector at top of constant-flux layer
	844	REAL(wp) :: energy !< kinetic energy
	845	REAL(wp) :: umax !< maximum of u
	846	REAL(wp) :: uv_total !< horizontal wind speed
	847	REAL(wp) :: vmax !< maximum of v
[1]	848
	849	!
	850	!-- Output
	851	IF ( myid == 0 ) THEN
	852	!
	853	!-- If necessary, write header
	854	IF ( .NOT. run_control_header_1d ) THEN
[184]	855	CALL check_open( 15 )
[1]	856	WRITE ( 15, 100 )
	857	run_control_header_1d = .TRUE.
	858	ENDIF
	859
	860	!
	861	!-- Compute control quantities
	862	!-- grid level nzp is excluded due to mirror boundary condition
[1353]	863	umax = 0.0_wp; vmax = 0.0_wp; energy = 0.0_wp
[1]	864	DO k = nzb+1, nzt+1
	865	umax = MAX( ABS( umax ), ABS( u1d(k) ) )
	866	vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
[1353]	867	energy = energy + 0.5_wp * ( u1d(k)2 + v1d(k)2 )
[1]	868	ENDDO
[1322]	869	energy = energy / REAL( nzt - nzb + 1, KIND=wp )
[1]	870
	871	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
[1691]	872	IF ( ABS( v1d(nzb+1) ) < 1.0E-5_wp ) THEN
[1346]	873	alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) )
[1]	874	ELSE
	875	alpha = ACOS( u1d(nzb+1) / uv_total )
[1353]	876	IF ( v1d(nzb+1) <= 0.0_wp ) alpha = 2.0_wp * pi - alpha
[1]	877	ENDIF
[1353]	878	alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
[1]	879
	880	WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &
	881	dt_1d, umax, vmax, us1d, alpha, energy
	882	!
	883	!-- Write buffer contents to disc immediately
[1808]	884	FLUSH( 15 )
[1]	885
	886	ENDIF
	887
	888	!
	889	!-- formats
[2299]	890	100 FORMAT (///'1D run control output:'/ &
[1]	891	&'------------------------------'// &
	892	&'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ &
	893	&'-------------------------------------------------------------')
[1697]	894	101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
[1]	895
	896
	897	END SUBROUTINE run_control_1d
	898
	899
	900
	901	!------------------------------------------------------------------------------!
	902	! Description:
	903	! ------------
[1682]	904	!> Compute the time step w.r.t. the diffusion criterion
[1]	905	!------------------------------------------------------------------------------!
[1682]	906
	907	SUBROUTINE timestep_1d
[1]	908
	909	IMPLICIT NONE
	910
[2338]	911	INTEGER(iwp) :: k !< loop index
[1]	912
[2338]	913	REAL(wp) :: dt_diff !< time step accorind to diffusion criterion
	914	REAL(wp) :: fac !< factor of criterion
	915	REAL(wp) :: value !< auxiliary variable
[1]	916
	917	!
	918	!-- Compute the currently feasible time step according to the diffusion
	919	!-- criterion. At nzb+1 the half grid length is used.
[1353]	920	fac = 0.35_wp
[1]	921	dt_diff = dt_max_1d
	922	DO k = nzb+2, nzt
[1353]	923	value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20_wp )
[1]	924	dt_diff = MIN( value, dt_diff )
	925	ENDDO
[1353]	926	value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20_wp )
[1]	927	dt_1d = MIN( value, dt_diff )
	928
	929	!
	930	!-- Set flag when the time step becomes too small
[1353]	931	IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) ) THEN
[1]	932	stop_dt_1d = .TRUE.
[254]	933
	934	WRITE( message_string, * ) 'timestep has exceeded the lower limit &', &
	935	'dt_1d = ',dt_1d,' s simulation stopped!'
	936	CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
	937
[1]	938	ENDIF
	939
	940	END SUBROUTINE timestep_1d
	941
	942
	943
	944	!------------------------------------------------------------------------------!
	945	! Description:
	946	! ------------
[1682]	947	!> List output of profiles from the 1D-model
[1]	948	!------------------------------------------------------------------------------!
[1682]	949
	950	SUBROUTINE print_1d_model
[1]	951
	952	IMPLICIT NONE
	953
[2338]	954	INTEGER(iwp) :: k !< loop parameter
[1]	955
[2338]	956	LOGICAL, SAVE :: write_first = .TRUE. !< flag for writing header
[1]	957
	958
	959	IF ( myid == 0 ) THEN
	960	!
	961	!-- Open list output file for profiles from the 1D-model
	962	CALL check_open( 17 )
	963
	964	!
	965	!-- Write Header
[2338]	966	IF ( write_first ) THEN
	967	WRITE ( 17, 100 ) TRIM( run_description_header )
	968	write_first = .FALSE.
	969	ENDIF
[1]	970
	971	!
	972	!-- Write the values
[2338]	973	WRITE ( 17, 104 ) TRIM( simulated_time_chr )
	974	WRITE ( 17, 101 )
[1]	975	WRITE ( 17, 102 )
	976	WRITE ( 17, 101 )
	977	DO k = nzt+1, nzb, -1
	978	WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
	979	rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
	980	ENDDO
	981	WRITE ( 17, 101 )
	982	WRITE ( 17, 102 )
	983	WRITE ( 17, 101 )
	984
	985	!
	986	!-- Write buffer contents to disc immediately
[1808]	987	FLUSH( 17 )
[1]	988
	989	ENDIF
	990
	991	!
	992	!-- Formats
[2338]	993	100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
	994	104 FORMAT (//'# Time: ',A)
	995	101 FORMAT ('#',79('-'))
	996	102 FORMAT ('# k zu u v pt e rif Km Kh ', &
[1]	997	'l zu k')
	998	103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
	999	1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
	1000
	1001
	1002	END SUBROUTINE print_1d_model
[2338]	1003
	1004
	1005	END MODULE

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |