Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |