FMS 2025.01.02-dev
Flexible Modeling System
Loading...
Searching...
No Matches
astronomy.F90
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @defgroup astronomy_mod astronomy_mod
20!> @ingroup astronomy
21!> @brief Provides astronomical variables for use
22!! by other modules within fms. The only currently used interface is
23!! for determination of astronomical values needed by the shortwave
24!! radiation packages.
25!> @author Fei Liu
26
27!> @addtogroup astronomy_mod
28!> @{
29module astronomy_mod
30
31
32 use fms_mod, only: fms_init, &
33 mpp_pe, mpp_root_pe, stdlog, &
34 write_version_number, &
35 check_nml_error, error_mesg, &
36 fatal, note, warning
37 use time_manager_mod, only: time_type, set_time, get_time, &
38 get_date_julian, set_date_julian, &
41 operator(-), operator(+), &
42 operator( // ), operator(<)
43 use fms_io_utils_mod, only: get_data_type_string
44 use constants_mod, only: constants_init, pi
45 use mpp_mod, only: input_nml_file
46 use platform_mod, only: r4_kind, r8_kind, i4_kind, i8_kind
47
48 !--------------------------------------------------------------------
49
50 implicit none
51 private
52
53 !---------------------------------------------------------------------
54 !----------- version number for this module --------------------------
55
56! Include variable "version" to be written to log file.
57#include<file_version.h>
58
59
60 !---------------------------------------------------------------------
61 !------- interfaces --------
62
63 public &
69
71 module procedure set_orbital_parameters_r4, set_orbital_parameters_r8
72 end interface set_orbital_parameters
73
75 module procedure get_orbital_parameters_r4, get_orbital_parameters_r8
76 end interface get_orbital_parameters
77
78 !> @}
79
80 !> @brief Calculates solar information for the given location(lat & lon) and time
81 !!
82 !> ~~~~~~~~~~{.f90}
83 !! call diurnal_solar (lat, lon, time, cosz, fracday, rrsun, dt_time)
84 !! call diurnal_solar (lat, lon, gmt, time_since_ae, cosz, fracday, rrsun, dt)
85 !! ~~~~~~~~~~
86 !!
87 !! The first option (used in conjunction with time_manager_mod)
88 !! generates the real variables gmt and time_since_ae from the
89 !! time_type input, and then calls diurnal_solar with these real inputs.
90 !!
91 !! The time of day is set by
92 !! ~~~~~~~~~~{.f90}
93 !! real, intent(in) :: gmt
94 !! ~~~~~~~~~~
95 !! The time of year is set by
96 !! ~~~~~~~~~~{.f90}
97 !! real, intent(in) :: time_since_ae
98 !! ~~~~~~~~~~
99 !! with time_type input, both of these are extracted from
100 !! ~~~~~~~~~~{.f90}
101 !! type(time_type), intent(in) :: time
102 !! ~~~~~~~~~~
103 !!
104 !! Separate routines exist within this interface for scalar,
105 !! 1D or 2D input and output fields:
106 !!
107 !! ~~~~~~~~~~{.f90}
108 !! real, intent(in), dimension(:,:) :: lat, lon
109 !! real, intent(in), dimension(:) :: lat, lon
110 !! real, intent(in) :: lat, lon
111 !!
112 !! real, intent(out), dimension(:,:) :: cosz, fracday
113 !! real, intent(out), dimension(:) :: cosz, fracday
114 !! real, intent(out) :: cosz, fracday
115 !! ~~~~~~~~~~
116 !!
117 !! One may also average the output fields over the time interval
118 !! between gmt and gmt + dt by including the optional argument dt (or
119 !! dt_time). dt is measured in radians and must be less than pi
120 !! (1/2 day). This average is computed analytically, and should be
121 !! exact except for the fact that changes in earth-sun distance over
122 !! the time interval dt are ignored. In the context of a diurnal GCM,
123 !! this option should always be employed to insure that the total flux
124 !! at the top of the atmosphere is not modified by time truncation error.
125 !!
126 !! ~~~~~~~~~~{.f90}
127 !! real, intent(in), optional :: dt
128 !! type(time_type), optional :: dt_time
129 !! ~~~~~~~~~~
130 !!
131 !! @param [in] <lat> Latitudes of model grid points [radians]
132 !! @param [in] <lon> Longitudes of model grid points [radians]
133 !! @param [in] <gmt> Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi [radians]
134 !! @param [in] <time_since_ae> Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
135 !! @param [in] <time> Time at which astronomical values are desired (time_type variable) [seconds, days]
136 !! @param [out] <cosz> Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
137 !! @param [out] <fracday> Daylight fraction of time interval [dimensionless]
138 !! @param [out] <rrsun> Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
139 !! (a):(a/r)**2 [dimensionless]
140 !! @param [in] <dt> OPTIONAL: Time interval after gmt over which the astronomical variables are to be
141 !! averaged. this produces averaged output rather than instantaneous. [radians], (1 day = 2 * pi)
142 !! @param [in] <dt_time> OPTIONAL: Time interval after gmt over which the astronomical variables are to be
143 !! averaged. this produces averaged output rather than instantaneous. time_type,
144 !! [days, seconds]
145 !! @param [in] <allow_negative_cosz> Allow negative values for cosz?
146 !! @param [out] <half_day_out> half_day_out
147 !> @ingroup astronomy_mod
149 module procedure diurnal_solar_2d_r4, diurnal_solar_2d_r8
150 module procedure diurnal_solar_1d_r4, diurnal_solar_1d_r8
151 module procedure diurnal_solar_0d_r4, diurnal_solar_0d_r8
152 module procedure diurnal_solar_cal_2d_r4, diurnal_solar_cal_2d_r8
153 module procedure diurnal_solar_cal_1d_r4, diurnal_solar_cal_1d_r8
154 module procedure diurnal_solar_cal_0d_r4, diurnal_solar_cal_0d_r8
155 end interface diurnal_solar
156
157 !> @brief Calculates the daily mean solar information for a given time and latitude.
158 !!
159 !> ~~~~~~~~~~{.f90}
160 !! call daily_mean_solar (lat, time, cosz, fracday, rrsun)
161 !! call daily_mean_solar (lat, time_since_ae, cosz, fracday, rrsun)
162 !! call daily_mean_solar (lat, time, cosz, solar)
163 !! call daily_mean_solar (lat, time_since_ae, cosz, solar)
164 !! ~~~~~~~~~~
165 !!
166 !! The first option (used in conjunction with time_manager_mod)
167 !! generates the real variable time_since_ae from the time_type
168 !! input time, and then calls daily_mean_solar with this real input
169 !! (option 2). The third and fourth options correspond to the first
170 !! and second and are used with then spectral 2-layer model, where
171 !! only cosz and solar are desired as output. These routines generate
172 !! dummy arguments and then call option 2, where the calculation is done.
173 !!
174 !! The time of year is set by
175 !! ~~~~~~~~~~{.f90}
176 !! real, intent(in) :: time_since_ae
177 !! ~~~~~~~~~~
178 !! With time_type input, the time of year is extracted from
179 !! ~~~~~~~~~~{.f90}
180 !! type(time_type), intent(in) :: time
181 !! ~~~~~~~~~~
182 !!
183 !! Separate routines exist within this interface for scalar,
184 !! 1D or 2D input and output fields:
185 !!
186 !! ~~~~~~~~~~{.f90}
187 !! real, intent(in), dimension(:,:) :: lat
188 !! real, intent(in), dimension(:) :: lat
189 !! real, intent(in) :: lat
190 !!
191 !! real, intent(out), dimension(:,:) :: cosz, fracday
192 !! real, intent(out), dimension(:) :: cosz, fracday
193 !! real, intent(out) :: cosz, fracday
194 !! ~~~~~~~~~~
195 !!
196 !! @param [in] <lat> Latitudes of model grid points [radians]
197 !! @param [in] <time_since_ae> Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
198 !! @param [in] <time> Time at which astronomical values are desired (time_type variable) [seconds, days]
199 !! @param [out] <cosz> Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
200 !! @param [out] <fracday> Daylight fraction of time interval [dimensionless]
201 !! @param [out] <rrsun> Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
202 !! (a):(a/r)**2 [dimensionless]
203 !! @param [out] <solar> shortwave flux factor: cosine of zenith angle * daylight fraction /
204 !! (earth-sun distance squared) [dimensionless]
205 !> @ingroup astronomy_mod
207 module procedure daily_mean_solar_2d_r4, daily_mean_solar_2d_r8
208 module procedure daily_mean_solar_1d_r4, daily_mean_solar_1d_r8
209 module procedure daily_mean_solar_2level_r4, daily_mean_solar_2level_r8
210 module procedure daily_mean_solar_0d_r4, daily_mean_solar_0d_r8
211 module procedure daily_mean_solar_cal_2d_r4, daily_mean_solar_cal_2d_r8
212 module procedure daily_mean_solar_cal_1d_r4, daily_mean_solar_cal_1d_r8
213 module procedure daily_mean_solar_cal_2level_r4, daily_mean_solar_cal_2level_r8
214 module procedure daily_mean_solar_cal_0d_r4, daily_mean_solar_cal_0d_r8
215 end interface daily_mean_solar
216
217 !> Calculates the annual mean of solar information for a given latitude and time.
218 !!
219 !> ~~~~~~~~~~{.f90}
220 !! call annual_mean_solar (js, je, lat, cosz, solar, fracday, rrsun)
221 !! call annual_mean_solar (lat, cosz, solar)
222 !! ~~~~~~~~~~
223 !!
224 !! The second interface above is used by the spectral 2-layer model,
225 !! which requires only cosz and solar as output arguments, and which
226 !! makes this call during the initialization phase of the model.
227 !! Separate routines exist within this interface for 1D or 2D input
228 !! and output fields:
229 !!
230 !! ~~~~~~~~~~{.f90}
231 !! real, intent(in), dimension(:,:) :: lat
232 !! real, intent(in), dimension(:) :: lat
233 !!
234 !! real, intent(out), dimension(:,:) :: cosz, solar, fracday
235 !! real, intent(out), dimension(:) :: cosz, solar, fracday
236 !! ~~~~~~~~~~
237 !!
238 !! @param [in] <jst> Starting subdomain j indices of data in the physics wiondow being integrated
239 !! @param [in] <jnd> Ending subdomain j indices of data in the physics wiondow being integrated
240 !! @param [in] <lat> Latitudes of model grid points [radians]
241 !! @param [out] <cosz> cosz is the average over the year of the cosine of an effective zenith angle
242 !! that would produce the correct daily solar flux if the sun were fixed at that
243 !! single position for the period of daylight on the given day. in this average,
244 !! the daily mean effective cosz is weighted by the daily mean solar flux. [dimensionless]
245 !! @param [out] <solar> Normalized solar flux, averaged over the year, equal to the product of
246 !! fracday*cosz*rrsun [dimensionless]
247 !! @param [out] <fracday> Daylight fraction calculated so as to make the average flux (solar) equal to the
248 !! product of the flux-weighted avg cosz * this fracday * assumed annual mean avg
249 !! Earth-Sun distance of 1.0. [dimensionless]
250 !! @param [out] <rrsun> Annual mean Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
251 !! (a):(a/r)**2 [dimensionless]
252 !> @ingroup astronomy_mod
254 module procedure annual_mean_solar_2d_r4, annual_mean_solar_2d_r8
255 module procedure annual_mean_solar_1d_r4, annual_mean_solar_1d_r8
256 module procedure annual_mean_solar_2level_r4, annual_mean_solar_2level_r8
257 end interface annual_mean_solar
258
259 !> Gets the length of year for current calendar
260 !!
261 !> ~~~~~~~~~~{.f90}
262 !! call get_period (period)
263 !! ~~~~~~~~~~
264 !!
265 !! Separate routines exist within this interface for integer
266 !! and time_type output:
267 !!
268 !! ~~~~~~~~~~{.f90}
269 !! integer, intent(out) :: period
270 !! type(time_type), intent(out) :: period
271 !! ~~~~~~~~~~
272 !!
273 !! @param [out] <period_out> Length of year for calendar in use
274 !> @ingroup astronomy_mod
275 interface get_period
277 end interface
278
279 !> Sets the length of a year for the calendar in use
280 !!
281 !> ~~~~~~~~~~{.f90}
282 !! call set_period (period_in)
283 !! ~~~~~~~~~~
284 !!
285 !! Separate routines exist within this interface for integer
286 !! and time_type output:
287 !!
288 !! ~~~~~~~~~~{.f90}
289 !! integer, intent(out) :: period_in
290 !! type(time_type), intent(out) :: period_in
291 !! ~~~~~~~~~~
292 !!
293 !! @param [in] <period_in> Length of year for calendar in use
294 !> @ingroup astronomy_mod
295 interface set_period
297 end interface
298
299
300 private &
301 orbit, & ! Called from astronomy_init and set_orbital_parameters
302 r_inv_squared, & ! Called from diurnal_solar, daily_mean_solar and orbit
303 angle, declination, half_day ! called from diurnal_solar and daily_mean_solar
304 ! half_day, orbital_time, & ! called from diurnal_solar and daily_mean_solar
305 ! universal_time ! called from diurnal_solar:
306
308 module procedure r_inv_squared_r4, r_inv_squared_r8
309 end interface r_inv_squared
310
311 interface angle
312 module procedure angle_r4, angle_r8
313 end interface angle
314
315 interface declination
316 module procedure declination_r4, declination_r8
317 end interface declination
318
319 !> Private interface for internal use by dirunal_solar and daily_mean_solar.
320 !!
321 !> Example usage:
322 !! ~~~~~~~~~~{.f90}
323 !! half_day (latitude, dec) result (h)
324 !! ~~~~~~~~~~
325 !!
326 !! Separate routines exist within this interface for scalar,
327 !! or 2D input and output fields:
328 !!
329 !! ~~~~~~~~~~{.f90}
330 !! real, intent(in), dimension(:,:) :: latitude
331 !! real, intent(in) :: latitude
332 !!
333 !! real, dimension(size(latitude,1),size(latitude,2)) :: h
334 !! real :: h
335 !! ~~~~~~~~~~
336 !!
337 !! @param [in] <latitude> Latitudes of model grid points [radians]
338 !! @param [in] <dec> Solar declination [radians]
339 !! @param [out] <h> Half of the length of daylight at the given latitude and orbital position (dec); value
340 !! ranges between 0 (all darkness) and pi (all daylight) [dimensionless]
341 !> @ingroup astronomy_mod
342 interface half_day
343 module procedure half_day_2d_r4, half_day_2d_r8
344 module procedure half_day_0d_r4, half_day_0d_r8
345 end interface half_day
346
347
348!> @addtogroup astronomy_mod
349!> @{
350
351!---------------------------------------------------------------------
352!-------- namelist ---------
353
354real(r8_kind) :: ecc = 0.01671_r8_kind !< Eccentricity of Earth's orbit [dimensionless]
355real(r8_kind) :: obliq = 23.439_r8_kind !< Obliquity [degrees]
356real(r8_kind) :: per = 102.932_r8_kind !< Longitude of perihelion with respect
357 !! to autumnal equinox in NH [degrees]
358integer :: period = 0 !< Specified length of year [seconds];
359 !! must be specified to override default
360 !! value given by length_of_year in
361 !! time_manager_mod
362integer :: day_ae = 23 !< Day of specified autumnal equinox
363integer :: month_ae = 9 !< Month of specified autumnal equinox
364integer :: year_ae = 1998 !< Year of specified autumnal equinox
365integer :: hour_ae = 5 !< Hour of specified autumnal equinox
366integer :: minute_ae = 37 !< Minute of specified autumnal equinox
367integer :: second_ae = 0 !< Second of specified autumnal equinox
368integer :: num_angles = 3600 !< Number of intervals into which the year
369 !! is divided to compute orbital positions
370
371
372namelist /astronomy_nml/ ecc, obliq, per, period, &
376
377!--------------------------------------------------------------------
378!------ public data ----------
379
380
381!--------------------------------------------------------------------
382!------ private data ----------
383
384type(time_type) :: autumnal_eq_ref !< time_type variable containing
385 !! specified time of reference
386 !! NH autumnal equinox
387
388type(time_type) :: period_time_type !< time_type variable containing
389 !! period of one orbit
390
391real(r8_kind), dimension(:), allocatable :: orb_angle !< table of orbital positions (0 to
392 !! 2*pi) as a function of time used
393 !! to find actual orbital position
394 !! via interpolation
395
396real(r8_kind) :: seconds_per_day = 86400.0_r8_kind !< seconds in a day
397real(r8_kind) :: deg_to_rad !< conversion from degrees to radians
398real(r8_kind) :: twopi !< 2 *PI
399logical :: module_is_initialized=.false. !< has the module been initialized ?
400
401real(r8_kind), dimension(:,:), allocatable :: &
402 cosz_ann, & !< annual mean cos of zenith angle
403 solar_ann, & !< annual mean solar factor
404 fracday_ann !< annual mean daylight fraction
405real(r8_kind) :: rrsun_ann !< annual mean earth-sun distance
406logical :: annual_mean_calculated=.false. !< have the annual mean values been calculated?
407integer :: num_pts = 0 !< count of grid_boxes for which
408 !! annual mean astronomy values have
409 !! been calculated
410integer :: total_pts !< number of grid boxes owned by the processor
411
412!--------------------------------------------------------------------
413
414
415
416contains
417
418
419
420!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
421!
422! PUBLIC SUBROUTINES
423!
424!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425
426!> @brief astronomy_init is the constructor for astronomy_mod.
427!!
428!! @throw FATAL, "astronomy_mod ecc must be between 0 and 0.99"
429!! @throw FATAL, "astronomy_mod obliquity must be between -90 and 90 degrees"
430!! @throw FATAL, "astronomy_mod perihelion must be between 0 and 360 degrees"
431subroutine astronomy_init (latb, lonb)
432
433class(*), dimension(:,:), intent(in), optional :: latb !< 2d array of model latitudes at cell corners [radians]
434class(*), dimension(:,:), intent(in), optional :: lonb !< 2d array of model longitudes at cell corners [radians]
435logical :: is_valid
436
437!-------------------------------------------------------------------
438! local variables:
439!-------------------------------------------------------------------
440integer :: iunit, ierr, io, seconds, days, jd, id
441character(len=17) :: err_str
442
443!-------------------------------------------------------------------
444! if module has already been initialized, exit.
445!-------------------------------------------------------------------
446 if (module_is_initialized) return
447
448!-------------------------------------------------------------------
449!> This routine will:
450!> Verify that modules used by this module have been initialized.
451!-------------------------------------------------------------------
452 call fms_init
454 call constants_init
455
456!-----------------------------------------------------------------------
457!> Read namelist.
458!-----------------------------------------------------------------------
459 read (input_nml_file, astronomy_nml, iostat=io)
460 ierr = check_nml_error(io,'astronomy_nml')
461!---------------------------------------------------------------------
462!> Write version number and namelist to logfile.
463!---------------------------------------------------------------------
464 call write_version_number("ASTRONOMY_MOD", version)
465 if (mpp_pe() == mpp_root_pe() ) then
466 iunit = stdlog()
467 write (iunit, nml=astronomy_nml)
468 endif
469!--------------------------------------------------------------------
470!> Be sure input values are within valid ranges.
471!---------------------------------------------------------------------
472 if (ecc < 0.0_r8_kind .or. ecc > 0.99_r8_kind) &
473 call error_mesg ('astronomy_mod', &
474 'ecc must be between 0 and 0.99', fatal)
475 if (obliq < -90.0_r8_kind .or. obliq > 90.0_r8_kind) &
476 call error_mesg ('astronomy_mod', &
477 'obliquity must be between -90 and 90 degrees', fatal)
478 if (per < 0.0_r8_kind .or. per > 360.0_r8_kind) &
479 call error_mesg ('astronomy_mod', &
480 'perihelion must be between 0 and 360 degrees', fatal)
481
482!----------------------------------------------------------------------
483!> Set up time-type variable defining specified time of autumnal equinox.
484!----------------------------------------------------------------------
487
488!---------------------------------------------------------------------
489!> Set up time-type variable defining length of year.
490!----------------------------------------------------------------------
491 if (period == 0) then
493 call get_time (period_time_type, seconds, days)
494 period = int(seconds_per_day) * days + seconds
495 else
497 endif
498
499!---------------------------------------------------------------------
500!> Define useful module variables.
501!----------------------------------------------------------------------
502 twopi = 2.0_r8_kind * pi
503 deg_to_rad = twopi/360.0_r8_kind
504
505!---------------------------------------------------------------------
506!> Call orbit to define table of orbital angles as function of
507!! orbital time.
508!----------------------------------------------------------------------
509! wfc moved here from orbit
510 allocate ( orb_angle(0:num_angles) )
511 call orbit
512
513!--------------------------------------------------------------------
514!> If annual mean radiation is desired, then latb will be present.
515!! allocate arrays to hold the needed astronomical factors. define
516!! the total number of points that the processor is responsible for.
517!--------------------------------------------------------------------
518 ! check that no invalid types (integers or characters) are given as optional arg
519
520 is_valid = .false.
521 if (present(latb) .and. present(lonb)) then
522 select type (latb)
523 type is (real(r4_kind))
524 select type (lonb)
525 type is (real(r4_kind))
526 is_valid = .true.
527 class default
528 call get_data_type_string(lonb, err_str)
529 call error_mesg('astronomy_mod', 'kind mismatch, argument latb is real(r4_kind) but lonb has type: '// &
530 err_str, fatal)
531 end select
532 type is (real(r8_kind))
533 select type (lonb)
534 type is (real(r8_kind))
535 is_valid = .true.
536 class default
537 call get_data_type_string(lonb, err_str)
538 call error_mesg('astronomy_mod', 'kind mismatch, argument latb is real(r8_kind) but lonb has type: '//&
539 err_str, fatal)
540 end select
541 end select
542 if( is_valid ) then
543 jd = size(latb,2) - 1
544 id = size(lonb,1) - 1
545 allocate (cosz_ann(id, jd))
546 allocate (solar_ann(id, jd))
547 allocate (fracday_ann(id, jd))
548 total_pts = jd*id
549 else
550 call error_mesg('astronomy_mod', 'latb has unsupported kind size.' // &
551 'latb and lonb should both be real(r4_kind) or real(r8_kind)', fatal)
552 end if
553 elseif ( (present(latb) .and. .not. present(lonb)) .or. (present(lonb) .and. .not. present(latb)) ) then
554 call error_mesg ('astronomy_mod', 'lat and lon must both be present', fatal)
555 endif
556
557!---------------------------------------------------------------------
558!> Mark the module as initialized.
559!---------------------------------------------------------------------
561
562!---------------------------------------------------------------------
563
564end subroutine astronomy_init
565
566
567!> @brief get_period_integer returns the length of the year as an
568!! integer number of seconds.
569!!
570!! @throw FATAL, "astronomy_mod module has not been initialized"
571subroutine get_period_integer (period_out)
572
573integer, intent(out) :: period_out !< Length of year [seconds]
574
575!--------------------------------------------------------------------
576! local variables:
577
578integer :: seconds, days
579
580!---------------------------------------------------------------------
581! exit if module has not been initialized.
582!---------------------------------------------------------------------
583 if (.not. module_is_initialized) &
584 call error_mesg ( 'astronomy_mod', ' module has not been initialized', fatal)
585
586!--------------------------------------------------------------------
587! define length of year in seconds.
588!--------------------------------------------------------------------
589 call get_time (period_time_type, seconds, days)
590 period_out = int(seconds_per_day) * days + seconds
591
592
593end subroutine get_period_integer
594
595!> @brief get_period_time_type returns the length of the year as a time_type
596!! variable.
597!!
598!! @throw FATAL, "astronomy_mod module has not been initialized"
599subroutine get_period_time_type (period_out)
600
601type(time_type), intent(inout) :: period_out !< Length of year as time_type variable
602
603!---------------------------------------------------------------------
604! exit if module has not been initialized.
605!---------------------------------------------------------------------
606 if (.not. module_is_initialized) &
607 call error_mesg ('astronomy_mod', 'module has not been initialized', fatal)
608
609!--------------------------------------------------------------------
610! define length of year as a time_type variable.
611!--------------------------------------------------------------------
612 period_out = period_time_type
613
614end subroutine get_period_time_type
615
616!> @brief set_period_integer saves as the input length of the year (an
617!! integer) in a time_type module variable.
618!!
619!! @throw FATAL, "astronomy_mod module has not been initialized"
620subroutine set_period_integer (period_in)
621
622integer, intent(in) :: period_in !< Length of year as a time_type
623
624!---------------------------------------------------------------------
625! exit if module has not been initialized.
626!---------------------------------------------------------------------
627 if (.not. module_is_initialized) &
628 call error_mesg ('astronomy_mod', 'module has not been initialized', fatal)
629
630!---------------------------------------------------------------------
631! define time_type variable defining the length of year from the
632! input value (integer seconds).
633!---------------------------------------------------------------------
634 period_time_type = set_time(period_in, 0)
635
636end subroutine set_period_integer
637
638
639!> @brief Set_period_time_type saves the length of the year (input as a
640!! time_type variable) into a time_type module variable.
641!!
642!! @throw FATAL, "astronomy_mod module has not been initialized"
643subroutine set_period_time_type(period_in)
644
645type(time_type), intent(in) :: period_in
646
647!---------------------------------------------------------------------
648! exit if module has not been initialized.
649!---------------------------------------------------------------------
650 if (.not. module_is_initialized) &
651 call error_mesg ('astronomy_mod', 'module has not been initialized', fatal)
652
653!---------------------------------------------------------------------
654! define time_type variable defining the length of year from the
655! input value (time_type).
656!---------------------------------------------------------------------
657 period_time_type = period_in
658
659
660end subroutine set_period_time_type
661
662!> @brief set_ref_date_of_ae provides a means of specifying the reference
663!! date of the NH autumnal equinox for a particular year.
664!!
665!! @details set_ref_date_of_ae provides a means of specifying the reference
666!! date of the NH autumnal equinox for a particular year. It is only
667!! used if calls are made to the calandar versions of the routines
668!! diurnal_solar and daily_mean_solar. If the NOLEAP calendar is
669!! used, then the date of autumnal equinox will be the same every
670!! year. If JULIAN is used, then the date of autumnal equinox will
671!! return to the same value every 4th year.
672!!
673!! @param [in] <day_in> Day of reference autumnal equinox
674!! @param [in] <month_in> Month of reference autumnal equinox
675!! @param [in] <year_in> Year of reference autumnal equinox
676!! @param [out] <second_in> OPTIONAL: Second of reference autumnal equinox
677!! @param [out] <minute_in> OPTIONAL: Minute of reference autumnal equinox
678!! @param [out] <hour_in> OPTIONAL: Hour of reference autumnal equinox
679!!
680!! @throw FATAL, "astronomy_mod module has not been initialized"
681
682subroutine set_ref_date_of_ae (day_in,month_in,year_in, &
683 second_in,minute_in,hour_in)
684
685integer, intent(in) :: day_in, month_in, year_in
686integer, intent(in), optional :: second_in, minute_in, hour_in
687
688!---------------------------------------------------------------------
689! exit if module has not been initialized.
690!---------------------------------------------------------------------
691 if (.not. module_is_initialized) &
692 call error_mesg ('astronomy_mod', 'module has not been initialized', fatal)
693
694!--------------------------------------------------------------------
695! save the input time of ae specification into a time_type module
696! variable autumnal_eq_ref.
697!--------------------------------------------------------------------
698 day_ae = day_in
699 month_ae = month_in
700 year_ae = year_in
701
702 if (present(second_in)) then
703 second_ae = second_in
704 minute_ae = minute_in
705 hour_ae = hour_in
706 else
707 second_ae = 0
708 minute_ae = 0
709 hour_ae = 0
710 endif
711
714
715!---------------------------------------------------------------------
716
717
718end subroutine set_ref_date_of_ae
719
720
721!> @brief get_ref_date_of_ae retrieves the reference date of the autumnal
722!! equinox as integer variables.
723!!
724!! @throw FATAL, "astronomy_mod module has not been initialized"
725!!
726!! @param [out] <day_out> Day of reference autumnal equinox
727!! @param [out] <month_out> Month of reference autumnal equinox
728!! @param [out] <year_out> Year of reference autumnal equinox
729!! @param [out] <second_out> Second of reference autumnal equinox
730!! @param [out] <minute_out> Minute of reference autumnal equinox
731!! @param [out] <hour_out> Hour of reference autumnal equinox
732subroutine get_ref_date_of_ae (day_out,month_out,year_out,&
733 second_out,minute_out,hour_out)
734
735integer, intent(out) :: day_out, month_out, year_out, &
736 second_out, minute_out, hour_out
737
738!---------------------------------------------------------------------
739! exit if module has not been initialized.
740!---------------------------------------------------------------------
741 if (.not. module_is_initialized) &
742 call error_mesg ('astronomy_mod', 'module has not been initialized', fatal)
743
744!---------------------------------------------------------------------
745! fill the output fields with the proper module data.
746!---------------------------------------------------------------------
747 day_out = day_ae
748 month_out = month_ae
749 year_out = year_ae
750 second_out = second_ae
751 minute_out = minute_ae
752 hour_out = hour_ae
753
754
755end subroutine get_ref_date_of_ae
756
757!> @brief astronomy_end is the destructor for astronomy_mod.
759
760 !----------------------------------------------------------------------
761 !> check if the module has been initialized.
762 !----------------------------------------------------------------------
763 if (.not. module_is_initialized) return
764 ! call error_mesg ( 'astronomy_mod', &
765 ! ' module has not been initialized', FATAL)
766
767 !----------------------------------------------------------------------
768 !> deallocate module variables.
769 !----------------------------------------------------------------------
770 if (allocated(orb_angle)) deallocate (orb_angle)
771 if (allocated(cosz_ann) ) then
772 deallocate (cosz_ann)
773 deallocate (fracday_ann)
774 deallocate (solar_ann)
775 endif
776
777 !----------------------------------------------------------------------
778 !> mark the module as uninitialized.
779 !----------------------------------------------------------------------
780 module_is_initialized = .false.
781
782end subroutine astronomy_end
783
784!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
785!
786! PRIVATE SUBROUTINES
787!
788!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
789
790!> @brief Orbit computes and stores a table of value of orbital angles as a
791!! function of orbital time (both the angle and time are zero at
792!! autumnal equinox in the NH, and range from 0 to 2*pi).
793subroutine orbit
794!---------------------------------------------------------------------
795! local variables
796
797integer :: n
798real(kind=r8_kind) :: d1, d2, d3, d4, d5, dt, norm
799
800!--------------------------------------------------------------------
801!> allocate the orbital angle array, sized by the namelist parameter
802!! num_angles, defining the annual cycle resolution of the earth's
803!! orbit. define some constants to be used.
804!--------------------------------------------------------------------
805! wfc moving to astronomy_init
806! allocate ( orb_angle(0:num_angles) )
807orb_angle(0) = 0.0_r8_kind
808dt = twopi/real(num_angles, r8_kind)
809norm = sqrt(1.0_r8_kind - ecc**2)
810dt = dt*norm
811
812!---------------------------------------------------------------------
813!> define the orbital angle at each of the num_angles locations in
814!! the orbit.
815!---------------------------------------------------------------------
816 do n = 1, num_angles
817 d1 = dt*r_inv_squared(orb_angle(n-1))
818 d2 = dt*r_inv_squared(orb_angle(n-1) + 0.5_r8_kind * d1)
819 d3 = dt*r_inv_squared(orb_angle(n-1) + 0.5_r8_kind * d2)
820 d4 = dt*r_inv_squared(orb_angle(n-1) + d3)
821 d5 = d1/6.0_r8_kind + d2/3.0_r8_kind + d3/3.0_r8_kind + d4/6.0_r8_kind
822 orb_angle(n) = orb_angle(n-1) + d5
823 end do
824
825end subroutine orbit
826
827
828!> @brief Orbital time returns the time (1 year = 2*pi) since autumnal
829!! equinox
830!!
831!! @details Orbital time returns the time (1 year = 2*pi) since autumnal
832!! equinox; autumnal_eq_ref is a module variable of time_type and
833!! will have been defined by default or by a call to
834!! set_ref_date_of_ae; length_of_year is available through the time
835!! manager and is set at the value approriate for the calandar being used
836function orbital_time(time) result(t)
837
838type(time_type), intent(in) :: time !< time (1 year = 2*pi) since autumnal equinox
839real(kind=r8_kind) :: t
840
842 t = twopi*(t - real(floor(t), r8_kind))
843 if (time < autumnal_eq_ref) t = twopi - t
844
845end function orbital_time
846
847
848!> @brief universal_time returns the time of day at longitude = 0.0
849!! (1 day = 2*pi)
850function universal_time(time) result(t)
851
852type(time_type), intent(in) :: time !< Time (1 year = 2*pi) since autumnal equinox
853real(kind=r8_kind) :: t
854
855 !--------------------------------------------------------------------
856 ! local variables
857 !--------------------------------------------------------------------
858 integer :: seconds, days
859
860 call get_time(time, seconds, days)
861 t = twopi * real(seconds, r8_kind)/real(seconds_per_day, r8_kind)
862
863end function universal_time
864
865#include "astronomy_r4.fh"
866#include "astronomy_r8.fh"
867
868end module astronomy_mod
869!> @}
870! close documentation grouping
subroutine get_period_time_type(period_out)
get_period_time_type returns the length of the year as a time_type variable.
real(r8_kind) ecc
Eccentricity of Earth's orbit [dimensionless].
subroutine set_period_time_type(period_in)
Set_period_time_type saves the length of the year (input as a time_type variable) into a time_type mo...
integer second_ae
Second of specified autumnal equinox.
subroutine, public astronomy_end
astronomy_end is the destructor for astronomy_mod.
type(time_type) autumnal_eq_ref
time_type variable containing specified time of reference NH autumnal equinox
subroutine, private orbit
Orbit computes and stores a table of value of orbital angles as a function of orbital time (both the ...
subroutine, public astronomy_init(latb, lonb)
astronomy_init is the constructor for astronomy_mod.
integer minute_ae
Minute of specified autumnal equinox.
integer period
Specified length of year [seconds]; must be specified to override default value given by length_of_ye...
real(r8_kind), dimension(:), allocatable orb_angle
table of orbital positions (0 to 2*pi) as a function of time used to find actual orbital position via...
real(r8_kind), dimension(:,:), allocatable fracday_ann
annual mean daylight fraction
integer total_pts
number of grid boxes owned by the processor
real(kind=r8_kind) function, public universal_time(time)
universal_time returns the time of day at longitude = 0.0 (1 day = 2*pi)
subroutine, public set_ref_date_of_ae(day_in, month_in, year_in, second_in, minute_in, hour_in)
set_ref_date_of_ae provides a means of specifying the reference date of the NH autumnal equinox for a...
integer num_angles
Number of intervals into which the year is divided to compute orbital positions.
real(r8_kind), dimension(:,:), allocatable cosz_ann
annual mean cos of zenith angle
integer num_pts
count of grid_boxes for which annual mean astronomy values have been calculated
subroutine set_period_integer(period_in)
set_period_integer saves as the input length of the year (an integer) in a time_type module variable.
integer year_ae
Year of specified autumnal equinox.
real(r8_kind) obliq
Obliquity [degrees].
real(kind=r8_kind) function, public orbital_time(time)
Orbital time returns the time (1 year = 2*pi) since autumnal equinox.
integer month_ae
Month of specified autumnal equinox.
subroutine, public get_ref_date_of_ae(day_out, month_out, year_out, second_out, minute_out, hour_out)
get_ref_date_of_ae retrieves the reference date of the autumnal equinox as integer variables.
real(r8_kind) seconds_per_day
seconds in a day
real(r8_kind) rrsun_ann
annual mean earth-sun distance
real(r8_kind), dimension(:,:), allocatable solar_ann
annual mean solar factor
integer day_ae
Day of specified autumnal equinox.
subroutine get_period_integer(period_out)
get_period_integer returns the length of the year as an integer number of seconds.
real(r8_kind) deg_to_rad
conversion from degrees to radians
real(r8_kind) per
Longitude of perihelion with respect to autumnal equinox in NH [degrees].
integer hour_ae
Hour of specified autumnal equinox.
logical module_is_initialized
has the module been initialized ?
real(r8_kind) twopi
2 *PI
type(time_type) period_time_type
time_type variable containing period of one orbit
logical annual_mean_calculated
have the annual mean values been calculated?
Calculates the annual mean of solar information for a given latitude and time.
Calculates the daily mean solar information for a given time and latitude.
Calculates solar information for the given location(lat & lon) and time.
Gets the length of year for current calendar.
Private interface for internal use by dirunal_solar and daily_mean_solar.
Sets the length of a year for the calendar in use.
type(time_type) function, public length_of_year()
Returns the mean length of the year in the default calendar setting.
subroutine, public get_time(time, seconds, days, ticks, err_msg)
Returns days and seconds ( < 86400 ) corresponding to a time. err_msg should be checked for any error...
subroutine, public time_manager_init()
Initialization routine. Writes the version information to the log file.
Given an input date in year, month, days, etc., creates a time_type that represents this time interva...
Given some number of seconds and days, returns the corresponding time_type.
Type to represent amounts of time. Implemented as seconds and days to allow for larger intervals.