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

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