MONC
Functions/Subroutines
solar_position_angle_mod Module Reference

Functions/Subroutines

subroutine solar_pos_calculation (socrates_opt, socrates_derived_fields)
 
subroutine solar_angle_calculation (socrates_opt, socrates_derived_fields)
 

Function/Subroutine Documentation

◆ solar_angle_calculation()

subroutine solar_position_angle_mod::solar_angle_calculation ( type (str_socrates_options), intent(in)  socrates_opt,
type (str_socrates_derived_fields), intent(inout)  socrates_derived_fields 
)

Definition at line 69 of file solar_position_angle_mod.F90.

70 
71  type (str_socrates_options), intent(in) :: socrates_opt
72  type (str_socrates_derived_fields), intent(inout) :: socrates_derived_fields
73 
74  real(kind=default_precision), parameter :: &
75  degrees_to_radians = 0.017453292, & ! conversion factor
76  seconds_in_hour = 3600.0
77 
78  real(kind=default_precision) :: &
79  sinlat, & ! sin of latitude in radians
80  longit ! longitude in radians
81 
82  real(kind=default_precision) :: &
83  twopi, & ! 2*pi
84  s2r, & ! seconds-to-radians converter
85  sinsin, & ! products of the sines and of the cosines
86  coscos, & ! of solar declination and of latitude.
87  hld, & ! half-length of the day in radians (equal
88  ! to the hour-angle of sunset, and minus
89  coshld, & ! the hour-angle of sunrise) its cosine.
90  hat, & ! local hour angle at the start time.
91  omegab, & ! beginning and end of the timestep and
92  omegae, & ! of the period over which cosz is
93  omega1, & ! integrated, and sunset - all measured in
94  omega2, & ! radians after local sunrise, not from
95  omegas, & ! local noon as the true hour angle is.
96  difsin, & ! a difference-of-sines intermediate value
97  diftim, & ! and the corresponding time period
98  ! the start-time and length of the timestep (time_radians & dt_radians)
99  ! converted to radians after midday gmt, or equivalently, hour
100  ! angle of the sun on the greenwich meridian.
101  time_radians, &
102  dt_radians
103 
104  twopi = 2. * pi
105  s2r = pi / 43200.
106 
107  sinlat = sin(socrates_opt%latitude*degrees_to_radians)
108  longit = socrates_opt%longitude*degrees_to_radians
109 
110  time_radians = &
111  (socrates_opt%rad_time_hours*seconds_in_hour) * s2r - pi
112  dt_radians = socrates_derived_fields%dt_secs * s2r
113 
114  ! original code from LEM looped over k, which is number of points
115  ! in the branch for version 1.0 of MONC, number of points will
116  ! always be 1, so removed the loop
117  hld = 0. ! logically unnecessary
118  ! statement without which the cray compiler will not vectorize this code
119  sinsin = socrates_derived_fields%sindec * sinlat
120  coscos = sqrt( (1.0- socrates_derived_fields%sindec**2.0) * &
121  (1.0- sinlat**2.0) )
122  coshld = sinsin / coscos
123  if (coshld.lt.-1.) then ! perpetual night
124  socrates_derived_fields%fraction_lit = 0.
125  socrates_derived_fields%cosz = 0.
126  else
127  hat = longit + time_radians ! (3.2.2)
128  if (coshld.gt.1.) then ! perpetual day - hour
129  omega1 = hat ! angles for (3.2.3) are
130  omega2 = hat + dt_radians ! start & end of timestep
131  else ! at this latitude some
132  ! points are sunlit, some not. different ones need different treatment.
133  hld = acos(-coshld) ! (3.2.4)
134  ! the logic seems simplest if one takes all "times" - actually hour
135  ! angles - relative to sunrise (or sunset), but they must be kept in the
136  ! range 0 to 2pi for the tests on their orders to work.
137  omegab = hat + hld
138  if (omegab.lt.0.) omegab = omegab + twopi
139  if (omegab.ge.twopi) omegab = omegab - twopi
140  if (omegab.ge.twopi) omegab = omegab - twopi
141  ! line repeated - otherwise could have failure if
142  ! longitudes w are > pi rather than < 0.
143  omegae = omegab + dt_radians
144  if (omegae.gt.twopi) omegae = omegae - twopi
145  omegas = 2. * hld
146  ! now that the start-time, end-time and sunset are set in terms of hour
147  ! angle, can set the two hour-angles for (3.2.3). the simple cases are
148  ! start-to-end-of-timestep, start-to-sunset, sunrise-to-end and sunrise-
149  ! -to-sunset, but two other cases exist and need special treatment.
150  if (omegab.le.omegas .or. omegab.lt.omegae) then
151  omega1 = omegab - hld
152  else
153  omega1 = - hld
154  endif
155  if (omegae.le.omegas) then
156  omega2 = omegae - hld
157  else
158  omega2 = omegas - hld
159  endif
160  if (omegae.gt.omegab.and.omegab.gt.omegas) omega2=omega1
161  ! put in an arbitrary marker for the case when the sun does not rise
162  ! during the timestep (though it is up elsewhere at this latitude).
163  ! (cannot set cosz & lit within the else ( coshld < 1 ) block
164  ! because 3.2.3 is done outside this block.)
165  endif ! this finishes the else (perpetual day) block
166  difsin = sin(omega2) - sin(omega1) ! begin (3.2.3)
167  diftim = omega2 - omega1
168  ! next, deal with the case where the sun sets and then rises again
169  ! within the timestep. there the integration has actually been done
170  ! backwards over the night, and the resulting negative difsin and diftim
171  ! must be combined with positive values representing the whole of the
172  ! timestep to get the right answer, which combines contributions from
173  ! the two separate daylit periods. a simple analytic expression for the
174  ! total sun throughout the day is used. (this could of course be used
175  ! alone at points where the sun rises and then sets within the timestep)
176  if (diftim.lt.0.) then
177  difsin = difsin + 2. * sqrt(1.-coshld**2)
178  diftim = diftim + 2. * hld
179  endif
180  if (diftim.eq.0.) then
181  ! pick up the arbitrary marker for night points at a partly-lit latitude
182  socrates_derived_fields%cosz = 0.
183  socrates_derived_fields%fraction_lit = 0.
184  else
185  socrates_derived_fields%cosz = difsin*coscos/diftim + sinsin ! (3.2.3)
186  socrates_derived_fields%fraction_lit = diftim / dt_radians
187  endif
188  endif ! this finishes the else (perpetual night) block
189 

