C Copyright(c) 2006, Space Science and Engineering Center, UW-Madison C Refer to "McIDAS Software Acquisition and Distribution Policies" C in the file mcidas/data/license.txt C *** $Id: nvxgeos.dlm,v 1.14 2024/08/19 21:32:27 tommyj Exp $ *** C**************************************************************** C AUTHOR : Oscar Alonso Lasheras C COPYRIGHT GMV S.A. 2006 C PROPERTY OF GMV S.A.; ALL RIGHTS RESERVED C C PROJECT : S3GICast C FILE NAME : nvxgeos.dlm C LANGUAGE : FORTRAN-77 C TYPE : Funcion auxliar para creacion de navegacion GEOS C DESCRIPTION : Navega los productos MPEF y OSIS, ademas de los productos C de los servidores FSDX y SAFN-HDF5 en McIDAS C C**************************************************************** C HISTORY C C DATE VERSION AUTHOR REASONS C MAY 2006 1.0 OOAL Creation C C**************************************************************** FUNCTION NVXINI(IFUNC,IPARMS) C C IPARMS: Navigation Block: C 1: GEOS C 2: LOFF C 3: COFF C 4: LFAC C 5: CFAC C 6: Proyection Longitude C 7: Base Resolution C 8: 1/2 pixel offset flag=1 (HIMAWARI) C 9: Sweep angle axis: Y (0, Meteosat), or X (1, GOES) DIMENSION IPARMS(*) REAL*4 PLAXH, PLON INTEGER LOFF,COFF,LFAC,CFAC,BRES,PFLAG CHARACTER*4 CLIT COMMON/GEOS/ITYPE COMMON/NVPARAM/LOFF,COFF,LFAC,CFAC,PLON,PLAXH,BRES,PFLAG PLAXH=0.0 ITYPE=0 NVXINI=0 IF (IFUNC.EQ.1) THEN IF (IPARMS(1).NE.LIT('GEOS')) THEN NVXINI=-1 RETURN ENDIF LOFF=IPARMS(2) COFF=IPARMS(3) LFAC=IPARMS(4) CFAC=IPARMS(5) PLON=IPARMS(6)/10. BRES=IPARMS(7) if( BRES.LE.0 ) BRES = 1 PFLAG=IPARMS(8) ELSE IF (IFUNC.EQ.2) THEN IF(INDEX(CLIT(IPARMS(1)),'XY').NE.0) ITYPE=1 IF(INDEX(CLIT(IPARMS(1)),'LL').NE.0) ITYPE=2 ENDIF RETURN END C ----------------------------------------------------------- FUNCTION NVXSAE(XLIN,XELE,XDUM,XLAT,XLON,Z) COMMON/GEOS/ITYPE INTEGER STATUS NVXSAE=0 STATUS = IMG_TO_LL(XLIN,XELE,XLAT,XLON) XLON=-XLON !Return longitudes in McIDAS representation (W+) IF(STATUS.NE.0) THEN NVXSAE=-1 RETURN ENDIF IF(ITYPE.EQ.1) THEN A=XLAT B=XLON CALL NLLXYZ(A,B,XLAT,XLON,Z) ENDIF RETURN END C ----------------------------------------------------------- FUNCTION NVXEAS(XLAT,XLON,Z,XLIN,XELE,XDUM) COMMON/GEOS/ITYPE NVXEAS=0 A=XLAT B=-XLON !Convert longitudes to standard representation (E+) IF(ITYPE.EQ.1) THEN CALL NXYZLL(XLAT,XLON,Z,A,B) B=-B ENDIF NVXEAS = LL_TO_IMG(A,B,XLIN,XELE) RETURN END C ----------------------------------------------------------- FUNCTION NVXOPT(IFUNC,XIN,XOUT) REAL*4 XIN(*),XOUT(*) CHARACTER*4 CLIT,CFUNC REAL*4 PLAXH, PLON INTEGER LOFF,COFF,LFAC,CFAC,BRES,PFLAG COMMON/NVPARAM/LOFF,COFF,LFAC,CFAC,PLON,PLAXH,BRES,PFLAG DATA LASDAY/-1/, LASTIM/-1/ NVXOPT=0 CFUNC=CLIT(IFUNC) IF(CFUNC.EQ.'SPOS') THEN XOUT(1)=0. XOUT(2)=-PLON ELSE IF(CFUNC(1:3) .EQ. 'HGT') THEN PLAXH = XIN(1) ELSE IF (CFUNC(1:3).eq.'ANG') THEN JDAY = IROUND(XIN(1)) JTIME = M0ITIME(XIN(2)) FLAT = XIN(3) FLON = XIN(4) IF (JDAY.NE.LASDAY.OR.JTIME.NE.LASTIM) THEN CALL SOLARP(JDAY,JTIME,GHA,DEC,XLAT,XLON) LASDAY = JDAY LASTIM = JTIME ENDIF CALL GEOSANG(JDAY,JTIME,FLAT,FLON,GHA,DEC, & XOUT(1),XOUT(2),XOUT(3)) NVXOPT = 0 ELSE NVXOPT=-1 ENDIF RETURN END C ----------------------------------------------------------- INTEGER FUNCTION LL_TO_IMG (XLAT,XLON,RLIN,RELE) C C S/R gives line and element in GEOS projectios for C a specified latitude and longitude C The subroutine uses the navigation parameters COFF, LOFF, C CFAC, LFAC & PLON provided in the GEOS navigation block C C Inputs: C XLAT,XLON (REAL) : latitude and longitude of selected point C latitude North is positive, C longitude East is positive C ================ C C Outputs: C XLIN, XELE (REAL) : line and element number C Output is -999. if specified XLAT/XLON is C not within image frame C IMPLICIT NONE INTEGER RETURN_CODE REAL RLIN, RELE REAL XLAT, XLON, XLIN, XELE REAL C_LAT, COSC_LAT, RN, R1, R2, R3, RL REAL X,Y REAL LAT,LON,SPLON REAL AD2, BD, CD, DELTA2, HALFSOM, R_EQ2, R_POL2 REAL PI PARAMETER (PI = 3.141592654) REAL*4 PLAXH, PLON INTEGER LOFF,COFF,LFAC,CFAC,BRES,PFLAG COMMON/NVPARAM/LOFF,COFF,LFAC,CFAC,PLON,PLAXH,BRES,PFLAG RETURN_CODE = -1 C --- Coordinates are computed accroding EUMETSAT's LRIT/HRIT Global Spec C --- Doc No: CGMS 03 c --- Coordinates are converted to Radians LAT = XLAT*PI/180. LON = XLON*PI/180.0 SPLON = PLON * PI/180.0 c --- Intermediate data C_LAT=ATAN(0.993243*TAN(LAT)) COSC_LAT=COS(C_LAT) R_POL2= (6356.7523+PLAXH)*(6356.7523+PLAXH) R_EQ2 = (6378.1370+PLAXH)*(6378.1370+PLAXH) RL=(6356.7523+PLAXH)/ * (SQRT(1-((R_EQ2-R_POL2)/R_EQ2)*COSC_LAT*COSC_LAT)) R1=42164.537-RL*COSC_LAT*COS(LON-SPLON) R2=-RL*COSC_LAT*SIN(LON-SPLON) R3=RL*SIN(C_LAT) RN=SQRT(R1*R1+R2*R2+R3*R3) c --- Compute variables useful to check if pixel is visible AD2 = R1*R1 + R2*R2 + R3*R3*R_EQ2 / R_POL2 BD = 42164.537*R1 CD = 42164.537*42164.537 - R_EQ2 DELTA2 = BD*BD-AD2*CD HALFSOM = BD*RN/AD2 IF ((DELTA2 .GE. 0.) .AND. (RN .LE. HALFSOM)) THEN c ------- Intermediate coordinates X = ATAN(-R2/R1) Y = ASIN(-R3/RN) X = X * 180./PI Y = Y * 180./PI XELE = COFF/10. + X / 2**16 * CFAC/10. XLIN = LOFF/10. + Y / 2**16 * LFAC/10. C --- Check is 1/2 pixel offset flag is set if (PFLAG .eq. 1) then RLIN = (XLIN*BRES)- ( (BRES-1)/2. ) RELE = (XELE*BRES)- ( (BRES-1)/2. ) else RLIN = (XLIN*BRES)- ( (BRES-1) ) RELE = (XELE*BRES)- ( (BRES-1) ) endif RETURN_CODE = 0 ELSE RLIN=-999 RELE=-999 END IF LL_TO_IMG = RETURN_CODE RETURN END INTEGER FUNCTION IMG_TO_LL (RLIN,RELE,XLAT,XLON) C C S/R gives latitude and longitude for C a specified line and element in GEOS projection C The subroutine uses the navigation parameters COFF, LOFF, C CFAC, LFAC & PLON provided in the GEOS navigation block C C Inputs: C RLIN, RELE (REAL) : line and element number C C Outputs: C XLAT,XLON (REAL) : latitude and longitude of selected point C latitude North is positive, C longitude East is positive C output is -999. if line/element is off the disk C IMPLICIT NONE REAL RLIN, RELE REAL XLAT, XLON, XLIN, XELE REAL X,Y REAL S1, S2, S3, SXY, SN, SD, SDD REAL AUX, AUX2 REAL COSX, COSY, SINX, SINY REAL R_EQ2, CD REAL PI PARAMETER ( PI = 3.14159265 ) REAL*4 PLAXH, PLON INTEGER LOFF,COFF,LFAC,CFAC,BRES,PFLAG COMMON/NVPARAM/LOFF,COFF,LFAC,CFAC,PLON,PLAXH,BRES,PFLAG C --- Coordinates are computed accroding EUMETSAT's LRIT/HRIT Global Spec C --- Doc No: CGMS 03 C --- USE BRES TO ADJUST THE COORDINATES if(PFLAG .eq. 1) then XLIN = (RLIN + ( (BRES-1)/2.) ) / BRES XELE = (RELE + ( (BRES-1)/2.) ) / BRES else XLIN = (RLIN + ( (BRES-1)) ) / BRES XELE = (RELE + ( (BRES-1)) ) / BRES endif c --- Intermediate coordinates X = (XELE - COFF/10.) * 2**16 / (CFAC/10.) Y = (XLIN - LOFF/10.) * 2**16 / (LFAC/10.) X = X * PI/180. Y = Y * PI/180. c --- Intermediate data COSX=COS(X) COSY=COS(Y) SINX=SIN(X) SINY=SIN(Y) R_EQ2 = (6378.1370+PLAXH)*(6378.1370+PLAXH) CD = 42164.537*42164.537 - R_EQ2 c --- CD=1737121856 when PLAXH=0.0 AUX=42164.537*COSX*COSY AUX2=COSY*COSY+1.006803*SINY*SINY SDD=AUX*AUX-AUX2*CD IF(SDD.LT.0) THEN XLAT=-999 XLON=-999 IMG_TO_LL = -1 RETURN END IF SD=SQRT(SDD) SN=(AUX-SD)/AUX2 S1=42164.537 - SN*COSX*COSY S2=SN*SINX*COSY S3= -SN*SINY SXY=SQRT(S1*S1+S2*S2) c --- Computation XLON = ATAN(S2/S1) XLON = XLON * 180./PI + PLON XLAT = ATAN(1.006803*S3/SXY)* 180./PI c --- Longitudes in [-180,180] IF(XLON .GT. 180) XLON = XLON - 360. IF(XLON .LT. -180) XLON = XLON + 360. IMG_TO_LL = 0 RETURN END SUBROUTINE GEOSANG(JDAY,JTIME,XLAT,XLON,GHA,DEC,SATANG,SUNANG,RELA &NG) C C $ SUBROUTINE ANGLES(JDAY,JTIME,XLAT,XLON,GHA,DEC,SATANG,SUNANG,RELANG) C $ ANGLES - computes zenith angles of sun and satellite and relative C $ azimuth angle (DAS) C $ INPUT: C $ JDAY = (I) picture day (YYDDD) C $ JTIME = (I) picture start time C $ XLAT = (R) latitude of point C $ XLON = (R) longitude of point C $ GHA = (R) Greenwich hour angle of sun C $ DEC = (R) declination of sun C $ OUTPUT: C $ SATANG = (R) zenith angle of satellite C $ SUNANG = (R) zenith angle of sun C $ RELANG = (R) relative angle C $$ ANGLES = COMPUTATION, NAVIGATION C ANGLES MOSHER 1074 WINLIB ZENITH ANGLES TO SAT,SUN,AND REL AZIMUTH AN C C INCLUDE 'ELCONS.INC' C=========================== DIELCONS ============================= C $ (JR) C $ THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE C $ C $ DESCRIPTION: C $ C $ C $$ DIELCONS = INCLUDE DOUBLE PRECISION PI PARAMETER (PI=3.141592653589793D0) DOUBLE PRECISION DEG PARAMETER (DEG=180.D0/PI) DOUBLE PRECISION RAD PARAMETER (RAD=PI/180.D0) C DEGREES TO RADIANS CONVERSION PI/180 DOUBLE PRECISION NOMORB PARAMETER (NOMORB=42164.537D0) C NOMINAL RADIAL DISTANCE OF SATELLITE (KM) DOUBLE PRECISION AE PARAMETER (AE=6378.1370D0) C EARTH EQUATORIAL RADIUS (KM) DOUBLE PRECISION FER PARAMETER (FER=1.D0-(6356.7523D0/AE)) C EARTH FLATTENING COEFFICIENT = 1-(BE/AE) REAL AEBE2 PARAMETER (AEBE2=1.D0/(1.D0-FER)**2) REAL AEBE3 PARAMETER (AEBE3=AEBE2-1.) REAL AEBE4 PARAMETER (AEBE4=(1.D0-FER)**4-1.) C COMMON VARIABLES C DOUBLE PRECISION XS(3) C NORMALIZED S/C POSITION IN ECEF COORDINATES REAL*4 XPLAT, XPLON, SNLT, CSLT, CSLN, SNLN REAL*4 PLAXH, PLON INTEGER LOFF,COFF,LFAC,CFAC,BRES,PFLAG COMMON/NVPARAM/LOFF,COFF,LFAC,CFAC,PLON,PLAXH,BRES,PFLAG DATA IDAY/0/ DATA R/6371.221/ RDPDG=PI/180.0 C XPLAT = 0.0 XPLON = PLON SNLT=SIN(XPLAT*RDPDG) CSLT=COS(XPLAT*RDPDG) CSLN=COS(XPLON*RDPDG) SNLN=SIN(XPLON*RDPDG) XS(1)=NOMORB*CSLT*CSLN/AE XS(2)=NOMORB*CSLT*SNLN/AE XS(3)=NOMORB*SNLT/AE IF(IDAY.EQ.JDAY)GO TO 1 IDAY=JDAY INORB=0 1 PICTIM=FTIME(JTIME) C C DETERMINE SATELLITE POSITION C XSAT = XS(1) * AE YSAT = XS(2) * AE ZSAT = XS(3) * AE HEIGHT=SQRT(XSAT**2+YSAT**2+ZSAT**2) YLAT=RDPDG*XLAT YLAT=GEOLAT(YLAT,1) YLON=RDPDG*XLON SLAT=SIN(YLAT) CLAT=COS(YLAT) SLON=SIN(YLON) CLON=COS(YLON) XSAM=R*CLAT*CLON YSAM=R*CLAT*SLON ZSAM=R*SLAT C C DETERMINE ZENITH ANGLE OF SUN C SNLG=-PICTIM*PI/12.0-RDPDG*GHA SNDC=RDPDG*DEC COSDEC=COS(SNDC) US=COS(SNLG)*COSDEC VS=SIN(SNLG)*COSDEC WS=SIN(SNDC) SUNANG=ACOS((US*XSAM+VS*YSAM+WS*ZSAM)/R)/RDPDG C C DETERMINE ZENITH ANGLE OF SATELLITE C XVEC=XSAT-XSAM YVEC=YSAT-YSAM ZVEC=ZSAT-ZSAM XFACT=SQRT(XVEC**2+YVEC**2+ZVEC**2) SATANG=ACOS((XVEC*XSAM+YVEC*YSAM+ZVEC*ZSAM)/(R*XFACT))/RDPDG C C DETERMINE RELATIVE ANGLE C X1=CLAT*CLON Y1=CLAT*SLON Z1=SLAT X2=SLON Y2=-CLON X3=-SLAT*CLON Y3=-SLAT*SLON Z3=CLAT XC1=US-X1 YC1=VS-Y1 ZC1=WS-Z1 XC2=XSAT/HEIGHT-X1 YC2=YSAT/HEIGHT-Y1 ZC2=ZSAT/HEIGHT-Z1 XAN1=XC1*X3+YC1*Y3+ZC1*Z3 XAN2=XC2*X3+YC2*Y3+ZC2*Z3 YAN1=XC1*X2+YC1*Y2 YAN2=XC2*X2+YC2*Y2 XAN3=XAN1*XAN2+YAN1*YAN2 YAN3=-YAN1*XAN2+XAN1*YAN2 RELANG=ATAN2(YAN3,XAN3)/RDPDG RELANG=ABS(RELANG) RETURN END