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

Last change on this file since 3395 was 3337, checked in by kanani, 6 years ago

reintegrate branch resler to trunk

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