◆ solar_pos_calculation()

subroutine solar_position_angle_mod::solar_pos_calculation ( type (str_socrates_options), intent(in)  socrates_opt,
type (str_socrates_derived_fields), intent(inout)  socrates_derived_fields 
)

Definition at line 13 of file solar_position_angle_mod.F90.

14  ! -------------------------------------------------------------------
15  ! Calculations of the earth's orbit described in UMDP 23. Using
16  ! the day of the year and the orbital "constants" (which vary over
17  ! "Milankovitch" timescales) it calculates the sin of the solar
18  ! declination and the inverse-square scaling factor for the solar
19  ! "constant".
20  ! -----------------------------------------------------------------
21 
22  type (str_socrates_options), intent(in) :: socrates_opt
23  type (str_socrates_derived_fields), intent(inout) :: socrates_derived_fields
24 
25  real(kind=default_precision) :: &
26  gamma_rad, e, tau0, sinobl, e1, e2, e3, e4, diny_360, &
27  tau1_360, tau1_365
28  real(kind=default_precision) :: diny_365 ! number of days in year
29  real(kind=default_precision) :: m, v ! mean & true anomaly
30 
31  gamma_rad = 1.352631
32  e=.0167
33  tau0 = 2.5 ! true date of perihelion
34  sinobl = 0.397789 ! sin (obliquity)
35  e1 = e * (2.0-0.25*e*e) ! coefficients for 3.1.2
36  e2 = 1.25 * e*e ! coefficients for 3.1.2
37  e3 = e*e*e * 13.0/12.0
38  e4 = ((1.0+e*e*0.5)/(1.0-e*e))**2.0 ! constant for 3.1.4
39  diny_360 = 360. ! number of days in year
40  tau1_360 = tau0*diny_360/365.25+0.71+.5
41  tau1_365 = tau0+.5
42 
43  if (.not. socrates_opt%l_360) then
44  if (mod(socrates_opt%rad_year,4.0) .eq. 0 .and. & ! is this a leap year?
45  (mod(socrates_opt%rad_year,400.0) .eq. 0 .or. &
46  mod(socrates_opt%rad_year,100.0) .ne. 0)) then
47  diny_365 = 366.
48  else
49  diny_365 = 365.
50  end if
51  end if
52 
53  ! tau1 is modified so as to include the conversion of day-ordinal into
54  ! fractional-number-of-days-into-the-year-at-12-z-on-this-day.
55 
56  if (socrates_opt%l_360) then
57  m = (2*pi) * (socrates_opt%rad_day-tau1_360) / diny_360 ! eq 3.1.1
58  else
59  m = (2*pi) * (socrates_opt%rad_day-tau1_365) / diny_365 ! eq 3.1.1
60  end if
61  v = m + e1*sin(m) + e2*sin(2.*m) + e3*sin(3.*m) ! eq 3.1.2
62 
63  ! Fields required by monc-socrates coupling
64  socrates_derived_fields%scs = e4 * ( 1. + e * cos(v) ) **2 ! eq 3.1.4
65  socrates_derived_fields%sindec = sinobl * sin(v - gamma_rad)
66 
datadefn_mod::default_precision
integer, parameter, public default_precision
MPI communication type which we use for the prognostic and calculation data.
Definition: datadefn.F90:17
Generated by doxygen 1.8.20