Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |