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

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

bugfix in indoor model, code layout change in urban surface model, indoor model integrated in spinup mechanism

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