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

Last change on this file since 3697 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 22.2 KB
Line 
1!> @file time_integration_spinup.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration_spinup.f90 3655 2019-01-07 16:51:22Z hellstea $
27! Removed call to calculation of near air (10 cm) potential temperature (now in
28! surface layer fluxes)
29!
30! 3467 2018-10-30 19:05:21Z suehring
31! call to material_heat_model now with check if spinup runs (rvtils)
32!
33! 3337 2018-10-12 15:17:09Z kanani
34! (from branch resler)
35! Add pt1 initialization
36!
37! 3274 2018-09-24 15:42:55Z knoop
38! Modularization of all bulk cloud physics code components
39!
40! 3241 2018-09-12 15:02:00Z raasch
41! unused variables removed
42!
43! 2983 2018-04-18 10:43:40Z suehring
44! Revise limitation of wall-adjacent velocity.
45!
46! 2934 2018-03-26 19:13:22Z suehring
47! Synchronize parent and child models after spinup.
48!
49! 2881 2018-03-13 16:24:40Z maronga
50! Added flag for switching on/off calculation of soil moisture
51!
52! 2818 2018-02-19 16:42:36Z maronga
53! Velocity components near walls/ground are now set to the profiles stored in
54! u_init and v_init. Activated soil moisture calculation during spinup.
55!
56! 2782 2018-02-02 11:51:10Z maronga
57! Bugfix and re-activation of homogeneous setting of velocity components
58! during spinup
59!
60! 2758 2018-01-17 12:55:21Z suehring
61! Comment out homogeneous setting of wind velocity as this will lead to zero
62! friction velocity and cause problems in MOST relationships.
63!
64! 2728 2018-01-09 07:03:53Z maronga
65! Set velocity componenets to homogeneous values during spinup
66!
67! 2724 2018-01-05 12:12:38Z maronga
68! Use dt_spinup for all active components during spinup
69!
70! 2723 2018-01-05 09:27:03Z maronga
71! Bugfix: array rad_sw_in no longer exists and is thus removed from RUN_CONTROL
72! output.
73! Added output of XY and 3D data during spinup.
74! Bugfix: time step in LSM and USM was set to dt_3d instead of dt_spinup
75!
76! 2718 2018-01-02 08:49:38Z maronga
77! Corrected "Former revisions" section
78!
79! 2696 2017-12-14 17:12:51Z kanani
80! Change in file header (GPL part)
81! Added radiation interactions (moved from USM) (MS)
82!
83! 2544 2017-10-13 18:09:32Z maronga
84! Date and time quantities are now read from date_and_time_mod
85!
86! 2299 2017-06-29 10:14:38Z maronga
87! Call of soil model adjusted to avoid prognostic equation for soil moisture
88! during spinup.
89! Better representation of diurnal cycle of near-surface temperature.
90! Excluded prognostic equation for soil moisture during spinup.
91! Added output of run control data for spinup.
92!
93! 2297 2017-06-28 14:35:57Z scharf
94! bugfixes
95!
96! 2296 2017-06-28 07:53:56Z maronga
97! Initial revision
98!
99!
100! Description:
101! ------------
102!> Integration in time of the non-atmospheric model components such as land
103!> surface model and urban surface model
104!------------------------------------------------------------------------------!
105 SUBROUTINE time_integration_spinup
106 
107    USE arrays_3d,                                                             &
108        ONLY:  pt, pt_p, u, u_init, v, v_init
109
110    USE control_parameters,                                                    &
111        ONLY:  averaging_interval_pr, calc_soil_moisture_during_spinup,        &
112               constant_diffusion, constant_flux_layer, coupling_start_time,   &
113               data_output_during_spinup, dopr_n, do_sum,                      &
114               dt_averaging_input_pr, dt_dopr, dt_dots, dt_do2d_xy, dt_do3d,   &
115               dt_spinup, dt_3d, humidity, intermediate_timestep_count,        &
116               intermediate_timestep_count_max, land_surface,                  &
117               simulated_time, simulated_time_chr, skip_time_dopr,             &
118               skip_time_do2d_xy, skip_time_do3d, spinup_pt_amplitude,         &
119               spinup_pt_mean, spinup_time, timestep_count, time_dopr,         &
120               time_dopr_av, time_dots, time_do2d_xy, time_do3d,               &
121               time_run_control, time_since_reference_point, urban_surface
122
123    USE cpulog,                                                                &
124        ONLY:  cpu_log, log_point, log_point_s
125
126    USE date_and_time_mod,                                                     &
127        ONLY: day_of_year_init, time_utc_init
128
129    USE indices,                                                               &
130        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
131
132
133    USE land_surface_model_mod,                                                &
134        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
135    USE pegrid
136
137    USE pmc_interface,                                                         &
138        ONLY:  nested_run
139
140    USE kinds
141
142    USE radiation_model_mod,                                                   &
143        ONLY:  force_radiation_call, radiation, radiation_control,             &
144               radiation_interaction, radiation_interactions, time_radiation
145
146    USE statistics,                                                            &
147        ONLY:  flow_statistics_called
148
149    USE surface_layer_fluxes_mod,                                              &
150        ONLY:  surface_layer_fluxes
151
152    USE surface_mod,                                                           &
153        ONLY :  surf_lsm_h, surf_lsm_v, surf_usm_h,    &
154                surf_usm_v
155
156    USE urban_surface_mod,                                                     &
157        ONLY:  usm_material_heat_model, usm_material_model,                    &
158               usm_surface_energy_balance, usm_swap_timelevel,                 &
159               usm_green_heat_model
160
161
162
163
164    IMPLICIT NONE
165
166    CHARACTER (LEN=9) ::  time_to_string          !<
167 
168    INTEGER(iwp) ::  i !< running index
169    INTEGER(iwp) ::  j !< running index
170    INTEGER(iwp) ::  k !< running index
171    INTEGER(iwp) ::  l !< running index
172    INTEGER(iwp) ::  m !< running index
173
174    INTEGER(iwp) :: current_timestep_number_spinup = 0  !< number if timestep during spinup
175 
176    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
177
178    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
179    REAL(wp) ::  dt_save     !< temporary storage for time step
180                 
181    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save  !< temporary storage of temperature
182    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_save   !< temporary storage of u wind component
183    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_save   !< temporary storage of v wind component
184
185
186!
187!-- Save 3D arrays because they are to be changed for spinup purpose
188    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
189    ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
190    ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
191
192    CALL exchange_horiz( pt, nbgp )   
193    CALL exchange_horiz( u,  nbgp ) 
194    CALL exchange_horiz( v,  nbgp ) 
195 
196    pt_save = pt
197    u_save  = u
198    v_save  = v
199
200!
201!-- Set the same wall-adjacent velocity to all grid points. The sign of the
202!-- original velocity field must be preserved because the surface schemes crash
203!-- otherwise. The precise reason is still unknown. A minimum velocity of 0.1
204!-- m/s is used to maintain turbulent transfer at the surface.
205    IF ( land_surface )  THEN
206       DO  m = 1, surf_lsm_h%ns
207          i   = surf_lsm_h%i(m)           
208          j   = surf_lsm_h%j(m)
209          k   = surf_lsm_h%k(m)
210          u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
211          v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
212       ENDDO
213
214       DO  l = 0, 3
215          DO  m = 1, surf_lsm_v(l)%ns
216             i   = surf_lsm_v(l)%i(m)           
217             j   = surf_lsm_v(l)%j(m)
218             k   = surf_lsm_v(l)%k(m)
219             u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
220             v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
221          ENDDO
222       ENDDO
223    ENDIF
224
225    IF ( urban_surface )  THEN
226       DO  m = 1, surf_usm_h%ns
227          i   = surf_usm_h%i(m)           
228          j   = surf_usm_h%j(m)
229          k   = surf_usm_h%k(m)
230          u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
231          v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
232       ENDDO
233
234       DO  l = 0, 3
235          DO  m = 1, surf_usm_v(l)%ns
236             i   = surf_usm_v(l)%i(m)           
237             j   = surf_usm_v(l)%j(m)
238             k   = surf_usm_v(l)%k(m)
239             u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
240             v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
241          ENDDO
242       ENDDO
243    ENDIF
244
245    CALL exchange_horiz( u,  nbgp )
246    CALL exchange_horiz( v,  nbgp )
247
248    dt_save = dt_3d
249    dt_3d   = dt_spinup
250
251    CALL location_message( 'starting spinup-sequence', .TRUE. )
252!
253!-- Start of the time loop
254    DO  WHILE ( simulated_time < spinup_time )
255
256       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
257   
258!
259!--    Start of intermediate step loop
260       intermediate_timestep_count = 0
261       DO  WHILE ( intermediate_timestep_count < &
262                   intermediate_timestep_count_max )
263
264          intermediate_timestep_count = intermediate_timestep_count + 1
265
266!
267!--       Set the steering factors for the prognostic equations which depend
268!--       on the timestep scheme
269          CALL timestep_scheme_steering
270
271
272!
273!--       Estimate a near-surface air temperature based on the position of the
274!--       sun and user input about mean temperature and amplitude. The time is
275!--       shifted by one hour to simulate a lag between air temperature and
276!--       incoming radiation
277          pt_spinup = spinup_pt_mean + spinup_pt_amplitude                     &
278             * solar_angle (time_utc_init + time_since_reference_point - 3600.0)
279
280!
281!--       Map air temperature to all grid points in the vicinity of a surface
282!--       element
283          IF ( land_surface )  THEN
284             DO  m = 1, surf_lsm_h%ns
285                i   = surf_lsm_h%i(m)           
286                j   = surf_lsm_h%j(m)
287                k   = surf_lsm_h%k(m)
288                pt(k,j,i) = pt_spinup
289             ENDDO
290
291             DO  l = 0, 3
292                DO  m = 1, surf_lsm_v(l)%ns
293                   i   = surf_lsm_v(l)%i(m)           
294                   j   = surf_lsm_v(l)%j(m)
295                   k   = surf_lsm_v(l)%k(m)
296                   pt(k,j,i) = pt_spinup
297                ENDDO
298             ENDDO
299          ENDIF
300
301          IF ( urban_surface )  THEN
302             DO  m = 1, surf_usm_h%ns
303                i   = surf_usm_h%i(m)           
304                j   = surf_usm_h%j(m)
305                k   = surf_usm_h%k(m)
306                pt(k,j,i) = pt_spinup
307                !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
308                surf_usm_h%pt1 = pt_spinup
309                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310             ENDDO
311
312             DO  l = 0, 3
313                DO  m = 1, surf_usm_v(l)%ns
314                   i   = surf_usm_v(l)%i(m)           
315                   j   = surf_usm_v(l)%j(m)
316                   k   = surf_usm_v(l)%k(m)
317                   pt(k,j,i) = pt_spinup
318                   !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
319                   surf_usm_v(l)%pt1 = pt_spinup
320                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
321                ENDDO
322             ENDDO
323          ENDIF
324
325          CALL exchange_horiz( pt,  nbgp )   
326
327
328!
329!--       Swap the time levels in preparation for the next time step.
330          timestep_count = timestep_count + 1
331     
332          IF ( land_surface )  THEN
333              CALL lsm_swap_timelevel ( 0 )
334          ENDIF
335
336          IF ( urban_surface )  THEN
337             CALL usm_swap_timelevel ( 0 )
338          ENDIF
339
340          IF ( land_surface )  THEN
341             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
342          ENDIF
343
344          IF ( urban_surface )  THEN
345             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
346          ENDIF
347         
348!
349!--       If required, compute virtual potential temperature
350          IF ( humidity )  THEN
351             CALL compute_vpt 
352          ENDIF 
353
354!
355!--       Compute the diffusion quantities
356          IF ( .NOT. constant_diffusion )  THEN
357
358!
359!--          First the vertical (and horizontal) fluxes in the surface
360!--          (constant flux) layer are computed
361             IF ( constant_flux_layer )  THEN
362                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
363                CALL surface_layer_fluxes
364                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
365             ENDIF
366
367!
368!--          If required, solve the energy balance for the surface and run soil
369!--          model. Call for horizontal as well as vertical surfaces.
370!--          The prognostic equation for soil moisure is switched off
371             IF ( land_surface )  THEN
372
373                CALL cpu_log( log_point(54), 'land_surface', 'start' )
374!
375!--             Call for horizontal upward-facing surfaces
376                CALL lsm_energy_balance( .TRUE., -1 )
377                CALL lsm_soil_model( .TRUE., -1, calc_soil_moisture_during_spinup )
378!
379!--             Call for northward-facing surfaces
380                CALL lsm_energy_balance( .FALSE., 0 )
381                CALL lsm_soil_model( .FALSE., 0, calc_soil_moisture_during_spinup )
382!
383!--             Call for southward-facing surfaces
384                CALL lsm_energy_balance( .FALSE., 1 )
385                CALL lsm_soil_model( .FALSE., 1, calc_soil_moisture_during_spinup )
386!
387!--             Call for eastward-facing surfaces
388                CALL lsm_energy_balance( .FALSE., 2 )
389                CALL lsm_soil_model( .FALSE., 2, calc_soil_moisture_during_spinup )
390!
391!--             Call for westward-facing surfaces
392                CALL lsm_energy_balance( .FALSE., 3 )
393                CALL lsm_soil_model( .FALSE., 3, calc_soil_moisture_during_spinup )
394               
395                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
396             ENDIF
397
398!
399!--          If required, solve the energy balance for urban surfaces and run
400!--          the material heat model
401             IF (urban_surface) THEN
402                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
403                CALL usm_surface_energy_balance( .TRUE. )
404                IF ( usm_material_model )  THEN
405                   CALL usm_green_heat_model
406                   CALL usm_material_heat_model( .TRUE. )
407                ENDIF
408               
409                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
410             ENDIF
411
412          ENDIF
413
414!
415!--       If required, calculate radiative fluxes and heating rates
416          IF ( radiation .AND. intermediate_timestep_count                     &
417               == intermediate_timestep_count_max )  THEN
418
419               time_radiation = time_radiation + dt_3d
420
421             IF ( time_radiation >= dt_3d .OR. force_radiation_call )          &
422             THEN
423
424                CALL cpu_log( log_point(50), 'radiation', 'start' )
425
426                IF ( .NOT. force_radiation_call )  THEN
427                   time_radiation = time_radiation - dt_3d
428                ENDIF
429
430                CALL radiation_control
431
432                CALL cpu_log( log_point(50), 'radiation', 'stop' )
433
434                IF ( radiation_interactions )  THEN
435                   CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
436                   CALL radiation_interaction
437                   CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
438                ENDIF
439             ENDIF
440          ENDIF
441
442       ENDDO   ! Intermediate step loop
443
444!
445!--    Increase simulation time and output times
446       current_timestep_number_spinup = current_timestep_number_spinup + 1
447       simulated_time             = simulated_time   + dt_3d
448       simulated_time_chr         = time_to_string( simulated_time )
449       time_since_reference_point = simulated_time - coupling_start_time
450
451       IF ( data_output_during_spinup )  THEN
452          IF ( simulated_time >= skip_time_do2d_xy )  THEN
453             time_do2d_xy       = time_do2d_xy     + dt_3d
454          ENDIF
455          IF ( simulated_time >= skip_time_do3d    )  THEN
456             time_do3d          = time_do3d        + dt_3d
457          ENDIF
458          time_dots          = time_dots        + dt_3d
459          IF ( simulated_time >= skip_time_dopr )  THEN
460             time_dopr       = time_dopr        + dt_3d
461          ENDIF
462          time_run_control   = time_run_control + dt_3d
463
464!
465!--       Carry out statistical analysis and output at the requested output times.
466!--       The MOD function is used for calculating the output time counters (like
467!--       time_dopr) in order to regard a possible decrease of the output time
468!--       interval in case of restart runs
469
470!
471!--       Set a flag indicating that so far no statistics have been created
472!--       for this time step
473          flow_statistics_called = .FALSE.
474
475!
476!--       If required, call flow_statistics for averaging in time
477          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
478             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
479             simulated_time >= skip_time_dopr )  THEN
480             time_dopr_av = time_dopr_av + dt_3d
481             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
482                do_sum = .TRUE.
483                time_dopr_av = MOD( time_dopr_av,                              &
484                               MAX( dt_averaging_input_pr, dt_3d ) )
485             ENDIF
486          ENDIF
487          IF ( do_sum )  CALL flow_statistics
488
489!
490!--       Output of profiles
491          IF ( time_dopr >= dt_dopr )  THEN
492             IF ( dopr_n /= 0 )  CALL data_output_profiles
493             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
494             time_dopr_av = 0.0_wp    ! due to averaging (see above)
495          ENDIF
496
497!
498!--       Output of time series
499          IF ( time_dots >= dt_dots )  THEN
500             CALL data_output_tseries
501             time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
502          ENDIF
503
504!
505!--       2d-data output (cross-sections)
506          IF ( time_do2d_xy >= dt_do2d_xy )  THEN
507             CALL data_output_2d( 'xy', 0 )
508             time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
509          ENDIF
510
511!
512!--       3d-data output (volume data)
513          IF ( time_do3d >= dt_do3d )  THEN
514             CALL data_output_3d( 0 )
515             time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
516          ENDIF
517
518
519       ENDIF
520
521!
522!--    Computation and output of run control parameters.
523!--    This is also done whenever perturbations have been imposed
524!        IF ( time_run_control >= dt_run_control  .OR.                           &
525!             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
526!        THEN
527!           CALL run_control
528!           IF ( time_run_control >= dt_run_control )  THEN
529!              time_run_control = MOD( time_run_control,                         &
530!                                      MAX( dt_run_control, dt_3d ) )
531!           ENDIF
532!        ENDIF
533
534       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
535
536
537!
538!--    Run control output
539       IF ( myid == 0 )  THEN
540!
541!--       If necessary, write header
542          IF ( .NOT. run_control_header_spinup )  THEN
543             CALL check_open( 15 )
544             WRITE ( 15, 100 )
545             run_control_header_spinup = .TRUE.
546          ENDIF
547!
548!--       Write some general information about the spinup in run control file
549          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_3d, pt_spinup
550!
551!--       Write buffer contents to disc immediately
552          FLUSH( 15 )
553       ENDIF
554
555
556
557    ENDDO   ! time loop
558
559!
560!-- Write back saved arrays to the 3D arrays
561    pt   = pt_save
562    pt_p = pt_save
563    u    = u_save
564    v    = v_save
565
566!
567!-- Reset time step
568    dt_3d = dt_save
569
570    DEALLOCATE(pt_save)
571    DEALLOCATE(u_save)
572    DEALLOCATE(v_save)
573
574#if defined( __parallel )
575    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
576#endif
577
578    CALL location_message( 'finished spinup-sequence', .TRUE. )
579
580
581!
582!-- Formats
583100 FORMAT (///'Spinup control output:'/  &
584            '--------------------------------'// &
585            'ITER.  HH:MM:SS    DT   PT(z_MO)'/   &
586            '--------------------------------')
587101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
588
589 CONTAINS
590
591!
592!-- Returns the cosine of the solar zenith angle at a given time. This routine
593!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
594    FUNCTION solar_angle( local_time ) 
595
596       USE basic_constants_and_equations_mod,                                  &
597       ONLY:  pi
598     
599       USE kinds
600
601       USE radiation_model_mod,                                                &
602           ONLY:  decl_1, decl_2, decl_3, lat, lon
603
604       IMPLICIT NONE
605
606
607       REAL(wp) ::  solar_angle     !< cosine of the solar zenith angle
608
609       REAL(wp) ::  day             !< day of the year
610       REAL(wp) ::  declination     !< solar declination angle
611       REAL(wp) ::  hour_angle      !< solar hour angle
612       REAL(wp) ::  time_utc        !< current time in UTC
613       REAL(wp), INTENT(IN) :: local_time
614!
615!--    Calculate current day and time based on the initial values and simulation
616!--    time
617
618       day = day_of_year_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
619       time_utc = MOD(local_time, 86400.0_wp)
620
621
622!
623!--    Calculate solar declination and hour angle   
624       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
625       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
626
627!
628!--    Calculate cosine of solar zenith angle
629       solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
630                     * COS(hour_angle)
631
632
633    END FUNCTION solar_angle
634
635
636 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.