Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |