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

Last change on this file since 2881 was 2881, checked in by maronga, 4 years ago

bugfix in land surface model and new option in spinup

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