Home

Context Navigation

time_integration_spinup.f90 @ 4829

Last change on this file since 4829 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

Rev	Line
[2296]	1	!> @file time_integration_spinup.f90
[4540]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[2296]	4	!
[4540]	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.
[2296]	8	!
[4540]	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.
[2296]	12	!
[4540]	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/>.
[2296]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4540]	17	!--------------------------------------------------------------------------------------------------!
[2296]	18	!
[4540]	19	!
[2296]	20	! Current revisions:
[4540]	21	! -----------------
[4631]	22	!
[4668]	23	!
[2296]	24	! Former revisions:
	25	! -----------------
	26	! $Id: time_integration_spinup.f90 4828 2021-01-05 11:21:41Z suehring $
[4802]	27	! Bugfix in time-control of indoor model at the very first time step
	28	!
	29	! 4750 2020-10-16 14:27:48Z suehring
[4750]	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
[4687]	34	! Indoor model is now available during spinup
	35	!
	36	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	37	! Implementation of downward facing USM and LSM surfaces
	38	!
	39	! 4669 2020-09-09 13:43:47Z pavelkrc
[4669]	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
[4668]	45	! Improve debug messages during timestepping
[4655]	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
[4635]	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
[4631]	54	! Set pt1 attribute only element-wise
	55	!
	56	! 4540 2020-05-18 15:23:29Z raasch
[4540]	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	!
[4457]	62	! 4444 2020-03-05 15:59:50Z raasch
[4540]	63	! Bugfix: cpp-directives for serial mode added
	64	!
[4444]	65	! 4360 2020-01-07 11:25:50Z suehring
[4331]	66	! Enable output of diagnostic quantities, e.g. 2-m temperature
[4540]	67	!
[4331]	68	! 4227 2019-09-10 18:04:34Z gronemeier
[4540]	69	! Implement new palm_date_time_mod
	70	!
[4227]	71	! 4223 2019-09-10 09:20:47Z gronemeier
[4182]	72	! Corrected "Former revisions" section
[4540]	73	!
[4182]	74	! 4064 2019-07-01 05:33:33Z gronemeier
[4064]	75	! Moved call to radiation module out of intermediate time loop
[4540]	76	!
[4064]	77	! 4023 2019-06-12 13:20:01Z maronga
[4023]	78	! Time stamps are now negative in run control output
[4540]	79	!
[4023]	80	! 3885 2019-04-11 11:29:34Z kanani
[4540]	81	! Changes related to global restructuring of location messages and introduction of additional debug
	82	! messages
	83	!
[3885]	84	! 3766 2019-02-26 16:23:41Z raasch
[4540]	85	! Unused variable removed
	86	!
[3766]	87	! 3719 2019-02-06 13:10:18Z kanani
[4540]	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	!
[3719]	92	! 3655 2019-01-07 16:51:22Z knoop
[4540]	93	! Removed call to calculation of near air (10 cm) potential temperature (now in surface layer fluxes)
	94	!
[4182]	95	! 2296 2017-06-28 07:53:56Z maronga
	96	! Initial revision
[2296]	97	!
[4182]	98	!
[2296]	99	! Description:
	100	! ------------
[4540]	101	!> Integration in time of the non-atmospheric model components such as land surface model and urban
	102	!> surface model
	103	!--------------------------------------------------------------------------------------------------!
[2296]	104	SUBROUTINE time_integration_spinup
	105
[4540]	106	USE arrays_3d, &
	107	ONLY: pt, &
	108	pt_p, &
	109	u, &
	110	u_init, &
	111	v, &
	112	v_init
[2296]	113
[4540]	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, &
[4667]	120	debug_output_timestep, &
	121	debug_string, &
[4540]	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, &
[4687]	132	indoor_model, &
[4540]	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, &
[4654]	144	surface_output, &
[4540]	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
[2296]	154
[4540]	155	USE cpulog, &
	156	ONLY: cpu_log, &
	157	log_point_s
	158
	159	USE diagnostic_output_quantities_mod, &
[4331]	160	ONLY: doq_calculate
	161
[4540]	162	USE exchange_horiz_mod, &
[4457]	163	ONLY: exchange_horiz
	164
[4540]	165	USE indices, &
	166	ONLY: nbgp, &
	167	nzb, &
	168	nzt, &
	169	nysg, &
	170	nyng, &
	171	nxlg, &
	172	nxrg
[2296]	173
[4687]	174	USE indoor_model_mod, &
[4750]	175	ONLY: dt_indoor, &
	176	im_main_heatcool, &
	177	indoor_during_spinup, &
	178	time_indoor
[4687]	179
[4540]	180	USE land_surface_model_mod, &
	181	ONLY: lsm_energy_balance, &
	182	lsm_swap_timelevel
[4227]	183
[2934]	184	USE pegrid
[2296]	185
[4444]	186	#if defined( __parallel )
[4540]	187	USE pmc_interface, &
[2934]	188	ONLY: nested_run
[4444]	189	#endif
[2934]	190
[2296]	191	USE kinds
	192
[4540]	193	USE palm_date_time_mod, &
	194	ONLY: get_date_time, &
	195	seconds_per_hour
[4227]	196
[4540]	197	USE radiation_model_mod, &
[4635]	198	ONLY: dt_radiation, &
	199	force_radiation_call, &
[4540]	200	radiation, &
	201	radiation_control, &
	202	radiation_interaction, &
	203	radiation_interactions, &
	204	time_radiation
[2296]	205
[4540]	206	USE statistics, &
[2296]	207	ONLY: flow_statistics_called
	208
[4654]	209	USE surface_data_output_mod, &
	210	ONLY: dt_dosurf, &
	211	skip_time_dosurf, &
	212	surface_data_output, &
	213	time_dosurf
	214
[4540]	215	USE surface_layer_fluxes_mod, &
[2296]	216	ONLY: surface_layer_fluxes
	217
[4540]	218	USE surface_mod, &
	219	ONLY : surf_lsm_h, &
	220	surf_lsm_v, surf_usm_h, &
[2296]	221	surf_usm_v
	222
[4540]	223	USE urban_surface_mod, &
[4671]	224	ONLY: usm_energy_balance, &
	225	usm_swap_timelevel
[2296]	226
	227
	228
	229
	230	IMPLICIT NONE
	231
[4540]	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 !<
[2299]	235
	236
[4540]	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
[4227]	251	REAL(wp) :: pt_spinup !< temporary storage of temperature
	252	REAL(wp) :: second_of_day !< second of the day
[4540]	253
[2728]	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
[2296]	257
[2728]	258
	259	!
	260	!-- Save 3D arrays because they are to be changed for spinup purpose
[2296]	261	ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[2728]	262	ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	263	ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[2296]	264
[4540]	265	CALL exchange_horiz( pt, nbgp )
	266	CALL exchange_horiz( u, nbgp )
	267	CALL exchange_horiz( v, nbgp )
	268
[2299]	269	pt_save = pt
[2728]	270	u_save = u
	271	v_save = v
[2296]	272
[2728]	273	!
[4540]	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.
[2782]	277	IF ( land_surface ) THEN
[4671]	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
[2782]	286	ENDDO
[2728]	287
[2782]	288	DO l = 0, 3
	289	DO m = 1, surf_lsm_v(l)%ns
[4540]	290	i = surf_lsm_v(l)%i(m)
[2782]	291	j = surf_lsm_v(l)%j(m)
	292	k = surf_lsm_v(l)%k(m)
[4540]	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)
[2782]	295	ENDDO
	296	ENDDO
	297	ENDIF
	298
	299	IF ( urban_surface ) THEN
[4671]	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
[2782]	308	ENDDO
	309
	310	DO l = 0, 3
	311	DO m = 1, surf_usm_v(l)%ns
