Home

Context Navigation

source: palm/trunk/SOURCE/time_integration_spinup.f90 @ 2297

Last change on this file since 2297 was 2297, checked in by scharf, 7 years ago
bugfixes
Property svn:keywords set to `Id`
File size: 13.4 KB

Line
1	!> @file time_integration_spinup.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
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.
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	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: time_integration_spinup.f90 2297 2017-06-28 14:35:57Z scharf $
27	! bugfixes
28	!
29	! 2296 2017-06-28 07:53:56Z maronga
30	! Initial revision
31	!
32	!
33	!
34	!
35	! Description:
36	! ------------
37	!> Integration in time of the non-atmospheric model components such as land
38	!> surface model and urban surface model
39	!------------------------------------------------------------------------------!
40	SUBROUTINE time_integration_spinup
41
42	USE arrays_3d, &
43	ONLY: pt, pt_p
44
45	USE control_parameters, &
46	ONLY: averaging_interval_pr, constant_diffusion, constant_flux_layer, &
47	coupling_start_time, current_timestep_number, &
48	data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
49	dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control, &
50	dt_spinup, humidity, intermediate_timestep_count, &
51	intermediate_timestep_count_max, land_surface, &
52	nr_timesteps_this_run, simulated_time, simulated_time_chr, &
53	skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean, &
54	spinup_time, timestep_count, timestep_scheme, time_dopr, &
55	time_dopr_av, time_dots, time_run_control, &
56	time_since_reference_point, urban_surface
57
58	USE constants, &
59	ONLY: pi
60
61	USE cpulog, &
62	ONLY: cpu_log, log_point, log_point_s
63
64	USE indices, &
65	ONLY: nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
66
67
68	USE land_surface_model_mod, &
69	ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel, &
70	skip_time_do_lsm
71
72	USE pegrid
73
74	USE kinds
75
76	USE radiation_model_mod, &
77	ONLY: dt_radiation, force_radiation_call, radiation, &
78	radiation_control, skip_time_do_radiation, time_radiation
79
80	USE statistics, &
81	ONLY: flow_statistics_called
82
83	USE surface_layer_fluxes_mod, &
84	ONLY: surface_layer_fluxes
85
86	USE surface_mod, &
87	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
88	surf_usm_v
89
90	USE urban_surface_mod, &
91	ONLY: usm_material_heat_model, usm_material_model, &
92	usm_radiation, usm_surface_energy_balance, usm_swap_timelevel
93
94
95
96
97	IMPLICIT NONE
98
99	CHARACTER (LEN=9) :: time_to_string !<
100
101	INTEGER(iwp) :: i
102	INTEGER(iwp) :: j
103	INTEGER(iwp) :: k
104	INTEGER(iwp) :: l
105	INTEGER(iwp) :: m
106
107	REAL(wp) :: pt_spinup !< temporary storage of temperature
108
109	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_save !< temporary storage of temperature
110
111	REAL(wp) :: day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
112
113	ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
114
115	CALL exchange_horiz( pt_p, nbgp )
116	pt_save = pt_p
117
118	CALL location_message( 'starting spinup-sequence', .TRUE. )
119	!
120	!-- Start of the time loop
121	DO WHILE ( simulated_time < spinup_time )
122
123	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
124
125	!
126	!-- Start of intermediate step loop
127	intermediate_timestep_count = 0
128	DO WHILE ( intermediate_timestep_count < &
129	intermediate_timestep_count_max )
130
131	intermediate_timestep_count = intermediate_timestep_count + 1
132
133	!
134	!-- Set the steering factors for the prognostic equations which depend
135	!-- on the timestep scheme
136	CALL timestep_scheme_steering
137
138
139	!!!! Set new values of pt_p here
140	pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
141
142	IF ( land_surface ) THEN
143	DO m = 1, surf_lsm_h%ns
144	i = surf_lsm_h%i(m)
145	j = surf_lsm_h%j(m)
146	k = surf_lsm_h%k(m)
147	pt_p(k,j,i) = pt_spinup
148	ENDDO
149
150	DO l = 0, 3
151	DO m = 1, surf_lsm_v(l)%ns
152	i = surf_lsm_v(l)%i(m)
153	j = surf_lsm_v(l)%j(m)
154	k = surf_lsm_v(l)%k(m)
155	pt_p(k,j,i) = pt_spinup
156	ENDDO
157	ENDDO
158	ENDIF
159
160	IF ( urban_surface ) THEN
161	DO m = 1, surf_usm_h%ns
162	i = surf_usm_h%i(m)
163	j = surf_usm_h%j(m)
164	k = surf_usm_h%k(m)
165	pt_p(k,j,i) = pt_spinup
166	ENDDO
167
168	DO l = 0, 3
169	DO m = 1, surf_usm_v(l)%ns
170	i = surf_usm_v(l)%i(m)
171	j = surf_usm_v(l)%j(m)
172	k = surf_usm_v(l)%k(m)
173	pt_p(k,j,i) = pt_spinup
174	ENDDO
175	ENDDO
176	ENDIF
177
178	!
179	!-- Swap the time levels in preparation for the next time step.
180	timestep_count = timestep_count + 1
181
182	IF ( land_surface ) THEN
183	CALL lsm_swap_timelevel ( 0 )
184	ENDIF
185
186	IF ( urban_surface ) THEN
187	CALL usm_swap_timelevel ( 0 )
188	ENDIF
189
190	IF ( land_surface ) THEN
191	CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
192	ENDIF
193
194	IF ( urban_surface ) THEN
195	CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
196	ENDIF
197
198	!
199	!-- If required, compute virtual potential temperature
200	IF ( humidity ) THEN
201	CALL compute_vpt
202	ENDIF
203
204	!
205	!-- Compute the diffusion quantities
206	IF ( .NOT. constant_diffusion ) THEN
207
208	!
209	!-- First the vertical (and horizontal) fluxes in the surface
210	!-- (constant flux) layer are computed
211	IF ( constant_flux_layer ) THEN
212	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
213	CALL surface_layer_fluxes
214	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
215	ENDIF
216
217	!
218	!-- If required, solve the energy balance for the surface and run soil
219	!-- model. Call for horizontal as well as vertical surfaces
220	IF ( land_surface .AND. simulated_time > skip_time_do_lsm) THEN
221
222	CALL cpu_log( log_point(54), 'land_surface', 'start' )
223	!
224	!-- Call for horizontal upward-facing surfaces
225	CALL lsm_energy_balance( .TRUE., -1 )
226	CALL lsm_soil_model( .TRUE., -1 )
227	!
228	!-- Call for northward-facing surfaces
229	CALL lsm_energy_balance( .FALSE., 0 )
230	CALL lsm_soil_model( .FALSE., 0 )
231	!
232	!-- Call for southward-facing surfaces
233	CALL lsm_energy_balance( .FALSE., 1 )
234	CALL lsm_soil_model( .FALSE., 1 )
235	!
236	!-- Call for eastward-facing surfaces
237	CALL lsm_energy_balance( .FALSE., 2 )
238	CALL lsm_soil_model( .FALSE., 2 )
239	!
240	!-- Call for westward-facing surfaces
241	CALL lsm_energy_balance( .FALSE., 3 )
242	CALL lsm_soil_model( .FALSE., 3 )
243
244	CALL cpu_log( log_point(54), 'land_surface', 'stop' )
245	ENDIF
246
247	!
248	!-- If required, solve the energy balance for urban surfaces and run
249	!-- the material heat model
250	IF (urban_surface) THEN
251	CALL cpu_log( log_point(74), 'urban_surface', 'start' )
252	CALL usm_surface_energy_balance
253	IF ( usm_material_model ) THEN
254	CALL usm_material_heat_model
255	ENDIF
256	CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
257	ENDIF
258
259	ENDIF
260
261	!
262	!-- If required, calculate radiative fluxes and heating rates
263	IF ( radiation .AND. intermediate_timestep_count &
264	== intermediate_timestep_count_max .AND. simulated_time > &
265	skip_time_do_radiation ) THEN
266
267	time_radiation = time_radiation + dt_spinup
268
269	IF ( time_radiation >= dt_radiation .OR. force_radiation_call ) &
270	THEN
271
272	CALL cpu_log( log_point(50), 'radiation', 'start' )
273
274	IF ( .NOT. force_radiation_call ) THEN
275	time_radiation = time_radiation - dt_radiation
276	ENDIF
277
278	CALL radiation_control
279
280	CALL cpu_log( log_point(50), 'radiation', 'stop' )
281
282	IF (urban_surface) THEN
283	CALL cpu_log( log_point(75), 'usm_radiation', 'start' )
284	CALL usm_radiation
285	CALL cpu_log( log_point(75), 'usm_radiation', 'stop' )
286	ENDIF
287	ENDIF
288	ENDIF
289
290	ENDDO ! Intermediate step loop
291
292	!
293	!-- Increase simulation time and output times
294	nr_timesteps_this_run = nr_timesteps_this_run + 1
295	current_timestep_number = current_timestep_number + 1
296	simulated_time = simulated_time + dt_spinup
297	simulated_time_chr = time_to_string( simulated_time )
298	time_since_reference_point = simulated_time - coupling_start_time
299
300	IF ( data_output_during_spinup ) THEN
301	time_dots = time_dots + dt_spinup
302	IF ( simulated_time >= skip_time_dopr ) THEN
303	time_dopr = time_dopr + dt_spinup
304	ENDIF
305	time_run_control = time_run_control + dt_spinup
306
307	!
308	!-- Carry out statistical analysis and output at the requested output times.
309	!-- The MOD function is used for calculating the output time counters (like
310	!-- time_dopr) in order to regard a possible decrease of the output time
311	!-- interval in case of restart runs
312
313	!
314	!-- Set a flag indicating that so far no statistics have been created
315	!-- for this time step
316	flow_statistics_called = .FALSE.
317
318	!
319	!-- If required, call flow_statistics for averaging in time
320	IF ( averaging_interval_pr /= 0.0_wp .AND. &
321	( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. &
322	simulated_time >= skip_time_dopr ) THEN
323	time_dopr_av = time_dopr_av + dt_spinup
324	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
325	do_sum = .TRUE.
326	time_dopr_av = MOD( time_dopr_av, &
327	MAX( dt_averaging_input_pr, dt_spinup ) )
328	ENDIF
329	ENDIF
330	IF ( do_sum ) CALL flow_statistics
331
332	!
333	!-- Output of profiles
334	IF ( time_dopr >= dt_dopr ) THEN
335	IF ( dopr_n /= 0 ) CALL data_output_profiles
336	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
337	time_dopr_av = 0.0_wp ! due to averaging (see above)
338	ENDIF
339
340	!
341	!-- Output of time series
342	IF ( time_dots >= dt_dots ) THEN
343	CALL data_output_tseries
344	time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
345	ENDIF
346
347	ENDIF
348
349	!
350	!-- Computation and output of run control parameters.
351	!-- This is also done whenever perturbations have been imposed
352	IF ( time_run_control >= dt_run_control .OR. &
353	timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) &
354	THEN
355	CALL run_control
356	IF ( time_run_control >= dt_run_control ) THEN
357	time_run_control = MOD( time_run_control, &
358	MAX( dt_run_control, dt_spinup ) )
359	ENDIF
360	ENDIF
361
362	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
363
364	IF ( myid == 0 ) THEN
365	PRINT*, time_since_reference_point
366	ENDIF
367
368	ENDDO ! time loop
369
370	!
371	!-- Write back saved temperature to the 3D arrays
372	pt(:,:,:) = pt_save
373	pt_p(:,:,:) = pt_save
374
375	DEALLOCATE(pt_save)
376
377	CALL location_message( 'finished time-stepping spinup', .TRUE. )
378
379	END SUBROUTINE time_integration_spinup

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |