Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |