Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |