Home

Context Navigation

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

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

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 4687 2020-09-21 19:40:16Z suehring $
[4687]	27	! Indoor model is now available during spinup
	28	!
	29	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	30	! Implementation of downward facing USM and LSM surfaces
	31	!
	32	! 4669 2020-09-09 13:43:47Z pavelkrc
[4669]	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
[4668]	38	! Improve debug messages during timestepping
[4655]	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
[4635]	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
[4631]	47	! Set pt1 attribute only element-wise
	48	!
	49	! 4540 2020-05-18 15:23:29Z raasch
[4540]	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	!
[4457]	55	! 4444 2020-03-05 15:59:50Z raasch
[4540]	56	! Bugfix: cpp-directives for serial mode added
	57	!
[4444]	58	! 4360 2020-01-07 11:25:50Z suehring
[4331]	59	! Enable output of diagnostic quantities, e.g. 2-m temperature
[4540]	60	!
[4331]	61	! 4227 2019-09-10 18:04:34Z gronemeier
[4540]	62	! Implement new palm_date_time_mod
	63	!
[4227]	64	! 4223 2019-09-10 09:20:47Z gronemeier
[4182]	65	! Corrected "Former revisions" section
[4540]	66	!
[4182]	67	! 4064 2019-07-01 05:33:33Z gronemeier
[4064]	68	! Moved call to radiation module out of intermediate time loop
[4540]	69	!
[4064]	70	! 4023 2019-06-12 13:20:01Z maronga
[4023]	71	! Time stamps are now negative in run control output
[4540]	72	!
[4023]	73	! 3885 2019-04-11 11:29:34Z kanani
[4540]	74	! Changes related to global restructuring of location messages and introduction of additional debug
	75	! messages
	76	!
[3885]	77	! 3766 2019-02-26 16:23:41Z raasch
[4540]	78	! Unused variable removed
	79	!
[3766]	80	! 3719 2019-02-06 13:10:18Z kanani
[4540]	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	!
[3719]	85	! 3655 2019-01-07 16:51:22Z knoop
[4540]	86	! Removed call to calculation of near air (10 cm) potential temperature (now in surface layer fluxes)
	87	!
[4182]	88	! 2296 2017-06-28 07:53:56Z maronga
	89	! Initial revision
[2296]	90	!
[4182]	91	!
[2296]	92	! Description:
	93	! ------------
[4540]	94	!> Integration in time of the non-atmospheric model components such as land surface model and urban
	95	!> surface model
	96	!--------------------------------------------------------------------------------------------------!
[2296]	97	SUBROUTINE time_integration_spinup
	98
[4540]	99	USE arrays_3d, &
	100	ONLY: pt, &
	101	pt_p, &
	102	u, &
	103	u_init, &
	104	v, &
	105	v_init
[2296]	106
[4540]	107	USE control_parameters, &
	108	ONLY: averaging_interval_pr, &
	109	constant_diffusion, &
	110	constant_flux_layer, &
	111	coupling_start_time, &
[4687]	112	current_timestep_number, &
[4540]	113	data_output_during_spinup, &
[4667]	114	debug_output_timestep, &
	115	debug_string, &
[4540]	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, &
[4687]	126	indoor_model, &
[4540]	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, &
[4654]	138	surface_output, &
[4540]	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
[2296]	148
[4540]	149	USE cpulog, &
	150	ONLY: cpu_log, &
	151	log_point_s
	152
	153	USE diagnostic_output_quantities_mod, &
[4331]	154	ONLY: doq_calculate
	155
[4540]	156	USE exchange_horiz_mod, &
[4457]	157	ONLY: exchange_horiz
	158
[4540]	159	USE indices, &
	160	ONLY: nbgp, &
	161	nzb, &
	162	nzt, &
	163	nysg, &
	164	nyng, &
	165	nxlg, &
	166	nxrg
[2296]	167
[4687]	168	USE indoor_model_mod, &
	169	ONLY: dt_indoor, im_main_heatcool, time_indoor
	170
[4540]	171	USE land_surface_model_mod, &
	172	ONLY: lsm_energy_balance, &
	173	lsm_swap_timelevel
[4227]	174
[2934]	175	USE pegrid
[2296]	176
[4444]	177	#if defined( __parallel )
[4540]	178	USE pmc_interface, &
[2934]	179	ONLY: nested_run
[4444]	180	#endif
[2934]	181
[2296]	182	USE kinds
	183
[4540]	184	USE palm_date_time_mod, &
	185	ONLY: get_date_time, &
	186	seconds_per_hour
