Home

Context Navigation

time_integration_spinup.f90 @ 4867

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

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

Context Navigation

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

Download in other formats: