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

Last change on this file since 4230 was 4227, checked in by gronemeier, 5 years ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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