Home

Context Navigation

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

Last change on this file since 2774 was 2758, checked in by suehring, 7 years ago
Bugfix, no initial masking of wind velocity at first prognostic grid level in case of land- or urban-surface spin-up.
Property svn:keywords set to `Id`
File size: 19.6 KB

Line
1	!> @file time_integration_spinup.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: time_integration_spinup.f90 2758 2018-01-17 12:55:21Z knoop $
27	! Comment out homogeneous setting of wind velocity as this will lead to zero
28	! friction velocity and cause problems in MOST relationships.
29	!
30	! 2728 2018-01-09 07:03:53Z maronga
31	! Set velocity componenets to homogeneous values during spinup
32	!
33	! 2724 2018-01-05 12:12:38Z maronga
34	! Use dt_spinup for all active components during spinup
35	!
36	! 2723 2018-01-05 09:27:03Z maronga
37	! Bugfix: array rad_sw_in no longer exists and is thus removed from RUN_CONTROL
38	! output.
39	! Added output of XY and 3D data during spinup.
40	! Bugfix: time step in LSM and USM was set to dt_3d instead of dt_spinup
41	!
42	! 2718 2018-01-02 08:49:38Z maronga
43	! Corrected "Former revisions" section
44	!
45	! 2696 2017-12-14 17:12:51Z kanani
46	! Change in file header (GPL part)
47	! Added radiation interactions (moved from USM) (MS)
48	!
49	! 2544 2017-10-13 18:09:32Z maronga
50	! Date and time quantities are now read from date_and_time_mod
51	!
52	! 2299 2017-06-29 10:14:38Z maronga
53	! Call of soil model adjusted to avoid prognostic equation for soil moisture
54	! during spinup.
55	! Better representation of diurnal cycle of near-surface temperature.
56	! Excluded prognostic equation for soil moisture during spinup.
57	! Added output of run control data for spinup.
58	!
59	! 2297 2017-06-28 14:35:57Z scharf
60	! bugfixes
61	!
62	! 2296 2017-06-28 07:53:56Z maronga
63	! Initial revision
64	!
65	!
66	! Description:
67	! ------------
68	!> Integration in time of the non-atmospheric model components such as land
69	!> surface model and urban surface model
70	!------------------------------------------------------------------------------!
71	SUBROUTINE time_integration_spinup
72
73	USE arrays_3d, &
74	ONLY: pt, pt_p, u, v
75
76	USE control_parameters, &
77	ONLY: averaging_interval_pr, constant_diffusion, constant_flux_layer, &
78	coupling_start_time, current_timestep_number, &
79	data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
80	dt_averaging_input_pr, dt_dopr, dt_dots, dt_do2d_xy, dt_do3d, &
81	dt_run_control, dt_spinup, dt_3d, humidity, &
82	intermediate_timestep_count, &
83	intermediate_timestep_count_max, land_surface, &
84	simulated_time, simulated_time_chr, &
85	skip_time_dopr, skip_time_do2d_xy, skip_time_do3d, spinup, &
86	spinup_pt_amplitude, spinup_pt_mean, spinup_time, &
87	timestep_count, timestep_scheme, time_dopr, time_dopr_av, &
88	time_dots, time_do2d_xy, time_do3d, time_run_control, &
89	time_since_reference_point, ug_surface, vg_surface, urban_surface
90
91	USE constants, &
92	ONLY: pi
93
94	USE cpulog, &
95	ONLY: cpu_log, log_point, log_point_s
96
97	USE date_and_time_mod, &
98	ONLY: day_of_year_init, time_utc_init
99
100	USE indices, &
101	ONLY: nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
102
103
104	USE land_surface_model_mod, &
105	ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
106
107	USE pegrid, &
108	ONLY: myid
109
110	USE kinds
111
112	USE radiation_model_mod, &
113	ONLY: force_radiation_call, radiation, &
114	radiation_control, rad_sw_in, time_radiation, &
115	radiation_interaction, radiation_interactions
116
117	USE statistics, &
118	ONLY: flow_statistics_called
119
120	USE surface_layer_fluxes_mod, &
121	ONLY: surface_layer_fluxes
122
123	USE surface_mod, &
124	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
125	surf_usm_v
126
127	USE urban_surface_mod, &
128	ONLY: usm_material_heat_model, usm_material_model, &
129	usm_surface_energy_balance, usm_swap_timelevel, &
130	usm_green_heat_model, usm_temperature_near_surface
131
132
133
134
135	IMPLICIT NONE
136
137	CHARACTER (LEN=9) :: time_to_string !<
138
139	INTEGER(iwp) :: i !< running index
140	INTEGER(iwp) :: j !< running index
141	INTEGER(iwp) :: k !< running index
142	INTEGER(iwp) :: l !< running index
143	INTEGER(iwp) :: m !< running index
144
145	INTEGER(iwp) :: current_timestep_number_spinup = 0 !< number if timestep during spinup
146
147	LOGICAL :: run_control_header_spinup = .FALSE. !< flag parameter for steering whether the header information must be output
148
149	REAL(wp) :: pt_spinup !< temporary storage of temperature
150	REAL(wp) :: dt_save !< temporary storage for time step
151
152	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_save !< temporary storage of temperature
153	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_save !< temporary storage of u wind component
154	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_save !< temporary storage of v wind component
155
156
157	!
158	!-- Save 3D arrays because they are to be changed for spinup purpose
159	ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
160	ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
161	ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
162
163	CALL exchange_horiz( pt, nbgp )
164	CALL exchange_horiz( u, nbgp )
165	CALL exchange_horiz( v, nbgp )
166
167	pt_save = pt
168	u_save = u
169	v_save = v
170
171	!
172	!-- Set the same wall-adjacent velocity to all grid points. The sign of the
173	!-- original velocity field must be preserved because the surface schemes crash
174	!-- otherwise. The precise reason is still unknown. A minimum velocity of 0.1
175	!-- m/s is used to maintain turbulent transfer at the surface.
176	! u = SIGN(1.0_wp,u) * MAX(ug_surface,0.1_wp)
177	! v = SIGN(1.0_wp,v) * MAX(vg_surface,0.1_wp)
178
179	dt_save = dt_3d
180	dt_3d = dt_spinup
181
182	CALL location_message( 'starting spinup-sequence', .TRUE. )
183	!
184	!-- Start of the time loop
185	DO WHILE ( simulated_time < spinup_time )
186
187	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
188
189	!
190	!-- Start of intermediate step loop
191	intermediate_timestep_count = 0
192	DO WHILE ( intermediate_timestep_count < &
193	intermediate_timestep_count_max )
194
195	intermediate_timestep_count = intermediate_timestep_count + 1
196
197	!
198	!-- Set the steering factors for the prognostic equations which depend
199	!-- on the timestep scheme
200	CALL timestep_scheme_steering
201
202
203	!
204	!-- Estimate a near-surface air temperature based on the position of the
205	!-- sun and user input about mean temperature and amplitude. The time is
206	!-- shifted by one hour to simulate a lag between air temperature and
207	!-- incoming radiation
208	pt_spinup = spinup_pt_mean + spinup_pt_amplitude &
209	* solar_angle (time_utc_init + time_since_reference_point - 3600.0)
210
211	!
212	!-- Map air temperature to all grid points in the vicinity of a surface
213	!-- element
214	IF ( land_surface ) THEN
215	DO m = 1, surf_lsm_h%ns
216	i = surf_lsm_h%i(m)
217	j = surf_lsm_h%j(m)
218	k = surf_lsm_h%k(m)
219	pt(k,j,i) = pt_spinup
220	ENDDO
221
222	DO l = 0, 3
223	DO m = 1, surf_lsm_v(l)%ns
224	i = surf_lsm_v(l)%i(m)
225	j = surf_lsm_v(l)%j(m)
226	k = surf_lsm_v(l)%k(m)
227	pt(k,j,i) = pt_spinup
228	ENDDO
229	ENDDO
230	ENDIF
231
232	IF ( urban_surface ) THEN
233	DO m = 1, surf_usm_h%ns
234	i = surf_usm_h%i(m)
235	j = surf_usm_h%j(m)
236	k = surf_usm_h%k(m)
237	pt(k,j,i) = pt_spinup
238	ENDDO
239
240	DO l = 0, 3
241	DO m = 1, surf_usm_v(l)%ns
242	i = surf_usm_v(l)%i(m)
243	j = surf_usm_v(l)%j(m)
244	k = surf_usm_v(l)%k(m)
245	pt(k,j,i) = pt_spinup
246	ENDDO
247	ENDDO
248	ENDIF
249
250	!
251	!-- Swap the time levels in preparation for the next time step.
252	timestep_count = timestep_count + 1
253
254	IF ( land_surface ) THEN
255	CALL lsm_swap_timelevel ( 0 )
256	ENDIF
257
258	IF ( urban_surface ) THEN
259	CALL usm_swap_timelevel ( 0 )
260	ENDIF
261
262	IF ( land_surface ) THEN
263	CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
264	ENDIF
265
266	IF ( urban_surface ) THEN
267	CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
268	ENDIF
269
270	!
271	!-- If required, compute virtual potential temperature
272	IF ( humidity ) THEN
273	CALL compute_vpt
274	ENDIF
275
276	!
277	!-- Compute the diffusion quantities
278	IF ( .NOT. constant_diffusion ) THEN
279
280	!
281	!-- First the vertical (and horizontal) fluxes in the surface
282	!-- (constant flux) layer are computed
283	IF ( constant_flux_layer ) THEN
284	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
285	CALL surface_layer_fluxes
286	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
287	ENDIF
288
289	!
290	!-- If required, solve the energy balance for the surface and run soil
291	!-- model. Call for horizontal as well as vertical surfaces.
292	!-- The prognostic equation for soil moisure is switched off
293	IF ( land_surface ) THEN
294
295	CALL cpu_log( log_point(54), 'land_surface', 'start' )
296	!
297	!-- Call for horizontal upward-facing surfaces
298	CALL lsm_energy_balance( .TRUE., -1 )
299	CALL lsm_soil_model( .TRUE., -1, .FALSE. )
300	!
301	!-- Call for northward-facing surfaces
302	CALL lsm_energy_balance( .FALSE., 0 )
303	CALL lsm_soil_model( .FALSE., 0, .FALSE. )
304	!
305	!-- Call for southward-facing surfaces
306	CALL lsm_energy_balance( .FALSE., 1 )
307	CALL lsm_soil_model( .FALSE., 1, .FALSE. )
308	!
309	!-- Call for eastward-facing surfaces
310	CALL lsm_energy_balance( .FALSE., 2 )
311	CALL lsm_soil_model( .FALSE., 2, .FALSE. )
312	!
313	!-- Call for westward-facing surfaces
314	CALL lsm_energy_balance( .FALSE., 3 )
315	CALL lsm_soil_model( .FALSE., 3, .FALSE. )
316
317	CALL cpu_log( log_point(54), 'land_surface', 'stop' )
318	ENDIF
319
320	!
321	!-- If required, solve the energy balance for urban surfaces and run
322	!-- the material heat model
323	IF (urban_surface) THEN
324	CALL cpu_log( log_point(74), 'urban_surface', 'start' )
325	CALL usm_surface_energy_balance
326	IF ( usm_material_model ) THEN
327	CALL usm_green_heat_model
328	CALL usm_material_heat_model
329	ENDIF
330	IF ( urban_surface ) THEN
331	CALL usm_temperature_near_surface
332	ENDIF
333	CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
334	ENDIF
335
336	ENDIF
337
338	!
339	!-- If required, calculate radiative fluxes and heating rates
340	IF ( radiation .AND. intermediate_timestep_count &
341	== intermediate_timestep_count_max ) THEN
342
343	time_radiation = time_radiation + dt_3d
344
345	IF ( time_radiation >= dt_3d .OR. force_radiation_call ) &
346	THEN
347
348	CALL cpu_log( log_point(50), 'radiation', 'start' )
349
350	IF ( .NOT. force_radiation_call ) THEN
351	time_radiation = time_radiation - dt_3d
352	ENDIF
353
354	CALL radiation_control
355
356	CALL cpu_log( log_point(50), 'radiation', 'stop' )
357
358	IF ( radiation_interactions ) THEN
359	CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
360	CALL radiation_interaction
361	CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
362	ENDIF
363	ENDIF
364	ENDIF
365
366	ENDDO ! Intermediate step loop
367
368	!
369	!-- Increase simulation time and output times
370	current_timestep_number_spinup = current_timestep_number_spinup + 1
371	simulated_time = simulated_time + dt_3d
372	simulated_time_chr = time_to_string( simulated_time )
373	time_since_reference_point = simulated_time - coupling_start_time
374
375	IF ( data_output_during_spinup ) THEN
376	IF ( simulated_time >= skip_time_do2d_xy ) THEN
377	time_do2d_xy = time_do2d_xy + dt_3d
378	ENDIF
379	IF ( simulated_time >= skip_time_do3d ) THEN
380	time_do3d = time_do3d + dt_3d
381	ENDIF
382	time_dots = time_dots + dt_3d
383	IF ( simulated_time >= skip_time_dopr ) THEN
384	time_dopr = time_dopr + dt_3d
385	ENDIF
386	time_run_control = time_run_control + dt_3d
387
388	!
389	!-- Carry out statistical analysis and output at the requested output times.
390	!-- The MOD function is used for calculating the output time counters (like
391	!-- time_dopr) in order to regard a possible decrease of the output time
392	!-- interval in case of restart runs
393
394	!
395	!-- Set a flag indicating that so far no statistics have been created
396	!-- for this time step
397	flow_statistics_called = .FALSE.
398
399	!
400	!-- If required, call flow_statistics for averaging in time
401	IF ( averaging_interval_pr /= 0.0_wp .AND. &
402	( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. &
403	simulated_time >= skip_time_dopr ) THEN
404	time_dopr_av = time_dopr_av + dt_3d
405	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
406	do_sum = .TRUE.
407	time_dopr_av = MOD( time_dopr_av, &
408	MAX( dt_averaging_input_pr, dt_3d ) )
409	ENDIF
410	ENDIF
411	IF ( do_sum ) CALL flow_statistics
412
413	!
414	!-- Output of profiles
415	IF ( time_dopr >= dt_dopr ) THEN
416	IF ( dopr_n /= 0 ) CALL data_output_profiles
417	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
418	time_dopr_av = 0.0_wp ! due to averaging (see above)
419	ENDIF
420
421	!
422	!-- Output of time series
423	IF ( time_dots >= dt_dots ) THEN
424	CALL data_output_tseries
425	time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
426	ENDIF
427
428	!
429	!-- 2d-data output (cross-sections)
430	IF ( time_do2d_xy >= dt_do2d_xy ) THEN
431	CALL data_output_2d( 'xy', 0 )
432	time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
433	ENDIF
434
435	!
436	!-- 3d-data output (volume data)
437	IF ( time_do3d >= dt_do3d ) THEN
438	CALL data_output_3d( 0 )
439	time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
440	ENDIF
441
442
443	ENDIF
444
445	!
446	!-- Computation and output of run control parameters.
447	!-- This is also done whenever perturbations have been imposed
448	! IF ( time_run_control >= dt_run_control .OR. &
449	! timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) &
450	! THEN
451	! CALL run_control
452	! IF ( time_run_control >= dt_run_control ) THEN
453	! time_run_control = MOD( time_run_control, &
454	! MAX( dt_run_control, dt_3d ) )
455	! ENDIF
456	! ENDIF
457
458	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
459
460
461	!
462	!-- Run control output
463	IF ( myid == 0 ) THEN
464	!
465	!-- If necessary, write header
466	IF ( .NOT. run_control_header_spinup ) THEN
467	CALL check_open( 15 )
468	WRITE ( 15, 100 )
469	run_control_header_spinup = .TRUE.
470	ENDIF
471	!
472	!-- Write some general information about the spinup in run control file
473	WRITE ( 15, 101 ) current_timestep_number_spinup, simulated_time_chr, dt_3d, pt_spinup
474	!
475	!-- Write buffer contents to disc immediately
476	FLUSH( 15 )
477	ENDIF
478
479
480
481	ENDDO ! time loop
482
483	!
484	!-- Write back saved arrays to the 3D arrays
485	pt = pt_save
486	pt_p = pt_save
487	u = u_save
488	v = v_save
489
490	!
491	!-- Reset time step
492	dt_3d = dt_save
493
494	DEALLOCATE(pt_save)
495	DEALLOCATE(u_save)
496	DEALLOCATE(v_save)
497
498	CALL location_message( 'finished spinup-sequence', .TRUE. )
499
500
501	!
502	!-- Formats
503	100 FORMAT (///'Spinup control output:'/ &
504	'--------------------------------'// &
505	'ITER. HH:MM:SS DT PT(z_MO)'/ &
506	'--------------------------------')
507	101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
508
509	CONTAINS
510
511	!
512	!-- Returns the cosine of the solar zenith angle at a given time. This routine
513	!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
514	FUNCTION solar_angle( local_time )
515
516	USE constants, &
517	ONLY: pi
518
519	USE kinds
520
521	USE radiation_model_mod, &
522	ONLY: decl_1, decl_2, decl_3, lat, lon
523
524	IMPLICIT NONE
525
526
527	REAL(wp) :: solar_angle !< cosine of the solar zenith angle
528
529	REAL(wp) :: day !< day of the year
530	REAL(wp) :: declination !< solar declination angle
531	REAL(wp) :: hour_angle !< solar hour angle
532	REAL(wp) :: time_utc !< current time in UTC
533	REAL(wp), INTENT(IN) :: local_time
534	!
535	!-- Calculate current day and time based on the initial values and simulation
536	!-- time
537
538	day = day_of_year_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
539	time_utc = MOD(local_time, 86400.0_wp)
540
541
542	!
543	!-- Calculate solar declination and hour angle
544	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
545	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
546
547	!
548	!-- Calculate cosine of solar zenith angle
549	solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
550	* COS(hour_angle)
551
552
553	END FUNCTION solar_angle
554
555
556	END SUBROUTINE time_integration_spinup

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |