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

Last change on this file since 4470 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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