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

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

bugfix/adjusment in spinup mechanism

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