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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

  • Property svn:keywords set to Id
File size: 17.1 KB
Line 
1!> @file time_integration_spinup.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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 2716 2017-12-29 16:35:59Z kanani $
27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
31! Added radiation interactions (moved from USM) (MS)
32!
33! 2544 2017-10-13 18:09:32Z maronga
34! Date and time quantities are now read from date_and_time_mod
35!
36! 2299 2017-06-29 10:14:38Z maronga
37! Call of soil model adjusted to avoid prognostic equation for soil moisture
38! during spinup.
39! Better representation of diurnal cycle of near-surface temperature.
40! Excluded prognostic equation for soil moisture during spinup.
41! Added output of run control data for spinup.
42!
43! 2297 2017-06-28 14:35:57Z scharf
44! bugfixes
45!
46! 2296 2017-06-28 07:53:56Z maronga
47! Initial revision
48!
49!
50! Description:
51! ------------
52!> Integration in time of the non-atmospheric model components such as land
53!> surface model and urban surface model
54!------------------------------------------------------------------------------!
55 SUBROUTINE time_integration_spinup
56 
57    USE arrays_3d,                                                             &
58        ONLY:  pt, pt_p
59
60    USE control_parameters,                                                    &
61        ONLY:  averaging_interval_pr, constant_diffusion, constant_flux_layer, &
62               coupling_start_time, current_timestep_number,                   &
63               data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
64               dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control,        &
65               dt_spinup, humidity, intermediate_timestep_count,               &
66               intermediate_timestep_count_max, land_surface,                  &
67               simulated_time, simulated_time_chr,      &
68               skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
69               spinup_time, timestep_count, timestep_scheme, time_dopr,        &
70               time_dopr_av, time_dots, time_run_control,                      &
71               time_since_reference_point, urban_surface
72
73    USE constants,                                                             &
74        ONLY:  pi
75
76    USE cpulog,                                                                &
77        ONLY:  cpu_log, log_point, log_point_s
78
79    USE date_and_time_mod,                                                     &
80        ONLY: day_of_year_init, time_utc_init
81
82    USE indices,                                                               &
83        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
84
85
86    USE land_surface_model_mod,                                                &
87        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
88
89    USE pegrid,                                                                &
90        ONLY:  myid
91
92    USE kinds
93
94    USE radiation_model_mod,                                                   &
95        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
96               radiation_control, rad_sw_in, time_radiation,                   &
97               radiation_interaction, radiation_interactions
98
99    USE statistics,                                                            &
100        ONLY:  flow_statistics_called
101
102    USE surface_layer_fluxes_mod,                                              &
103        ONLY:  surface_layer_fluxes
104
105    USE surface_mod,                                                           &
106        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
107                surf_usm_v
108
109    USE urban_surface_mod,                                                     &
110        ONLY:  usm_material_heat_model, usm_material_model,                    &
111               usm_surface_energy_balance, usm_swap_timelevel,                 &
112               usm_green_heat_model, usm_temperature_near_surface
113
114
115
116
117    IMPLICIT NONE
118
119    CHARACTER (LEN=9) ::  time_to_string          !<
120 
121    INTEGER(iwp) ::  i !< running index
122    INTEGER(iwp) ::  j !< running index
123    INTEGER(iwp) ::  k !< running index
124    INTEGER(iwp) ::  l !< running index
125    INTEGER(iwp) ::  m !< running index
126
127    INTEGER(iwp) :: current_timestep_number_spinup = 0  !< number if timestep during spinup
128 
129    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
130
131    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
132                 
133    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
134
135    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
136
137    CALL exchange_horiz( pt, nbgp )   
138    pt_save = pt
139
140    CALL location_message( 'starting spinup-sequence', .TRUE. )
141!
142!-- Start of the time loop
143    DO  WHILE ( simulated_time < spinup_time )
144
145       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
146   
147!
148!--    Start of intermediate step loop
149       intermediate_timestep_count = 0
150       DO  WHILE ( intermediate_timestep_count < &
151                   intermediate_timestep_count_max )
152
153          intermediate_timestep_count = intermediate_timestep_count + 1
154
155!
156!--       Set the steering factors for the prognostic equations which depend
157!--       on the timestep scheme
158          CALL timestep_scheme_steering
159
160
161!
162!--       Estimate a near-surface air temperature based on the position of the
163!--       sun and user input about mean temperature and amplitude. The time is
164!--       shifted by one hour to simulate a lag between air temperature and
165!--       incoming radiation
166          pt_spinup = spinup_pt_mean + spinup_pt_amplitude                     &
167             * solar_angle (time_utc_init + time_since_reference_point - 3600.0)
168
169!
170!--       Map air temperature to all grid points in the vicinity of a surface
171!--       element
172          IF ( land_surface )  THEN
173             DO  m = 1, surf_lsm_h%ns
174                i   = surf_lsm_h%i(m)           
175                j   = surf_lsm_h%j(m)
176                k   = surf_lsm_h%k(m)
177                pt(k,j,i) = pt_spinup
178             ENDDO
179
180             DO  l = 0, 3
181                DO  m = 1, surf_lsm_v(l)%ns
182                   i   = surf_lsm_v(l)%i(m)           
183                   j   = surf_lsm_v(l)%j(m)
184                   k   = surf_lsm_v(l)%k(m)
185                   pt(k,j,i) = pt_spinup
186                ENDDO
187             ENDDO
188          ENDIF
189
190          IF ( urban_surface )  THEN
191             DO  m = 1, surf_usm_h%ns
192                i   = surf_usm_h%i(m)           
193                j   = surf_usm_h%j(m)
194                k   = surf_usm_h%k(m)
195                pt(k,j,i) = pt_spinup
196             ENDDO
197
198             DO  l = 0, 3
199                DO  m = 1, surf_usm_v(l)%ns
200                   i   = surf_usm_v(l)%i(m)           
201                   j   = surf_usm_v(l)%j(m)
202                   k   = surf_usm_v(l)%k(m)
203                   pt(k,j,i) = pt_spinup
204                ENDDO
205             ENDDO
206          ENDIF
207
208!
209!--       Swap the time levels in preparation for the next time step.
210          timestep_count = timestep_count + 1
211     
212          IF ( land_surface )  THEN
213              CALL lsm_swap_timelevel ( 0 )
214          ENDIF
215
216          IF ( urban_surface )  THEN
217             CALL usm_swap_timelevel ( 0 )
218          ENDIF
219
220          IF ( land_surface )  THEN
221             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
222          ENDIF
223
224          IF ( urban_surface )  THEN
225             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
226          ENDIF
227         
228!
229!--       If required, compute virtual potential temperature
230          IF ( humidity )  THEN
231             CALL compute_vpt 
232          ENDIF 
233
234!
235!--       Compute the diffusion quantities
236          IF ( .NOT. constant_diffusion )  THEN
237
238!
239!--          First the vertical (and horizontal) fluxes in the surface
240!--          (constant flux) layer are computed
241             IF ( constant_flux_layer )  THEN
242                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
243                CALL surface_layer_fluxes
244                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
245             ENDIF
246
247!
248!--          If required, solve the energy balance for the surface and run soil
249!--          model. Call for horizontal as well as vertical surfaces.
250!--          The prognostic equation for soil moisure is switched off
251             IF ( land_surface )  THEN
252
253                CALL cpu_log( log_point(54), 'land_surface', 'start' )
254!
255!--             Call for horizontal upward-facing surfaces
256                CALL lsm_energy_balance( .TRUE., -1 )
257                CALL lsm_soil_model( .TRUE., -1, .FALSE. )
258!
259!--             Call for northward-facing surfaces
260                CALL lsm_energy_balance( .FALSE., 0 )
261                CALL lsm_soil_model( .FALSE., 0, .FALSE. )
262!
263!--             Call for southward-facing surfaces
264                CALL lsm_energy_balance( .FALSE., 1 )
265                CALL lsm_soil_model( .FALSE., 1, .FALSE. )
266!
267!--             Call for eastward-facing surfaces
268                CALL lsm_energy_balance( .FALSE., 2 )
269                CALL lsm_soil_model( .FALSE., 2, .FALSE. )
270!
271!--             Call for westward-facing surfaces
272                CALL lsm_energy_balance( .FALSE., 3 )
273                CALL lsm_soil_model( .FALSE., 3, .FALSE. )
274
275                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
276             ENDIF
277
278!
279!--          If required, solve the energy balance for urban surfaces and run
280!--          the material heat model
281             IF (urban_surface) THEN
282                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
283                CALL usm_surface_energy_balance
284                IF ( usm_material_model )  THEN
285                   CALL usm_green_heat_model
286                   CALL usm_material_heat_model
287                ENDIF
288                IF ( urban_surface ) THEN
289                   CALL usm_temperature_near_surface
290                ENDIF
291                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
292             ENDIF
293
294          ENDIF
295
296!
297!--       If required, calculate radiative fluxes and heating rates
298          IF ( radiation .AND. intermediate_timestep_count                     &
299               == intermediate_timestep_count_max )  THEN
300
301               time_radiation = time_radiation + dt_spinup
302
303             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
304             THEN
305
306                CALL cpu_log( log_point(50), 'radiation', 'start' )
307
308                IF ( .NOT. force_radiation_call )  THEN
309                   time_radiation = time_radiation - dt_radiation
310                ENDIF
311
312                CALL radiation_control
313
314                CALL cpu_log( log_point(50), 'radiation', 'stop' )
315
316                IF ( radiation_interactions )  THEN
317                   CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
318                   CALL radiation_interaction
319                   CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
320                ENDIF
321             ENDIF
322          ENDIF
323
324       ENDDO   ! Intermediate step loop
325
326!
327!--    Increase simulation time and output times
328       current_timestep_number_spinup = current_timestep_number_spinup + 1
329       simulated_time             = simulated_time   + dt_spinup
330       simulated_time_chr         = time_to_string( simulated_time )
331       time_since_reference_point = simulated_time - coupling_start_time
332
333       IF ( data_output_during_spinup )  THEN
334          time_dots          = time_dots        + dt_spinup
335          IF ( simulated_time >= skip_time_dopr )  THEN
336             time_dopr       = time_dopr        + dt_spinup
337          ENDIF
338          time_run_control   = time_run_control + dt_spinup
339
340!
341!--       Carry out statistical analysis and output at the requested output times.
342!--       The MOD function is used for calculating the output time counters (like
343!--       time_dopr) in order to regard a possible decrease of the output time
344!--       interval in case of restart runs
345
346!
347!--       Set a flag indicating that so far no statistics have been created
348!--       for this time step
349          flow_statistics_called = .FALSE.
350
351!
352!--       If required, call flow_statistics for averaging in time
353          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
354             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
355             simulated_time >= skip_time_dopr )  THEN
356             time_dopr_av = time_dopr_av + dt_spinup
357             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
358                do_sum = .TRUE.
359                time_dopr_av = MOD( time_dopr_av,                              &
360                               MAX( dt_averaging_input_pr, dt_spinup ) )
361             ENDIF
362          ENDIF
363          IF ( do_sum )  CALL flow_statistics
364
365!
366!--       Output of profiles
367          IF ( time_dopr >= dt_dopr )  THEN
368             IF ( dopr_n /= 0 )  CALL data_output_profiles
369             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
370             time_dopr_av = 0.0_wp    ! due to averaging (see above)
371          ENDIF
372
373!
374!--       Output of time series
375          IF ( time_dots >= dt_dots )  THEN
376             CALL data_output_tseries
377             time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
378          ENDIF
379
380       ENDIF
381
382!
383!--    Computation and output of run control parameters.
384!--    This is also done whenever perturbations have been imposed
385!        IF ( time_run_control >= dt_run_control  .OR.                           &
386!             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
387!        THEN
388!           CALL run_control
389!           IF ( time_run_control >= dt_run_control )  THEN
390!              time_run_control = MOD( time_run_control,                         &
391!                                      MAX( dt_run_control, dt_spinup ) )
392!           ENDIF
393!        ENDIF
394
395       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
396
397
398!
399!--    Run control output
400       IF ( myid == 0 )  THEN
401!
402!--       If necessary, write header
403          IF ( .NOT. run_control_header_spinup )  THEN
404             CALL check_open( 15 )
405             WRITE ( 15, 100 )
406             run_control_header_spinup = .TRUE.
407          ENDIF
408!
409!--       Write some general information about the spinup in run control file
410          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_spinup, pt_spinup, rad_sw_in(0,nysg,nxlg)
411!
412!--       Write buffer contents to disc immediately
413          FLUSH( 15 )
414       ENDIF
415
416
417
418    ENDDO   ! time loop
419
420!
421!-- Write back saved temperature to the 3D arrays
422    pt(:,:,:)   = pt_save
423    pt_p(:,:,:) = pt_save
424
425    DEALLOCATE(pt_save)
426
427    CALL location_message( 'finished spinup-sequence', .TRUE. )
428
429
430!
431!-- Formats
432100 FORMAT (///'Spinup control output:'/  &
433            '----------------------------------------'// &
434            'ITER.  HH:MM:SS    DT   PT(z_MO)     SWD'/   &
435            '----------------------------------------')
436101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
437
438 CONTAINS
439
440!
441!-- Returns the cosine of the solar zenith angle at a given time. This routine
442!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
443    FUNCTION solar_angle( local_time ) 
444
445       USE constants,                                                          &
446       ONLY:  pi
447     
448       USE kinds
449
450       USE radiation_model_mod,                                                &
451           ONLY:  decl_1, decl_2, decl_3, lat, lon
452
453       IMPLICIT NONE
454
455
456       REAL(wp) ::  solar_angle     !< cosine of the solar zenith angle
457
458       REAL(wp) ::  day             !< day of the year
459       REAL(wp) ::  declination     !< solar declination angle
460       REAL(wp) ::  hour_angle      !< solar hour angle
461       REAL(wp) ::  time_utc        !< current time in UTC
462       REAL(wp), INTENT(IN) :: local_time
463!
464!--    Calculate current day and time based on the initial values and simulation
465!--    time
466
467       day = day_of_year_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
468       time_utc = MOD(local_time, 86400.0_wp)
469
470
471!
472!--    Calculate solar declination and hour angle   
473       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
474       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
475
476!
477!--    Calculate cosine of solar zenith angle
478       solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
479                     * COS(hour_angle)
480
481
482    END FUNCTION solar_angle
483
484
485 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.