[4227]	187
[4540]	188	USE radiation_model_mod, &
[4635]	189	ONLY: dt_radiation, &
	190	force_radiation_call, &
[4540]	191	radiation, &
	192	radiation_control, &
	193	radiation_interaction, &
	194	radiation_interactions, &
	195	time_radiation
[2296]	196
[4540]	197	USE statistics, &
[2296]	198	ONLY: flow_statistics_called
	199
[4654]	200	USE surface_data_output_mod, &
	201	ONLY: dt_dosurf, &
	202	skip_time_dosurf, &
	203	surface_data_output, &
	204	time_dosurf
	205
[4540]	206	USE surface_layer_fluxes_mod, &
[2296]	207	ONLY: surface_layer_fluxes
	208
[4540]	209	USE surface_mod, &
	210	ONLY : surf_lsm_h, &
	211	surf_lsm_v, surf_usm_h, &
[2296]	212	surf_usm_v
	213
[4540]	214	USE urban_surface_mod, &
[4671]	215	ONLY: usm_energy_balance, &
	216	usm_swap_timelevel
[2296]	217
	218
	219
	220
	221	IMPLICIT NONE
	222
[4540]	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 !<
[2299]	226
	227
[4540]	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
[4227]	242	REAL(wp) :: pt_spinup !< temporary storage of temperature
	243	REAL(wp) :: second_of_day !< second of the day
[4540]	244
[2728]	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
[2296]	248
[2728]	249
	250	!
	251	!-- Save 3D arrays because they are to be changed for spinup purpose
[2296]	252	ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[2728]	253	ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	254	ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[2296]	255
[4540]	256	CALL exchange_horiz( pt, nbgp )
	257	CALL exchange_horiz( u, nbgp )
	258	CALL exchange_horiz( v, nbgp )
	259
[2299]	260	pt_save = pt
[2728]	261	u_save = u
	262	v_save = v
[2296]	263
[2728]	264	!
[4540]	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.
[2782]	268	IF ( land_surface ) THEN
[4671]	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
[2782]	277	ENDDO
[2728]	278
[2782]	279	DO l = 0, 3
	280	DO m = 1, surf_lsm_v(l)%ns
[4540]	281	i = surf_lsm_v(l)%i(m)
[2782]	282	j = surf_lsm_v(l)%j(m)
	283	k = surf_lsm_v(l)%k(m)
[4540]	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)
[2782]	286	ENDDO
	287	ENDDO
	288	ENDIF
	289
	290	IF ( urban_surface ) THEN
[4671]	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
[2782]	299	ENDDO
	300
	301	DO l = 0, 3
	302	DO m = 1, surf_usm_v(l)%ns
[4540]	303	i = surf_usm_v(l)%i(m)
[2782]	304	j = surf_usm_v(l)%j(m)
	305	k = surf_usm_v(l)%k(m)
[4540]	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)
[2782]	308	ENDDO
	309	ENDDO
	310	ENDIF
	311
[4540]	312	CALL exchange_horiz( u, nbgp )
	313	CALL exchange_horiz( v, nbgp )
[2818]	314
[2723]	315	dt_save = dt_3d
	316	dt_3d = dt_spinup
	317
[3885]	318	CALL location_message( 'wall/soil spinup time-stepping', 'start' )
[2296]	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' )
[4540]	324
[4667]	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
[2296]	330	!
	331	!-- Start of intermediate step loop
	332	intermediate_timestep_count = 0
[4540]	333	DO WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
[2296]	334
	335	intermediate_timestep_count = intermediate_timestep_count + 1
	336
	337	!
[4540]	338	!-- Set the steering factors for the prognostic equations which depend on the timestep scheme
[2296]	339	CALL timestep_scheme_steering
	340
	341
[2299]	342	!
[4540]	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 )
[4227]	348
[4540]	349	pt_spinup = spinup_pt_mean + spinup_pt_amplitude * &
	350	solar_angle( day_of_year, second_of_day )
[2296]	351
[2299]	352	!
[4540]	353	!-- Map air temperature to all grid points in the vicinity of a surface element
[2296]	354	IF ( land_surface ) THEN
[4671]	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
[2296]	362	ENDDO
	363
	364	DO l = 0, 3
	365	DO m = 1, surf_lsm_v(l)%ns
[4540]	366	i = surf_lsm_v(l)%i(m)
[2296]	367	j = surf_lsm_v(l)%j(m)
	368	k = surf_lsm_v(l)%k(m)
[2299]	369	pt(k,j,i) = pt_spinup
[2296]	370	ENDDO
	371	ENDDO
	372	ENDIF
	373
	374	IF ( urban_surface ) THEN
[4671]	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
[2296]	385	ENDDO
	386
	387	DO l = 0, 3
	388	DO m = 1, surf_usm_v(l)%ns
[4540]	389	i = surf_usm_v(l)%i(m)
[2296]	390	j = surf_usm_v(l)%j(m)
	391	k = surf_usm_v(l)%k(m)
[2299]	392	pt(k,j,i) = pt_spinup
[3337]	393	!!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
[4631]	394	surf_usm_v(l)%pt1(m) = pt_spinup
[3337]	395	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2296]	396	ENDDO
	397	ENDDO
	398	ENDIF
	399
[4540]	400	CALL exchange_horiz( pt, nbgp )
[2818]	401
	402
[2296]	403	!
	404	!-- Swap the time levels in preparation for the next time step.
	405	timestep_count = timestep_count + 1
[4540]	406
[2296]	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
[4540]	416	CALL lsm_swap_timelevel ( MOD( timestep_count, 2 ) )
[2296]	417	ENDIF
	418
	419	IF ( urban_surface ) THEN
[4540]	420	CALL usm_swap_timelevel ( MOD( timestep_count, 2 ) )
[2296]	421	ENDIF
[4540]	422
[2296]	423	!
[4540]	424	!-- If required, compute virtual potential temperature
	425	IF ( humidity ) THEN
	426	CALL compute_vpt
	427	ENDIF
[2296]	428
	429	!
	430	!-- Compute the diffusion quantities
	431	IF ( .NOT. constant_diffusion ) THEN
	432
	433	!
[4540]	434	!-- First the vertical (and horizontal) fluxes in the surface (constant flux) layer are
	435	!-- computed
[2296]	436	IF ( constant_flux_layer ) THEN
	437	CALL surface_layer_fluxes
	438	ENDIF
	439
	440	!
[4540]	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
[2299]	444	IF ( land_surface ) THEN
[2296]	445
[4671]	446	CALL lsm_energy_balance( .TRUE. )
[3719]	447
[2296]	448	ENDIF
	449
	450	!
[4540]	451	!-- If required, solve the energy balance for urban surfaces and run the material heat model
[2296]	452	IF (urban_surface) THEN
[3719]	453
[4671]	454	CALL usm_energy_balance( .TRUE. )
[3719]	455
[2296]	456	ENDIF
	457
	458	ENDIF
	459
[4064]	460	ENDDO ! Intermediate step loop
	461
[2296]	462	!
[4064]	463	!-- If required, calculate radiative fluxes and heating rates
	464	IF ( radiation ) THEN
[2296]	465
[4635]	466	time_radiation = time_radiation + dt_3d
[2296]	467
[4635]	468	IF ( time_radiation >= dt_radiation .OR. force_radiation_call ) THEN
[2296]	469
[4064]	470	IF ( .NOT. force_radiation_call ) THEN
[4635]	471	time_radiation = time_radiation - dt_radiation
[4064]	472	ENDIF
[2296]	473
[4064]	474	CALL radiation_control
[2296]	475
[4064]	476	IF ( radiation_interactions ) THEN
	477	CALL radiation_interaction
[2296]	478	ENDIF
[4669]	479	!
	480	!-- Reset forcing of radiation call
	481	force_radiation_call = .FALSE.
	482
[2296]	483	ENDIF
[4064]	484	ENDIF
[2296]	485
	486	!
[4687]	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	!
[2296]	506	!-- Increase simulation time and output times
[2299]	507	current_timestep_number_spinup = current_timestep_number_spinup + 1
[4540]	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
[4023]	513	IF ( time_since_reference_point < 0.0_wp ) THEN
	514	sign_chr = '-'
	515	ELSE
	516	sign_chr = ' '
	517	ENDIF
[4540]	518
	519
[2296]	520	IF ( data_output_during_spinup ) THEN
[2723]	521	IF ( simulated_time >= skip_time_do2d_xy ) THEN
[4540]	522	time_do2d_xy = time_do2d_xy + dt_3d
[2723]	523	ENDIF
	524	IF ( simulated_time >= skip_time_do3d ) THEN
[4540]	525	time_do3d = time_do3d + dt_3d
[2723]	526	ENDIF
[4540]	527	time_dots = time_dots + dt_3d
[2296]	528	IF ( simulated_time >= skip_time_dopr ) THEN
[4540]	529	time_dopr = time_dopr + dt_3d
[2296]	530	ENDIF
[4654]	531	IF ( surface_output ) THEN
	532	IF ( simulated_time >= skip_time_dosurf ) THEN
	533	time_dosurf = time_dosurf + dt_3d
	534	ENDIF
	535	ENDIF
