cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 31-May-2007 JPDM c Routine to get satellite/sun angles if McIDAS doesn't work c Converted from Australian MTSAT code to F77 syntax real function get_satellite_za_navpix( longitude, latitude, & satellite_long, satellite_lat ) IMPLICIT NONE real longitude real latitude real satellite_long real satellite_lat real*8 semi_major_axis parameter( semi_major_axis = 6378.137d0 ) real*8 f_minus parameter( f_minus = 289.257222101d0 ) real*8 e_squared parameter( e_squared = (2./f_minus) - (1./(f_minus*f_minus)) ) real*8 r parameter( r = 42200.d0 ) real*8 h parameter( h = 0.d0 ) ! geodetic height of point of interest ! since we are doing SST this is the sea ! height which we'll approximate by 0. real*8 M_PI parameter( M_PI = 3.14159265358979323846d0 ) real*8 DEG_TO_RAD parameter( DEG_TO_RAD = M_PI / 180.d0 ) real*8 lambda real*8 phi real*8 lambda_s real*8 phi_s real*8 W real*8 value real*8 N real*8 xp, yp, zp real*8 xs, ys, zs real*8 input(3) real*8 output(3) real*8 nu c if( latitude .lt. 10 )then c print *,longitude,latitude,satellite_long,satellite_lat c endif lambda = longitude*DEG_TO_RAD phi = latitude*DEG_TO_RAD lambda_s = satellite_long*DEG_TO_RAD phi_s = satellite_lat*DEG_TO_RAD value = sin(phi) W = sqrt(1 - e_squared * value * value) ! sqrt( 1 - e^2 * sin^2(lat)) N = semi_major_axis/W value = (N+h) xp = value*cos(lambda)*cos(phi) yp = value*sin(lambda)*cos(phi) zp = (N*(1-e_squared)+h)*sin(phi) xs = r*cos(lambda_s) ys = r*sin(lambda_s) zs = 0.d0 input(1) = xs - xp input(2) = ys - yp input(3) = zs - zp call rotate_matrix_navpix( lambda, phi, input, output ) c! alpha = atan2(output(1),output(2)) nu = atan2(output(3), & sqrt(output(1)*output(1)+output(2)*output(2))) nu = nu / DEG_TO_RAD get_satellite_za_navpix = 90. - nu return end function subroutine rotate_matrix_navpix( lambda, phi, input, output ) IMPLICIT NONE real*8 lambda real*8 phi real*8 input(3) real*8 output(3) integer i, j real*8 R(3,3) C if( 3 .ne. size(input,DIM=1) )then C write(*,'(/,''ERROR: input array not 1x3'')') C stop C endif C if( 3 .ne. size(output,DIM=1) )then C write(*,'(/,''ERROR: output array not 1x3'')') C stop C endif c ! Rotation matrix R(1,1) = -sin(lambda) R(2,1) = cos(lambda) R(3,1) = 0. R(1,2) = -sin(phi)*cos(lambda) R(2,2) = -sin(phi)*sin(lambda) R(3,2) = cos(phi) R(1,3) = cos(phi)*cos(lambda) R(2,3) = cos(phi)*sin(lambda) R(3,3) = sin(phI) do i=1,3 output(i) = 0. do j=1,3 output(i) = output(i) + R(j,i)*input(j) end do end do end subroutine c---------------------------------------------------------------------------------------- real function get_sun_angle_navpix( longitude, latitude, year, & month, day, hours ) IMPLICIT NONE real longitude real latitude integer year integer month integer day real hours c Parameters real*8 M_PI parameter( M_PI = 3.14159265358979323846d0 ) real*8 M_PI_2 parameter( M_PI_2 = M_PI/2.0d0) real*8 DEGREES_TO_RADIANS parameter( DEGREES_TO_RADIANS = M_PI / 180.d0 ) integer yr integer mnth real*8 a real*8 b real*8 c real*8 e real*8 f real*8 numberOfDays real*8 numberOfDaysUT real*8 numberCenturies real*8 numberJulianCenturies real*8 siderialTime0Greenwich integer nHours real*8 localSiderialTime real*8 solarMeanLong real*8 solarMeanAnomoly real*8 solarCentre real*8 solarTrueLong real*8 numPI real*8 solarObliquity real*8 sunRA real*8 sunDec real*8 hourAngle real*8 altitude c ! Number of days and centuries from epoch J2000.0 c ! Number of days from J2000.0 Jan 1st c ! Get JD for this date if( 1 .eq. Month .or. 2 .eq. Month )then yr = Year - 1 mnth = Month+12 else yr = Year; mnth = Month; endif a = (INT(yr/100.))*1. b = (INT(a/4.))*1. c = 2.-a+b e = (INT(365.25*(yr+4716)))*1. f = (INT(30.6001*(mnth+1.)))*1. c ! This is number of days since 2000 1 jan 0h UT numberOfDays = c + e + f + Day - 1524.5 - 2451545.5; numberJulianCenturies = numberOfDays / 36525.; c ! Siderial time c ! Taken from http://aa.usno.navy.mil/faq/docs/GAST.html siderialTime0Greenwich = 6.697374558 + &0.06570982441908*numberOfDays + & 1.00273790935*Hours nHours = INT(siderialTime0Greenwich/24.) c ! Get to range 0 - 24 siderialTime0Greenwich = siderialTime0Greenwich - nHours*24. localSiderialTime = siderialTime0Greenwich + longitude/15. if( 0 .gt. localSiderialTime )then localSiderialTime = localSiderialTime + 24. else if( 24. .lt. localSiderialTime )then localSiderialTime = localSiderialTime - 24. endif c ! Solar coordinates numberOfDaysUT = numberOfDays + Hours/24. numberCenturies = numberOfDaysUT / 36525. solarMeanLong = 280.466 + 36000.770 * numberCenturies solarMeanAnomoly = 357.529 + 35999.050 * numberCenturies solarCentre = (1.915 - 0.005*numberCenturies)* &sin(solarMeanAnomoly*DEGREES_TO_RADIANS) + &0.020*sin(2*solarMeanAnomoly*DEGREES_TO_RADIANS) c ! Get true ecliptic longitude solarTrueLong = solarMeanLong + solarCentre c ! Make sure we're in 0 - 360. numPI = INT(solarTrueLong / 360.) if( 0 .le. solarTrueLong )then solarTrueLong = solarTrueLong - numPI * 360. else solarTrueLong = (numPI+1) * 360. - solarTrueLong endif solarTrueLong = solarTrueLong * DEGREES_TO_RADIANS solarObliquity = (23.439 - 0.013 * numberCenturies)* & DEGREES_TO_RADIANS sunRA = tan(solarTrueLong) * cos(solarObliquity) c ! Make sure we're in the right quadrant if( 0. .le. solarTrueLong .and. M_PI_2 .ge. solarTrueLong )then sunRA = atan(sunRA) else if ( M_PI_2 .lt. solarTrueLong .and. & 1.5*M_PI .ge. solarTrueLong )then sunRA = M_PI + atan(sunRA) else if ( 1.5*M_PI .le. solarTrueLong .and. &2.*M_PI .ge. solarTrueLong )then sunRA = 2*M_PI + atan(sunRA) endif sunDec = sin(sunRA) * sin(solarObliquity) hourAngle = localSiderialTime*15. - sunRA/DEGREES_TO_RADIANS altitude = sin(latitude*DEGREES_TO_RADIANS)*sunDec + &cos(latitude*DEGREES_TO_RADIANS)*cos(asin(sunDec))* &cos(hourAngle*DEGREES_TO_RADIANS); c ! Sun angle in degrees get_sun_angle_navpix = asin(altitude) c! if (get_sun_angle .le. 0.0) then c! get_sun_angle = M_PI+get_sun_angle c! get_sun_angle = get_sun_angle - M_PI c! else if (get_sun_angle .gt. 0.0)then c! get_sun_angle = get_sun_angle - M_PI c! get_sun_angle = get_sun_angle + M_PI c! endif get_sun_angle_navpix = SNGL((M_PI/2-get_sun_angle_navpix)/ & DEGREES_TO_RADIANS) end function c---------------------------------------------------------------------------------------- subroutine day_number_to_date_navpix( year, day_num, month, day ) integer year integer day_num integer month integer day integer i integer Month_Day_Number(12) data Month_Day_Number/0 , 31 , 59 , 90 , 120 , 151 , 181 , 212 , & 243 , 273 , 304 , 334/ integer Month_Day_Number_Leap(12) data Month_Day_Number_Leap/0 , 31 , 60 , 91 , 121 , 152 , 182 , & 213 , 244 , 274 , 305 , 335 / month = 0 if( 0 .eq. MOD(Year,4) )then do i=12,1,-1 if( day_num .gt. Month_Day_Number_Leap(i) )then month = i goto 101 endif end do 101 continue if( 0 .eq. month )then write(*,'(''ERROR: day_number_to_day: Month not found'')') stop endif day = day_num - Month_Day_Number_Leap(month) else do i=12,1,-1 if( day_num .gt. Month_Day_Number(i) )then month = i goto 102 endif end do 102 continue if( 0 .eq. month )then write(*,'(''ERROR: day_number_to_day: Month not found'')') stop endif day = day_num - Month_Day_Number(month) endif end subroutine