!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 !-------------------------------------------------------------------- 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: astronomy.F90,v 17.0 2009/07/21 03:18:20 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !--------------------------------------------------------------------- !------- 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_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. !----------------------------------------------------------------------- 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 !--------------------------------------------------------------------- ! 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(jd)) allocate (solar_ann(jd)) allocate (fracday_ann(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_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(js:je) = cosz(1,:) solar_ann(js:je) = solar(1,:) fracday_ann(js:je) = fracday(1,:) 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 do i=1, size(lat,1) cosz(i,:) = cosz_ann(js:je) solar(i,:) = solar_ann(js:je) fracday(i,:) = fracday_ann(js:je) end do 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(jst:jnd) solar(:) = solar_ann(jst:jnd) fracday(:) = fracday_ann(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) & 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