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

Last change on this file since 4869 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 29.4 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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration_spinup.f90 4828 2021-01-05 11:21:41Z schwenkel $
27! Bugfix in time-control of indoor model at the very first time step
28!
29! 4750 2020-10-16 14:27:48Z suehring
30! - bugfix, call hourly-based indoor model only once per hour, not every timestep.
31! - optionally switch-off/on the indoor model during spinup
32!
33! 4687 2020-09-21 19:40:16Z maronga
34! Indoor model is now available during spinup
35!
36! 4671 2020-09-09 20:27:58Z pavelkrc
37! Implementation of downward facing USM and LSM surfaces
38!
39! 4669 2020-09-09 13:43:47Z pavelkrc
40! - Fix missing call of radiation after spinup
41! - Fix calculation of force_radiation_call
42! - Fix calculation of radiation times
43!
44! 4668 2020-09-09 13:00:16Z pavelkrc
45! Improve debug messages during timestepping
46!
47! 4654 2020-08-28 13:47:23Z pavelkrc
48! Add possibility to output surface data during spinup (author: J. Resler)
49!
50! 4635 2020-08-05 16:44:15Z suehring
51! Bugfix in calling radiation_control, was called every dt_spinup rather than only dt_radiation
52!
53! 4631 2020-08-03 13:48:14Z suehring
54! Set pt1 attribute only element-wise
55!
56! 4540 2020-05-18 15:23:29Z raasch
57! File re-formatted to follow the PALM coding standard
58!
59! 4457 2020-03-11 14:20:43Z raasch
60! Use statement for exchange horiz added
61!
62! 4444 2020-03-05 15:59:50Z raasch
63! Bugfix: cpp-directives for serial mode added
64!
65! 4360 2020-01-07 11:25:50Z suehring
66! Enable output of diagnostic quantities, e.g. 2-m temperature
67!
68! 4227 2019-09-10 18:04:34Z gronemeier
69! Implement new palm_date_time_mod
70!
71! 4223 2019-09-10 09:20:47Z gronemeier
72! Corrected "Former revisions" section
73!
74! 4064 2019-07-01 05:33:33Z gronemeier
75! Moved call to radiation module out of intermediate time loop
76!
77! 4023 2019-06-12 13:20:01Z maronga
78! Time stamps are now negative in run control output
79!
80! 3885 2019-04-11 11:29:34Z kanani
81! Changes related to global restructuring of location messages and introduction of additional debug
82! messages
83!
84! 3766 2019-02-26 16:23:41Z raasch
85! Unused variable removed
86!
87! 3719 2019-02-06 13:10:18Z kanani
88! Removed log_point(19,54,74,50,75), since they count together with same log points in
89! time_integration, impossible to separate the contributions. Instead, the entire spinup gets an
90! individual log_point in palm.f90
91!
92! 3655 2019-01-07 16:51:22Z knoop
93! Removed call to calculation of near air (10 cm) potential temperature (now in surface layer fluxes)
94!
95! 2296 2017-06-28 07:53:56Z maronga
96! Initial revision
97!
98!
99! Description:
100! ------------
101!> Integration in time of the non-atmospheric model components such as land surface model and urban
102!> surface model
103!--------------------------------------------------------------------------------------------------!
104 SUBROUTINE time_integration_spinup
105
106    USE arrays_3d,                                                                                 &
107        ONLY:  pt,                                                                                 &
108               pt_p,                                                                               &
109               u,                                                                                  &
110               u_init,                                                                             &
111               v,                                                                                  &
112               v_init
113
114    USE control_parameters,                                                                        &
115        ONLY:  averaging_interval_pr,                                                              &
116               constant_diffusion,                                                                 &
117               constant_flux_layer,                                                                &
118               coupling_start_time,                                                                &
119               data_output_during_spinup,                                                          &
120               debug_output_timestep,                                                              &
121               debug_string,                                                                       &
122               dopr_n,                                                                             &
123               do_sum,                                                                             &
124               dt_averaging_input_pr,                                                              &
125               dt_dopr,                                                                            &
126               dt_dots,                                                                            &
127               dt_do2d_xy,                                                                         &
128               dt_do3d,                                                                            &
129               dt_spinup,                                                                          &
130               dt_3d,                                                                              &
131               humidity,                                                                           &
132               indoor_model,                                                                       &
133               intermediate_timestep_count,                                                        &
134               intermediate_timestep_count_max,                                                    &
135               land_surface,                                                                       &
136               simulated_time,                                                                     &
137               simulated_time_chr,                                                                 &
138               skip_time_dopr,                                                                     &
139               skip_time_do2d_xy,                                                                  &
140               skip_time_do3d,                                                                     &
141               spinup_pt_amplitude,                                                                &
142               spinup_pt_mean,                                                                     &
143               spinup_time,                                                                        &
144               surface_output,                                                                     &
145               timestep_count,                                                                     &
146               time_dopr,                                                                          &
147               time_dopr_av,                                                                       &
148               time_dots,                                                                          &
149               time_do2d_xy,                                                                       &
150               time_do3d,                                                                          &
151               time_run_control,                                                                   &
152               time_since_reference_point,                                                         &
153               urban_surface
154
155    USE cpulog,                                                                                    &
156        ONLY:  cpu_log,                                                                            &
157               log_point_s
158
159    USE diagnostic_output_quantities_mod,                                                          &
160        ONLY:  doq_calculate
161
162    USE exchange_horiz_mod,                                                                        &
163        ONLY:  exchange_horiz
164
165    USE indices,                                                                                   &
166        ONLY:  nbgp,                                                                               &
167               nzb,                                                                                &
168               nzt,                                                                                &
169               nysg,                                                                               &
170               nyng,                                                                               &
171               nxlg,                                                                               &
172               nxrg
173
174    USE indoor_model_mod,                                                                          &
175        ONLY:  dt_indoor,                                                                          &
176               im_main_heatcool,                                                                   &
177               indoor_during_spinup,                                                               &
178               time_indoor
179
180    USE land_surface_model_mod,                                                                    &
181        ONLY:  lsm_energy_balance,                                                                 &
182               lsm_swap_timelevel
183
184    USE pegrid
185
186#if defined( __parallel )
187    USE pmc_interface,                                                                             &
188        ONLY:  nested_run
189#endif
190
191    USE kinds
192
193    USE palm_date_time_mod,                                                                        &
194        ONLY:  get_date_time,                                                                      &
195               seconds_per_hour
196
197    USE radiation_model_mod,                                                                       &
198        ONLY:  dt_radiation,                                                                       &
199               force_radiation_call,                                                               &
200               radiation,                                                                          &
201               radiation_control,                                                                  &
202               radiation_interaction,                                                              &
203               radiation_interactions,                                                             &
204               time_radiation
205
206    USE statistics,                                                                                &
207        ONLY:  flow_statistics_called
208
209    USE surface_data_output_mod,                                                                   &
210        ONLY:  dt_dosurf,                                                                          &
211               skip_time_dosurf,                                                                   &
212               surface_data_output,                                                                &
213               time_dosurf
214
215    USE surface_layer_fluxes_mod,                                                                  &
216        ONLY:  surface_layer_fluxes
217
218    USE surface_mod,                                                                               &
219        ONLY :  surf_lsm_h,                                                                        &
220                surf_lsm_v, surf_usm_h,                                                            &
221                surf_usm_v
222
223    USE urban_surface_mod,                                                                         &
224        ONLY:  usm_energy_balance,                                                                 &
225               usm_swap_timelevel
226
227
228
229
230    IMPLICIT NONE
231
232    CHARACTER(LEN=1) ::  sign_chr                        !< String containing '-' or ' '
233    CHARACTER(LEN=9) ::  time_since_reference_point_chr  !< time since reference point, i.e., negative during spinup
234    CHARACTER(LEN=9) ::  time_to_string                  !<
235
236
237    INTEGER(iwp) ::  current_timestep_number_spinup = 0  !< number if timestep during spinup
238    INTEGER(iwp) ::  day_of_year                         !< day of the year
239
240    INTEGER(iwp) ::  i  !< running index
241    INTEGER(iwp) ::  j  !< running index
242    INTEGER(iwp) ::  k  !< running index
243    INTEGER(iwp) ::  l  !< running index
244    INTEGER(iwp) ::  m  !< running index
245
246
247    LOGICAL ::  run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
248
249
250    REAL(wp) ::  dt_save        !< temporary storage for time step
251    REAL(wp) ::  pt_spinup      !< temporary storage of temperature
252    REAL(wp) ::  second_of_day  !< second of the day
253
254    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save  !< temporary storage of temperature
255    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_save   !< temporary storage of u wind component
256    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_save   !< temporary storage of v wind component
257
258
259!
260!-- Save 3D arrays because they are to be changed for spinup purpose
261    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
262    ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
263    ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
264
265    CALL exchange_horiz( pt, nbgp )
266    CALL exchange_horiz( u,  nbgp )
267    CALL exchange_horiz( v,  nbgp )
268
269    pt_save = pt
270    u_save  = u
271    v_save  = v
272
273!
274!-- Set the same wall-adjacent velocity to all grid points. The sign of the original velocity field
275!-- must be preserved because the surface schemes crash otherwise. The precise reason is still
276!-- unknown. A minimum velocity of 0.1 m/s is used to maintain turbulent transfer at the surface.
277    IF ( land_surface )  THEN
278       DO  l = 0, 1
279          DO  m = 1, surf_lsm_h(l)%ns
280             i   = surf_lsm_h(l)%i(m)
281             j   = surf_lsm_h(l)%j(m)
282             k   = surf_lsm_h(l)%k(m)
283             u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp)
284             v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp)
285          ENDDO
286       ENDDO
287
288       DO  l = 0, 3
289          DO  m = 1, surf_lsm_v(l)%ns
290             i   = surf_lsm_v(l)%i(m)
291             j   = surf_lsm_v(l)%j(m)
292             k   = surf_lsm_v(l)%k(m)
293             u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp)
294             v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp)
295          ENDDO
296       ENDDO
297    ENDIF
298
299    IF ( urban_surface )  THEN
300       DO  l = 0, 1
301          DO  m = 1, surf_usm_h(l)%ns
302             i   = surf_usm_h(l)%i(m)
303             j   = surf_usm_h(l)%j(m)
304             k   = surf_usm_h(l)%k(m)
305             u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp)
306             v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp)
307          ENDDO
308       ENDDO
309
310       DO  l = 0, 3
311          DO  m = 1, surf_usm_v(l)%ns
312             i   = surf_usm_v(l)%i(m)
313             j   = surf_usm_v(l)%j(m)
314             k   = surf_usm_v(l)%k(m)
315             u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp)
316             v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp)
317          ENDDO
318       ENDDO
319    ENDIF
320
321    CALL exchange_horiz( u, nbgp )
322    CALL exchange_horiz( v, nbgp )
323
324    dt_save = dt_3d
325    dt_3d   = dt_spinup
326
327    CALL location_message( 'wall/soil spinup time-stepping', 'start' )
328!
329!-- Start of the time loop
330    DO  WHILE ( simulated_time < spinup_time )
331
332       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
333
334       IF ( debug_output_timestep )  THEN
335           WRITE( debug_string, * ) 'time_integration_spinup', simulated_time
336           CALL debug_message( debug_string, 'start' )
337       ENDIF
338
339!
340!--    Start of intermediate step loop
341       intermediate_timestep_count = 0
342       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
343
344          intermediate_timestep_count = intermediate_timestep_count + 1
345
346!
347!--       Set the steering factors for the prognostic equations which depend on the timestep scheme
348          CALL timestep_scheme_steering
349
350
351!
352!--       Estimate a near-surface air temperature based on the position of the sun and user input
353!--       about mean temperature and amplitude. The time is shifted by one hour to simulate a lag
354!--       between air temperature and incoming radiation.
355          CALL get_date_time( simulated_time - spinup_time - seconds_per_hour,                     &
356                              day_of_year = day_of_year, second_of_day = second_of_day )
357
358          pt_spinup = spinup_pt_mean + spinup_pt_amplitude *                                       &
359                      solar_angle( day_of_year, second_of_day )
360
361!
362!--       Map air temperature to all grid points in the vicinity of a surface element
363          IF ( land_surface )  THEN
364             DO  l = 0, 1
365                DO  m = 1, surf_lsm_h(l)%ns
366                   i   = surf_lsm_h(l)%i(m)
367                   j   = surf_lsm_h(l)%j(m)
368                   k   = surf_lsm_h(l)%k(m)
369                   pt(k,j,i) = pt_spinup
370                ENDDO
371             ENDDO
372
373             DO  l = 0, 3
374                DO  m = 1, surf_lsm_v(l)%ns
375                   i   = surf_lsm_v(l)%i(m)
376                   j   = surf_lsm_v(l)%j(m)
377                   k   = surf_lsm_v(l)%k(m)
378                   pt(k,j,i) = pt_spinup
379                ENDDO
380             ENDDO
381          ENDIF
382
383          IF ( urban_surface )  THEN
384             DO  l = 0, 1
385                DO  m = 1, surf_usm_h(l)%ns
386                   i   = surf_usm_h(l)%i(m)
387                   j   = surf_usm_h(l)%j(m)
388                   k   = surf_usm_h(l)%k(m)
389                   pt(k,j,i) = pt_spinup
390                   !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
391                   surf_usm_h(l)%pt1(m) = pt_spinup
392                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393                ENDDO
394             ENDDO
395
396             DO  l = 0, 3
397                DO  m = 1, surf_usm_v(l)%ns
398                   i   = surf_usm_v(l)%i(m)
399                   j   = surf_usm_v(l)%j(m)
400                   k   = surf_usm_v(l)%k(m)
401                   pt(k,j,i) = pt_spinup
402                   !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
403                   surf_usm_v(l)%pt1(m) = pt_spinup
404                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
405                ENDDO
406             ENDDO
407          ENDIF
408
409          CALL exchange_horiz( pt, nbgp )
410
411
412!
413!--       Swap the time levels in preparation for the next time step.
414          timestep_count = timestep_count + 1
415
416          IF ( land_surface )  THEN
417              CALL lsm_swap_timelevel ( 0 )
418          ENDIF
419
420          IF ( urban_surface )  THEN
421             CALL usm_swap_timelevel ( 0 )
422          ENDIF
423
424          IF ( land_surface )  THEN
425             CALL lsm_swap_timelevel ( MOD( timestep_count, 2 ) )
426          ENDIF
427
428          IF ( urban_surface )  THEN
429             CALL usm_swap_timelevel ( MOD( timestep_count, 2 ) )
430          ENDIF
431
432!
433!--       If required, compute virtual potential temperature
434          IF ( humidity )  THEN
435             CALL compute_vpt
436          ENDIF
437
438!
439!--       Compute the diffusion quantities
440          IF ( .NOT. constant_diffusion )  THEN
441
442!
443!--          First the vertical (and horizontal) fluxes in the surface (constant flux) layer are
444!--          computed
445             IF ( constant_flux_layer )  THEN
446                CALL surface_layer_fluxes
447             ENDIF
448
449!
450!--          If required, solve the energy balance for the surface and run soil model. Call for
451!--          horizontal as well as vertical surfaces. The prognostic equation for soil moisure is
452!--          switched off
453             IF ( land_surface )  THEN
454
455                CALL lsm_energy_balance( .TRUE. )
456
457             ENDIF
458
459!
460!--          If required, solve the energy balance for urban surfaces and run the material heat model
461             IF (urban_surface) THEN
462
463                CALL usm_energy_balance( .TRUE. )
464
465             ENDIF
466
467          ENDIF
468
469       ENDDO   ! Intermediate step loop
470
471!
472!--    If required, calculate radiative fluxes and heating rates
473       IF ( radiation )  THEN
474
475          time_radiation = time_radiation + dt_3d
476
477          IF ( time_radiation >= dt_radiation .OR. force_radiation_call )  THEN
478
479             IF ( .NOT. force_radiation_call )  THEN
480                time_radiation = time_radiation - dt_radiation
481             ENDIF
482
483             CALL radiation_control
484
485             IF ( radiation_interactions )  THEN
486                CALL radiation_interaction
487             ENDIF
488!
489!--          Reset forcing of radiation call
490             force_radiation_call = .FALSE.
491
492          ENDIF
493       ENDIF
494
495!
496!--    If required, calculate indoor temperature, waste heat, heat flux
497!--    through wall, etc.
498!--    dt_indoor steers the frequency of the indoor model calculations.
499!--    Note, at first timestep indoor model is called, in order to provide
500!--    a waste heat flux.
501       IF ( indoor_model  .AND.  indoor_during_spinup )  THEN
502
503          time_indoor = time_indoor + dt_3d
504
505          IF ( time_indoor >= dt_indoor  .OR.  current_timestep_number_spinup == 0 )  THEN
506
507             IF ( time_indoor >= dt_indoor )  time_indoor = time_indoor - dt_indoor
508
509             CALL im_main_heatcool
510
511          ENDIF
512       ENDIF
513
514!
515!--    Increase simulation time and output times
516       current_timestep_number_spinup = current_timestep_number_spinup + 1
517       simulated_time                 = simulated_time   + dt_3d
518       simulated_time_chr             = time_to_string( simulated_time )
519       time_since_reference_point     = simulated_time - coupling_start_time
520       time_since_reference_point_chr = time_to_string( ABS( time_since_reference_point ) )
521
522       IF ( time_since_reference_point < 0.0_wp )  THEN
523          sign_chr = '-'
524       ELSE
525          sign_chr = ' '
526       ENDIF
527
528
529       IF ( data_output_during_spinup )  THEN
530          IF ( simulated_time >= skip_time_do2d_xy )  THEN
531             time_do2d_xy      = time_do2d_xy     + dt_3d
532          ENDIF
533          IF ( simulated_time >= skip_time_do3d    )  THEN
534             time_do3d         = time_do3d        + dt_3d
535          ENDIF
536          time_dots            = time_dots        + dt_3d
537          IF ( simulated_time >= skip_time_dopr )  THEN
538             time_dopr         = time_dopr        + dt_3d
539          ENDIF
540          IF ( surface_output )  THEN
541             IF ( simulated_time >= skip_time_dosurf )  THEN
542                time_dosurf    = time_dosurf + dt_3d
543             ENDIF
544          ENDIF
545          time_run_control     = time_run_control + dt_3d
546
547!
548!--       Carry out statistical analysis and output at the requested output times.
549!--       The MOD function is used for calculating the output time counters (like time_dopr) in
550!--       order to regard a possible decrease of the output time interval in case of restart runs.
551
552!
553!--       Set a flag indicating that so far no statistics have been created for this time step
554          flow_statistics_called = .FALSE.
555
556!
557!--       If required, call flow_statistics for averaging in time
558          IF ( averaging_interval_pr /= 0.0_wp  .AND.                                              &
559               ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.                             &
560               simulated_time >= skip_time_dopr )                                                  &
561          THEN
562             time_dopr_av = time_dopr_av + dt_3d
563             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
564                do_sum = .TRUE.
565                time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
566             ENDIF
567          ENDIF
568          IF ( do_sum )  CALL flow_statistics
569
570!
571!--       Output of profiles
572          IF ( time_dopr >= dt_dopr )  THEN
573             IF ( dopr_n /= 0 )  CALL data_output_profiles
574             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
575             time_dopr_av = 0.0_wp    ! Due to averaging (see above)
576          ENDIF
577
578!
579!--       Output of time series
580          IF ( time_dots >= dt_dots )  THEN
581             CALL data_output_tseries
582             time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
583          ENDIF
584
585!
586!--       2d-data output (cross-sections)
587          IF ( time_do2d_xy >= dt_do2d_xy )  THEN
588             CALL doq_calculate
589             CALL data_output_2d( 'xy', 0 )
590             time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
591          ENDIF
592
593!
594!--       3d-data output (volume data)
595          IF ( time_do3d >= dt_do3d )  THEN
596             CALL doq_calculate
597             CALL data_output_3d( 0 )
598             time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
599          ENDIF
600!
601!--       Output of surface data
602          IF ( surface_output )  THEN
603             IF ( time_dosurf >= dt_dosurf )  THEN
604                CALL surface_data_output( 0 )
605                time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
606             ENDIF
607          ENDIF
608
609       ENDIF
610
611!
612!--    Computation and output of run control parameters. This is also done whenever perturbations
613!--    have been imposed
614!        IF ( time_run_control >= dt_run_control  .OR.                                              &
615!             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )  THEN
616!           CALL run_control
617!           IF ( time_run_control >= dt_run_control )  THEN
618!              time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d ) )
619!           ENDIF
620!        ENDIF
621
622       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
623
624
625!
626!--    Run control output
627       IF ( myid == 0 )  THEN
628!
629!--       If necessary, write header
630          IF ( .NOT. run_control_header_spinup )  THEN
631             CALL check_open( 15 )
632             WRITE ( 15, 100 )
633             run_control_header_spinup = .TRUE.
634          ENDIF
635!
636!--       Write some general information about the spinup in run control file
637          WRITE ( 15, 101 )  current_timestep_number_spinup, sign_chr,                             &
638                             time_since_reference_point_chr, dt_3d, pt_spinup
639!
640!--       Write buffer contents to disc immediately
641          FLUSH( 15 )
642       ENDIF
643
644       IF ( debug_output_timestep )  THEN
645           WRITE( debug_string, * ) 'time_integration_spinup', simulated_time
646           CALL debug_message( debug_string, 'end' )
647       ENDIF
648
649
650    ENDDO   ! Time loop
651
652!
653!-- Write back saved arrays to the 3D arrays
654    pt   = pt_save
655    pt_p = pt_save
656    u    = u_save
657    v    = v_save
658
659!
660!-- Reset time step
661    dt_3d = dt_save
662!-- Force radiation step in time zero
663!-- It is performed at the end of init_radiation in case of run without spinup
664    time_radiation = dt_radiation
665
666    DEALLOCATE(pt_save)
667    DEALLOCATE(u_save)
668    DEALLOCATE(v_save)
669
670#if defined( __parallel )
671    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
672#endif
673
674    CALL location_message( 'wall/soil spinup time-stepping', 'finished' )
675
676
677!
678!-- Formats
679100 FORMAT (///'Spinup control output:---------------------------------'//                         &
680            'ITER.   HH:MM:SS    DT   PT(z_MO)---------------------------------')
681101 FORMAT (I5,2X,A1,A9,1X,F6.2,3X,F6.2,2X,F6.2)
682
683 CONTAINS
684
685!
686!-- Returns the cosine of the solar zenith angle at a given time. This routine is similar to that
687!-- for calculation zenith (see radiation_model_mod.f90)
688    !> @todo Load function calc_zenith of radiation model instead of rewrite the function here.
689    FUNCTION solar_angle( day_of_year, second_of_day )
690
691       USE basic_constants_and_equations_mod,                                                      &
692           ONLY:  pi
693
694       USE kinds
695
696       USE radiation_model_mod,                                                                    &
697           ONLY:  decl_1,                                                                          &
698                  decl_2,                                                                          &
699                  decl_3,                                                                          &
700                  lat,                                                                             &
701                  lon
702
703       IMPLICIT NONE
704
705
706       INTEGER(iwp), INTENT(IN) ::  day_of_year  !< day of the year
707
708       REAL(wp)             ::  declination    !< solar declination angle
709       REAL(wp)             ::  hour_angle     !< solar hour angle
710       REAL(wp), INTENT(IN) ::  second_of_day  !< current time of the day in UTC
711       REAL(wp)             ::  solar_angle    !< cosine of the solar zenith angle
712!
713!--    Calculate solar declination and hour angle
714       declination = ASIN( decl_1 * SIN( decl_2 * REAL( day_of_year, KIND = wp) - decl_3 ) )
715       hour_angle  = 2.0_wp * pi * ( second_of_day / 86400.0_wp ) + lon - pi
716
717!
718!--    Calculate cosine of solar zenith angle
719       solar_angle = SIN( lat ) * SIN( declination ) + COS( lat ) * COS( declination ) *           &
720                     COS( hour_angle )
721
722    END FUNCTION solar_angle
723
724
725 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.