[4540]	536	time_run_control = time_run_control + dt_3d
[2296]	537
	538	!
	539	!-- Carry out statistical analysis and output at the requested output times.
[4540]	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.
[2296]	542
	543	!
[4540]	544	!-- Set a flag indicating that so far no statistics have been created for this time step
[2296]	545	flow_statistics_called = .FALSE.
	546
	547	!
	548	!-- If required, call flow_statistics for averaging in time
[4540]	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
[2723]	553	time_dopr_av = time_dopr_av + dt_3d
[2296]	554	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
	555	do_sum = .TRUE.
[4540]	556	time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
[2296]	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
[2723]	565	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
[4540]	566	time_dopr_av = 0.0_wp ! Due to averaging (see above)
[2296]	567	ENDIF
	568
	569	!
	570	!-- Output of time series
	571	IF ( time_dots >= dt_dots ) THEN
	572	CALL data_output_tseries
[2723]	573	time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
[2296]	574	ENDIF
	575
[2723]	576	!
	577	!-- 2d-data output (cross-sections)
	578	IF ( time_do2d_xy >= dt_do2d_xy ) THEN
[4331]	579	CALL doq_calculate
[2723]	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
[4331]	587	CALL doq_calculate
[2723]	588	CALL data_output_3d( 0 )
	589	time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
	590	ENDIF
[4654]	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
[2723]	599
[2296]	600	ENDIF
	601
	602	!
[4540]	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
[2299]	607	! CALL run_control
	608	! IF ( time_run_control >= dt_run_control ) THEN
[4540]	609	! time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d ) )
[2299]	610	! ENDIF
	611	! ENDIF
[2296]	612
	613	CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
	614
[2299]	615
	616	!
	617	!-- Run control output
[2296]	618	IF ( myid == 0 ) THEN
[2299]	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
[4540]	628	WRITE ( 15, 101 ) current_timestep_number_spinup, sign_chr, &
	629	time_since_reference_point_chr, dt_3d, pt_spinup
[2299]	630	!
	631	!-- Write buffer contents to disc immediately
	632	FLUSH( 15 )
[2296]	633	ENDIF
	634
[4667]	635	IF ( debug_output_timestep ) THEN
	636	WRITE( debug_string, * ) 'time_integration_spinup', simulated_time
	637	CALL debug_message( debug_string, 'end' )
	638	ENDIF
[2299]	639
	640
[4540]	641	ENDDO ! Time loop
[2296]	642
	643	!
[2728]	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
[2296]	649
[2723]	650	!
	651	!-- Reset time step
	652	dt_3d = dt_save
[4669]	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
[2723]	656
[2296]	657	DEALLOCATE(pt_save)
[2728]	658	DEALLOCATE(u_save)
	659	DEALLOCATE(v_save)
[2296]	660
[2934]	661	#if defined( __parallel )
	662	IF ( nested_run ) CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
	663	#endif
	664
[3885]	665	CALL location_message( 'wall/soil spinup time-stepping', 'finished' )
[2296]	666
[2299]	667
	668	!
	669	!-- Formats
[4540]	670	100 FORMAT (///'Spinup control output:---------------------------------'// &
	671	'ITER. HH:MM:SS DT PT(z_MO)---------------------------------')
[4023]	672	101 FORMAT (I5,2X,A1,A9,1X,F6.2,3X,F6.2,2X,F6.2)
[2299]	673
	674	CONTAINS
	675
	676	!
[4540]	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 )
[2299]	681
[4540]	682	USE basic_constants_and_equations_mod, &
[4331]	683	ONLY: pi
[4540]	684
[2299]	685	USE kinds
	686
[4540]	687	USE radiation_model_mod, &
	688	ONLY: decl_1, &
	689	decl_2, &
	690	decl_3, &
	691	lat, &
	692	lon
[2299]	693
[4540]	694	IMPLICIT NONE
[2299]	695
	696
[4227]	697	INTEGER(iwp), INTENT(IN) :: day_of_year !< day of the year
[2299]	698
[4540]	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
[2299]	703	!
[4540]	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
[2299]	707
	708	!
	709	!-- Calculate cosine of solar zenith angle
[4540]	710	solar_angle = SIN( lat ) * SIN( declination ) + COS( lat ) * COS( declination ) * &
	711	COS( hour_angle )
[2299]	712
	713	END FUNCTION solar_angle
	714
	715
[2296]	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 |