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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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