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

Last change on this file since 2956 was 2934, checked in by suehring, 7 years ago

Synchronize parent and child model after initialization and spinup phase; Check for consistent setting of spinup times in parent and child model; remove obsolete masking of tendency arrays during initialization

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