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

Last change on this file since 2297 was 2297, checked in by scharf, 7 years ago

bugfixes

  • Property svn:keywords set to Id
File size: 13.4 KB
Line 
1!> @file time_integration_spinup.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration_spinup.f90 2297 2017-06-28 14:35:57Z scharf $
27! bugfixes
28!
29! 2296 2017-06-28 07:53:56Z maronga
30! Initial revision
31!
32!
33!
34!
35! Description:
36! ------------
37!> Integration in time of the non-atmospheric model components such as land
38!> surface model and urban surface model
39!------------------------------------------------------------------------------!
40 SUBROUTINE time_integration_spinup
41 
42    USE arrays_3d,                                                             &
43        ONLY:  pt, pt_p
44
45    USE control_parameters,                                                    &
46        ONLY:  averaging_interval_pr, constant_diffusion, constant_flux_layer, &
47               coupling_start_time, current_timestep_number,                   &
48               data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
49               dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control,        &
50               dt_spinup, humidity, intermediate_timestep_count,               &
51               intermediate_timestep_count_max, land_surface,                  &
52               nr_timesteps_this_run, simulated_time, simulated_time_chr,      &
53               skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
54               spinup_time, timestep_count, timestep_scheme, time_dopr,        &
55               time_dopr_av, time_dots, time_run_control,                      &
56               time_since_reference_point, urban_surface
57
58    USE constants,                                                             &
59        ONLY:  pi
60
61    USE cpulog,                                                                &
62        ONLY:  cpu_log, log_point, log_point_s
63
64    USE indices,                                                               &
65        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
66
67
68    USE land_surface_model_mod,                                                &
69        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
70               skip_time_do_lsm
71
72    USE pegrid
73
74    USE kinds
75
76    USE radiation_model_mod,                                                   &
77        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
78               radiation_control, skip_time_do_radiation, time_radiation
79
80    USE statistics,                                                            &
81        ONLY:  flow_statistics_called
82
83    USE surface_layer_fluxes_mod,                                              &
84        ONLY:  surface_layer_fluxes
85
86    USE surface_mod,                                                           &
87        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
88                surf_usm_v
89
90    USE urban_surface_mod,                                                     &
91        ONLY:  usm_material_heat_model, usm_material_model,                    &
92               usm_radiation, usm_surface_energy_balance, usm_swap_timelevel
93
94
95
96
97    IMPLICIT NONE
98
99    CHARACTER (LEN=9) ::  time_to_string          !<
100 
101    INTEGER(iwp) :: i
102    INTEGER(iwp) :: j
103    INTEGER(iwp) :: k
104    INTEGER(iwp) :: l
105    INTEGER(iwp) :: m
106 
107    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
108                 
109    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
110
111    REAL(wp) ::  day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
112
113    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
114
115    CALL exchange_horiz( pt_p, nbgp )   
116    pt_save = pt_p
117
118    CALL location_message( 'starting spinup-sequence', .TRUE. )
119!
120!-- Start of the time loop
121    DO  WHILE ( simulated_time < spinup_time )
122
123       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
124   
125!
126!--    Start of intermediate step loop
127       intermediate_timestep_count = 0
128       DO  WHILE ( intermediate_timestep_count < &
129                   intermediate_timestep_count_max )
130
131          intermediate_timestep_count = intermediate_timestep_count + 1
132
133!
134!--       Set the steering factors for the prognostic equations which depend
135!--       on the timestep scheme
136          CALL timestep_scheme_steering
137
138
139!!!!      Set new values of pt_p here
140          pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
141
142          IF ( land_surface )  THEN
143             DO  m = 1, surf_lsm_h%ns
144                i   = surf_lsm_h%i(m)           
145                j   = surf_lsm_h%j(m)
146                k   = surf_lsm_h%k(m)
147                pt_p(k,j,i) = pt_spinup
148             ENDDO
149
150             DO  l = 0, 3
151                DO  m = 1, surf_lsm_v(l)%ns
152                   i   = surf_lsm_v(l)%i(m)           
153                   j   = surf_lsm_v(l)%j(m)
154                   k   = surf_lsm_v(l)%k(m)
155                   pt_p(k,j,i) = pt_spinup
156                ENDDO
157             ENDDO
158          ENDIF
159
160          IF ( urban_surface )  THEN
161             DO  m = 1, surf_usm_h%ns
162                i   = surf_usm_h%i(m)           
163                j   = surf_usm_h%j(m)
164                k   = surf_usm_h%k(m)
165                pt_p(k,j,i) = pt_spinup
166             ENDDO
167
168             DO  l = 0, 3
169                DO  m = 1, surf_usm_v(l)%ns
170                   i   = surf_usm_v(l)%i(m)           
171                   j   = surf_usm_v(l)%j(m)
172                   k   = surf_usm_v(l)%k(m)
173                   pt_p(k,j,i) = pt_spinup
174                ENDDO
175             ENDDO
176          ENDIF
177
178!
179!--       Swap the time levels in preparation for the next time step.
180          timestep_count = timestep_count + 1
181     
182          IF ( land_surface )  THEN
183              CALL lsm_swap_timelevel ( 0 )
184          ENDIF
185
186          IF ( urban_surface )  THEN
187             CALL usm_swap_timelevel ( 0 )
188          ENDIF
189
190          IF ( land_surface )  THEN
191             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
192          ENDIF
193
194          IF ( urban_surface )  THEN
195             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
196          ENDIF
197         
198!
199!--       If required, compute virtual potential temperature
200          IF ( humidity )  THEN
201             CALL compute_vpt 
202          ENDIF 
203
204!
205!--       Compute the diffusion quantities
206          IF ( .NOT. constant_diffusion )  THEN
207
208!
209!--          First the vertical (and horizontal) fluxes in the surface
210!--          (constant flux) layer are computed
211             IF ( constant_flux_layer )  THEN
212                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
213                CALL surface_layer_fluxes
214                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
215             ENDIF
216
217!
218!--          If required, solve the energy balance for the surface and run soil
219!--          model. Call for horizontal as well as vertical surfaces
220             IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
221
222                CALL cpu_log( log_point(54), 'land_surface', 'start' )
223!
224!--             Call for horizontal upward-facing surfaces
225                CALL lsm_energy_balance( .TRUE., -1 )
226                CALL lsm_soil_model( .TRUE., -1 )
227!
228!--             Call for northward-facing surfaces
229                CALL lsm_energy_balance( .FALSE., 0 )
230                CALL lsm_soil_model( .FALSE., 0 )
231!
232!--             Call for southward-facing surfaces
233                CALL lsm_energy_balance( .FALSE., 1 )
234                CALL lsm_soil_model( .FALSE., 1 )
235!
236!--             Call for eastward-facing surfaces
237                CALL lsm_energy_balance( .FALSE., 2 )
238                CALL lsm_soil_model( .FALSE., 2 )
239!
240!--             Call for westward-facing surfaces
241                CALL lsm_energy_balance( .FALSE., 3 )
242                CALL lsm_soil_model( .FALSE., 3 )
243
244                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
245             ENDIF
246
247!
248!--          If required, solve the energy balance for urban surfaces and run
249!--          the material heat model
250             IF (urban_surface) THEN
251                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
252                CALL usm_surface_energy_balance
253                IF ( usm_material_model )  THEN
254                   CALL usm_material_heat_model
255                ENDIF
256                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
257             ENDIF
258
259          ENDIF
260
261!
262!--       If required, calculate radiative fluxes and heating rates
263          IF ( radiation .AND. intermediate_timestep_count                     &
264               == intermediate_timestep_count_max .AND. simulated_time >       &
265               skip_time_do_radiation )  THEN
266
267               time_radiation = time_radiation + dt_spinup
268
269             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
270             THEN
271
272                CALL cpu_log( log_point(50), 'radiation', 'start' )
273
274                IF ( .NOT. force_radiation_call )  THEN
275                   time_radiation = time_radiation - dt_radiation
276                ENDIF
277
278                CALL radiation_control
279
280                CALL cpu_log( log_point(50), 'radiation', 'stop' )
281
282                IF (urban_surface)  THEN
283                   CALL cpu_log( log_point(75), 'usm_radiation', 'start' )
284                   CALL usm_radiation
285                   CALL cpu_log( log_point(75), 'usm_radiation', 'stop' )
286                ENDIF
287             ENDIF
288          ENDIF
289
290       ENDDO   ! Intermediate step loop
291
292!
293!--    Increase simulation time and output times
294       nr_timesteps_this_run      = nr_timesteps_this_run + 1
295       current_timestep_number    = current_timestep_number + 1
296       simulated_time             = simulated_time   + dt_spinup
297       simulated_time_chr         = time_to_string( simulated_time )
298       time_since_reference_point = simulated_time - coupling_start_time
299
300       IF ( data_output_during_spinup )  THEN
301          time_dots          = time_dots        + dt_spinup
302          IF ( simulated_time >= skip_time_dopr )  THEN
303             time_dopr       = time_dopr        + dt_spinup
304          ENDIF
305          time_run_control   = time_run_control + dt_spinup
306
307!
308!--       Carry out statistical analysis and output at the requested output times.
309!--       The MOD function is used for calculating the output time counters (like
310!--       time_dopr) in order to regard a possible decrease of the output time
311!--       interval in case of restart runs
312
313!
314!--       Set a flag indicating that so far no statistics have been created
315!--       for this time step
316          flow_statistics_called = .FALSE.
317
318!
319!--       If required, call flow_statistics for averaging in time
320          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
321             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
322             simulated_time >= skip_time_dopr )  THEN
323             time_dopr_av = time_dopr_av + dt_spinup
324             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
325                do_sum = .TRUE.
326                time_dopr_av = MOD( time_dopr_av,                              &
327                               MAX( dt_averaging_input_pr, dt_spinup ) )
328             ENDIF
329          ENDIF
330          IF ( do_sum )  CALL flow_statistics
331
332!
333!--       Output of profiles
334          IF ( time_dopr >= dt_dopr )  THEN
335             IF ( dopr_n /= 0 )  CALL data_output_profiles
336             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
337             time_dopr_av = 0.0_wp    ! due to averaging (see above)
338          ENDIF
339
340!
341!--       Output of time series
342          IF ( time_dots >= dt_dots )  THEN
343             CALL data_output_tseries
344             time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
345          ENDIF
346
347       ENDIF
348
349!
350!--    Computation and output of run control parameters.
351!--    This is also done whenever perturbations have been imposed
352       IF ( time_run_control >= dt_run_control  .OR.                           &
353            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
354       THEN
355          CALL run_control
356          IF ( time_run_control >= dt_run_control )  THEN
357             time_run_control = MOD( time_run_control,                         &
358                                     MAX( dt_run_control, dt_spinup ) )
359          ENDIF
360       ENDIF
361
362       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
363
364       IF ( myid == 0 )  THEN
365          PRINT*, time_since_reference_point
366       ENDIF
367
368    ENDDO   ! time loop
369
370!
371!-- Write back saved temperature to the 3D arrays
372    pt(:,:,:)   = pt_save
373    pt_p(:,:,:) = pt_save
374
375    DEALLOCATE(pt_save)
376
377    CALL location_message( 'finished time-stepping spinup', .TRUE. )
378
379 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.