Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |