!------------------------------------------------------------------------------ ! Global Modeling and Assimilation Office (GMAO) ! ! Goddard Earth Observing System (GEOS) ! ! MAPL Component ! !------------------------------------------------------------------------------ ! #include "unused_dummy.H" ! !> !### MODULE: `MAPL_NominalOrbitsMod` ! ! Author: GMAO SI-Team ! ! The module `MAPL_NominalOrbitsMod` provides necessary functions for generating ! ground tracks, and corresponding masks for polar orbiter satellites. ! !#### History !- 07Jul2009 Albayrak Initial implementation. !- 30Jul2009 Albayrak Beta version 2 implemantation ! MODULE MAPL_NominalOrbitsMod use MAPL_Constants IMPLICIT NONE PRIVATE ! !PUBLIC MEMBER FUNCTIONS: PUBLIC Orbits_Track ! Generate satellite ground tracks PUBLIC Orbits_Swath ! Generate satellite ground tracks' mask PUBLIC Orbits_Track0 ! Information provider for track ! PUBLIC Orbits_Mask0 ! Information provider for mask ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! Internal constants (not public!) ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) ! 1- Satellite names INTEGER, PARAMETER :: Aqua = 1 INTEGER, PARAMETER :: Calipso = 2 INTEGER, PARAMETER :: CloudSat = 3 INTEGER, PARAMETER :: Aura = 4 INTEGER, PARAMETER :: Terra = 5 ! 2- Extrapolation Coef INTEGER, PARAMETER :: Num_sat = 5 ! Number of satellites INTEGER, PARAMETER :: Num_coef = 9 ! extrapolation coef REAL(dp), DIMENSION (Num_sat*Num_coef), PARAMETER :: Coefs = (/ & !Aqua 91.499324355103511, -3.203889864874015, 0.000273159230172, & 91.499293811827172, -0.377148321262191, 0.000091480174841, & 91.499384350660918, 1.280726447439133, -0.001159920580549, & 91.499573289741733, -3.297445861946248, 0.000239609802876, & !Calipso 91.499584024805969, -0.460273580159880, 0.000125365994414, & 91.499608453633982, 1.205360260436871, -0.001192777696086, & 91.499480902088152, -3.285668060865482, 0.000203038327945, & !CloudSat 91.499455508570492, -0.448712225445062, 0.000158569625164, & 91.499464053621409, 1.216826349594128, -0.001182646395663, & 91.499391954502713, 2.573431689964800, 0.000258633680290, & !Aura 91.499394689118887, -0.873061358923764, 0.000110177145433, & 91.499402404413175, 0.792159351805744, -0.001228074893676, & 91.499598756142376, -1.925320207499946, -0.000582299168678, & !Terra 91.499618148202927, -1.413146107480887, 0.000030982760855, & 91.499493908174912, -3.023526958545179, -0.001275012087734 & /) ! Earth/Time related constants REAL(dp), PARAMETER :: g1p = 0.0172027912; ! increase in the lon of GR per day (rad/day) ! g2p need to be calculated (optimized) version in matlab values are copied here ! Rotational rate of the earth rad per day; this is calculated for ! each sat. Think of a parameter .... REAL(dp), DIMENSION(Num_sat), PARAMETER :: g2p = (/ & 6.283200, & !Aqua 6.283295, & !Calipso 6.283200, & !CloudSat 6.283185, & !Aura 6.283190 & !Terra /) INTEGER, PARAMETER :: Num_normcoef = 3 ! extrapolation coef ! Learning parameters REAL(dp), DIMENSION(1:Num_sat*3), PARAMETER :: Norm_factor = (/ & 3.426270004720844e+06, 5.456731200281643e+06, 6.292978192496602e+06, & !Aqua 3.627175963236969e+06, 5.325345682183397e+06, 6.292978192496602e+06, & !Calipso 3.622642340324375e+06, 5.328452165105113e+06, 6.292978192496602e+06, & !CloudSat 3.615314416474894e+06, 5.333592634678283e+06, 6.292978192496602e+06, & !Aura 1.922132316029709e+06, 6.149792370758657e+06, 6.292978192496602e+06 & !Terra /) ! INTEGER, PARAMETER :: Day2year = 81 ! First day of the year for learning ! INTEGER, PARAMETER :: RefDate = 20090323 ! this is date learning started ! INTEGER, PARAMETER :: RefTime = 000000 INTEGER, DIMENSION(1:Num_sat), PARAMETER :: Day2year = (/ 281, 281, 281, 281, 281/) ! First day of the year for learning INTEGER, DIMENSION(1:Num_sat), PARAMETER :: RefDate = (/20091009, 20091009,20091009,20091009,20091009/) ! this is date learning started INTEGER, DIMENSION(1:Num_sat), PARAMETER :: RefTime = (/000000,000000,000000,000000,210000/) ! 3-Other Coef REAL(dp), PARAMETER :: myeps = 0.0000001 !for zero checks ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! END OF INTERNAL CONSTANTS ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! Definition of overloaded functions ........... interface linspace module procedure linspace1 module procedure linspace2 end interface interface Orbits_Track0 module procedure Orbits_Track1 module procedure Orbits_Track2 end interface ! ERROR codes ! 0 - OK ! 3 - ALLOCATE (memory) problem ! 90 - Swath width is not given correctly ! 99 - Satellite name is not correct CONTAINS !------------------------------------------------------------------------- !> ! The subroutine `Orbits_Track` calculates the satellite ground track during ! a given time interval. These calculations are performed using pre-calculated ! coefficients using combination of coordinate transformations, optimization ! techniques and mathematical modeling techniques. ! SUBROUTINE Orbits_Track(lons, lats, Sat_name, nymd, nhms, deltat, rc) IMPLICIT NONE ! !INPUT PARAMETERS: character(len=*), intent(in) :: Sat_name ! Satellite name integer, intent(in) :: nymd(2) ! Beginning/ending date: YYYYMMDD integer, intent(in) :: nhms(2) ! Beginning/ending time: HHMMSS integer, intent(in) :: deltat ! Time step [secs] ! !OUTPUT PARAMETERS: real(dp), pointer :: lons(:) ! Ground track longitudes [degrees] real(dp), pointer :: lats(:) ! Ground track latitudes [degrees] integer, intent(out):: rc ! Error code = 0 all is well ! = 3 memory allocation error REAL(dp) :: time_day ! ! Other parameters! REAL(dp), DIMENSION(:), ALLOCATABLE :: obstime real(dp) :: sim_start, sim_end real(dp) :: s_fraction2day, e_fraction2day INTEGER :: ts, las, los integer :: Sat integer :: ierr1, ierr2, ierr3 REAL(dp) :: distlat !lat is equal dis INTEGER :: say !()()()()()() ! ............................................................................ rc = 0 ! Initiate error code to 0 CALL satname2int(Sat_name, Sat, rc) if ( rc /= 0 ) return CALL get_time(Sat, nymd(1), nhms(1), s_fraction2day, sim_start) ! handle time CALL get_time(Sat, nymd(2), nhms(2), e_fraction2day, sim_end) time_day = deltat/(24.0*60.0*60.0) !()()()() CALL built_vecsize(sim_start, sim_end, time_day, say) ! ()()()()() ALLOCATE(lons(1:say), stat=ierr1) ! deallocate in the driver ALLOCATE(lats(1:say), stat=ierr2) ! deallocate in the driver ALLOCATE(obstime(1:say),stat=ierr3) ! deallocate in the driver if ( ierr1 /= 0 .or. ierr2 /= 0 .or. ierr3 /= 0 ) then rc = 3 ! print*, rc return endif distlat = 1 !for flagauto 0 option it is not used. later change to optional CALL find_waypointslat(Sat, deltat,sim_start,sim_end,distlat, & obstime, lats, lons, ts, las, los, rc ) if (rc>0) return END SUBROUTINE Orbits_Track ! SUBROUTINE Orbits_Track0(Sat_name) ! IMPLICIT NONE ! ! character(len=*), intent(in) :: Sat_name ! Satellite name ! integer :: Sat ! ! Sat = satname2int(Sat_name) ! Get satellite coef ! print*, " " ! print*, "----------------------------------------------------------" ! print*, "satellite name: ", Sat_name, " and carresponding Sat number: ", Sat ! print*, "Learning is performed on:", RefDate, " at:", 000000, " hour/min/sec" ! print*, "This is the ", Day2year, " th day of the year" ! ! print*, "----------------------------------------------------------" ! print*, " " ! ! END SUBROUTINE Orbits_Track0 ! subroutine orb_getsat(sat,sat_name,rc) ! integer, intent(in) :: sat ! integer, intent(out) :: rc ! character(len=*), intent(out) :: sat_name ! ! if ( sat == AQUA ) then ! sat_name = "Aqua" ! else if ( sat == CALIPSO ) then ! sat_name = "Calipso" ! else if ( sat == CLOUDSAT ) then ! sat_name = "Cloudsat" ! else if ( sat == AURA ) then ! sat_name = "Aura" ! else ! rc = 99 ! return ! end if ! rc = 0 ! end Subroutine satname2int(name, satnum, rc) IMPLICIT NONE character(len=*), intent(in) :: name ! Satellite name integer, intent(out) :: satnum ! Satellite number integer, intent(out) :: rc rc = 0 if ((name.EQ."AQUA").or.(name.EQ."aqua").or.(name.EQ."Aqua")) then satnum = 1 elseif ((name.EQ." calipso").or.(name.EQ."CALIPSO").or.(name.EQ."Calipso")) then satnum = 2 elseif ((name.EQ." cloudsat").or.(name.EQ."CLOUDSAT").or.(name.EQ."CloudSat")) then satnum = 3 elseif ((name.EQ." aura").or.(name.EQ."AURA").or.(name.EQ."Aura")) then satnum = 4 elseif ((name.EQ." terra").or.(name.EQ."TERRA").or.(name.EQ."Terra")) then satnum = 5 else rc =99 return endif end Subroutine satname2int SUBROUTINE Orbits_Track1(name, Satcoef_vec1, Normcoef_vec1) ! See interface IMPLICIT NONE character(len=*), intent(in) :: name ! Satellite name real(dp), intent(out), dimension(1:Num_coef) :: Satcoef_vec1 real(dp), intent(out), dimension(1:Num_normcoef) :: Normcoef_vec1 integer :: start_ind, end_ind integer :: Normcoef_ind integer :: iSat integer :: rc ! ................................... ! Num_sat, Num_coef needs to be defined at the begining of the module ! Num_coef is 9 for each sat... ! Also Norm factor coef. need to be retrived. ! Here locations of these vectors are calculated ! ................................... CALL satname2int(name, iSat, rc) if ( rc /= 0 ) return start_ind = (iSat - 1) * Num_coef + 1 !CALL calc_ind end_ind = start_ind + (Num_coef-1) Satcoef_vec1(1:9) = Coefs(start_ind:end_ind) ! x,y,z ! % ________________Mult by Norm fator_________________________ !Normcoef_ind = (iSat-1)*Num_sat+1 Normcoef_ind = iSat*Num_normcoef-(Num_normcoef-1) ! print*, Normcoef_ind ! print*, start_ind, end_ind Normcoef_vec1 = (/Norm_factor(Normcoef_ind), & Norm_factor(Normcoef_ind+1), & Norm_factor(Normcoef_ind+2)/) END SUBROUTINE Orbits_Track1 SUBROUTINE Orbits_Track2(iSat, Satcoef_vec, Normcoef_vec) ! See interface IMPLICIT NONE integer, intent(in) :: iSat ! this is the sat number ex: Aqua=1 real(dp), intent(out), dimension(1:Num_coef) :: Satcoef_vec real(dp), intent(out), dimension(1:Num_normcoef) :: Normcoef_vec integer :: start_ind, end_ind integer :: Normcoef_ind ! ................................... ! Num_sat, Num_coef needs to be defined at the begining of the module ! Num_coef is 9 for each sat... ! Also Norm factor coef. need to be retrived. ! Here locations of these vectors are calculated ! ................................... start_ind = (iSat - 1) * Num_coef + 1 !CALL calc_ind end_ind = start_ind + (Num_coef-1) Satcoef_vec(1:9) = Coefs(start_ind:end_ind) ! x,y,z ! % ________________Mult by Norm fator_________________________ !Normcoef_ind = (iSat-1)*Num_sat+1 Normcoef_ind = iSat*Num_normcoef-(Num_normcoef-1) ! print*, Normcoef_ind ! print*, start_ind, end_ind Normcoef_vec = (/Norm_factor(Normcoef_ind), & Norm_factor(Normcoef_ind+1), & Norm_factor(Normcoef_ind+2)/) END SUBROUTINE Orbits_Track2 ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE Orbits_Swath (slons, slats, Sat_name, nymd, nhms, deltat, & SwathWidth, rc, wrap) ! SUBROUTINE Orbits_Swath (slons, slats, Sat_name, nymd, nhms, deltat, SwathWidth, wrapon, rc) IMPLICIT NONE ! !INPUT PARAMETERS: real(dp), pointer :: slons(:,:) !Track longitude [degree] real(dp), pointer :: slats(:,:) !Track latitude [degree] character(len=*),intent(in) :: Sat_name !Satellite name Aqua etc ... integer, intent(in) :: nymd(2) !Beginning/ending date: YYYYMMDD integer, intent(in) :: nhms(2) !Beginning/ending time: HHMMSS REAL(dp), INTENT(IN) :: SwathWidth(:) !Right, left swath [km] integer, intent(in) :: deltat ! Time step [secs] LOGICAL, intent(in), optional:: wrap ! if true then wrap to -180 to 180 !LOGICAL :: wrapon ! if true then wrap to -180 to 180 ! !OUTPUT PARAMETERS: integer, intent(out) :: rc ! Other parameters! real(dp), pointer :: lons(:) !Ground Track longitude [degree] real(dp), pointer :: lats(:) !Ground Track latitude [degree] REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: latlonshift real(dp) :: temp1 integer, parameter :: iii = 3 !this is used for the number of swaths !here is 3 indicating mid left right integer :: Sat ! Satellite number integer :: ierr1 logical :: wrapon if ( present(wrap) ) then wrapon = wrap else wrapon = .true. end if rc = 0 ! Initiate error code to 0 CALL Orbits_Track(lons, lats, Sat_name, nymd, nhms, deltat, rc) !rc will return 3 if ( rc /= 0 ) return CALL satname2int(Sat_name, Sat, rc) ! name to sat number rc will return 99 if ( rc /= 0 ) return ! Check if the Swath width is given correctly if ( SwathWidth(1)<0.AND.SwathWidth(2) < 0 ) then rc = 90 return else if (SwathWidth(1)<0.AND.SwathWidth(2) >= 0) then temp1 = SwathWidth(2) + SwathWidth(1) if (temp1<=0) then rc = 90 return end if else if (SwathWidth(2)<0.AND.SwathWidth(1) >= 0) then temp1 = SwathWidth(1) + SwathWidth(2) if (temp1<=0) then rc = 90 return end if end if ! 2- Find sweep points if (allocated(latlonshift)) deallocate(latlonshift) ! this is main ALLOCATE(latlonshift(iii, SIZE(lats) ,2), stat = ierr1 ) ! either is ok size of lat or lon if ( ierr1 /= 0 ) then rc = 3 ! return endif latlonshift = 0 !CALL find_sweeppoints(iii,lat_l(1:las),lon_l(1:los), real(swath(1),dp), & ! real(swath(2),dp), latlonshift) CALL find_sweeppoints(iii, lats, lons, -SwathWidth(1), SwathWidth(2), wrapon, latlonshift) ALLOCATE(slats(iii, SIZE(lats)), stat = ierr1 ) ! either is ok size of lat or lon ALLOCATE(slons(iii, SIZE(lons)), stat = ierr1 ) ! either is ok size of lat or lon slats(1:iii, 1:SIZE(lats)) = latlonshift(1:iii, 1:SIZE(lats) ,1) slons(1:iii, 1:SIZE(lons)) = latlonshift(1:iii, 1:SIZE(lons) ,2) slats(2,1:SIZE(lats)) = lats(1:SIZE(lats)) ! slats(1:3, :) 1- right 2- original ground track 3-left slons(2, 1:SIZE(lons)) = lons(1:SIZE(lons)) DEALLOCATE(latlonshift, stat = ierr1) if (ierr1 /= 0) then rc = 3 ! return endif END SUBROUTINE Orbits_Swath ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !------------------------------------------------------------------------- ! ----------------- ! INTERNAL ROUTINES ! ----------------- !------------------------------------------------------------------------- Subroutine do_gridcorrections(lon_vec_centers1, lat_vec_centers1, & lon_vec_centers, lat_vec_centers, rc) IMPLICIT NONE REAL(dp), DIMENSION(:), INTENT(IN) :: lon_vec_centers1, lat_vec_centers1 REAL(dp), DIMENSION(:), INTENT(INOUT) :: lon_vec_centers, lat_vec_centers integer, intent(inout) :: rc ! one need to have data center analysis example lat has to be in between ! 90 to -90, lon has to be in between -180 to 180 if (size(lon_vec_centers1) < 30) then print*, " you need to have minimum of 30 members in lon centers" stop end if if (size(lat_vec_centers1) < 30) then print*, " you need to have minimum of 30 members in lat centers" stop end if lon_vec_centers = lon_vec_centers1 if ( (lon_vec_centers1(1)>180) .OR. (lon_vec_centers1(size(lon_vec_centers1))>180) ) then lon_vec_centers(1:size(lon_vec_centers1)) = lon_vec_centers1(1:size(lon_vec_centers1)) - 180.00 end if if (lon_vec_centers1(1) > 0) then CALL updownlon(lon_vec_centers1, lon_vec_centers, rc) if(rc>0) return end if if (lat_vec_centers1(1) < 0) then CALL updownlon(lat_vec_centers1, lat_vec_centers, rc) if (rc>0) return else lat_vec_centers = lat_vec_centers1 end if end subroutine do_gridcorrections SUBROUTINE updownlon(vec_centers, upsidedown, rc) IMPLICIT NONE REAL(dp), DIMENSION(:), INTENT(IN) :: vec_centers REAL(dp), DIMENSION(:), INTENT(OUT) :: upsidedown integer, intent(inout) :: rc REAL(dp), DIMENSION(:), allocatable :: tempcenters INTEGER :: i, ierr1 allocate( tempcenters(1:size(vec_centers)), stat=ierr1) if (ierr1 /= 0) then rc = 3 ! return endif tempcenters = vec_centers DO i=1, size(vec_centers) upsidedown(size(vec_centers)+1-i) = tempcenters(i) END DO deallocate(tempcenters, stat=ierr1) if (ierr1 /= 0) then rc = 3 ! return endif END SUBROUTINE updownlon !------------------------------------------------------------------------------ !> ! The subroutine `find_waypointslat` ! finds waypoints for the possible max distance on the equator line. ! If flagauto = 0 then equal deltat data is created. Otherwise optimum ! number of minimum data is created (nonlinear). ! SUBROUTINE find_waypointslat(iSat, time_sec, sim_start,sim_end, distlat, & t1_l, lat_l, lon_l, ts, las, los, rc) ! sat name - time flag - interval second IMPLICIT NONE ! !ARGUMENTS: ! INTEGER, INTENT(IN) :: iSat, time_sec REAL(dp), INTENT(IN) :: sim_start, sim_end REAL(dp), OPTIONAL, INTENT(IN) :: distlat REAL(dp), DIMENSION(:), INTENT(INOUT):: t1_l, lat_l, lon_l ! data matrices (dm) INTEGER, INTENT(OUT) :: ts, las, los ! size of dm !!! INTEGER, INTENT(OUT):: time integer, intent(inout):: rc REAL(dp), DIMENSION(:), ALLOCATABLE :: t1 REAL(dp), DIMENSION(:,:),ALLOCATABLE :: lat_estwayp, long_estwayp ! temp arrays REAL(dp) :: latest, lonest REAL(dp) :: time_day INTEGER :: say, i INTEGER :: ierr1, ierr2, ierr3 ! _UNUSED_DUMMY(distlat) IF (time_sec.LE.0.0) THEN print*, "Error: time is negative. For now it is not allowed" STOP ENDIF time_day = time_sec/(24.0*60.0*60.0) CALL built_vecsize(sim_start, sim_end, time_day, say) ALLOCATE (t1(1:say), stat = ierr1) ALLOCATE (lat_estwayp(1:say,1), stat = ierr2) ALLOCATE (long_estwayp(1:say,1), stat = ierr3) if ( ierr1 /= 0 .or. ierr2 /= 0 .or. ierr3 /= 0 ) then rc = 3 ! return endif CALL built_vec(sim_start, sim_end, time_day, t1) DO i=1,size(t1) CALL get_latlon(iSat, t1(i), latest, lonest) lat_estwayp(i,1) = latest long_estwayp(i,1) = lonest END DO ! i ! % shift long (-pi) long_estwayp = (long_estwayp-180.0); ts = size(t1) las = size(lat_estwayp,1) los = size(long_estwayp,1) t1_l(1:ts) = t1 lat_l(1:las) = lat_estwayp(:,1) lon_l(1:los) = long_estwayp(:,1) DEALLOCATE (t1, stat=ierr1) DEALLOCATE (lat_estwayp, stat=ierr2) DEALLOCATE (long_estwayp, stat=ierr3) if ( ierr1 /= 0 .or. ierr2 /= 0 .or. ierr3 /= 0 ) then rc = 3 ! return endif END SUBROUTINE find_waypointslat !------------------------------------------------------------------------------ !> ! Swap points are the points that gives the satellite instrument view ! points that is different than the original track. For now, calculation is basix. ! `r`, and `l` give how far sweep points will go vertically from the orbit points. ! !#### See Also !- get_recon !- get_azimuth ! SUBROUTINE find_sweeppoints(iii, latwayp, longwayp, l, r, wrapon, latlonshift) IMPLICIT NONE ! REAL(dp), DIMENSION(:), INTENT(IN) :: latwayp, longwayp REAL(dp), INTENT(IN) :: l, r ! left right swap space INTEGER, INTENT(IN) :: iii REAL(dp), DIMENSION(:,:,:), INTENT(INOUT):: latlonshift LOGICAL, intent(in) :: wrapon ! if true then wrap to -180 to 180 REAL(dp), DIMENSION(:), ALLOCATABLE :: myvec, myvec_deg REAL(dp), DIMENSION(:), ALLOCATABLE :: latshift1, lonshift1 REAL(dp) :: az, latout1, lonout1 INTEGER :: say, i, count1 ! ALLOCATE(myvec(iii), myvec_deg(iii)) if (r>l) then CALL linspace(l,r,iii,myvec) else CALL linspace(r,l,iii,myvec) end if do i=1, size(myvec) myvec_deg(i) = MAPL_KM_PER_DEG * myvec(i) end do DO say = 1,size(myvec) ALLOCATE(latshift1(1:size(longwayp)),lonshift1(1:size(longwayp))) count1 = 1 DO i=1,size(longwayp)-1 az = get_azimuth(latwayp(i),longwayp(i),latwayp(i+1),longwayp(i+1), wrapon) CALL get_reckon(latwayp(i),longwayp(i), az+90.0, & myvec_deg(say), latout1, lonout1, wrapon) latshift1(count1) = latout1 lonshift1(count1) = lonout1 count1 = count1 + 1 END DO ! i increased by 1 CALL get_reckon(latwayp(i),longwayp(i), az+90.0, & myvec_deg(say), latout1, lonout1, wrapon) latshift1(count1) = latout1 lonshift1(count1) = lonout1 latlonshift(say,:,1) = latshift1(:) latlonshift(say,:,2) = lonshift1(:) DEALLOCATE(latshift1, lonshift1) END DO DEALLOCATE(myvec, myvec_deg) END SUBROUTINE find_sweeppoints !------------------------------------------------------------------------------ !> ! The subroutine `get_latlon` ! calls three important functions to calculate lat lon values of ! a satellite for a given time. Those functions are described in the ! following subsections. ! SUBROUTINE get_latlon(iSat, t1, lat_estwayp, long_estwayp) IMPLICIT NONE REAL(dp), INTENT(IN) :: t1 INTEGER, INTENT(IN) :: iSat REAL(dp), INTENT(INOUT) :: lat_estwayp, long_estwayp REAL(dp), DIMENSION(3,1) :: ECI_est, ECEF_est REAL(dp) :: lat, lon, alt, fraction2day ! CALL get_estimateECI(iSat, t1, ECI_est) fraction2day = get_fraction(t1) CALL ECI2ECEF(iSat,fraction2day, ECI_est, ECEF_est) CALL ECEF2LLA(ECEF_est(1,1), ECEF_est(2,1), & ECEF_est(3,1), lat, lon, alt) lat_estwayp = rad2deg(lat) long_estwayp = rad2deg(lon) END SUBROUTINE SUBROUTINE get_estimateECI(iSat, ntime_day, ECI_est) IMPLICIT NONE INTEGER, INTENT(IN) :: iSat REAL(dp), INTENT(IN) :: ntime_day REAL(dp), DIMENSION(:,:), INTENT(OUT) :: ECI_est REAL(dp), DIMENSION(Num_coef) :: myvec REAL(dp), DIMENSION(Num_normcoef) :: normvec REAL(dp) :: xdata_est, ydata_est, zdata_est CALL Orbits_Track0(iSat, myvec, normvec) ! myvec has 9 sat coef ! normvec has normalization coef ! % _________________________Estimates_________________________ xdata_est = sin(myvec(1) * ntime_day + myvec(2)) + myvec(3); ydata_est = sin(myvec(4) * ntime_day + myvec(5)) + myvec(6); zdata_est = sin(myvec(7) * ntime_day + myvec(8)) + myvec(9); ! % ________________Mult by Norm fator_________________________ xdata_est = xdata_est * normvec(1) ! max(xECI); ydata_est = ydata_est * normvec(2) ! max(yECI); zdata_est = zdata_est * normvec(3) ! max(zECI); ECI_est(1,1) = xdata_est ECI_est(2,1) = ydata_est ECI_est(3,1) = zdata_est END SUBROUTINE get_estimateECI SUBROUTINE ECI2ECEF(iSat,fraction2day, ECI_est, ECEF_est) IMPLICIT NONE INTEGER, INTENT(IN) :: iSat REAL(dp), DIMENSION(:,:), INTENT(IN) :: ECI_est REAL(dp), DIMENSION(:,:), INTENT(OUT) :: ECEF_est REAL(dp), DIMENSION(3,3) :: RotM ! rotation matrix from In2ECEF INTEGER :: t1 REAL(dp) :: Goft, a1, a2 REAL(dp) :: fraction2day associate(g0 => MAPL_OMEGA_R8) t1 = Day2year(iSat) Goft = g0 + g1p*t1 + g2p(iSat) * fraction2day a1 = cos(Goft); a2 = sin(Goft); RotM = 0.0 RotM(1,1) = a1 RotM(1,2) = a2 RotM(2,1) = -a2 RotM(2,2) = a1 RotM(3,3) = 1 ECEF_est = MATMUL(RotM, ECI_est) end associate END SUBROUTINE ECI2ECEF SUBROUTINE get_time(iSat, nymd, nhms, fraction2day, ntime_day) IMPLICIT NONE INTEGER, INTENT(IN) :: iSat, nymd, nhms REAL(dp), INTENT(OUT) :: fraction2day, ntime_day INTEGER :: Year, Month, Day, Hour, Minute, Seconds Year = int(nymd / 10000 ) Month = mod ( nymd, 10000 ) / 100 Day = mod ( nymd, 100 ) Hour = int(nhms / 10000 ) Minute = mod ( nhms, 10000 ) / 100 Seconds = mod ( nhms, 100 ) fraction2day = (Hour*60*60+Minute*60+Seconds)/(24.0*60.0*60.0) !fraction !calculated from seconds ntime_day=(-ODS_Julian(RefDate(iSat))+ODS_Julian(nymd))+fraction2day END SUBROUTINE get_time SUBROUTINE ECEF2LLA(x, y, z, lat, lon, alt) ! % ECEF3LLA - convert earth-centered earth-fixed (ECEF) ! % cartesian coordinates to latitude, longitude, ! % and altitude ! % ! % USAGE: ! % [lat,lon,alt] = ecef2lla(x,y,z) ! % ! % lat = geodetic latitude (radians) ! % lon = longitude (radians) ! % alt = height above WGS84 ellipsoid (m) ! % x = ECEF X-coordinate (m) ! % y = ECEF Y-coordinate (m) ! % z = ECEF Z-coordinate (m) ! % ! % Notes: (1) This function assumes the WGS84 model. ! % (2) Latitude is customary geodetic (not geocentric). ! % (3) Inputs may be scalars, vectors, or matrices of the same ! % size and shape. Outputs will have that same size and shape. ! % (4) Tested but no warranty; use at your own risk. ! % (5) Michael Kleder, April 2006 REAL(dp), INTENT(IN) :: x,y,z REAL(dp), INTENT(OUT) :: lat,lon,alt ! % WGS84 ellipsoid constants: REAL(dp) :: b, ep, p, th, N associate(a => MAPL_EARTH_SEMIMAJOR_AXIS, e => MAPL_EARTH_ECCENTRICITY, pi => MAPL_PI_R8) b = sqrt(a**2 * (1-e**2)) ep = sqrt((a**2-b**2)/b**2) p = sqrt(x**2+y**2) th = atan2(a*z,b*p) lon = atan2(y,x) lat = atan2((z+ep**2*b*sin(th)**3),(p-e**2*a*cos(th)**3)) N = a/sqrt(1-e**2*sin(lat)**2) alt = p/cos(lat)-N !% !return lon in range [0,2*pi) !lon = mod(lon,2*pi) ! does not the same as matlab lon = lon-floor(lon/(2*pi))*(2*pi) ! % correct for numerical instability in altitude near exact poles: ! % (after this correction, error is about 2 millimeters, which is about ! % the same as the numerical precision of the overall function) !k= (abs(x)<1.AND.abs(y)<1 ) !print *, k !!alt(k) = abs(z(k))-b end associate END SUBROUTINE ECEF2LLA integer function ODS_Julian ( CalDate ) implicit NONE integer CalDate ! Calendar date in the format YYYYMMDD ! where YYYY is the year, MM is the ! month and DD is the day. A negative ! number implies that the year is B.C. ! Other variables ! --------------- integer Year integer Month integer Day integer iGreg ! Gregorian Calendar adopted Oct 12, 1582 parameter ( iGreg = 15 + 31 * ( 10 + 12 * 1582 ) ) integer JulDay integer jy, jm, ja Year = CalDate / 10000 Month = mod ( CalDate, 10000 ) / 100 Day = mod ( CalDate, 100 ) ! Change year 0 to year 1 if ( Year .eq. 0 ) Year = 1 ! Account for the nonexisting year 0 if ( Year .lt. 0 ) Year = Year + 1 if ( Month .gt. 2 ) then jy = Year jm = Month + 1 else jy = Year - 1 jm = Month + 13 endif JulDay = int ( 365.25 * jy ) & + int ( 30.6001 * jm ) & + Day + 1720995 ! Test whether to change to Gregorian Celendar if ( Day + 31 * ( Month + 12 * Year ) .ge. iGreg) then ja = int ( 0.01 * jy ) Julday = JulDay + 2 - ja + int ( 0.25 * ja ) endif ODS_Julian = JulDay return end function ODS_Julian ! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE linspace1(d1,d2,n,y) ! from d1 to d2, n by n -> y IMPLICIT NONE REAL(dp), INTENT(IN) :: d1, d2, n INTEGER :: n2, i REAL(dp), DIMENSION(:), INTENT(INOUT) :: y n2 = n DO i = 0, n2-2 y(i+1) = i END DO y(1:n2-1) = (d1+y(1:n2-1)*(d2-d1)/(floor(n)-1)) ! d2); y(n2) = d2 END SUBROUTINE linspace1 SUBROUTINE linspace2(d1,d2,n,y) IMPLICIT NONE REAL(dp), INTENT(IN) :: d1, d2 INTEGER, INTENT(IN) :: n INTEGER :: n2, i REAL(dp), DIMENSION(:), INTENT(INOUT) :: y n2 = n DO i = 0, n2-2 y(i+1) = i END DO y(1:n2-1) = (d1+y(1:n2-1)*(d2-d1)/(n-1)) ! d2); y(n2) = d2 END SUBROUTINE linspace2 SUBROUTINE build_array(s,e,dp, mat) IMPLICIT NONE REAL, DIMENSION(:) :: mat REAL :: s,dp ! start, end, delta increase, matrix INTEGER :: k, e mat(1) = s do k=2,e mat(k) = mat(k-1) + dp end do END SUBROUTINE build_array REAL(dp) FUNCTION deg2rad(angle) IMPLICIT NONE REAL(dp), INTENT(IN) :: angle deg2rad = MAPL_DEGREES_TO_RADIANS_R8 * angle END FUNCTION deg2rad REAL(dp) FUNCTION rad2deg(rad) IMPLICIT NONE REAL(dp), INTENT(IN) :: rad rad2deg = MAPL_RADIANS_TO_DEGREES * rad END FUNCTION rad2deg !------------------------------------------------------------------------------ !> ! The function `get_distance` computes great circle distance and azimuth ! Lat - lon in degrees. ! REAL(dp) FUNCTION get_distance(lt1, ln1, lt2, ln2) IMPLICIT NONE REAL(dp) :: lat1, lat2, lon1, lon2 REAL(dp) :: lt1, lt2, ln1, ln2 REAL(dp) :: a, r lat1 = deg2rad(lt1) lon1 = deg2rad(ln1) lat2 = deg2rad(lt2) lon2 = deg2rad(ln2) r=1.0 a = sin((lat2-lat1)/2.0)**2.0 + cos(lat1) * cos(lat2) * & sin((lon2-lon1)/2.0)**2.0; get_distance = r * 2.0 * atan2(sqrt(a),sqrt(1.0 - a)); get_distance = rad2deg(get_distance) END FUNCTION get_distance REAL(dp) FUNCTION get_azimuth(lt1, ln1, lt2, ln2, wrapon) ! need to be ! Computes great circle distance and azimuth ! Lat - lon in degrees IMPLICIT NONE REAL(dp) :: lat1, lat2, lon1, lon2 REAL(dp) :: lt1, lt2, ln1, ln2 REAL(dp) :: az, epsilone, saz LOGICAL, intent(in) :: wrapon ! if true then wrap to -180 to 180 lat1 = deg2rad(lt1) lon1 = deg2rad(ln1) lat2 = deg2rad(lt2) lon2 = deg2rad(ln2) ! Inputs LAT1, LON1, LAT2, LON2 are in units of radians. az = atan2(cos(lat2) * sin(lon2-lon1), cos(lat1) * sin(lat2)- & sin(lat1) * cos(lat2) * cos(lon2-lon1)) ! Azimuths are undefined at the poles, so we choose a convention: zero at ! the north pole and pi at the south pole. associate(pi => MAPL_PI_R8) if (lat1 <= -pi/2.0) az = 0.0 if (lat2 >= pi/2.0) az = 0.0 if (lat2 <= -pi/2.0) az = pi if (lat1 >= pi/2.0) az = pi epsilone = 1.74E-8 if (wrapon) then if ( az>-epsilone .AND. az < epsilone) then saz = 0 else saz = az/abs(az) end if az = pi*((abs(az)/pi) - 2*ceiling(((abs(az)/pi)-1)/2)) * saz end if if ( az < -epsilone) then ! change to 0 to 2pi az = az + 2*pi end if end associate ! Reset near-zero points if (az<0.0) az=0 get_azimuth = rad2deg(az) END FUNCTION get_azimuth SUBROUTINE get_reckon(phi1, lambda1, az1, rng1, phi, lambda, wrapon) ! phi0 - > lat, lambda0 -> lon ! calculates a position (LATOUT, LONOUT) at a given range RNG and azimuth ! AZ along a great circle from a starting point defined by LAT and LON. ! LAT and LON are in degrees. The range is in degrees of arc length on a ! sphere. The input azimuth is in degrees, measured clockwise from due ! north. Translated from MatLab IMPLICIT NONE REAL(dp), INTENT(IN) :: phi1, lambda1, az1, rng1 REAL(dp), INTENT(INOUT) :: phi, lambda REAL(dp) :: phi0, lambda0 REAL(dp) :: epsilone, az, saz, rng LOGICAL, intent(in) :: wrapon ! if true then wrap to -180 to 180 epsilone = 10*1.74E-8 phi0 = deg2rad(phi1) lambda0 = deg2rad(lambda1) az = deg2rad(az1) rng = deg2rad(rng1) associate(pi => MAPL_PI_R8) if (phi0 >= pi/2-epsilone) az = pi ! starting at north pole if (phi0 <= epsilone-pi/2) az = 0 ! starting at south pole ! Calculate coordinates of great circle end point using spherical trig. phi = asin(sin(phi0)*cos(rng) + cos(phi0)*sin(rng)*cos(az)) lambda = lambda0 + atan2( sin(rng)*sin(az), & cos(phi0)*cos(rng) - sin(phi0)*sin(rng)*cos(az) ) if (wrapon) then if ( lambda>-epsilone .AND. lambda < epsilone) then saz = 0 else saz = lambda/abs(lambda) end if lambda = pi*((abs(lambda)/pi) - 2*ceiling(((abs(lambda)/pi)-1)/2)) * saz end if end associate lambda = rad2deg(lambda) phi = rad2deg(phi) END SUBROUTINE get_reckon subroutine built_vecsize(sim_start, sim_end, time_day, say) implicit none real(dp), intent(in) :: sim_start, sim_end, time_day integer :: say real(dp) :: temp1 say = 1 temp1 = sim_start !print*, "inside simend-----", sim_end, temp1, time_day !WRITE(*,"(1F10.2)") sim_end do while (sim_end.GE.temp1) temp1 = temp1 + time_day if (temp1.LE.sim_end) then say = say + 1 end if enddo end subroutine built_vecsize subroutine built_vec(sim_start, sim_end, time_day, t1) implicit none real(dp), intent(in) :: sim_start, sim_end, time_day integer :: say real(dp) :: temp1 real(dp), dimension(:) :: t1 say = 1 temp1 = sim_start ! garante that it lesim start !time t1(1) = sim_start do while (sim_end.GE.temp1) temp1 = temp1 + time_day say = say + 1 if (say>size(t1)) then exit end if t1(say) = temp1 enddo end subroutine built_vec !------------------------------------------------------------------------------ !> ! The function `find_LTGE` finds minimum possible long value for carresponding ! lat interval. ! SUBROUTINE find_LTGE(distlon, lat_vec_centers, vec_lat, long_limits) IMPLICIT NONE INTEGER :: mys, i, k REAL(dp), POINTER, DIMENSION(:) :: long_limits !INTENT(INOUT) REAL(dp), DIMENSION(:),INTENT(IN)::lat_vec_centers, distlon, vec_lat REAL(dp), DIMENSION(:), ALLOCATABLE :: long_limitstemp mys = size(lat_vec_centers) ALLOCATE(long_limitstemp(1:mys)) do i=1,size(vec_lat) - 1 do k=1,mys if ( (lat_vec_centers(k).LT.vec_lat(i)).AND. & (lat_vec_centers(k).GE.vec_lat(i+1)) ) then long_limitstemp(i) = MAPL_DEG_PER_KM * distlon(k) exit end if end do end do long_limitstemp(i)=long_limitstemp(i-1) ALLOCATE(long_limits(1:i)) long_limits = long_limitstemp(1:i) DEALLOCATE(long_limitstemp) END SUBROUTINE find_LTGE SUBROUTINE redifine_timeint(kk, longind, ti, iSat, & lat_new, lon_new, t1_new, carylat, carylon, caryt) IMPLICIT NONE REAL(dp), DIMENSION(:), INTENT(INOUT) :: lat_new, lon_new, t1_new INTEGER, INTENT(INOUT) :: carylat, carylon, caryt REAL(dp), DIMENSION(:), INTENT(INOUT) :: longind INTEGER, INTENT(IN) :: kk INTEGER, INTENT(IN) :: ti REAL(dp), DIMENSION(:), ALLOCATABLE :: t1 REAL(dp), DIMENSION(:), ALLOCATABLE :: t1_f, lat_final, lon_final REAL(dp), DIMENSION(:,:), ALLOCATABLE :: latestwp, longestwp REAL(dp), DIMENSION(:), ALLOCATABLE :: t_tempo REAL(dp), DIMENSION(:,:), ALLOCATABLE:: lat_wayptemp, long_wayptemp INTEGER :: temp1, temp2, temp4, iSat, say, ii, temp_size REAL(dp) :: simtimestart_temp, simtimeend_temp, ti_day REAL(dp) :: lat_estwayp, lon_estwayp ALLOCATE(latestwp(1:carylat,1), longestwp(1:carylon,1), t1(1:caryt)) latestwp(:,1) = lat_new(1:carylat) longestwp(:,1) = lon_new(1:carylon) t1(:) = t1_new(1:caryt) ti_day = REAL(ti)/(24.0*60.0*60.0) temp1 = longind(kk) temp2 = longind(kk)+1 if (temp2 <= size(latestwp,1) ) then simtimestart_temp = t1(temp1) simtimeend_temp = t1(temp2) CALL built_vecsize(simtimestart_temp, simtimeend_temp, & ti_day, say) ALLOCATE(t_tempo(1:say)) CALL built_vec(simtimestart_temp,simtimeend_temp, & ti_day, t_tempo) ALLOCATE(lat_wayptemp(1:size(t_tempo),1) , & long_wayptemp(1:size(t_tempo),1) ) do ii=1,size(t_tempo) ! Calculate carresponding lat lon points CALL get_latlon(iSat, t_tempo(ii), lat_estwayp, lon_estwayp) lat_wayptemp(ii,1) = lat_estwayp long_wayptemp(ii,1) = lon_estwayp end do long_wayptemp = (long_wayptemp-180) temp_size = size(t1)+size(t_tempo)-2 ALLOCATE(t1_f(1:temp_size)) ALLOCATE(lat_final(1:temp_size)) ALLOCATE(lon_final(1:temp_size)) CALL add_newdata(t1, temp1, t_tempo, t1_f) CALL add_newdata(latestwp(:,1), temp1, lat_wayptemp(:,1), lat_final) CALL add_newdata(longestwp(:,1), temp1, long_wayptemp(:,1), lon_final) temp4 = size(lat_wayptemp,1)-2; longind(kk+1:size(longind)) = longind(kk+1:size(longind)) + temp4; carylat = size(lat_final) carylon = size(lon_final) caryt = size(t1_f) lat_new(1:carylat) = lat_final lon_new(1:carylon) = lon_final t1_new(1:caryt) = t1_f DEALLOCATE(t1_f, lat_final, lon_final ) DEALLOCATE(lat_wayptemp,long_wayptemp) DEALLOCATE (t_tempo) end if DEALLOCATE(latestwp, longestwp, t1) END SUBROUTINE redifine_timeint SUBROUTINE add_newdata(t1, temp1, t_tempo, rtime) IMPLICIT NONE REAL(dp), DIMENSION(:), INTENT(IN) :: t1, t_tempo REAL(dp), DIMENSION(:), INTENT(INOUT) :: rtime INTEGER, INTENT(IN) :: temp1 INTEGER :: i, ii, k ii = 1 do i=1, temp1-1 rtime(ii) = t1(i) ii=ii+1 end do do k=1, size(t_tempo) rtime(ii) = t_tempo(k) ii=ii+1 end do do i=temp1+2, size(t1) rtime(ii) = t1(i) ii=ii+1 end do END SUBROUTINE add_newdata REAL(dp) FUNCTION find_delta_t(iSat, kk, long_limits, longind,t1,time_day) IMPLICIT NONE REAL(dp) :: ti, ti_day INTEGER, INTENT(IN) :: kk REAL(dp), INTENT(IN) :: time_day REAL(dp), DIMENSION(:), INTENT(IN) :: longind REAL(dp), INTENT(IN) :: long_limits REAL(dp), DIMENSION(:), INTENT(IN) :: t1 INTEGER, INTENT(IN) :: iSat REAL(dp), DIMENSION(:), ALLOCATABLE :: t_tempo REAL(dp), DIMENSION(:,:), ALLOCATABLE :: lat_wayptemp, long_wayptemp INTEGER :: flag_add, temp1, temp2, say, ii REAL(dp) :: simtimestart_temp, simtimeend_temp REAL(dp) :: lat_estwayp, lon_estwayp, dista flag_add = 1 ti = time_day do while (flag_add.EQ.1) ti = ceiling(ti/2.0); !create new way points ti_day = ti/(24*60*60) temp1 = longind(kk) temp2 = longind(kk)+1 simtimestart_temp = t1(temp1) simtimeend_temp = t1(temp2) CALL built_vecsize(simtimestart_temp, simtimeend_temp, & ti_day, say) ALLOCATE(t_tempo(1:say)) CALL built_vec(simtimestart_temp,simtimeend_temp, & ti_day, t_tempo) ALLOCATE(lat_wayptemp(1:size(t_tempo),1) , & long_wayptemp(1:size(t_tempo),1) ) DO ii=1,size(t_tempo) ! Cal caresponding lat lon CALL get_latlon(iSat, t_tempo(ii), lat_estwayp, lon_estwayp) lat_wayptemp(ii,1) = lat_estwayp long_wayptemp(ii,1) = lon_estwayp END DO long_wayptemp = (long_wayptemp-180) flag_add = 0 dista = MAPL_DEG_PER_KM * get_distance( lat_wayptemp(1,1),long_wayptemp(1,1), & lat_wayptemp(2,1),long_wayptemp(2,1)) if (dista.GT.long_limits) then flag_add = 1 end if DEALLOCATE( lat_wayptemp,long_wayptemp ) DEALLOCATE (t_tempo) end do find_delta_t = ti END FUNCTION find_delta_t REAL(dp) FUNCTION get_fraction(t1) IMPLICIT NONE REAL(dp), INTENT(IN) :: t1 REAL(dp) :: fraction2day, temp3, temp4 if (t1.EQ.0.0) then fraction2day = 0.0 else if (t1.lt.1.0 .AND. t1.gt.0) then fraction2day = t1 else temp3 = floor(t1) temp4 = t1 if (temp3 == temp4) then fraction2day = temp3 ! either is ok else fraction2day = MOD(temp4,temp3) end if end if get_fraction = fraction2day END FUNCTION get_fraction SUBROUTINE find_LTGE_sc(myarray, num1, num2, myindex) IMPLICIT NONE INTEGER :: k, say INTEGER, INTENT(IN) :: num1, num2 REAL(dp), POINTER, DIMENSION(:) :: myindex !INTENT(INOUT) REAL(dp), DIMENSION(:), INTENT(IN) :: myarray !REAL(dp), DIMENSION(:), ALLOCATABLE :: myindextemp REAL(dp), DIMENSION(:), ALLOCATABLE :: myindextemp say = 0 ALLOCATE(myindextemp(1:size(myarray))) myindextemp = 0.0 do k=1, size(myarray) if ( (myarray(k).LT.num1).AND. & (myarray(k).GE.num2) ) then say = say + 1 myindextemp(say) = k end if end do if (say.EQ.0) then ALLOCATE(myindex(1:1)) myindex(1) = -999999 else ALLOCATE(myindex(1:say)) myindex(1:say) = myindextemp(1:say) end if DEALLOCATE(myindextemp) END SUBROUTINE find_LTGE_sc SUBROUTINE find_LEGT_sc(myarray, num1, num2, myindex) IMPLICIT NONE INTEGER :: k, say INTEGER, INTENT(IN) :: num1, num2 REAL(dp), POINTER, DIMENSION(:) :: myindex !INTENT(INOUT) REAL(dp), DIMENSION(:), INTENT(IN) :: myarray REAL(dp), DIMENSION(:), ALLOCATABLE :: myindextemp say = 0 ALLOCATE(myindextemp(1:size(myarray))) myindextemp = 0.0 do k=1, size(myarray) if ( (myarray(k).LE.num1).AND. & (myarray(k).GT.num2) ) then say = say + 1 myindextemp(say) = k end if end do if (say.EQ.0) then ALLOCATE(myindex(1:1)) myindex(1) = -999999.0 else ALLOCATE(myindex(1:say)) myindex = myindextemp(1:say) end if DEALLOCATE(myindextemp) END SUBROUTINE find_LEGT_sc SUBROUTINE ss(ar) real(dp), pointer, dimension(:) :: ar allocate(ar(1:2)) ar = (/ 3, 4 /) return END SUBROUTINE ss SUBROUTINE orbits(lat, lon, iSat, nymd, nhms) ! for now time has to be after IMPLICIT NONE INTEGER, INTENT(IN) :: iSat INTEGER, INTENT(IN) :: nymd, nhms REAL(dp), INTENT(OUT) :: lat,lon REAL(dp), DIMENSION(3,1) :: ECI_est, ECEF_est REAL(dp) :: alt REAL(dp) :: fraction2day, ntime_day CALL get_time(iSat, nymd, nhms, fraction2day, ntime_day) CALL get_estimateECI(iSat, ntime_day, ECI_est) CALL ECI2ECEF(iSat,fraction2day, ECI_est, ECEF_est) CALL ECEF2LLA(ECEF_est(1,1), ECEF_est(2,1), & ECEF_est(3,1), lat, lon, alt) END SUBROUTINE orbits END MODULE MAPL_NominalOrbitsMod