[4540]	312	i = surf_usm_v(l)%i(m)
[2782]	313	j = surf_usm_v(l)%j(m)
	314	k = surf_usm_v(l)%k(m)
[4540]	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)
[2782]	317	ENDDO
	318	ENDDO
	319	ENDIF
	320
[4540]	321	CALL exchange_horiz( u, nbgp )
	322	CALL exchange_horiz( v, nbgp )
[2818]	323
[2723]	324	dt_save = dt_3d
	325	dt_3d = dt_spinup
	326
[3885]	327	CALL location_message( 'wall/soil spinup time-stepping', 'start' )
[2296]	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' )
[4540]	333
[4667]	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
[2296]	339	!
	340	!-- Start of intermediate step loop
	341	intermediate_timestep_count = 0
[4540]	342	DO WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
[2296]	343
	344	intermediate_timestep_count = intermediate_timestep_count + 1
	345
	346	!
[4540]	347	!-- Set the steering factors for the prognostic equations which depend on the timestep scheme
[2296]	348	CALL timestep_scheme_steering
	349
	350
[2299]	351	!
[4540]	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 )
[4227]	357
[4540]	358	pt_spinup = spinup_pt_mean + spinup_pt_amplitude * &
	359	solar_angle( day_of_year, second_of_day )
[2296]	360
[2299]	361	!
[4540]	362	!-- Map air temperature to all grid points in the vicinity of a surface element
[2296]	363	IF ( land_surface ) THEN
[4671]	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
[2296]	371	ENDDO
	372
	373	DO l = 0, 3
	374	DO m = 1, surf_lsm_v(l)%ns
[4540]	375	i = surf_lsm_v(l)%i(m)
[2296]	376	j = surf_lsm_v(l)%j(m)
	377	k = surf_lsm_v(l)%k(m)
[2299]	378	pt(k,j,i) = pt_spinup
[2296]	379	ENDDO
	380	ENDDO
	381	ENDIF
	382
	383	IF ( urban_surface ) THEN
[4671]	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
[2296]	394	ENDDO
	395
	396	DO l = 0, 3
	397	DO m = 1, surf_usm_v(l)%ns
[4540]	398	i = surf_usm_v(l)%i(m)
[2296]	399	j = surf_usm_v(l)%j(m)
	400	k = surf_usm_v(l)%k(m)
[2299]	401	pt(k,j,i) = pt_spinup
[3337]	402	!!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
[4631]	403	surf_usm_v(l)%pt1(m) = pt_spinup
[3337]	404	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2296]	405	ENDDO
	406	ENDDO
	407	ENDIF
	408
[4540]	409	CALL exchange_horiz( pt, nbgp )
[2818]	410
	411
[2296]	412	!
	413	!-- Swap the time levels in preparation for the next time step.
	414	timestep_count = timestep_count + 1
[4540]	415
[2296]	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
[4540]	425	CALL lsm_swap_timelevel ( MOD( timestep_count, 2 ) )
[2296]	426	ENDIF
	427
	428	IF ( urban_surface ) THEN
[4540]	429	CALL usm_swap_timelevel ( MOD( timestep_count, 2 ) )
[2296]	430	ENDIF
[4540]	431
[2296]	432	!
[4540]	433	!-- If required, compute virtual potential temperature
	434	IF ( humidity ) THEN
	435	CALL compute_vpt
	436	ENDIF
[2296]	437
	438	!
	439	!-- Compute the diffusion quantities
	440	IF ( .NOT. constant_diffusion ) THEN
	441
	442	!
[4540]	443	!-- First the vertical (and horizontal) fluxes in the surface (constant flux) layer are
	444	!-- computed
[2296]	445	IF ( constant_flux_layer ) THEN
	446	CALL surface_layer_fluxes
	447	ENDIF
	448
	449	!
[4540]	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
[2299]	453	IF ( land_surface ) THEN
[2296]	454
[4671]	455	CALL lsm_energy_balance( .TRUE. )
[3719]	456
[2296]	457	ENDIF
	458
	459	!
[4540]	460	!-- If required, solve the energy balance for urban surfaces and run the material heat model
[2296]	461	IF (urban_surface) THEN
[3719]	462
[4671]	463	CALL usm_energy_balance( .TRUE. )
[3719]	464
[2296]	465	ENDIF
	466
	467	ENDIF
	468
[4064]	469	ENDDO ! Intermediate step loop
	470
[2296]	471	!
[4064]	472	!-- If required, calculate radiative fluxes and heating rates
	473	IF ( radiation ) THEN
[2296]	474
[4635]	475	time_radiation = time_radiation + dt_3d
[2296]	476
[4635]	477	IF ( time_radiation >= dt_radiation .OR. force_radiation_call ) THEN
[2296]	478
[4064]	479	IF ( .NOT. force_radiation_call ) THEN
[4635]	480	time_radiation = time_radiation - dt_radiation
[4064]	481	ENDIF
[2296]	482
[4064]	483	CALL radiation_control
[2296]	484
[4064]	485	IF ( radiation_interactions ) THEN
	486	CALL radiation_interaction
[2296]	487	ENDIF
[4669]	488	!
	489	!-- Reset forcing of radiation call
	490	force_radiation_call = .FALSE.
	491
[2296]	492	ENDIF
[4064]	493	ENDIF
[2296]	494
	495	!
[4687]	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.
[4750]	501	IF ( indoor_model .AND. indoor_during_spinup ) THEN
[4687]	502
	503	time_indoor = time_indoor + dt_3d
	504
[4750]	505	IF ( time_indoor >= dt_indoor .OR. current_timestep_number_spinup == 0 ) THEN
[4687]	506
[4802]	507	IF ( time_indoor >= dt_indoor ) time_indoor = time_indoor - dt_indoor
[4687]	508
	509	CALL im_main_heatcool
	510
	511	ENDIF
	512	ENDIF
	513
	514	!
[2296]	515	!-- Increase simulation time and output times
[2299]	516	current_timestep_number_spinup = current_timestep_number_spinup + 1
[4540]	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
[4023]	522	IF ( time_since_reference_point < 0.0_wp ) THEN
	523	sign_chr = '-'
	524	ELSE
	525	sign_chr = ' '
	526	ENDIF
[4540]	527
	528
[2296]	529	IF ( data_output_during_spinup ) THEN
[2723]	530	IF ( simulated_time >= skip_time_do2d_xy ) THEN
[4540]	531	time_do2d_xy = time_do2d_xy + dt_3d
[2723]	532	ENDIF
	533	IF ( simulated_time >= skip_time_do3d ) THEN
[4540]	534	time_do3d = time_do3d + dt_3d
[2723]	535	ENDIF
[4540]	536	time_dots = time_dots + dt_3d
[2296]	537	IF ( simulated_time >= skip_time_dopr ) THEN
[4540]	538	time_dopr = time_dopr + dt_3d
[2296]	539	ENDIF
[4654]	540	IF ( surface_output ) THEN
	541	IF ( simulated_time >= skip_time_dosurf ) THEN
	542	time_dosurf = time_dosurf + dt_3d
	543	ENDIF
	544	ENDIF
[4540]	545	time_run_control = time_run_control + dt_3d
[2296]	546
	547	!
	548	!-- Carry out statistical analysis and output at the requested output times.
[4540]	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.
[2296]	551
	552	!
[4540]	553	!-- Set a flag indicating that so far no statistics have been created for this time step
[2296]	554	flow_statistics_called = .FALSE.
	555
	556	!
	557	!-- If required, call flow_statistics for averaging in time
[4540]	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
[2723]	562	time_dopr_av = time_dopr_av + dt_3d
[2296]	563	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
	564	do_sum = .TRUE.
[4540]	565	time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
[2296]	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
[2723]	574	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
[4540]	575	time_dopr_av = 0.0_wp ! Due to averaging (see above)
[2296]	576	ENDIF
	577
	578	!
	579	!-- Output of time series
	580	IF ( time_dots >= dt_dots ) THEN
	581	CALL data_output_tseries
[2723]	582	time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
[2296]	583	ENDIF
	584
[2723]	585	!
	586	!-- 2d-data output (cross-sections)
	587	IF ( time_do2d_xy >= dt_do2d_xy ) THEN
[4331]	588	CALL doq_calculate
[2723]	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
[4331]	596	CALL doq_calculate
[2723]	597	CALL data_output_3d( 0 )
	598	time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
	599	ENDIF
[4654]	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
[2723]	608
[2296]	609	ENDIF
	610
	611	!
[4540]	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
[2299]	616	! CALL run_control
	617	! IF ( time_run_control >= dt_run_control ) THEN
[4540]	618	! time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d ) )
[2299]	619	! ENDIF
	620	! ENDIF
[2296]	621
	622	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
	623
[2299]	624
	625	!
	626	!-- Run control output
[2296]	627	IF ( myid == 0 ) THEN
[2299]	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
[4540]	637	WRITE ( 15, 101 ) current_timestep_number_spinup, sign_chr, &
	638	time_since_reference_point_chr, dt_3d, pt_spinup
[2299]	639	!
	640	!-- Write buffer contents to disc immediately
	641	FLUSH( 15 )
[2296]	642	ENDIF
	643
[4667]	644	IF ( debug_output_timestep ) THEN
	645	WRITE( debug_string, * ) 'time_integration_spinup', simulated_time
	646	CALL debug_message( debug_string, 'end' )
	647	ENDIF
[2299]	648
	649
[4540]	650	ENDDO ! Time loop
[2296]	651
	652	!
[2728]	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
[2296]	658
[2723]	659	!
	660	!-- Reset time step
	661	dt_3d = dt_save
[4669]	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
[2723]	665
[2296]	666	DEALLOCATE(pt_save)
[2728]	667	DEALLOCATE(u_save)
	668	DEALLOCATE(v_save)
[2296]	669
[2934]	670	#if defined( __parallel )
	671	IF ( nested_run ) CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
	672	#endif
	673
[3885]	674	CALL location_message( 'wall/soil spinup time-stepping', 'finished' )
[2296]	675
[2299]	676
	677	!
	678	!-- Formats
[4540]	679	100 FORMAT (///'Spinup control output:---------------------------------'// &
	680	'ITER. HH:MM:SS DT PT(z_MO)---------------------------------')
[4023]	681	101 FORMAT (I5,2X,A1,A9,1X,F6.2,3X,F6.2,2X,F6.2)
[2299]	682
	683	CONTAINS
	684
	685	!
[4540]	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 )
[2299]	690
[4540]	691	USE basic_constants_and_equations_mod, &
[4331]	692	ONLY: pi
[4540]	693
[2299]	694	USE kinds
	695
[4540]	696	USE radiation_model_mod, &
	697	ONLY: decl_1, &
	698	decl_2, &
	699	decl_3, &
	700	lat, &
	701	lon
[2299]	702
[4540]	703	IMPLICIT NONE
[2299]	704
	705
[4227]	706	INTEGER(iwp), INTENT(IN) :: day_of_year !< day of the year
[2299]	707
[4540]	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
[2299]	712	!
[4540]	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
[2299]	716
	717	!
	718	!-- Calculate cosine of solar zenith angle
[4540]	719	solar_angle = SIN( lat ) * SIN( declination ) + COS( lat ) * COS( declination ) * &
	720	COS( hour_angle )
[2299]	721
	722	END FUNCTION solar_angle
	723
	724
[2296]	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 @ 4829

Download in other formats: