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

Last change on this file since 4598 was 4540, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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