Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |