module astronomy_mod ! ! fil ! ! ! ! ! ! astronomy_mod provides astronomical variables for use ! by other modules within fms. the only currently used interface is ! for determination of astronomical values needed by the shortwave ! radiation packages. ! ! ! ! shared modules: use fms_mod, only: open_namelist_file, fms_init, & mpp_pe, mpp_root_pe, stdlog, & file_exist, write_version_number, & check_nml_error, error_mesg, & FATAL, NOTE, WARNING, close_file use time_manager_mod, only: time_type, set_time, get_time, & get_date_julian, set_date_julian, & set_date, length_of_year, & time_manager_init, & operator(-), operator(+), & operator( // ), operator(<) use constants_mod, only: constants_init, PI use mpp_mod, only: input_nml_file !-------------------------------------------------------------------- implicit none private !--------------------------------------------------------------------- ! astronomy_mod provides astronomical variables for use ! by other modules within fms. the only currently used interface is ! for determination of astronomical values needed by the shortwave ! radiation packages. !--------------------------------------------------------------------- !--------------------------------------------------------------------- !----------- version number for this module -------------------------- character(len=128) :: version = '$Id$' character(len=128) :: tagname = '$Name$' !--------------------------------------------------------------------- !------- interfaces -------- public & astronomy_init, get_period, set_period, & set_orbital_parameters, get_orbital_parameters, & set_ref_date_of_ae, get_ref_date_of_ae, & diurnal_solar, daily_mean_solar, annual_mean_solar, & astronomy_end, universal_time, orbital_time interface diurnal_solar module procedure diurnal_solar_2d module procedure diurnal_solar_1d module procedure diurnal_solar_0d module procedure diurnal_solar_cal_2d module procedure diurnal_solar_cal_1d module procedure diurnal_solar_cal_0d end interface interface daily_mean_solar module procedure daily_mean_solar_2d module procedure daily_mean_solar_1d module procedure daily_mean_solar_2level module procedure daily_mean_solar_0d module procedure daily_mean_solar_cal_2d module procedure daily_mean_solar_cal_1d module procedure daily_mean_solar_cal_2level module procedure daily_mean_solar_cal_2level_radwt module procedure daily_mean_solar_cal_0d end interface interface annual_mean_solar module procedure annual_mean_solar_2d module procedure annual_mean_solar_1d module procedure annual_mean_solar_2level end interface interface get_period module procedure get_period_time_type, get_period_integer end interface interface set_period module procedure set_period_time_type, set_period_integer end interface private & ! called from astronomy_init and set_orbital_parameters: orbit, & ! called from diurnal_solar, daily_mean_solar and orbit: r_inv_squared, & ! called from diurnal_solar and daily_mean_solar: angle, declination, & half_day ! half_day, orbital_time, & ! called from diurnal_solar: ! universal_time interface half_day module procedure half_day_2d, half_day_0d end interface !--------------------------------------------------------------------- !-------- namelist --------- real :: ecc = 0.01671 ! eccentricity of orbital ellipse ! [ non-dimensional ] real :: obliq = 23.439 ! obliquity [ degrees ] real :: per = 102.932 ! longitude of perihelion with respect ! to autumnal equinox in NH [ degrees ] integer :: period = 0 ! specified length of year [ seconds ] ; ! must be specified to override default ! value given by length_of_year in ! time_manager_mod integer :: day_ae = 23 ! day of specified autumnal equinox integer :: month_ae = 9 ! month of specified autumnal equinox integer :: year_ae = 1998 ! year of specified autumnal equinox integer :: hour_ae = 5 ! hour of specified autumnal equinox integer :: minute_ae = 37 ! minute of specified autumnal equinox integer :: second_ae = 0 ! second of specified autumnal equinox integer :: num_angles = 3600 ! number of intervals into which the year ! is divided to compute orbital positions namelist /astronomy_nml/ ecc, obliq, per, period, & year_ae, month_ae, day_ae, & hour_ae, minute_ae, second_ae, & num_angles !-------------------------------------------------------------------- !------ public data ---------- !-------------------------------------------------------------------- !------ private data ---------- type(time_type) :: autumnal_eq_ref ! time_type variable containing ! specified time of reference ! NH autumnal equinox type(time_type) :: period_time_type ! time_type variable containing ! period of one orbit real, dimension(:), allocatable :: & orb_angle ! table of orbital positions (0 to ! 2*pi) as a function of time used ! to find actual orbital position ! via interpolation real :: seconds_per_day=86400. ! seconds in a day real :: deg_to_rad ! conversion from degrees to radians real :: twopi ! 2 *PI logical :: module_is_initialized= & .false. ! has the module been initialized ? real, dimension(:,:), allocatable :: & cosz_ann, & ! annual mean cos of zenith angle solar_ann, & ! annual mean solar factor fracday_ann ! annual mean daylight fraction real :: rrsun_ann ! annual mean earth-sun distance logical :: annual_mean_calculated = & .false. ! have the annual mean values been ! calculated ? integer :: num_pts = 0 ! count of grid_boxes for which ! annual mean astronomy values have ! been calculated integer :: total_pts ! number of grid boxes owned by the ! processor !-------------------------------------------------------------------- !-------------------------------------------------------------------- contains !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PUBLIC SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !#################################################################### ! ! ! astronomy_init is the constructor for astronomy_mod. ! ! ! astronomy_init is the constructor for astronomy_mod. ! ! ! ! 2d array of model latitudes at cell corners [radians] ! ! ! 2d array of model longitudes at cell corners [radians] ! ! ! subroutine astronomy_init (latb, lonb) !------------------------------------------------------------------- ! astronomy_init is the constructor for astronomy_mod. !------------------------------------------------------------------- real, dimension(:,:), intent(in), optional :: latb real, dimension(:,:), intent(in), optional :: lonb !-------------------------------------------------------------------- ! intent(in) variables: ! ! latb 2d array of model latitudes at cell corners ! [ radians ] ! lonb 2d array of model longitudes at cell corners ! [ radians ] ! !-------------------------------------------------------------------- !------------------------------------------------------------------- ! local variables: integer :: unit, ierr, io, seconds, & days, jd, id !------------------------------------------------------------------- ! local variables: ! ! unit ! ierr ! io ! seconds ! days ! jd ! id ! !--------------------------------------------------------------------- !------------------------------------------------------------------- ! if module has already been initialized, exit. !------------------------------------------------------------------- if (module_is_initialized) return !------------------------------------------------------------------- ! verify that modules used by this module have been initialized. !------------------------------------------------------------------- call fms_init call time_manager_init call constants_init !----------------------------------------------------------------------- ! read namelist. !----------------------------------------------------------------------- #ifdef INTERNAL_FILE_NML read (input_nml_file, astronomy_nml, iostat=io) ierr = check_nml_error(io,'astronomy_nml') #else if ( file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=astronomy_nml, iostat=io, end=10) ierr = check_nml_error(io,'astronomy_nml') end do 10 call close_file (unit) endif #endif !--------------------------------------------------------------------- ! write version number and namelist to logfile. !--------------------------------------------------------------------- call write_version_number (version, tagname) if (mpp_pe() == mpp_root_pe() ) then unit = stdlog() write (unit, nml=astronomy_nml) endif !-------------------------------------------------------------------- ! be sure input values are within valid ranges. ! QUESTION : ARE THESE THE RIGHT LIMITS ??? !--------------------------------------------------------------------- if (ecc < 0.0 .or. ecc > 0.99) & call error_mesg ('astronomy_mod', & 'ecc must be between 0 and 0.99', FATAL) if (obliq < -90. .or. obliq > 90.) & call error_mesg ('astronomy_mod', & 'obliquity must be between -90 and 90 degrees', FATAL) if (per < 0.0 .or. per > 360.0) & call error_mesg ('astronomy_mod', & 'perihelion must be between 0 and 360 degrees', FATAL) !---------------------------------------------------------------------- ! set up time-type variable defining specified time of autumnal ! equinox. !---------------------------------------------------------------------- autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, & hour_ae,minute_ae,second_ae) !--------------------------------------------------------------------- ! set up time-type variable defining length of year. !---------------------------------------------------------------------- if (period == 0) then period_time_type = length_of_year() call get_time (period_time_type, seconds, days) period = seconds_per_day*days + seconds else period_time_type = set_time(period,0) endif !--------------------------------------------------------------------- ! define useful module variables. !---------------------------------------------------------------------- twopi = 2.*PI deg_to_rad = twopi/360. !--------------------------------------------------------------------- ! call orbit to define table of orbital angles as function of ! orbital time. !---------------------------------------------------------------------- ! wfc moved here from orbit allocate ( orb_angle(0:num_angles) ) call orbit !-------------------------------------------------------------------- ! if annual mean radiation is desired, then latb will be present. ! allocate arrays to hold the needed astronomical factors. define ! the total number of points that the processor is responsible for. !-------------------------------------------------------------------- if (present(latb)) then jd = size(latb,2) - 1 id = size(lonb,1) - 1 allocate (cosz_ann(id, jd)) allocate (solar_ann(id, jd)) allocate (fracday_ann(id, jd)) total_pts = jd*id endif !--------------------------------------------------------------------- ! mark the module as initialized. !--------------------------------------------------------------------- module_is_initialized=.true. !--------------------------------------------------------------------- end subroutine astronomy_init !################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE GET_PERIOD ! ! ! call get_period (period) ! ! separate routines exist within this interface for integer ! and time_type output: ! ! integer, intent(out) :: period ! OR ! type(time_type), intent(out) :: period ! !-------------------------------------------------------------------- ! ! intent(out) variable: ! ! period_out length of year for calendar in use ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! ! get_period_integer returns the length of the year as an integer ! number of seconds. ! ! ! get_period_integer returns the length of the year as an integer ! number of seconds. ! ! ! ! number of seconds as the length of the year ! ! ! subroutine get_period_integer (period_out) !-------------------------------------------------------------------- ! get_period_integer returns the length of the year as an integer ! number of seconds. !-------------------------------------------------------------------- integer, intent(out) :: period_out !-------------------------------------------------------------------- ! local variables: integer :: seconds, days !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !-------------------------------------------------------------------- ! define length of year in seconds. !-------------------------------------------------------------------- call get_time (period_time_type, seconds, days) period_out = seconds_per_day*days + seconds end subroutine get_period_integer !#################################################################### ! ! ! get_period_time_type returns the length of the year as a time_type ! variable. ! ! ! get_period_time_type returns the length of the year as a time_type ! variable. ! ! ! ! the length of the year as a time_type ! ! ! subroutine get_period_time_type (period_out) !-------------------------------------------------------------------- ! get_period_time_type returns the length of the year as a time_type ! variable. !-------------------------------------------------------------------- type(time_type), intent(inout) :: period_out !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !-------------------------------------------------------------------- ! define length of year as a time_type variable. !-------------------------------------------------------------------- period_out = period_time_type end subroutine get_period_time_type !##################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE GET_PERIOD ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE SET_PERIOD ! ! ! call set_period (period_in) ! ! separate routines exist within this interface for integer ! and time_type output: ! ! integer, intent(out) :: period_in ! OR ! type(time_type), intent(out) :: period_in ! !-------------------------------------------------------------------- ! ! intent(in) variable: ! ! period_in length of year for calendar in use ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! ! set_period_integer saves as the input length of the year (an ! integer) in a time_type module variable. ! ! ! set_period_integer saves as the input length of the year (an ! integer) in a time_type module variable. ! ! ! ! the length of the year as a time_type ! ! ! 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, intent(in) :: period_in !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !--------------------------------------------------------------------- ! define time_type variable defining the length of year from the ! input value (integer seconds). !--------------------------------------------------------------------- period_time_type = set_time(period_in, 0) end subroutine set_period_integer !##################################################################### 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 module variable. !-------------------------------------------------------------------- type(time_type), intent(in) :: period_in !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !--------------------------------------------------------------------- ! define time_type variable defining the length of year from the ! input value (time_type). !--------------------------------------------------------------------- period_time_type = period_in end subroutine set_period_time_type !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE SET_PERIOD ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !##################################################################### ! ! ! set_orbital_parameters saves the input values of eccentricity, ! obliquity and perihelion time as module variables for use by ! astronomy_mod. ! ! ! set_orbital_parameters saves the input values of eccentricity, ! obliquity and perihelion time as module variables for use by ! astronomy_mod. ! ! ! ! eccentricity of orbital ellipse ! ! ! obliquity fof orbital ellipse ! ! ! longitude of perihelion with respect to autumnal ! equinox in northern hemisphere ! ! ! subroutine set_orbital_parameters (ecc_in, obliq_in, per_in) !-------------------------------------------------------------------- ! set_orbital_parameters saves the input values of eccentricity, ! obliquity and perihelion time as module variables for use by ! astronomy_mod. !-------------------------------------------------------------------- real, intent(in) :: ecc_in real, intent(in) :: obliq_in real, intent(in) :: per_in !-------------------------------------------------------------------- ! ! intent(in) variables: ! ! ecc_in eccentricity of orbital ellipse ! [ non-dimensional ] ! obliq_in obliquity ! [ degrees ] ! per_in longitude of perihelion with respect to autumnal ! equinox in northern hemisphere ! [ degrees ] ! !------------------------------------------------------------------- !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !-------------------------------------------------------------------- ! be sure input values are within valid ranges. ! QUESTION : ARE THESE THE RIGHT LIMITS ??? !--------------------------------------------------------------------- if (ecc_in < 0.0 .or. ecc_in > 0.99) & call error_mesg ('astronomy_mod', & 'ecc must be between 0 and 0.99', FATAL) if (obliq_in < -90.0 .or. obliq > 90.0) & call error_mesg ('astronomy_mod', & 'obliquity must be between -90. and 90. degrees', FATAL) if (per_in < 0.0 .or. per_in > 360.0) & call error_mesg ('astronomy_mod', & 'perihelion must be between 0.0 and 360. degrees', FATAL) !--------------------------------------------------------------------- ! save input values into module variables. !--------------------------------------------------------------------- ecc = ecc_in obliq = obliq_in per = per_in !--------------------------------------------------------------------- ! call orbit to define table of orbital angles as function of ! orbital time using the input values of parameters just supplied. !---------------------------------------------------------------------- call orbit !---------------------------------------------------------------------- end subroutine set_orbital_parameters !#################################################################### ! ! ! get_orbital_parameters retrieves the orbital parameters for use ! by another module. ! ! ! get_orbital_parameters retrieves the orbital parameters for use ! by another module. ! ! ! ! eccentricity of orbital ellipse ! ! ! obliquity fof orbital ellipse ! ! ! longitude of perihelion with respect to autumnal ! equinox in northern hemisphere ! ! ! subroutine get_orbital_parameters (ecc_out, obliq_out, per_out) !------------------------------------------------------------------- ! get_orbital_parameters retrieves the orbital parameters for use ! by another module. !-------------------------------------------------------------------- real, intent(out) :: ecc_out, obliq_out, per_out !-------------------------------------------------------------------- ! ! intent(out) variables: ! ! ecc_out eccentricity of orbital ellipse ! [ non-dimensional ] ! obliq_out obliquity ! [ degrees ] ! per_out longitude of perihelion with respect to autumnal ! equinox in northern hemisphere ! [ degrees ] ! !------------------------------------------------------------------- !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !-------------------------------------------------------------------- ! fill the output arguments with the eccentricity, obliquity and ! perihelion angle. !-------------------------------------------------------------------- ecc_out = ecc obliq_out = obliq per_out = per end subroutine get_orbital_parameters !#################################################################### ! ! ! set_ref_date_of_ae provides a means of specifying the reference ! date of the NH autumnal equinox for a particular year. ! ! ! set_ref_date_of_ae provides a means of specifying the reference ! date of the NH autumnal equinox for a particular year. it is only ! used if calls are made to the calandar versions of the routines ! diurnal_solar and daily_mean_solar. if the NOLEAP calendar is ! used, then the date of autumnal equinox will be the same every ! year. if JULIAN is used, then the date of autumnal equinox will ! return to the same value every 4th year. ! ! ! ! day of reference autumnal equinox ! ! ! month of reference autumnal equinox ! ! ! year of reference autumnal equinox ! ! ! OPTIONAL: second of reference autumnal equinox ! ! ! OPTIONAL: minute of reference autumnal equinox ! ! ! OPTIONAL: hour of reference autumnal equinox ! ! ! subroutine 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 particular year. it is only ! used if calls are made to the calandar versions of the routines ! diurnal_solar and daily_mean_solar. if the NOLEAP calendar is ! used, then the date of autumnal equinox will be the same every ! year. if JULIAN is used, then the date of autumnal equinox will ! return to the same value every 4th year. !---------------------------------------------------------------------- integer, intent(in) :: day_in, month_in, year_in integer, intent(in), optional :: second_in, minute_in, hour_in !-------------------------------------------------------------------- ! ! intent(in) variables: ! ! day_in day of reference autumnal equinox ! [ non-dimensional ] ! month_in month of reference autumnal equinox ! [ non-dimensional ] ! year_in year of reference autumnal equinox ! [ non-dimensional ] ! ! intent(in), optional variables: ! ! second_in second of reference autumnal equinox ! [ non-dimensional ] ! minute_in minute of reference autumnal equinox ! [ non-dimensional ] ! hour_in hour of reference autumnal equinox ! [ non-dimensional ] ! !------------------------------------------------------------------- !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !-------------------------------------------------------------------- ! save the input time of ae specification into a time_type module ! variable autumnal_eq_ref. !-------------------------------------------------------------------- day_ae = day_in month_ae = month_in year_ae = year_in if (present(second_in)) then second_ae = second_in minute_ae = minute_in hour_ae = hour_in else second_ae = 0 minute_ae = 0 hour_ae = 0 endif autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, & hour_ae,minute_ae,second_ae) !--------------------------------------------------------------------- end subroutine set_ref_date_of_ae !#################################################################### ! ! ! get_ref_date_of_ae retrieves the reference date of the autumnal ! equinox as integer variables. ! ! ! get_ref_date_of_ae retrieves the reference date of the autumnal ! equinox as integer variables. ! ! ! ! day of reference autumnal equinox ! ! ! month of reference autumnal equinox ! ! ! year of reference autumnal equinox ! ! ! second of reference autumnal equinox ! ! ! minute of reference autumnal equinox ! ! ! hour of reference autumnal equinox ! ! ! subroutine 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. !--------------------------------------------------------------------- integer, intent(out) :: day_out, month_out, year_out, & second_out, minute_out, hour_out !-------------------------------------------------------------------- ! ! intent(out) variables: ! ! day_out day of reference autumnal equinox ! [ non-dimensional ] ! month_out month of reference autumnal equinox ! [ non-dimensional ] ! year_out year of reference autumnal equinox ! [ non-dimensional ] ! second_out second of reference autumnal equinox ! [ non-dimensional ] ! minute_out minute of reference autumnal equinox ! [ non-dimensional ] ! hour_out hour of reference autumnal equinox ! [ non-dimensional ] ! !------------------------------------------------------------------- !--------------------------------------------------------------------- ! exit if module has not been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized) & call error_mesg ( 'astronomy_mod', & ' module has not been initialized', FATAL) !--------------------------------------------------------------------- ! fill the output fields with the proper module data. !--------------------------------------------------------------------- day_out = day_ae month_out = month_ae year_out = year_ae second_out = second_ae minute_out = minute_ae hour_out = hour_ae end subroutine get_ref_date_of_ae !##################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE DIURNAL_SOLAR ! ! call diurnal_solar (lat, lon, time, cosz, fracday, rrsun, dt_time) ! or ! call diurnal_solar (lat, lon, gmt, time_since_ae, cosz, fracday, ! rrsun, dt) ! ! the first option (used in conjunction with time_manager_mod) ! generates the real variables gmt and time_since_ae from the ! time_type input, and then calls diurnal_solar with these real ! inputs. ! ! the time of day is set by ! real, intent(in) :: gmt ! the time of year is set by ! real, intent(in) :: time_since_ae ! with time_type input, both of these are extracted from ! type(time_type), intent(in) :: time ! ! ! separate routines exist within this interface for scalar, ! 1D or 2D input and output fields: ! ! real, intent(in), dimension(:,:) :: lat, lon ! OR real, intent(in), dimension(:) :: lat, lon ! OR real, intent(in) :: lat, lon ! ! real, intent(out), dimension(:,:) :: cosz, fracday ! OR real, intent(out), dimension(:) :: cosz, fracday ! OR real, intent(out) :: cosz, fracday ! ! one may also average the output fields over the time interval ! between gmt and gmt + dt by including the optional argument dt (or ! dt_time). dt is measured in radians and must be less than pi ! (1/2 day). this average is computed analytically, and should be ! exact except for the fact that changes in earth-sun distance over ! the time interval dt are ignored. in the context of a diurnal GCM, ! this option should always be employed to insure that the total flux ! at the top of the atmosphere is not modified by time truncation ! error. ! ! real, intent(in), optional :: dt ! type(time_type), optional :: dt_time ! (see test.90 for examples of the use of these types) ! !-------------------------------------------------------------------- ! ! intent(in) variables: ! ! lat latitudes of model grid points ! [ radians ] ! lon longitudes of model grid points ! [ radians ] ! gmt time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! [ radians ] ! time_since_ae time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! [ radians ] ! time time at which astronomical values are desired ! time_type variable [ seconds, days] ! ! ! intent(out) variables: ! ! cosz cosine of zenith angle, set to zero when entire ! period is in darkness ! [ dimensionless ] ! fracday daylight fraction of time interval ! [ dimensionless ] ! rrsun earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! [ dimensionless ] ! ! intent(in), optional variables: ! ! dt time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! [ radians ], (1 day = 2 * pi) ! dt_time as in dt, but dt_time is a time_type variable ! time_type [ days, seconds ] ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! ! diurnal_solar_2d returns 2d fields of cosine of zenith angle, ! daylight fraction and earth-sun distance at the specified lati- ! tudes, longitudes and time. these values may be instantaneous ! or averaged over a specified time interval. ! ! ! diurnal_solar_2d returns 2d fields of cosine of zenith angle, ! daylight fraction and earth-sun distance at the specified lati- ! tudes, longitudes and time. these values may be instantaneous ! or averaged over a specified time interval. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt, allow_negative_cosz, & half_day_out) !--------------------------------------------------------------------- ! diurnal_solar_2d returns 2d fields of cosine of zenith angle, ! daylight fraction and earth-sun distance at the specified lati- ! tudes, longitudes and time. these values may be instantaneous ! or averaged over a specified time interval. !--------------------------------------------------------------------- real, dimension(:,:), intent(in) :: lat, lon real, intent(in) :: gmt, time_since_ae real, dimension(:,:), intent(out) :: cosz, fracday real, intent(out) :: rrsun real, intent(in), optional :: dt logical, intent(in), optional :: allow_negative_cosz real, dimension(:,:), intent(out), optional :: half_day_out !--------------------------------------------------------------------- ! local variables real, dimension(size(lat,1),size(lat,2)) :: t, tt, h, aa, bb, & st, stt, sh real :: ang, dec logical :: Lallow_negative !--------------------------------------------------------------------- ! local variables ! ! t time of day with respect to local noon (2 pi = 1 day) ! [ radians ] ! tt end of averaging period [ radians ] ! h half of the daily period of daylight, centered at noon ! [ radians, -pi --> pi ] ! aa sin(lat) * sin(declination) ! bb cos(lat) * cos(declination) ! st sine of time of day ! stt sine of time of day at end of averaging period ! sh sine of half-day period ! ang position of earth in its orbit wrt autumnal equinox ! [ radians ] ! dec earth's declination [ radians ] ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! be sure the time in the annual cycle is legitimate. !--------------------------------------------------------------------- if (time_since_ae < 0.0 .or. time_since_ae > twopi) & call error_mesg('astronomy_mod', & 'time_since_ae not between 0 and 2pi', FATAL) !-------------------------------------------------------------------- ! be sure the time at longitude = 0.0 is legitimate. !--------------------------------------------------------------------- if (gmt < 0.0 .or. gmt > twopi) & call error_mesg('astronomy_mod', & 'gmt not between 0 and 2pi', FATAL) !--------------------------------------------------------------------- ! define the orbital angle (location in year), solar declination and ! earth sun distance factor. use functions contained in this module. !--------------------------------------------------------------------- ang = angle(time_since_ae) dec = declination(ang) rrsun = r_inv_squared(ang) !--------------------------------------------------------------------- ! define terms needed in the cosine zenith angle equation. !-------------------------------------------------------------------- aa = sin(lat)*sin(dec) bb = cos(lat)*cos(dec) !--------------------------------------------------------------------- ! define local time. force it to be between -pi and pi. !-------------------------------------------------------------------- t = gmt + lon - PI where(t >= PI) t = t - twopi where(t < -PI) t = t + twopi Lallow_negative = .false. if (present(allow_negative_cosz)) then if (allow_negative_cosz) Lallow_negative = .true. end if !--------------------------------------------------------------------- ! perform a time integration to obtain cosz and fracday if desired. ! output is valid over the period from t to t + dt. !-------------------------------------------------------------------- h = half_day (lat,dec) if ( present(half_day_out) ) then half_day_out = h end if if ( present(dt) ) then ! (perform time averaging) tt = t + dt st = sin(t) stt = sin(tt) sh = sin(h) cosz = 0.0 if (.not. Lallow_negative) then !------------------------------------------------------------------- ! case 1: entire averaging period is before sunrise. !------------------------------------------------------------------- where (t < -h .and. tt < -h) cosz = 0.0 !------------------------------------------------------------------- ! case 2: averaging period begins before sunrise, ends after sunrise ! but before sunset !------------------------------------------------------------------- where ( (tt+h) /= 0.0 .and. t < -h .and. abs(tt) <= h) & cosz = aa + bb*(stt + sh)/ (tt + h) !------------------------------------------------------------------- ! case 3: averaging period begins before sunrise, ends after sunset, ! but before the next sunrise. modify if averaging period extends ! past the next day's sunrise, but if averaging period is less than ! a half- day (pi) that circumstance will never occur. !------------------------------------------------------------------- where (t < -h .and. h /= 0.0 .and. h < tt) & cosz = aa + bb*( sh + sh)/(h+h) !------------------------------------------------------------------- ! case 4: averaging period begins after sunrise, ends before sunset. !------------------------------------------------------------------- where ( abs(t) <= h .and. abs(tt) <= h) & cosz = aa + bb*(stt - st)/ (tt - t) !------------------------------------------------------------------- ! case 5: averaging period begins after sunrise, ends after sunset. ! modify when averaging period extends past the next day's sunrise. !------------------------------------------------------------------- where ((h-t) /= 0.0 .and. abs(t) <= h .and. h < tt) & cosz = aa + bb*(sh - st)/(h-t) !------------------------------------------------------------------- ! case 6: averaging period begins after sunrise , ends after the ! next day's sunrise. note that this includes the case when the ! day length is one day (h = pi). !------------------------------------------------------------------- where (twopi - h < tt .and. (tt+h-twopi) /= 0.0 .and. t <= h ) & cosz = (cosz*(h - t) + (aa*(tt + h - twopi) + & bb*(stt + sh))) / ((h - t) + (tt + h - twopi)) !------------------------------------------------------------------- ! case 7: averaging period begins after sunset and ends before the ! next day's sunrise !------------------------------------------------------------------- where( h < t .and. twopi - h >= tt ) cosz = 0.0 !------------------------------------------------------------------- ! case 8: averaging period begins after sunset and ends after the ! next day's sunrise but before the next day's sunset. if the ! averaging period is less than a half-day (pi) the latter ! circumstance will never occur. !----------------------------------------------------------------- where( h < t .and. twopi - h < tt ) cosz = aa + bb*(stt + sh) / (tt + h - twopi) end where else cosz = aa + bb*(stt - st)/ (tt - t) end if !------------------------------------------------------------------- ! day fraction is the fraction of the averaging period contained ! within the (-h,h) period. !------------------------------------------------------------------- where (t < -h .and. tt < -h) fracday = 0.0 where (t < -h .and. abs(tt) <= h) fracday = (tt + h )/dt where (t < -h .and. h < tt) fracday = ( h + h )/dt where (abs(t) <= h .and. abs(tt) <= h) fracday = (tt - t )/dt where (abs(t) <= h .and. h < tt) fracday = ( h - t )/dt where ( h < t ) fracday = 0.0 where (twopi - h < tt) fracday = fracday + & (tt + h - & twopi)/dt !---------------------------------------------------------------------- ! if instantaneous values are desired, define cosz at time t. !---------------------------------------------------------------------- else ! (no time averaging) if (.not. Lallow_negative) then where (abs(t) < h) cosz = aa + bb*cos(t) fracday = 1.0 elsewhere cosz = 0.0 fracday = 0.0 end where else cosz = aa + bb*cos(t) where (abs(t) < h) fracday = 1.0 elsewhere fracday = 0.0 end where end if end if !---------------------------------------------------------------------- ! be sure that cosz is not negative. !---------------------------------------------------------------------- if (.not. Lallow_negative) then cosz = max(0.0, cosz) end if !-------------------------------------------------------------------- end subroutine diurnal_solar_2d !###################################################################### ! ! ! diurnal_solar_1d takes 1-d input fields, adds a second dimension ! and calls diurnal_solar_2d. on return, the 2d fields are returned ! to the original 1d fields. ! ! ! diurnal_solar_1d takes 1-d input fields, adds a second dimension ! and calls diurnal_solar_2d. on return, the 2d fields are returned ! to the original 1d fields. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_1d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt, allow_negative_cosz, & half_day_out) !-------------------------------------------------------------------- ! diurnal_solar_1d takes 1-d input fields, adds a second dimension ! and calls diurnal_solar_2d. on return, the 2d fields are returned ! to the original 1d fields. !---------------------------------------------------------------------- !--------------------------------------------------------------------- real, dimension(:), intent(in) :: lat, lon real, intent(in) :: gmt, time_since_ae real, dimension(:), intent(out) :: cosz, fracday real, intent(out) :: rrsun real, intent(in), optional :: dt logical, intent(in), optional :: allow_negative_cosz real, dimension(:), intent(out), optional :: half_day_out !--------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, & fracday_2d,halfday_2d !-------------------------------------------------------------------- ! define 2-d versions of input data arrays. !-------------------------------------------------------------------- lat_2d(:,1) = lat lon_2d(:,1) = lon !-------------------------------------------------------------------- ! call diurnal_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- ! if (present(dt)) then call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae,& cosz_2d, fracday_2d, rrsun, dt=dt, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) ! else ! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, & ! cosz_2d, fracday_2d, rrsun) ! endif !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- fracday = fracday_2d(:,1) cosz = cosz_2d (:,1) if (present(half_day_out)) then half_day_out = halfday_2d(:,1) end if end subroutine diurnal_solar_1d !##################################################################### ! ! ! diurnal_solar_0d takes scalar input fields, makes them into 2d ! arrays dimensioned (1,1), and calls diurnal_solar_2d. on return, ! the 2d fields are converted back to the desired scalar output. ! ! ! diurnal_solar_0d takes scalar input fields, makes them into 2d ! arrays dimensioned (1,1), and calls diurnal_solar_2d. on return, ! the 2d fields are converted back to the desired scalar output. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_0d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt, allow_negative_cosz, & half_day_out) !-------------------------------------------------------------------- ! diurnal_solar_0d takes scalar input fields, makes them into 2d ! arrays dimensioned (1,1), and calls diurnal_solar_2d. on return, ! the 2d fields are converted back to the desired scalar output. !---------------------------------------------------------------------- real, intent(in) :: lat, lon, gmt, time_since_ae real, intent(out) :: cosz, fracday, rrsun real, intent(in), optional :: dt logical,intent(in), optional :: allow_negative_cosz real, intent(out), optional :: half_day_out !-------------------------------------------------------------------- ! local variables: ! real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d !--------------------------------------------------------------------- ! create 2d arrays from the scalar input fields. !--------------------------------------------------------------------- lat_2d = lat lon_2d = lon !-------------------------------------------------------------------- ! call diurnal_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- ! if (present(dt)) then call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, & cosz_2d, fracday_2d, rrsun, dt=dt, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) ! else ! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, & ! cosz_2d, fracday_2d, rrsun) ! end if !------------------------------------------------------------------- ! place output fields into scalars for return to calling routine. !------------------------------------------------------------------- fracday = fracday_2d(1,1) cosz = cosz_2d(1,1) if (present(half_day_out)) then half_day_out = halfday_2d(1,1) end if end subroutine diurnal_solar_0d !#################################################################### ! ! ! diurnal_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! diurnal_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_cal_2d (lat, lon, time, cosz, fracday, & rrsun, dt_time, allow_negative_cosz, & half_day_out) !------------------------------------------------------------------- ! diurnal_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. !------------------------------------------------------------------- !------------------------------------------------------------------- real, dimension(:,:), intent(in) :: lat, lon type(time_type), intent(in) :: time real, dimension(:,:), intent(out) :: cosz, fracday real, intent(out) :: rrsun type(time_type), intent(in), optional :: dt_time logical, intent(in), optional :: allow_negative_cosz real, dimension(:,:), intent(out), optional :: half_day_out !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real :: dt real :: gmt, time_since_ae !--------------------------------------------------------------------- ! extract time of day (gmt) from time_type variable time with ! function universal_time. !--------------------------------------------------------------------- gmt = universal_time(time) !--------------------------------------------------------------------- ! extract the time of year (time_since_ae) from time_type variable ! time using the function orbital_time. !--------------------------------------------------------------------- time_since_ae = orbital_time(time) !--------------------------------------------------------------------- ! convert optional time_type variable dt_time (length of averaging ! period) to a real variable dt with the function universal_time. !--------------------------------------------------------------------- if (present(dt_time)) then dt = universal_time(dt_time) if (dt > PI) then call error_mesg ( 'astronomy_mod', & 'radiation time step must be no longer than 12 hrs', & FATAL) endif if (dt == 0.0) then call error_mesg ( 'astronomy_mod', & 'radiation time step must not be an integral & &number of days', FATAL) endif !-------------------------------------------------------------------- ! call diurnal_solar_2d to calculate astronomy fields, with or ! without the optional argument dt. !-------------------------------------------------------------------- call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt=dt, & allow_negative_cosz=allow_negative_cosz, & half_day_out=half_day_out) else call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, & allow_negative_cosz=allow_negative_cosz, & half_day_out=half_day_out) end if end subroutine diurnal_solar_cal_2d !##################################################################### ! ! ! diurnal_solar_cal_1d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! diurnal_solar_cal_1d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_cal_1d (lat, lon, time, cosz, fracday, & rrsun, dt_time, allow_negative_cosz, & half_day_out) !-------------------------------------------------------------------- real, dimension(:), intent(in) :: lat, lon type(time_type), intent(in) :: time real, dimension(:), intent(out) :: cosz, fracday real, intent(out) :: rrsun type(time_type), intent(in), optional :: dt_time logical, intent(in), optional :: allow_negative_cosz real, dimension(:), intent(out), optional :: half_day_out !-------------------------------------------------------------------- !------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, & fracday_2d, halfday_2d !-------------------------------------------------------------------- ! define 2-d versions of input data arrays. !-------------------------------------------------------------------- lat_2d(:,1) = lat lon_2d(:,1) = lon !-------------------------------------------------------------------- ! call diurnal_solar_cal_2d to convert the time_types to reals and ! then calculate the astronomy fields. !-------------------------------------------------------------------- if (present(dt_time)) then call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, & fracday_2d, rrsun, dt_time=dt_time, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) else call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, & fracday_2d, rrsun, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) end if !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- fracday = fracday_2d(:,1) cosz = cosz_2d (:,1) if (present(half_day_out)) then half_day_out = halfday_2d(:,1) end if end subroutine diurnal_solar_cal_1d !##################################################################### ! ! ! diurnal_solar_cal_0d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! diurnal_solar_cal_0d receives time_type inputs, converts ! them to real variables and then calls diurnal_solar_2d to ! compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! longitude of model grid points ! ! ! time of day at longitude 0.0; midnight = 0.0, ! one day = 2 * pi ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! OPTIONAL: time interval after gmt over which the astronomical ! variables are to be averaged. this produces averaged ! output rather than instantaneous. ! ! ! subroutine diurnal_solar_cal_0d (lat, lon, time, cosz, fracday, & rrsun, dt_time, allow_negative_cosz, & half_day_out) !--------------------------------------------------------------------- real, intent(in) :: lat, lon type(time_type), intent(in) :: time real, intent(out) :: cosz, fracday, rrsun type(time_type), intent(in), optional :: dt_time logical, intent(in), optional :: allow_negative_cosz real, intent(out), optional :: half_day_out !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d !-------------------------------------------------------------------- ! define 2-d versions of input data arrays. !-------------------------------------------------------------------- lat_2d = lat lon_2d = lon !-------------------------------------------------------------------- ! call diurnal_solar_cal_2d to convert the time_types to reals and ! then calculate the astronomy fields. !-------------------------------------------------------------------- if (present(dt_time)) then call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, & fracday_2d, rrsun, dt_time=dt_time, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) else call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, & fracday_2d, rrsun, & allow_negative_cosz=allow_negative_cosz, & half_day_out=halfday_2d) end if !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- fracday= fracday_2d(1,1) cosz = cosz_2d(1,1) if (present(half_day_out)) then half_day_out = halfday_2d(1,1) end if end subroutine diurnal_solar_cal_0d !#################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE DIURNAL_SOLAR ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE DAILY_MEAN_SOLAR ! ! ! call daily_mean_solar (lat, time, cosz, fracday, rrsun) ! or ! call daily_mean_solar (lat, time_since_ae, cosz, fracday, rrsun) ! or ! call daily_mean_solar (lat, time, cosz, solar) ! or ! call daily_mean_solar (lat, time_since_ae, cosz, solar) ! ! the first option (used in conjunction with time_manager_mod) ! generates the real variable time_since_ae from the time_type ! input time, and then calls daily_mean_solar with this real input ! (option 2). the third and fourth options correspond to the first ! and second and are used with then spectral 2-layer model, where ! only cosz and solar are desired as output. these routines generate ! dummy arguments and then call option 2, where the calculation is ! done. ! ! the time of year is set by ! real, intent(in) :: time_since_ae ! with time_type input, the time of year is extracted from ! type(time_type), intent(in) :: time ! ! ! separate routines exist within this interface for scalar, ! 1D or 2D input and output fields: ! ! real, intent(in), dimension(:,:) :: lat ! OR real, intent(in), dimension(:) :: lat ! OR real, intent(in) :: lat ! ! real, intent(out), dimension(:,:) :: cosz, fracday ! OR real, intent(out), dimension(:) :: cosz, fracday ! OR real, intent(out) :: cosz, fracday ! !-------------------------------------------------------------------- ! ! intent(in) variables: ! ! lat latitudes of model grid points ! [ radians ] ! time_since_ae time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! [ radians ] ! time time at which astronomical values are desired ! time_type variable [ seconds, days] ! ! ! intent(out) variables: ! ! cosz cosz is cosine of an effective zenith angle that ! would produce the correct daily solar flux if the ! sun were fixed at that position for the period of ! daylight. ! should one also, or instead, compute cosz weighted ! by the instantaneous flux, averaged over the day ?? ! [ dimensionless ] ! fracday daylight fraction of time interval ! [ dimensionless ] ! rrsun earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! [ dimensionless ] ! solar shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! [ dimensionless ] ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! ! daily_mean_solar_2d computes the daily mean astronomical ! parameters for the input points at latitude lat and time of year ! time_since_ae. ! ! ! daily_mean_solar_2d computes the daily mean astronomical ! parameters for the input points at latitude lat and time of year ! time_since_ae. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! 2-d array of half-day lengths at the latitudes ! ! ! the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. ! ! ! subroutine daily_mean_solar_2d (lat, time_since_ae, cosz, h_out, rr_out) !---------------------------------------------------------------------- ! daily_mean_solar_2d computes the daily mean astronomical ! parameters for the input points at latitude lat and time of year ! time_since_ae. ! !---------------------------------------------------------------------- !---------------------------------------------------------------------- real, dimension(:,:), intent(in) :: lat real, intent(in) :: time_since_ae real, dimension(:,:), intent(out) :: cosz, h_out real, intent(out) :: rr_out !---------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables real, dimension(size(lat,1),size(lat,2)) :: h real :: ang, dec, rr !-------------------------------------------------------------------- ! be sure the time in the annual cycle is legitimate. !--------------------------------------------------------------------- if (time_since_ae < 0.0 .or. time_since_ae > twopi) & call error_mesg('astronomy_mod', & 'time_since_ae not between 0 and 2pi', FATAL) !--------------------------------------------------------------------- ! define the orbital angle (location in year), solar declination, ! half-day length and earth sun distance factor. use functions ! contained in this module. !--------------------------------------------------------------------- ang = angle (time_since_ae) dec = declination(ang) h = half_day (lat, dec) rr = r_inv_squared (ang) !--------------------------------------------------------------------- ! where the entire day is dark, define cosz to be zero. otherwise ! use the standard formula. define the daylight fraction and earth- ! sun distance. !--------------------------------------------------------------------- where (h == 0.0) cosz = 0.0 elsewhere cosz = sin(lat)*sin(dec) + cos(lat)*cos(dec)*sin(h)/h end where h_out = h/PI rr_out = rr !-------------------------------------------------------------------- end subroutine daily_mean_solar_2d !##################################################################### ! ! ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. ! ! ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! 2-d array of half-day lengths at the latitudes ! ! ! the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. ! ! ! subroutine daily_mean_solar_1d (lat, time_since_ae, cosz, h_out, rr_out) !-------------------------------------------------------------------- ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. !---------------------------------------------------------------------- !---------------------------------------------------------------------- real, intent(in), dimension(:) :: lat real, intent(in) :: time_since_ae real, intent(out), dimension(size(lat(:))) :: cosz real, intent(out), dimension(size(lat(:))) :: h_out real, intent(out) :: rr_out !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat !-------------------------------------------------------------------- ! call daily_mean_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, & hout_2d, rr_out) !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- h_out = hout_2d(:,1) cosz = cosz_2d(:,1) end subroutine daily_mean_solar_1d !###################################################################### ! ! ! daily_mean_solar_2level takes 1-d input fields, adds a second ! dimension and calls daily_mean_solar_2d. on return, the 2d fields ! are returned to the original 1d fields. ! ! ! daily_mean_solar_2level takes 1-d input fields, adds a second ! dimension and calls daily_mean_solar_2d. on return, the 2d fields ! are returned to the original 1d fields. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! subroutine daily_mean_solar_2level (lat, time_since_ae, cosz, solar) !-------------------------------------------------------------------- ! daily_mean_solar_2level takes 1-d input fields, adds a second ! dimension and calls daily_mean_solar_2d. on return, the 2d fields ! are returned to the original 1d fields. !---------------------------------------------------------------------- !---------------------------------------------------------------------- real, intent(in), dimension(:) :: lat real, intent(in) :: time_since_ae real, intent(out), dimension(size(lat(:))) :: cosz, solar !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d real :: rr_out !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat !-------------------------------------------------------------------- ! call daily_mean_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, & hout_2d, rr_out) !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- solar = cosz_2d(:,1)*hout_2d(:,1)*rr_out cosz = cosz_2d(:,1) end subroutine daily_mean_solar_2level !#################################################################### ! ! ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. ! ! ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! 2-d array of half-day lengths at the latitudes ! ! ! the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. ! ! ! subroutine daily_mean_solar_0d (lat, time_since_ae, cosz, h_out, rr_out) !-------------------------------------------------------------------- ! daily_mean_solar_1d takes 1-d input fields, adds a second dimension ! and calls daily_mean_solar_2d. on return, the 2d fields are ! returned to the original 1d fields. !---------------------------------------------------------------------- real, intent(in) :: lat, time_since_ae real, intent(out) :: cosz, h_out, rr_out !-------------------------------------------------------------------- ! local variables real, dimension(1,1) :: lat_2d, cosz_2d, hout_2d !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d = lat !-------------------------------------------------------------------- ! call daily_mean_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, & hout_2d, rr_out) !------------------------------------------------------------------- ! return output fields to scalars for return to calling routine. !------------------------------------------------------------------- h_out = hout_2d(1,1) cosz = cosz_2d(1,1) end subroutine daily_mean_solar_0d !#################################################################### ! ! ! daily_mean_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. ! ! ! daily_mean_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! subroutine daily_mean_solar_cal_2d (lat, time, cosz, fracday, rrsun) !------------------------------------------------------------------- ! daily_mean_solar_cal_2d receives time_type inputs, converts ! them to real variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. !------------------------------------------------------------------- !------------------------------------------------------------------- real, dimension(:,:), intent(in) :: lat type(time_type), intent(in) :: time real, dimension(:,:), intent(out) :: cosz, fracday real, intent(out) :: rrsun !------------------------------------------------------------------- !------------------------------------------------------------------- ! local variables real :: time_since_ae !-------------------------------------------------------------------- ! be sure the time in the annual cycle is legitimate. !--------------------------------------------------------------------- time_since_ae = orbital_time(time) if (time_since_ae < 0.0 .or. time_since_ae > twopi) & call error_mesg ('astronomy_mod', & 'time_since_ae not between 0 and 2pi', FATAL) !-------------------------------------------------------------------- ! call daily_mean_solar_2d to calculate astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_2d (lat, time_since_ae, cosz, & fracday, rrsun) end subroutine daily_mean_solar_cal_2d !##################################################################### ! ! ! daily_mean_solar_cal_1d receives time_type inputs, converts ! them to real, 2d variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. ! ! ! daily_mean_solar_cal_1d receives time_type inputs, converts ! them to real, 2d variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! subroutine daily_mean_solar_cal_1d (lat, time, cosz, fracday, rrsun) !------------------------------------------------------------------- ! daily_mean_solar_cal_1d receives time_type inputs, converts ! them to real, 2d variables and then calls daily_mean_solar_2d to ! compute desired astronomical variables. !------------------------------------------------------------------- real, dimension(:), intent(in) :: lat type(time_type), intent(in) :: time real, dimension(:), intent(out) :: cosz, fracday real, intent(out) :: rrsun !--------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat !-------------------------------------------------------------------- ! call daily_mean_solar_cal_2d to convert the time_types to reals and ! then calculate the astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, & fracday_2d, rrsun) !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- fracday = fracday_2d(:,1) cosz = cosz_2d(:,1) end subroutine daily_mean_solar_cal_1d !################################################################### ! ! ! daily_mean_solar_cal_2level receives 1d arrays and time_type input, ! converts them to real, 2d variables and then calls ! daily_mean_solar_2d to compute desired astronomical variables. ! ! ! daily_mean_solar_cal_2level receives 1d arrays and time_type input, ! converts them to real, 2d variables and then calls ! daily_mean_solar_2d to compute desired astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! subroutine daily_mean_solar_cal_2level (lat, time, cosz, solar) !------------------------------------------------------------------- ! daily_mean_solar_cal_2level receives 1d arrays and time_type input, ! converts them to real, 2d variables and then calls ! daily_mean_solar_2d to compute desired astronomical variables. !------------------------------------------------------------------- real, dimension(:), intent(in) :: lat type(time_type), intent(in) :: time real, dimension(:), intent(out) :: cosz, solar !--------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d real :: rrsun !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat !-------------------------------------------------------------------- ! call daily_mean_solar_cal_2d to convert the time_types to reals and ! then calculate the astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, & fracday_2d, rrsun) !------------------------------------------------------------------- ! place output fields into 1-d arguments for return to ! calling routine. !------------------------------------------------------------------- solar = cosz_2d(:,1)*fracday_2d(:,1)*rrsun cosz = cosz_2d(:,1) end subroutine daily_mean_solar_cal_2level !################################################################### ! ! ! daily_mean_solar_cal_2level_radwt receives 1d arrays and time_type input. ! It uses code copied from daily_mean_solar_2d to compute desired astronomical variables. ! In addition to daily_mean solar, also radiation-weighted cosz ! is computed. ! ! ! daily_mean_solar_cal_2level receives 1d arrays and time_type input. ! It uses code copied from daily_mean_solar_2d to compute desired astronomical variables. ! In addition to daily_mean solar, also radiation-weighted cosz ! is computed. [Georg Feulner: Radiation-weighted daily average solar zenith angle, 27.11.13] ! Blowing-up to 2-D variables is avoided here. ! ! ! ! latitudes of model grid points ! ! ! time of year ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! radiation-weighted cosine of solar zenith angle ! ! ! subroutine daily_mean_solar_cal_2level_radwt (lat, time, cosz, solar, cosz_radwt) real, dimension(:), intent(in) :: lat type(time_type), intent(in) :: time real, dimension(:), intent(out) :: cosz, solar, cosz_radwt !--------------------------------------------------------------------- ! local variables real :: rrsun ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! [ dimensionless ] real, dimension(size(lat,1)) :: lat2 ! lat with slightly shifted poles to avoid singularity in tangens real, dimension(size(lat,1)) :: tan_lat, sin_lat, cos_lat real, dimension(size(lat,1)) :: cos_half_day real, dimension(size(lat,1)) :: HA_SS ! hour angle at sunset; same as half_day elsewhere in this module real, dimension(size(lat,1)) :: sin_HA_SS real, dimension(size(lat,1)) :: fracday real, dimension(size(lat,1)) :: factor_A real :: time_since_ae ! time of year; autumnal equinox = 0.0, ! one year = 2 * pi ! [ radians ] real :: ang ! orbital position relative to NH autumnal equinox [ radians ] real :: dec ! solar declination angle in earth's orbit real :: tan_dec, sin_dec, cos_dec real :: eps = 1.0E-05 integer :: i ! code from daily_mean_solar_cal_2d: convert the time_types to reals and ! then calculate the astronomy fields. time_since_ae = orbital_time(time) ! code from daily_mean_solar_2d: calculate the astronomy fields. ang = angle (time_since_ae) dec = declination(ang) rrsun = r_inv_squared(ang) ! code adapted from half_day_2d tan_dec = tan(dec) ! adjust latitude so that its tangent will be defined. ! TODO: compare to range 0.5*PI +/- eps instead of rounded 0.5*PI where (lat == 0.5*PI) lat2 = lat - eps elsewhere where (lat == -0.5*PI) lat2 = lat + eps elsewhere lat2 = lat end where end where tan_lat = tan(lat2) ! define the cosine of the half-day length. adjust for cases of ! all daylight or all night. cos_half_day = -tan_lat * tan_dec where (cos_half_day <= -1.0) HA_SS = PI elsewhere where (cos_half_day >= +1.0) HA_SS = 0.0 elsewhere HA_SS = acos(cos_half_day) end where end where !where(cos_half_day > -1.0 .and. cos_half_day < 1.0) HA_SS = acos(cos_half_day) ! code from daily_mean_solar_2d ! where the entire day is dark, define cosz to be zero. otherwise ! use the standard formula. define the daylight fraction and earth- ! sun distance. !write(*,*) 'HA_SS', HA_SS !write(*,*) 'lat', lat !write(*,*) 'dec', dec ! new code ! sin_lat(:)=sin(lat(:)) !write(*,*) 'sin_lat', sin_lat cos_lat(:)=cos(lat(:)) !write(*,*) 'cos_lat', cos_lat sin_dec=sin(dec) !write(*,*) 'sin_dec', sin_dec cos_dec=cos(dec) !write(*,*) 'cos_dec', cos_dec sin_HA_SS(:) = sin(HA_SS(:)) !write(*,*) 'sin_HA_SS', sin_HA_SS !write(*,*) 'sin(lat)*sin(dec)', sin_lat*sin_dec !write(*,*) 'cos(lat)*cos(dec)', cos_lat*cos_dec !write(*,*) 'daily_mean_solar_cal_2level_radwt: size lat ',size(lat), 'bounds', lbound(lat), ubound(lat) !write(*,*) ' size cosz ', size(cosz), 'bounds', lbound(cosz), ubound(cosz) !write(*,*) ' size HA_SS ', size(HA_SS), 'bounds', lbound(HA_SS), ubound(HA_SS) ! code from daily_mean_solar_2d where (HA_SS .eq. 0) ! (HA_SS .lt. eps .and. HA_SS .gt. -eps) cosz = 0.0 elsewhere ! cosz = sin(lat)*sin(dec) + cos(lat)*cos(dec)*sin(HA_SS)/HA_SS cosz = sin_lat*sin_dec + cos_lat*cos_dec*sin_HA_SS/HA_SS end where ! Sigh deeply. ifort 17.0.4 produces buggy code for the above where-construct ! when compiling with options -O3 -xCORE-AVX2 -qno-opt-dynamic-align , ! also with -fp-model precise . It works correctly with -fp-model consistent , ! or without vectorization. ! Symptoms: Error: floating invalid is raised in the elsewhere branch. ! The same problem occurs when the where-construct is replaced by an explicit do-loop. ! Also, when the write-statement inside the do-loop is un-commented, ! it apparently prevents the buggy optimization. !do i=1,size(cosz) ! if (HA_SS(i) .eq. 0) then ! cosz(i) = 0.0 ! else ! !write(*,*) i !, HA_SS(i), sin_HA_SS(i) ! !flush(*) ! cosz(i) = sin_lat(i)*sin_dec + cos_lat(i)*cos_dec*sin_HA_SS(i)/HA_SS(i) ! end if !end do !write(*,*) 'finished cosz', cosz fracday = HA_SS/PI ! code from daily_mean_solar_cal_2level solar = cosz*fracday*rrsun ! new code ! factor_A = 2 * sin_lat**2 * sin_dec**2 * HA_SS & + 4 * sin_lat * sin_dec * cos_lat * cos_dec * sin_HA_SS & + cos_lat**2 * cos_dec**2 * (HA_SS + sin(2*HA_SS)/2) ! from G.F. 27.11.13 where (cosz > 0) ! Note: <\cos\theta_\mbox{s}> in the equations of G.F. is equal to astronomy_mod HA_SS*cosz/PI cosz_radwt = factor_A / (2*cosz*HA_SS) ! from G.F. 28.11.13 elsewhere cosz_radwt = 0 end where end subroutine daily_mean_solar_cal_2level_radwt !################################################################### ! ! ! daily_mean_solar_cal_0d converts scalar input fields to real, ! 2d variables and then calls daily_mean_solar_2d to compute desired ! astronomical variables. ! ! ! daily_mean_solar_cal_0d converts scalar input fields to real, ! 2d variables and then calls daily_mean_solar_2d to compute desired ! astronomical variables. ! ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! subroutine daily_mean_solar_cal_0d (lat, time, cosz, fracday, rrsun) !------------------------------------------------------------------- ! daily_mean_solar_cal_0d converts scalar input fields to real, ! 2d variables and then calls daily_mean_solar_2d to compute desired ! astronomical variables. !------------------------------------------------------------------- !-------------------------------------------------------------------- real, intent(in) :: lat type(time_type), intent(in) :: time real, intent(out) :: cosz, fracday, rrsun !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables real, dimension(1,1) :: lat_2d, cosz_2d, fracday_2d !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d = lat !-------------------------------------------------------------------- ! call daily_mean_solar_cal_2d to convert the time_types to reals and ! then calculate the astronomy fields. !-------------------------------------------------------------------- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, & fracday_2d, rrsun) !------------------------------------------------------------------- ! place output fields into scalar arguments for return to ! calling routine. !------------------------------------------------------------------- fracday = fracday_2d(1,1) cosz = cosz_2d(1,1) end subroutine daily_mean_solar_cal_0d !#################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE DAILY_MEAN_SOLAR ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE ANNUAL_MEAN_SOLAR ! ! call annual_mean_solar (js, je, lat, cosz, solar, fracday, rrsun) ! or ! call annual_mean_solar (lat, cosz, solar) ! ! the second interface above is used by the spectral 2-layer model, ! which requires only cosz and solar as output arguments, and which ! makes this call during the initialization phase of the model. ! separate routines exist within this interface for 1D or 2D input ! and output fields: ! ! real, intent(in), dimension(:,:) :: lat ! OR real, intent(in), dimension(:) :: lat ! ! real, intent(out), dimension(:,:) :: cosz, solar, fracday ! OR real, intent(out), dimension(:) :: cosz, solar, fracday ! !--------------------------------------------------------------------- ! intent(in) variables: ! ! js, je starting/ ending subdomain j indices of data in ! the physics wiondow being integrated ! lat latitudes of model grid points ! [ radians ] ! ! ! intent(out) variables: ! ! cosz cosz is the average over the year of the cosine of ! an effective zenith angle that would produce the ! correct daily solar flux if the sun were fixed at ! that single position for the period of daylight on ! the given day. in this average, the daily mean ! effective cosz is weighted by the daily mean solar ! flux. ! [ dimensionless ] ! solar normalized solar flux, averaged over the year, ! equal to the product of fracday*cosz*rrsun ! [ dimensionless ] ! fracday daylight fraction calculated so as to make the ! average flux (solar) equal to the product of the ! flux-weighted avg cosz * this fracday * assumed ! annual mean avg earth-sun distance of 1.0. ! [ dimensionless ] ! rrsun annual mean earth-sun distance (r) relative to ! semi-major axis of orbital ellipse (a); assumed ! to be 1.0 ! [ dimensionless ] ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------------------------------------------------- ! ! ! annual_mean_solar_2d returns 2d fields of annual mean values of ! the cosine of zenith angle, daylight fraction and earth-sun ! distance at the specified latitude. ! ! ! annual_mean_solar_2d returns 2d fields of annual mean values of ! the cosine of zenith angle, daylight fraction and earth-sun ! distance at the specified latitude. ! ! ! ! Starting/ending index of latitude window ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! subroutine annual_mean_solar_2d (js, je, lat, cosz, solar, fracday, & rrsun) !--------------------------------------------------------------------- ! annual_mean_solar_2d returns 2d fields of annual mean values of ! the cosine of zenith angle, daylight fraction and earth-sun ! distance at the specified latitude. !--------------------------------------------------------------------- !-------------------------------------------------------------------- integer, intent(in) :: js, je real, dimension(:,:), intent(in) :: lat real, dimension(:,:), intent(out) :: solar, cosz, fracday real, intent(out) :: rrsun !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables real, dimension(size(lat,1),size(lat,2)) :: s,z real :: t integer :: n, i !-------------------------------------------------------------------- ! if the calculation has not yet been done, do it here. !-------------------------------------------------------------------- if (.not. annual_mean_calculated) then !---------------------------------------------------------------------- ! determine annual mean values of solar flux and product of cosz ! and solar flux by integrating the annual cycle in num_angles ! orbital increments. !---------------------------------------------------------------------- solar = 0.0 cosz = 0.0 do n =1, num_angles t = float(n-1)*twopi/float(num_angles) call daily_mean_solar(lat,t, z, fracday, rrsun) s = z*rrsun*fracday solar = solar + s cosz = cosz + z*s end do solar = solar/float(num_angles) cosz = cosz/float(num_angles) !-------------------------------------------------------------------- ! define the flux-weighted annual mean cosine of the zenith angle. !-------------------------------------------------------------------- where(solar.eq.0.0) cosz = 0.0 elsewhere cosz = cosz/solar end where !------------------------------------------------------------------- ! define avg fracday such as to make the avg flux (solar) equal to ! the product of the avg cosz * avg fracday * assumed mean avg ! radius of 1.0. it is unlikely that these avg fracday and avg rr ! will ever be used. !-------------------------------------------------------------------- where(solar .eq.0.0) fracday = 0.0 elsewhere fracday = solar/cosz end where rrsun = 1.00 !--------------------------------------------------------------------- ! save the values that have been calculated as module variables, if ! those variables are present; i.e., not the spectral 2-layer model. !--------------------------------------------------------------------- if (allocated (cosz_ann)) then cosz_ann = cosz solar_ann = solar fracday_ann = fracday rrsun_ann = rrsun !-------------------------------------------------------------------- ! increment the points computed counter. set flag to end execution ! once values have been calculated for all points owned by the ! processor. !-------------------------------------------------------------------- num_pts = num_pts + size(lat,1)*size(lat,2) if ( num_pts == total_pts) annual_mean_calculated = .true. endif !-------------------------------------------------------------------- ! if the calculation has been done, return the appropriate module ! variables. !-------------------------------------------------------------------- else if (allocated (cosz_ann)) then cosz = cosz_ann solar = solar_ann fracday = fracday_ann rrsun = rrsun_ann endif endif !---------------------------------------------------------------------- end subroutine annual_mean_solar_2d !##################################################################### ! ! ! annual_mean_solar_1d creates 2-d input fields from 1-d input fields ! and then calls annual_mean_solar_2d to obtain 2-d output fields ! which are then stored into 1-d fields for return to the calling ! subroutine. ! ! ! annual_mean_solar_1d creates 2-d input fields from 1-d input fields ! and then calls annual_mean_solar_2d to obtain 2-d output fields ! which are then stored into 1-d fields for return to the calling ! subroutine. ! ! ! ! Starting/ending index of latitude window ! ! ! latitudes of model grid points ! ! ! time of year; autumnal equinox = 0.0, one year = 2 * pi ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! daylight fraction of time interval ! ! ! earth-sun distance (r) relative to semi-major axis ! of orbital ellipse (a) : (a/r)**2 ! ! ! subroutine annual_mean_solar_1d (jst, jnd, lat, cosz, solar, & fracday, rrsun_out) !--------------------------------------------------------------------- ! annual_mean_solar_1d creates 2-d input fields from 1-d input fields ! and then calls annual_mean_solar_2d to obtain 2-d output fields ! which are then stored into 1-d fields for return to the calling ! subroutine. !--------------------------------------------------------------------- !--------------------------------------------------------------------- integer, intent(in) :: jst, jnd real, dimension(:), intent(in) :: lat(:) real, dimension(:), intent(out) :: cosz, solar, fracday real, intent(out) :: rrsun_out !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, solar_2d, cosz_2d, & fracday_2d real :: rrsun !-------------------------------------------------------------------- ! if the calculation has not been done, do it here. !-------------------------------------------------------------------- if ( .not. annual_mean_calculated) then !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat !-------------------------------------------------------------------- ! call annual_mean_solar_2d to calculate the astronomy fields. !-------------------------------------------------------------------- call annual_mean_solar_2d (jst, jnd, lat_2d, cosz_2d, & solar_2d, fracday_2d, rrsun) !------------------------------------------------------------------- ! place output fields into 1-D arrays for return to calling routine. !------------------------------------------------------------------- fracday = fracday_2d(:,1) rrsun_out = rrsun solar = solar_2d(:,1) cosz = cosz_2d(:,1) !-------------------------------------------------------------------- ! if the calculation has been done, simply return the module ! variables contain the results at the desired latitudes. !-------------------------------------------------------------------- else cosz(:) = cosz_ann(1,jst:jnd) solar(:) = solar_ann(1,jst:jnd) fracday(:) = fracday_ann(1,jst:jnd) rrsun = rrsun_ann endif end subroutine annual_mean_solar_1d !#################################################################### ! ! ! annual_mean_solar_2level creates 2-d input fields from 1-d input ! fields and then calls annual_mean_solar_2d to obtain 2-d output ! fields which are then stored into 1-d fields for return to the ! calling subroutine. this subroutine will be called during model ! initialization. ! ! ! annual_mean_solar_2level creates 2-d input fields from 1-d input ! fields and then calls annual_mean_solar_2d to obtain 2-d output ! fields which are then stored into 1-d fields for return to the ! calling subroutine. this subroutine will be called during model ! initialization. ! ! ! ! latitudes of model grid points ! ! ! cosine of solar zenith angle ! ! ! shortwave flux factor: cosine of zenith angle * ! daylight fraction / (earth-sun distance squared) ! ! ! subroutine annual_mean_solar_2level (lat, cosz, solar) !--------------------------------------------------------------------- ! annual_mean_solar_2level creates 2-d input fields from 1-d input ! fields and then calls annual_mean_solar_2d to obtain 2-d output ! fields which are then stored into 1-d fields for return to the ! calling subroutine. this subroutine will be called during model ! initialization. !--------------------------------------------------------------------- !--------------------------------------------------------------------- real, dimension(:), intent(in) :: lat real, dimension(:), intent(out) :: cosz, solar !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real, dimension(size(lat),1) :: lat_2d, solar_2d, cosz_2d, & fracday_2d integer :: jst, jnd real :: rrsun !-------------------------------------------------------------------- ! if the calculation has not been done, do it here. !-------------------------------------------------------------------- if ( .not. annual_mean_calculated) then !-------------------------------------------------------------------- ! define 2-d versions of input data array. !-------------------------------------------------------------------- lat_2d(:,1) = lat jst = 1 jnd = size(lat(:)) !-------------------------------------------------------------------- ! call annual_mean_solar_2d to calculate the astronomy fields. !-------------------------------------------------------------------- call annual_mean_solar_2d (jst, jnd, lat_2d, cosz_2d, & solar_2d, fracday_2d, rrsun) !------------------------------------------------------------------- ! place output fields into 1-D arrays for return to calling routine. !------------------------------------------------------------------- solar = solar_2d(:,1) cosz = cosz_2d(:,1) !-------------------------------------------------------------------- ! if the calculation has been done, print an error message since ! this subroutine should be called only once. !-------------------------------------------------------------------- else call error_mesg ('astronomy_mod', & 'annual_mean_solar_2level should be called only once', & FATAL) endif annual_mean_calculated = .true. end subroutine annual_mean_solar_2level !#################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE ANNUAL_MEAN_SOLAR ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !################################################################### ! ! ! astronomy_init is the destructor for astronomy_mod. ! ! ! astronomy_init is the destructor for astronomy_mod. ! ! ! ! subroutine astronomy_end !---------------------------------------------------------------------- ! astronomy_end is the destructor for astronomy_mod. !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! check if the module has been initialized. !---------------------------------------------------------------------- if (.not. module_is_initialized) return ! call error_mesg ( 'astronomy_mod', & ! ' module has not been initialized', FATAL) !---------------------------------------------------------------------- ! deallocate module variables. !---------------------------------------------------------------------- deallocate (orb_angle) if (allocated(cosz_ann) ) then deallocate (cosz_ann) deallocate (fracday_ann) deallocate (solar_ann) endif !---------------------------------------------------------------------- ! mark the module as uninitialized. !---------------------------------------------------------------------- module_is_initialized = .false. !--------------------------------------------------------------------- end subroutine astronomy_end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PRIVATE SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !#################################################################### ! ! ! orbit computes and stores a table of value of orbital angles as a ! function of orbital time (both the angle and time are zero at ! autumnal equinox in the NH, and range from 0 to 2*pi). ! ! ! orbit computes and stores a table of value of orbital angles as a ! function of orbital time (both the angle and time are zero at ! autumnal equinox in the NH, and range from 0 to 2*pi). ! ! ! ! subroutine orbit !--------------------------------------------------------------------- ! orbit computes and stores a table of value of orbital angles as a ! function of orbital time (both the angle and time are zero at ! autumnal equinox in the NH, and range from 0 to 2*pi). !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables integer :: n real :: d1, d2, d3, d4, d5, dt, norm !-------------------------------------------------------------------- ! allocate the orbital angle array, sized by the namelist parameter ! num_angles, defining the annual cycle resolution of the earth's ! orbit. define some constants to be used. !-------------------------------------------------------------------- ! wfc moving to astronomy_init ! allocate ( orb_angle(0:num_angles) ) orb_angle(0) = 0.0 dt = twopi/float(num_angles) norm = sqrt(1.0 - ecc**2) dt = dt*norm !--------------------------------------------------------------------- ! define the orbital angle at each of the num_angles locations in ! the orbit. !--------------------------------------------------------------------- do n = 1, num_angles d1 = dt*r_inv_squared(orb_angle(n-1)) d2 = dt*r_inv_squared(orb_angle(n-1)+0.5*d1) d3 = dt*r_inv_squared(orb_angle(n-1)+0.5*d2) d4 = dt*r_inv_squared(orb_angle(n-1)+d3) d5 = d1/6.0 + d2/3.0 +d3/3.0 +d4/6.0 orb_angle(n) = orb_angle(n-1) + d5 end do !------------------------------------------------------------------- end subroutine orbit !################################################################### ! ! ! r_inv_squared returns the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. ! ! ! r_inv_squared returns the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. ! ! ! ! angular position of earth in its orbit, relative to a ! value of 0.0 at the NH autumnal equinox, value between ! 0.0 and 2 * pi ! ! ! function r_inv_squared (ang) !-------------------------------------------------------------------- ! r_inv_squared returns the inverse of the square of the earth-sun ! distance relative to the mean distance at angle ang in the earth's ! orbit. !-------------------------------------------------------------------- !-------------------------------------------------------------------- real, intent(in) :: ang !-------------------------------------------------------------------- !--------------------------------------------------------------------- ! ! intent(in) variables: ! ! ang angular position of earth in its orbit, relative to a ! value of 0.0 at the NH autumnal equinox, value between ! 0.0 and 2 * pi ! [ radians ] ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real :: r_inv_squared, r, rad_per !--------------------------------------------------------------------- ! local variables: ! ! r_inv_squared the inverse of the square of the earth-sun ! distance relative to the mean distance ! [ dimensionless ] ! r earth-sun distance relative to mean distance ! [ dimensionless ] ! rad_per angular position of perihelion ! [ radians ] ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! define the earth-sun distance (r) and then return the inverse of ! its square (r_inv_squared) to the calling routine. !-------------------------------------------------------------------- rad_per = per*deg_to_rad r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) r_inv_squared = r**(-2) end function r_inv_squared !#################################################################### ! ! ! angle determines the position within the earth's orbit at time t ! in the year ( t = 0 at NH autumnal equinox.) by interpolating ! into the orbital position table. ! ! ! angle determines the position within the earth's orbit at time t ! in the year ( t = 0 at NH autumnal equinox.) by interpolating ! into the orbital position table. ! ! ! ! time of year (between 0 and 2*pi; t=0 at NH autumnal ! equinox ! ! ! function angle (t) !-------------------------------------------------------------------- ! angle determines the position within the earth's orbit at time t ! in the year ( t = 0 at NH autumnal equinox.) by interpolating ! into the orbital position table. !-------------------------------------------------------------------- !-------------------------------------------------------------------- real, intent(in) :: t !-------------------------------------------------------------------- !------------------------------------------------------------------- ! ! intent(in) variables: ! ! t time of year (between 0 and 2*pi; t=0 at NH autumnal ! equinox ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables real :: angle, norm_time, x integer :: int, int_1 !-------------------------------------------------------------------- ! local variables: ! ! angle orbital position relative to NH autumnal equinox ! [ radians ] ! norm_time index into orbital table corresponding to input time ! [ dimensionless ] ! x fractional distance between the orbital table entries ! bracketing the input time ! [ dimensionless ] ! int table index which is lower than actual position, but ! closest to it ! [ dimensionless ] ! int_1 next table index just larger than actual orbital ! position ! [ dimensionless ] ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! define orbital tables indices bracketing current orbital time ! (int and int_1). define table index distance between the lower ! table value (int) and the actual orbital time (x). define orbital ! position as being x of the way between int and int_1. renormalize ! angle to be within the range 0 to 2*pi. !-------------------------------------------------------------------- norm_time = t*float(num_angles)/twopi int = floor(norm_time) int = modulo(int,num_angles) int_1 = int+1 x = norm_time - floor(norm_time) angle = (1.0 -x)*orb_angle(int) + x*orb_angle(int_1) angle = modulo(angle, twopi) end function angle !#################################################################### ! ! ! declination returns the solar declination angle at orbital ! position ang in earth's orbit. ! ! ! declination returns the solar declination angle at orbital ! position ang in earth's orbit. ! ! ! ! solar orbital position ang in earth's orbit ! ! ! function declination (ang) !-------------------------------------------------------------------- ! declination returns the solar declination angle at orbital ! position ang in earth's orbit. !-------------------------------------------------------------------- !-------------------------------------------------------------------- real, intent(in) :: ang !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables real :: declination real :: rad_obliq, sin_dec !-------------------------------------------------------------------- ! local variables: ! ! declination solar declination angle ! [ radians ] ! rad_obliq obliquity of the ecliptic ! [ radians ] ! sin_dec sine of the solar declination ! [ dimensionless ] ! !-------------------------------------------------------------------- !--------------------------------------------------------------------- ! compute the solar declination. !--------------------------------------------------------------------- rad_obliq = obliq*deg_to_rad sin_dec = - sin(rad_obliq)*sin(ang) declination = asin(sin_dec) end function declination !##################################################################### !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE HALF_DAY ! ! half_day (latitude, dec) result (h) ! ! ! separate routines exist within this interface for scalar, ! or 2D input and output fields: ! ! real, intent(in), dimension(:,:) :: latitude ! OR real, intent(in) :: latitude ! ! real, dimension(size(latitude,1),size(latitude,2)) :: h ! OR real :: h ! !-------------------------------------------------------------------- ! ! intent(in) variables: ! ! latitude latitudes of model grid points ! [ radians ] ! dec solar declination ! [ radians ] ! ! intent(out) variables: ! ! h half of the length of daylight at the given ! latitude and orbital position (dec); value ! ranges between 0 (all darkness) and pi (all ! daylight) ! [ dimensionless ] ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! ! half_day_2d returns a 2-d array of half-day lengths at the ! latitudes and declination provided. ! ! ! half_day_2d returns a 2-d array of half-day lengths at the ! latitudes and declination provided. ! ! ! ! latitutde of view point ! ! ! solar declination angle at view point ! ! ! function half_day_2d (latitude, dec) result(h) !-------------------------------------------------------------------- ! half_day_2d returns a 2-d array of half-day lengths at the ! latitudes and declination provided. !-------------------------------------------------------------------- !--------------------------------------------------------------------- real, dimension(:,:), intent(in) :: latitude real, intent(in) :: dec real, dimension(size(latitude,1),size(latitude,2)) :: h !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables real, dimension (size(latitude,1),size(latitude,2)):: & cos_half_day, lat real :: tan_dec real :: eps = 1.0E-05 !--------------------------------------------------------------------- ! local variables ! ! cos_half_day cosine of half-day length ! [ dimensionless ] ! lat model latitude, adjusted so that it is never ! 0.5*pi or -0.5*pi ! tan_dec tangent of solar declination ! [ dimensionless ] ! eps small increment ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! define tangent of the declination. !-------------------------------------------------------------------- tan_dec = tan(dec) !-------------------------------------------------------------------- ! adjust latitude so that its tangent will be defined. !-------------------------------------------------------------------- lat = latitude where (latitude == 0.5*PI) lat= latitude - eps where (latitude == -0.5*PI) lat= latitude + eps !-------------------------------------------------------------------- ! define the cosine of the half-day length. adjust for cases of ! all daylight or all night. !-------------------------------------------------------------------- cos_half_day = -tan(lat)*tan_dec where (cos_half_day <= -1.0) h = PI where (cos_half_day >= +1.0) h = 0.0 where(cos_half_day > -1.0 .and. cos_half_day < 1.0) & h = acos(cos_half_day) end function half_day_2d !--------------------------------------------------------------- ! ! ! half_day_0d takes scalar input fields, makes them into 2-d fields ! dimensioned (1,1), and calls half_day_2d. on return, the 2-d ! fields are converted to the desired scalar output. ! ! ! half_day_0d takes scalar input fields, makes them into 2-d fields ! dimensioned (1,1), and calls half_day_2d. on return, the 2-d ! fields are converted to the desired scalar output. ! ! ! ! latitutde of view point ! ! ! solar declination angle at view point ! ! ! function half_day_0d(latitude, dec) result(h) !-------------------------------------------------------------------- ! half_day_0d takes scalar input fields, makes them into 2-d fields ! dimensioned (1,1), and calls half_day_2d. on return, the 2-d ! fields are converted to the desired scalar output. !-------------------------------------------------------------------- real, intent(in) :: latitude, dec real :: h !---------------------------------------------------------------------- ! local variables real, dimension(1,1) :: lat_2d, h_2d !--------------------------------------------------------------------- ! create 2d array from the input latitude field. !--------------------------------------------------------------------- lat_2d = latitude !--------------------------------------------------------------------- ! call half_day with the 2d arguments to calculate half-day length. !--------------------------------------------------------------------- h_2d = half_day (lat_2d, dec) !--------------------------------------------------------------------- ! create scalar from 2d array. !--------------------------------------------------------------------- h = h_2d(1,1) end function half_day_0d !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE HALF_DAY ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !#################################################################### ! ! ! orbital time returns the time (1 year = 2*pi) since autumnal ! equinox ! ! ! orbital time returns the time (1 year = 2*pi) since autumnal ! equinox; autumnal_eq_ref is a module variable of time_type and ! will have been defined by default or by a call to ! set_ref_date_of_ae; length_of_year is available through the time ! manager and is set at the value approriate for the calandar being ! used ! ! ! ! time (1 year = 2*pi) since autumnal equinox ! ! ! function orbital_time(time) result(t) !--------------------------------------------------------------------- ! orbital time returns the time (1 year = 2*pi) since autumnal ! equinox; autumnal_eq_ref is a module variable of time_type and ! will have been defined by default or by a call to ! set_ref_date_of_ae; length_of_year is available through the time ! manager and is set at the value approriate for the calandar being ! used !--------------------------------------------------------------------- type(time_type), intent(in) :: time real :: t !-------------------------------------------------------------------- t = real ( (time - autumnal_eq_ref)//period_time_type) t = twopi*(t - floor(t)) if (time < autumnal_eq_ref) t = twopi - t end function orbital_time !##################################################################### ! ! ! universal_time returns the time of day at longitude = 0.0 ! (1 day = 2*pi) ! ! ! universal_time returns the time of day at longitude = 0.0 ! (1 day = 2*pi) ! ! ! ! time (1 year = 2*pi) since autumnal equinox ! ! ! function universal_time(time) result(t) !-------------------------------------------------------------------- ! universal_time returns the time of day at longitude = 0.0 ! (1 day = 2*pi) !-------------------------------------------------------------------- type(time_type), intent(in) :: time real :: t !-------------------------------------------------------------------- ! local variables integer :: seconds, days call get_time (time, seconds, days) t = twopi*real(seconds)/seconds_per_day end function universal_time !##################################################################### end module astronomy_mod