Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |