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

Last change on this file since 4651 was 4635, checked in by suehring, 4 years ago

Bugfix in calling radiation_control, was called every dt_spinup rather than only dt_radiation

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