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

Last change on this file since 3795 was 3766, checked in by raasch, 5 years ago

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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