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

Last change on this file since 2544 was 2544, checked in by maronga, 6 years ago

introduced new module date_and_time_mod

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