Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |