Home

Context Navigation

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

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

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