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

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

improvements in spinup mechanism

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