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

Last change on this file since 4178 was 4064, checked in by gronemeier, 2 years ago

Moved call to radiation module out of intermediate time loop (time_integration, time_integration_spinup)

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