!> @file time_integration_spinup.f90 !------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2017 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: time_integration_spinup.f90 2297 2017-06-28 14:35:57Z scharf $ ! bugfixes ! ! 2296 2017-06-28 07:53:56Z maronga ! Initial revision ! ! ! ! ! Description: ! ------------ !> Integration in time of the non-atmospheric model components such as land !> surface model and urban surface model !------------------------------------------------------------------------------! SUBROUTINE time_integration_spinup USE arrays_3d, & ONLY: pt, pt_p USE control_parameters, & ONLY: averaging_interval_pr, constant_diffusion, constant_flux_layer, & coupling_start_time, current_timestep_number, & data_output_during_spinup, disturbance_created, dopr_n, do_sum, & dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control, & dt_spinup, humidity, intermediate_timestep_count, & intermediate_timestep_count_max, land_surface, & nr_timesteps_this_run, simulated_time, simulated_time_chr, & skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean, & spinup_time, timestep_count, timestep_scheme, time_dopr, & time_dopr_av, time_dots, time_run_control, & time_since_reference_point, urban_surface USE constants, & ONLY: pi USE cpulog, & ONLY: cpu_log, log_point, log_point_s USE indices, & ONLY: nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg USE land_surface_model_mod, & ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel, & skip_time_do_lsm USE pegrid USE kinds USE radiation_model_mod, & ONLY: dt_radiation, force_radiation_call, radiation, & radiation_control, skip_time_do_radiation, time_radiation USE statistics, & ONLY: flow_statistics_called USE surface_layer_fluxes_mod, & ONLY: surface_layer_fluxes USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v USE urban_surface_mod, & ONLY: usm_material_heat_model, usm_material_model, & usm_radiation, usm_surface_energy_balance, usm_swap_timelevel IMPLICIT NONE CHARACTER (LEN=9) :: time_to_string !< INTEGER(iwp) :: i INTEGER(iwp) :: j INTEGER(iwp) :: k INTEGER(iwp) :: l INTEGER(iwp) :: m REAL(wp) :: pt_spinup !< temporary storage of temperature REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_save !< temporary storage of temperature REAL(wp) :: day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL exchange_horiz( pt_p, nbgp ) pt_save = pt_p CALL location_message( 'starting spinup-sequence', .TRUE. ) ! !-- Start of the time loop DO WHILE ( simulated_time < spinup_time ) CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' ) ! !-- Start of intermediate step loop intermediate_timestep_count = 0 DO WHILE ( intermediate_timestep_count < & intermediate_timestep_count_max ) intermediate_timestep_count = intermediate_timestep_count + 1 ! !-- Set the steering factors for the prognostic equations which depend !-- on the timestep scheme CALL timestep_scheme_steering !!!! Set new values of pt_p here pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi) IF ( land_surface ) THEN DO m = 1, surf_lsm_h%ns i = surf_lsm_h%i(m) j = surf_lsm_h%j(m) k = surf_lsm_h%k(m) pt_p(k,j,i) = pt_spinup ENDDO DO l = 0, 3 DO m = 1, surf_lsm_v(l)%ns i = surf_lsm_v(l)%i(m) j = surf_lsm_v(l)%j(m) k = surf_lsm_v(l)%k(m) pt_p(k,j,i) = pt_spinup ENDDO ENDDO ENDIF IF ( urban_surface ) THEN DO m = 1, surf_usm_h%ns i = surf_usm_h%i(m) j = surf_usm_h%j(m) k = surf_usm_h%k(m) pt_p(k,j,i) = pt_spinup ENDDO DO l = 0, 3 DO m = 1, surf_usm_v(l)%ns i = surf_usm_v(l)%i(m) j = surf_usm_v(l)%j(m) k = surf_usm_v(l)%k(m) pt_p(k,j,i) = pt_spinup ENDDO ENDDO ENDIF ! !-- Swap the time levels in preparation for the next time step. timestep_count = timestep_count + 1 IF ( land_surface ) THEN CALL lsm_swap_timelevel ( 0 ) ENDIF IF ( urban_surface ) THEN CALL usm_swap_timelevel ( 0 ) ENDIF IF ( land_surface ) THEN CALL lsm_swap_timelevel ( MOD( timestep_count, 2) ) ENDIF IF ( urban_surface ) THEN CALL usm_swap_timelevel ( MOD( timestep_count, 2) ) ENDIF ! !-- If required, compute virtual potential temperature IF ( humidity ) THEN CALL compute_vpt ENDIF ! !-- Compute the diffusion quantities IF ( .NOT. constant_diffusion ) THEN ! !-- First the vertical (and horizontal) fluxes in the surface !-- (constant flux) layer are computed IF ( constant_flux_layer ) THEN CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' ) CALL surface_layer_fluxes CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' ) ENDIF ! !-- If required, solve the energy balance for the surface and run soil !-- model. Call for horizontal as well as vertical surfaces IF ( land_surface .AND. simulated_time > skip_time_do_lsm) THEN CALL cpu_log( log_point(54), 'land_surface', 'start' ) ! !-- Call for horizontal upward-facing surfaces CALL lsm_energy_balance( .TRUE., -1 ) CALL lsm_soil_model( .TRUE., -1 ) ! !-- Call for northward-facing surfaces CALL lsm_energy_balance( .FALSE., 0 ) CALL lsm_soil_model( .FALSE., 0 ) ! !-- Call for southward-facing surfaces CALL lsm_energy_balance( .FALSE., 1 ) CALL lsm_soil_model( .FALSE., 1 ) ! !-- Call for eastward-facing surfaces CALL lsm_energy_balance( .FALSE., 2 ) CALL lsm_soil_model( .FALSE., 2 ) ! !-- Call for westward-facing surfaces CALL lsm_energy_balance( .FALSE., 3 ) CALL lsm_soil_model( .FALSE., 3 ) CALL cpu_log( log_point(54), 'land_surface', 'stop' ) ENDIF ! !-- If required, solve the energy balance for urban surfaces and run !-- the material heat model IF (urban_surface) THEN CALL cpu_log( log_point(74), 'urban_surface', 'start' ) CALL usm_surface_energy_balance IF ( usm_material_model ) THEN CALL usm_material_heat_model ENDIF CALL cpu_log( log_point(74), 'urban_surface', 'stop' ) ENDIF ENDIF ! !-- If required, calculate radiative fluxes and heating rates IF ( radiation .AND. intermediate_timestep_count & == intermediate_timestep_count_max .AND. simulated_time > & skip_time_do_radiation ) THEN time_radiation = time_radiation + dt_spinup IF ( time_radiation >= dt_radiation .OR. force_radiation_call ) & THEN CALL cpu_log( log_point(50), 'radiation', 'start' ) IF ( .NOT. force_radiation_call ) THEN time_radiation = time_radiation - dt_radiation ENDIF CALL radiation_control CALL cpu_log( log_point(50), 'radiation', 'stop' ) IF (urban_surface) THEN CALL cpu_log( log_point(75), 'usm_radiation', 'start' ) CALL usm_radiation CALL cpu_log( log_point(75), 'usm_radiation', 'stop' ) ENDIF ENDIF ENDIF ENDDO ! Intermediate step loop ! !-- Increase simulation time and output times nr_timesteps_this_run = nr_timesteps_this_run + 1 current_timestep_number = current_timestep_number + 1 simulated_time = simulated_time + dt_spinup simulated_time_chr = time_to_string( simulated_time ) time_since_reference_point = simulated_time - coupling_start_time IF ( data_output_during_spinup ) THEN time_dots = time_dots + dt_spinup IF ( simulated_time >= skip_time_dopr ) THEN time_dopr = time_dopr + dt_spinup ENDIF time_run_control = time_run_control + dt_spinup ! !-- Carry out statistical analysis and output at the requested output times. !-- The MOD function is used for calculating the output time counters (like !-- time_dopr) in order to regard a possible decrease of the output time !-- interval in case of restart runs ! !-- Set a flag indicating that so far no statistics have been created !-- for this time step flow_statistics_called = .FALSE. ! !-- If required, call flow_statistics for averaging in time IF ( averaging_interval_pr /= 0.0_wp .AND. & ( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. & simulated_time >= skip_time_dopr ) THEN time_dopr_av = time_dopr_av + dt_spinup IF ( time_dopr_av >= dt_averaging_input_pr ) THEN do_sum = .TRUE. time_dopr_av = MOD( time_dopr_av, & MAX( dt_averaging_input_pr, dt_spinup ) ) ENDIF ENDIF IF ( do_sum ) CALL flow_statistics ! !-- Output of profiles IF ( time_dopr >= dt_dopr ) THEN IF ( dopr_n /= 0 ) CALL data_output_profiles time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) ) time_dopr_av = 0.0_wp ! due to averaging (see above) ENDIF ! !-- Output of time series IF ( time_dots >= dt_dots ) THEN CALL data_output_tseries time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) ) ENDIF ENDIF ! !-- Computation and output of run control parameters. !-- This is also done whenever perturbations have been imposed IF ( time_run_control >= dt_run_control .OR. & timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) & THEN CALL run_control IF ( time_run_control >= dt_run_control ) THEN time_run_control = MOD( time_run_control, & MAX( dt_run_control, dt_spinup ) ) ENDIF ENDIF CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' ) IF ( myid == 0 ) THEN PRINT*, time_since_reference_point ENDIF ENDDO ! time loop ! !-- Write back saved temperature to the 3D arrays pt(:,:,:) = pt_save pt_p(:,:,:) = pt_save DEALLOCATE(pt_save) CALL location_message( 'finished time-stepping spinup', .TRUE. ) END SUBROUTINE time_integration_spinup