MODULE onebyte_1band_rect_read ! ! This module reads the single banded mcidas area files in 1 byte ! rectalinear projection images ! ! Adapted from code that J. Kossin wrote while a post doc at CIRA (2002), now ! modified to read rectilinear projection instead of mercator ! ! code will return latitude,longitude,temperature of each pixel as well ! as date and time information in global variables. ! ! Written by John Knaff ! NOAA/NESDIS/StAR/RAMMB ! Fort Collin, CO ! ! Last Modified 01/10/2018 ! ! created onebyte_1band_rect_read module by modifying onebyte_1band_merc read ! module. Required using the SSEC software to provide naviagation for ! the rectilinear projection. ! NOTE: the offset to the Navigation directory had to be shortened by one ! word. (JAK) - 01/09/2018 ! NOTE: I did not change the global variable names from the mercator reader. ! This may cause problems if these two pieces of code are run together ! as they use the same global variable names. ! Documentation and clean-up (JAK) -01/10/2018 ! !********************************************************************* ! G L O B A L V A R I A B L E S: ! ! Description of Global Variables (those loaded in memory) ! ! TIME: ! fyear - fractional year (DOUBLE PRECISION) ! dtgj - date time group /year+year-day(julian date) (YYYYjjj) ! (INTEGER) ! dtgg - date time group /year+mo+day (YYYYmmdd) (INTEGER) ! time - time HHMMSS (INTEGER) ! ! IMAGE: ! nliny - number of lines (y/vertical) (INTEGER) ! nelex - number of elements per line (x/horizontal) (INTEGER) ! bt(,) - brightness temperature [K], array (nelex,nliny), (REAL) ! ! LOCATION: ! lat() - Latitude in fractional degrees ! lon() - Longitude in fractional degrees (0 to 360 convention) ! ! SATELLITE: ! isat - ssec satellite sensor number (a limited table ! is used in this program to map to ascii descriptions ! (INTEGER) ! sat - ascii description of the platform (20 characters) ! ich - channel (satellite specific) ! !********************************************************************* IMPLICIT NONE ! ! Global variables INTEGER :: dtgj ! julian date YYJJJ INTEGER :: dtgg ! gregorian date YYMMDD INTEGER :: time ! HHMMSS INTEGER :: nliny,nelex INTEGER :: isat ! ssec satellite sensor number INTEGER :: ich ! channel DOUBLE PRECISION :: fyear ! fractional year YYYY.(jday/365) REAL, ALLOCATABLE :: lat(:),lon(:) ! latitude (degree N) longitude (0-360E) REAL, ALLOCATABLE :: bt (:,:) ! brightness temperature [K] REAL:: STDLAT ! standard latitude for mercator projection [degrees] CHARACTER (LEN=20) :: sat ! satellite name associated with ss number CONTAINS LOGICAL FUNCTION rect_get(cfn) INTEGER :: i,nav_beg,nav_end,icount INTEGER :: luarea, indwrd INTEGER :: irl INTEGER :: ival,nchar,ielem,iline INTEGER :: inv,idi,m,n,iwest,icord INTEGER, DIMENSION(2) :: LL INTEGER, DIMENSION(64) :: areadir INTEGER (KIND=1) , ALLOCATABLE :: imgblk(:) INTEGER, ALLOCATABLE :: navdir(:) ! JFD added 12/30/2008 INTEGER, ALLOCATABLE :: caldir(:) INTEGER :: icb INTEGER:: iajday,iayear,iamon,iaday REAL, PARAMETER :: xmsng=-999.99 REAL :: tk, AAA REAL:: ul_y,ul_x,yl,xl,rlon,rlat REAL:: xdumy,xla,xlo,zed CHARACTER (LEN=*),INTENT(IN) :: cfn ! CHARACTER (LEN=4) :: CLIT LOGICAL :: lfn, swap ! Executable statements rect_get = .false. IF (ALLOCATED (lat) ) DEALLOCATE (lat ) IF (ALLOCATED (lon) ) DEALLOCATE (lon ) IF (ALLOCATED (navdir)) DEALLOCATE (navdir) IF (ALLOCATED (imgblk)) DEALLOCATE (imgblk) IF (ALLOCATED (bt) ) DEALLOCATE (bt ) luarea=50 INQUIRE (FILE=trim(cfn),exist=lfn) IF (lfn) THEN OPEN (UNIT=luarea,FILE=trim(cfn),ACCESS='direct',& STATUS='old',recl=256,err=900) ! determine endian of the machine CALL detend(indwrd) READ (luarea,rec=1)areadir CLOSE (luarea) ! swap bytes if necessary (match endian-ness between host and data) ! swap=.false. IF(areadir(2).EQ.67108864)THEN swap=.true. DO i=1,64 CALL swap4(areadir(i)) END DO END IF ! print*,'SWAP =',swap ! query file parameters ! IF (areadir(14) .NE. 1) THEN PRINT*, cfn PRINT*,'...is a multibanded image' GOTO 910 END IF IF (areadir(1).NE.0.OR.areadir(2).NE.4) THEN PRINT*, cfn PRINT*, '...is not a McIDAS area. Skipping file.' GOTO 910 END IF IF (areadir(11).NE.1) THEN PRINT*, cfn PRINT*, '...is not composed of 1 byte data and is not' PRINT*, 'from the expected IR archive. Skipping file.' GOTO 910 END IF IF(areadir(12).NE.areadir(13)) THEN PRINT*, 'In file ',cfn PRINT*, '...x-resolution is not equal to y-resolution.' PRINT*, 'Skipping file.' GOTO 910 END IF IF (MOD(areadir(35),4).NE.0) THEN PRINT*, 'In file ',cfn PRINT*, '...data were found between the area and' PRINT*, 'navigation blocks which are not composed' PRINT*, 'of words. Skipping file.' GOTO 910 END IF isat=areadir(3) IF (areadir(19).EQ. 1) ich=1 IF (areadir(19).EQ. 2) ich=2 IF (areadir(19).EQ. 4) ich=3 IF (areadir(19).EQ. 8) ich=4 IF (areadir(19).EQ. 16) ich=5 IF (areadir(19).EQ. 32) ich=6 IF (areadir(19).EQ. 64) ich=7 IF (areadir(19).EQ. 128) ich=8 IF (areadir(19).EQ. 256) ich=9 IF (areadir(19).EQ. 512) ich=10 IF (areadir(19).EQ. 1024) ich=11 IF (areadir(19).EQ. 2048) ich=12 IF (areadir(19).EQ. 4096) ich=13 IF (areadir(19).EQ. 8192) ich=14 IF (areadir(19).EQ. 16384) ich=15 IF (areadir(19).EQ. 32768) ich=16 sat='NOT in Table' IF (isat == 56) sat='Meteosat-5' IF (isat == 58) sat='Meteosat-7' IF (isat == 51) sat='MSG-1' IF (isat == 52) sat='MSG-2' IF (isat == 53) sat='MSG-3' IF (isat == 29) sat='GOES-5 IR' IF (isat == 31) sat='GOES-6 IR' IF (isat == 33) sat='GOES-7 IR' IF (isat == 70) sat='GOES-8 imager' IF (isat == 72) sat='GOES-9 imager' IF (isat == 74) sat='GOES-10 imager' IF (isat == 76) sat='GOES-11 imager' IF (isat == 78) sat='GOES-12 imager' IF (isat == 180) sat='GOES-13 imager' IF (isat == 182) sat='GOES-14 imager' IF (isat == 184) sat='GOES-15 imager' IF (isat == 186) sat='GOES-16 imager' IF (isat == 188) sat='GOES-17 imager' IF (isat == 85) sat='MTSAT-2' IF (isat == 84) sat='MTSAT-1R' IF (isat == 83) sat='GMS-5' IF (isat == 82) sat='GMS-4' ! ******************************************************************* ! Calculate total data file length (in bytes) to end of image block ! ! Note the AREA block is alway 64...defined above !irl=areadir(34)+areadir(9)*areadir(10) ! size of the dataset !inv= (areadir(34) - areadir(63))/4 - 64 ! size of the nav block !idi=irl - 64*4 -inv*4 ! size of the data block ! 12/30/2008 - Apparently, IMGREMAP puts a calibration ! block into the AREA ! file. The entries may all be -999, but its presence caused some ! problems. Rework the reading of the area file assuming ! the following ! structure for the file: ! Directory block - 256 bytes (64 words) ! Navigation block - areadir(63)-256 bytes ! Calibration block - areadir(34)-areadir(63) bytes ! Supplemental block - assumed size 0 ! Data block - areadir(9)*areadir(10)*4 bytes starting at ! areadir(34)+1 bytes !k*** !k areadir(34) for the isats version of the files -JAK 01/10/2018 !k*** ! ! So irl above is correct, but inv and idi need to be recomputed, ! and icb ! (size of the calibration block) added. !******************************************************************** IF (areadir(63) .EQ. 0) THEN ! No Calibration Block irl=areadir(34)+areadir(9)*areadir(10) ! size of the dataset inv= (areadir(34) - areadir(63))/4 - 64 ! size of the nav block idi=irl - 64*4 -inv*4 ! size of the data block ! Allocate array sizes ALLOCATE (navdir(inv)) ALLOCATE (imgblk(idi)) ALLOCATE (bt(areadir(10),areadir(9))) ALLOCATE (lat(areadir(9))) ALLOCATE (lon(areadir(10))) ! Reopen file OPEN (UNIT=luarea,FILE=cfn,ACCESS='direct',STATUS='old',recl=irl, & err=900) READ(luarea,rec=1)areadir,navdir,imgblk !k this is what is traditionally has worked !k nav_beg=areadir(35)/4+1 - 64 !k This is required to read the navigation directory correctly nav_beg=areadir(35)/4 - 64 nav_end=areadir(34)/4 - 64 ! print*,'areadir(35),/4 =',areadir(35),areadir(35)/4 ELSE ! there is a calibration block!! irl=areadir(34)+areadir(9)*areadir(10) ! size of the dataset inv= (areadir(63) - 256)/4 ! size of the nav block icb=(areadir(34)-areadir(63))/4 ! size of the cal block idi=areadir(9)*areadir(10) ! size of the data block ! Allocate array sizes ALLOCATE (navdir(inv)) ALLOCATE (caldir(icb)) ALLOCATE (imgblk(idi)) ALLOCATE (bt(areadir(10),areadir(9))) ALLOCATE (lat(areadir(9))) ALLOCATE (lon(areadir(10))) ! Reopen file OPEN (UNIT=luarea,FILE=cfn,ACCESS='direct',STATUS='old',recl=irl, & err=900) READ(luarea,rec=1)areadir,navdir,caldir,imgblk !k this is what is traditionally has worked !k nav_beg=areadir(35)/4+1 - 64 !k This is required to read the navigation directory correctly !k since isats version does not have calibaration (NOT TESTED) nav_beg=areadir(35)/4 - 64 nav_end=areadir(63)/4 - 64 END IF ! swap bytes in navigation block if necessary ! IF (swap) THEN DO i=nav_beg,nav_end CALL swap4(navdir(i)) END DO DO i=1,64 CALL swap4(areadir(i)) END DO END IF ! Image navigation: pixels, distances (km), lat/lon (deg) ! ! initialize navigation routines to accept a (lat,lon) pair and ! return (x,y) distance (km) from the top left corner of the image ! STDLAT=float(navdir(4))/10000.0 CALL nvxini(1,navdir,nav_beg,nav_end,cfn,indwrd) LL(1)=lit('LL ') CALL nvxini(2,LL,1,2,cfn,indwrd) ! return (lat,lon) pair given an (x,y) distance (km) ! CALL nvxini(1,navdir,nav_beg,nav_end,cfn,indwrd) icount=0 AAA=0.0 DO iline=1,areadir(9) DO ielem=1,areadir(10) icount=icount+1 CALL byttowrd(imgblk(icount),indwrd,ival) ! convert count value ival to degrees Kelvin; set to zero if ! count is zero ! IF (float(ival).EQ.0.)THEN tk=xmsng ELSE IF (ival.le.175)THEN tk=330.-0.5*float(ival) ELSE tk=418.-float(ival) END IF bt(ielem,iline)=tk AAA=AAA+tk END DO END DO ! Assign lines and elements (GLOBAL) ! nliny = areadir( 9) nelex = areadir(10) ! Create arrays for the latitude and longitude arrays (GLOBAL) ! ul_y=float(areadir(6)) ul_x=float(areadir(7)) DO iline=1,areadir(9) yl=ul_y + float(iline-1)*float(areadir(12)) CALL nvxsae(yl,ul_x,xdumy,xla,xlo,zed,rlat,rlon) LAT(iline)=rlat END DO DO ielem=1,areadir(10) xl=ul_x + float(ielem-1)*float(areadir(13)) CALL nvxsae(ul_y,xl,xdumy,xla,xlo,zed,rlat,rlon) rlon=-1.0*rlon IF (rlon .LT. 0.0)rlon=360.0 + rlon LON(ielem)=rlon END DO ! Calculate the various time variables (GLOBAL) ! dtgj=areadir(4)+1900000 iayear=dtgj/1000 iajday=dtgj-iayear*1000 CALL jdayi(iajday,iayear,iamon,iaday) dtgg=iayear*10000+iamon*100+iaday time=areadir(5) fyear=yearfrac(iayear,iajday,time) rect_get=.true. ELSE PRINT*, 'Failed to open ',TRIM (cfn) END IF RETURN 900 CONTINUE PRINT *, 'Error while opening file named ',TRIM(cfn) RETURN 910 CONTINUE RETURN END FUNCTION rect_get !*$ Name: !*$ clit - Returns the four bytes of an integer as a !*$ character string. !*$ !*$ Interface: !*$ character*4 function !*$ clit(integer l) !*$ !*$ Input: !*$ l - Integer that is copied to a character variable. !*$ !*$ Input and Output: !*$ none !*$ !*$ Output: !*$ none !*$ !*$ Return values: !*$ The four bytes of the input integer. !*$ !*$ Remarks: !*$ This function is designed to move the four bytes of an integer to a !*$ character variable. It enables programmers to avoid conversions !*$ implicit in replacement statements in which variables of different !*$ types are on opposite sides of the equal sign. !*$ !*$ Categories: !*$ converter !************************************************************************* ! Version used here is the same idea, but less complicated and requiring ! no additional/external calls (uses equivalence to make conversion) ! A old trick but a good trick ! John Knaff, NOAA/NESDIS !************************************************************************* CHARACTER (LEN=4) FUNCTION CLIT(L) IMPLICIT NONE INTEGER :: L ! --- local variables CHARACTER (LEN=4):: C INTEGER :: L1 EQUIVALENCE(L1,C) L1=L CLIT=C RETURN END FUNCTION CLIT SUBROUTINE swap4(in) ! ! swaps the byte order of a 4-byte integer (big endian<->little endian) ! IMPLICIT NONE INTEGER :: in,out INTEGER (KIND=1) :: bout(4),b ! BYTE :: bout(4),b ! some compilers/f77 EQUIVALENCE(out,bout) out=in b=bout(1) bout(1)=bout(4) bout(4)=b b=bout(2) bout(2)=bout(3) bout(3)=b in=out RETURN END SUBROUTINE swap4 SUBROUTINE detend(idum) ! ! Determines endian-ness of host machine. Returns index 1 for Little ! Endian, index 4 for Big Endian. The index specifies least significant ! byte position of 4-byte word. ! IMPLICIT NONE INTEGER:: i,j,idum INTEGER (KIND=1):: word(4) ! BYTE:: word(4) - some compilers/f77 EQUIVALENCE(word(1),i) DO j=1,4 word(j)=0 END DO i=1 IF(word(1).ne.0)then idum=1 ELSE idum=4 END IF RETURN END SUBROUTINE detend SUBROUTINE byttowrd(ibyt,ind,iwrd) ! ! Puts a 1-byte unsigned integer, with value between 0 and 255, into a ! 4-byte integer. Index ind contains information regarding the ! endian-ness of the host machine. The same thing might be accomplished ! by adding 256 to negative values of ibyt and putting the result into ! an integer, but I'm unsure about the portability, since it might ! require that the host machine represent negative numbers with two's ! complement binary sequences. This routine handles the steps ! explicitly. !************************************************************** IMPLICIT NONE INTEGER (KIND=1) :: ibyt,jbyt(4) ! BYTE :: ibyt,jbyt(4) ! CHARACTER (LEN=1) :: ibyt,jbyt(4) INTEGER:: ind,iwrd,jwrd EQUIVALENCE (jwrd,jbyt(1)) jwrd=0 jbyt(ind)=ibyt iwrd=jwrd ! iwrd=int(jbyt(1),kind=2) RETURN END SUBROUTINE byttowrd INTEGER*4 FUNCTION LIT(C) IMPLICIT NONE CHARACTER*4 C ! --- local variables INTEGER*4 L CHARACTER*4 C1 EQUIVALENCE (C1,L) C1=C LIT=L RETURN END FUNCTION LIT ! ! THIS IS SSEC PROPRIETARY SOFTWARE - ITS USE IS RESTRICTED. REAL FUNCTION FLALO(M) ! *** McIDAS Revision History *** ! 1 FLALO.FOR 19-Mar-90,21:49:24,`SSEC' PC-McIDAS ver 5.00 ! 2 FLALO.FOR 25-Sep-90,7:33:56,`SMG' First Release into COMmon ! 3 FLALO.FOR 1-Apr-94,1:11:14,`BARRYR' Add proprietary statement ! 4 FLALO.FOR 2-May-94,16:38:40,`USER' Released ! *** McIDAS Revision History *** ! $ FUNCTION FLALO(M) (BL) ! $ CONVERT PACKED INTEGER (SIGN DDD MM SS) LATITUDE-LONGITUDE TO REAL*4 ! $ M = (I) INPUT PACKED INTEGER (SIGN DDD MM SS) LATITUDE-LONGITUDE ! $$ FLALO = CONVERT, INTEGER, LATITUDE, LONGITUDE, REAL ! IMPLICIT NONE INTEGER :: m,n REAL :: x IF(M.LT.0)GO TO 1 N=M X=1.0 GO TO 2 1 N=-M X=-1.0 2 FLALO=FLOAT(N/10000)+FLOAT(MOD(N/100,100))/60.0+FLOAT(MOD(N,100))/ & 13600.0 FLALO=X*FLALO RETURN END FUNCTION FLALO ! ! THIS IS SSEC PROPRIETARY SOFTWARE - ITS USE IS RESTRICTED. SUBROUTINE LLCART(XLAT,XLON,X,Y,Z) ! *** McIDAS Revision History *** ! 1 LLCART.FOR 19-Mar-90,21:50:54,`SSEC' PC-McIDAS ver 5.00 ! 2 LLCART.FOR 25-Sep-90,7:33:58,`SMG' First Release into COMmon ! 3 LLCART.FOR 1-Apr-94,1:22:16,`BARRYR' Add proprietary statement ! 4 LLCART.FOR 2-May-94,16:57:16,`USER' Released ! *** McIDAS Revision History *** ! $ SUBROUTINE LLCART(XLAT, XLON, X,Y,Z) (DAS) ! $ CONVERT LAT, LON TO CARTESIAN CENTERED COORDS (X, Y, Z). ! $ XLAT = (R) INPUT LATITUDE IN DEGREES, NORTH + ! $ XLON = (R) INPUT LONGITUDE IN DEGREES, WEST + ! $ X, Y, Z = (R) OUTPUT COORDS IN SYSTEM WITH ORIGIN AT CENTER OF ! $ PLANET. POS X-AXIS PIERCES EQUATOR AT LON 0 DEG, POS Y-AXIS PIERCES ! $ EQUATOR AT LON 90 DEG, & POS Z-AXIS INTERSECTS THE NORTH POLE. ! $ (IN KM). ! $$ LLCART = COORDINATES,LATITUDE,LONGITUDE,NAVIGATION IMPLICIT REAL*8 (D) REAL*8 ASQ,BSQ,AB,RDPDG,ECC,ECCSQR REAL :: XLAT,XLON,RAD(*),x,y,z,tlat(*),tcos,ylat,ylon REAL :: snlt,cslt,csln,snln,tnlt,r,a INTEGER :: KWEST,KCORD,IWEST,ICORD DATA ASQ/40683833.48/,BSQ/40410330.18/,AB/40546851.22/ DATA ECC/.081992D0/,ECCSQR/6.72265D-3/ DATA RDPDG/1.74532925199D-02/ DATA KWEST/-1/ DATA KCORD/0/ GO TO 55 ENTRY LLOPT(DRAD,DECC,IWEST,ICORD) ! $ SUBROUTINE LLOPT(DRAD,DECC,IWEST,ICORD) ! $ LLOPT IS USED TO INPUT RADIUS & ECCENTRICITY OTHER THAN EARTH ! $ AND DETERMINE + LONGITUDE CONVETION & PLANETOCENTRIC/DETIC ! $ DRAD = (R*8) INPUT EQUATORIAL RADIUS ! $ DECC = (R*8) INPUT ECCENTRICITY ! $ IWEST = (I) INPUT >= 0 WEST POSITIVE, < 0 WEST NEGATIVE ! $ ICORD = (I) INPUT >= 0 PLANETODETIC, < 0 PLANETOCENTRIC ! $$ LLOPT = COORDINATES,LATITUDE,LONGITUDE,NAVIGATION ASQ=DRAD*DRAD ECC=DECC ECCSQR=ECC*ECC DPOLE=SQRT(ASQ*(1.D0-ECCSQR)) BSQ=DPOLE*DPOLE AB=DRAD*DPOLE IF(IWEST.LT.0) KWEST=1 IF(ICORD.LT.0) KCORD=-1 RETURN ENTRY LLOBL(TLAT,RAD) ! $ SUBROUTINE LLOBL(TLAT,RAD) (DAS) ! $ CALCULATE RADIUS AT LATITUDE TLAT ! $ TLAT = (R) INPUT LATITUDE ! $ RAD = (R) OUTPUT RADIUS IN KM ! $$ LLOBL = NAVIGATION,COORDINATES TCOS=COS(TLAT(1)*RDPDG) DDRAD=((1.D0-ECCSQR)*ASQ)/(1.D0-ECCSQR*TCOS*TCOS) RAD(1)=SQRT(DDRAD) RETURN 55 CONTINUE YLAT=RDPDG*XLAT !-----CONVERT TO GEOCENTRIC (SPHERICAL) LATITUDE IF KCORD >= 0 !CC YLAT=GEOLAT(YLAT,1) IF(KCORD.GE.0) YLAT=ATAN2(BSQ*SIN(YLAT),ASQ*COS(YLAT)) YLON=KWEST*RDPDG*XLON SNLT=SIN(YLAT) CSLT=COS(YLAT) CSLN=COS(YLON) SNLN=SIN(YLON) TNLT=(SNLT/CSLT)**2 R=AB*SQRT((1.0+TNLT)/(BSQ+ASQ*TNLT)) X=R*CSLT*CSLN Y=R*CSLT*SNLN Z=R*SNLT RETURN ENTRY CARTLL(X,Y,Z,XLAT,XLON) ! $ SUBROUTINE CARTLL(X, Y, Z, XLAT, XLON) (DALY 1978) ! $ CONVERT CARTESIAN CENTERED COORD (X, Y, Z) TO LAT, LON. ! $ X, Y, Z = (R) OUTPUT COORDS IN SYSTEM WITH ORIGIN AT CENTER OF ! $ PLANET. POS X-AXIS PIERCES EQUATOR AT LON 0 DEG, POS Y-AXIS PIERCES ! $ EQUATOR AT LON 90 DEG, & POS Z-AXIS INTERSECTS THE NORTH POLE. ! $ IN KM. ! $ XLAT = (R) INPUT LATITUDE IN DEGREES, NORTH + ! $ XLON = (R) INPUT LONGITUDE IN DEGREES, WEST + ! $$ CARTLL = COORDINATES,LATITUDE,LONGITUDE,NAVIGATION ! ! XLAT=100.0 XLON=200.0 IF(X.EQ.0..AND.Y.EQ.0..AND.Z.EQ.0.) GO TO 90 A=ATAN(Z/SQRT(X*X+Y*Y)) !-----CONVERT TO GEODETIC LATITUDE IF KCORD > 0 !CC XLAT=GEOLAT(ATAN(Z/SQRT(X*X+Y*Y)),2)/RDPDG IF(KCORD.GE.0) THEN XLAT=ATAN2(ASQ*SIN(A),BSQ*COS(A))/RDPDG ELSE XLAT=A/RDPDG ENDIF XLON=KWEST*ATAN2(Y,X)/RDPDG 90 RETURN END SUBROUTINE LLCART !k Rectilinear routines from Mcidas2017.2..... !k functions changed to subroutines - JAK SUBROUTINE NVXINI(IFUNC,IPARMS,nbeg,nend,fndum,indwrd) ! Copyright(c) 1998, Space Science and Engineering Center, UW-Madison ! Refer to "McIDAS Software Acquisition and Distribution Policies" ! in the file mcidas/data/license.txt ! *** $Id: nvxrect.dlm,v 1.13 1999/03/31 21:36:59 rodr Exp $ *** ! ***************************************************************** ! * NVXINI(IFUNC,IPARMS) -- Initializes this modules Nav routine data ! * IFUNC -- ??? ! * IPARMS -- array of Navblock entries as follows ! * Projection : RECT ! * WORD 1 -- 4 character projection type RECT MERC LAMB etc.. ! * 2 -- Image Row Number of the Center of the image ! * 3 -- Latitude(degrees)*10000 of Image Row (2) ! * 4 -- Image Column Number of the Center of the image ! * 5 -- Longitude *10000 of 4 ! * 6 -- Lat(degrees per line)*10000 ! * 7 -- Lon(degrees per element)*10000 ! * 8 -- radius of the planet in meters ! * 9 -- eccentricity of the planet * 1000000 ! * 10 -- coordinate type >=0 panetodetic <0 panetocentric ! * 11 -- longitude convention >=0 west positive <0 west neg. ! * 12 -- LAT power of 10 -- if 0 leave at 4 else this power ! * 13 -- LON power of 10 -- if 0 leave at 4 else this power ! * 14 -- Degress_per_line power ! * 15 -- Degrees_per_element power ! * 16 -- Degrees_radian power ! * 17 -- Degrees_eccentricity power ! * ** -- next for added to IMGREMAP but not currently used anyway ! * 21 -- ecor -- element number of upper left hand corner ! * 22 -- lcor -- line number of upper left hand corner ! * 23 -- esize -- size in elements of coverage (width) ! * 24 -- lsize -- size in lines of coverage (height) ! * !************************************************************************* ! converted to a subroutine - JAK for use in this module !************************************************************************* IMPLICIT NONE INTEGER:: IPARMS(*) REAL*8 :: DRAD,DECC INTEGER :: IPOWLAT INTEGER :: IPOWLON INTEGER :: IPOWDLIN ! Degrees Per Line Power INTEGER :: IPOWDELE ! Degrees Per Element Power INTEGER :: IPOWRAD ! Degrees Radian power INTEGER :: IPOWECC ! Eccentricity Power INTEGER :: nbeg,nend,indwrd,ifunc CHARACTER (LEN=80):: fndum ! CHARACTER (LEN=4):: CLIT REAL:: XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,r REAL:: RAD=.01745329 INTEGER:: ITYPE,IWEST,i COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST ! ,ULLON IF (IFUNC.EQ.1) THEN IF (IPARMS(nbeg+1).NE.LIT('RECT')) GO TO 900 ITYPE=1 XROW=IPARMS(nbeg+2) IPOWLAT=IPARMS(nbeg+12) IF (IPOWLAT .EQ. 0) IPOWLAT=4 ! default is 10000 (10^4) ZSLAT=IPARMS(nbeg+3)/10.**IPOWLAT ! REAL Latitude XCOL=IPARMS(nbeg+4) IPOWLON=IPARMS(nbeg+13) IF (IPOWLON .EQ. 0) IPOWLON=4 ZSLON=IPARMS(nbeg+5)/10.**IPOWLON ! REAL Longitude IPOWDLIN=IPARMS(nbeg+14) IF (IPOWDLIN .EQ. 0) IPOWDLIN=4 ZDLAT=IPARMS(nbeg+6)/10.**IPOWDLIN ! REAL Degrees_per_line_latitude IPOWDELE=IPARMS(nbeg+15) IF (IPOWDELE .EQ. 0) IPOWDELE=4 ZDLON=IPARMS(nbeg+7)/10.**IPOWDELE ! REAL Degrees_per_line_longitude IPOWRAD=IPARMS(nbeg+16) IF (IPOWRAD .EQ. 0) IPOWRAD=3 DRAD=IPARMS(nbeg+8)/10.D0**IPOWRAD! REAL Radius of the planet in meters R=DRAD IPOWECC=IPARMS(nbeg+17) IF (IPOWECC .EQ. 0) IPOWECC=6 DECC=IPARMS(nbeg+9)/10.D0**IPOWECC ! REAL Eccentricity IWEST=IPARMS(nbeg+11) ! West positive vs. West negative IF(IWEST.GE.0) IWEST=1 CALL LLOPT(DRAD,DECC,IWEST,IPARMS(nbeg+10))! Initialze LLCART code IF (XCOL.EQ.1) THEN! special case of XCOL not located at image center ZSLON=ZSLON-180.0*IWEST ! -- so assume it's the left edge(duh) ENDIF ELSE IF (IFUNC.EQ.2) THEN IF(INDEX(CLIT(IPARMS(nbeg+1)),'XY').NE.0) ITYPE=1 IF(INDEX(CLIT(IPARMS(nbeg+1)),'LL').NE.0) ITYPE=2 ENDIF !NVXINI=0 ! not a function RETURN 900 CONTINUE PRINT*, "SUBROUTINE NVXINI: Might not be a RECTILINEAR FILE" !NVXINI=-1 ! not a function RETURN END SUBROUTINE NVXINI !*************************************************************************** !* NVXSAE -- Line/Element to Lat/Lon * !*************************************************************************** SUBROUTINE NVXSAE(XLIN,XELE,XDUM,XLAT,XLON,Z,ylat,ylon) IMPLICIT NONE REAL:: xlin,xele,xdum,xlat,xlon,xldif,xedif REAL:: XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON INTEGER ITYPE,IWEST REAL ylat,ylon,z COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST !,ULLON XLDIF=XROW-XLIN if(xcol.eq.1) then XEDIF=IWEST*(XELE-XCOL) XLON=ZSLON+180*IWEST-XEDIF*ZDLON else XEDIF=IWEST*(XCOL-XELE) XLON=ZSLON+XEDIF*ZDLON endif XLAT=ZSLAT+XLDIF*ZDLAT IF(XLAT.GT.90. .OR. XLAT.LT.-90.) GO TO 900 ! -- All we want is to keep the XLON within 180 degrees ! -- of the center LON on both sides IF(XLON.GT.(ZSLON+180.0)) GO TO 900 IF(XLON.LT.(ZSLON-180.0)) GO TO 900 IF(XLON.LT.-180.0) THEN XLON=XLON+360.0 ELSE IF (XLON.GT.180) THEN XLON=XLON-360.0 ENDIF itype=1 IF(ITYPE.EQ.1) THEN YLAT=XLAT YLON=XLON CALL LLCART(YLAT,YLON,XLAT,XLON,Z) ENDIF ! NVXSAE=0 RETURN 900 CONTINUE ! NVXSAE=-1 PRINT*,"SUBROUTINE NVXSAE: may be a problem with lat/lon ranges" RETURN END SUBROUTINE NVXSAE !*************************************************************************** !* NVXSAE -- Lat/Lon to Line/Element * !*************************************************************************** INTEGER FUNCTION NVXEAS(ZLAT,ZLON,Z,XLIN,XELE,XDUM) IMPLICIT NONE CHARACTER*12 cfe REAL:: zlat,zlon,z,xlin,xele,xdum,xlat,xlon,x,y REAL:: XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON INTEGER:: ITYPE,IWEST COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST ! ,ULLON XLAT=ZLAT XLON=ZLON IF(ITYPE.EQ.1) THEN X=XLAT Y=XLON CALL CARTLL(X,Y,Z,XLAT,XLON) ENDIF ! -- Keep XLON within 180.0 degrees of ZSLON IF(XLON.GT.(ZSLON+180.0)) THEN XLON=XLON-360.0 ELSEIF(XLON.LT.ZSLON-180.0) THEN XLON=XLON+360.0 ENDIF XLIN=XROW-(XLAT-ZSLAT)/ZDLAT IF (XCOL.EQ.1) THEN ! -- Need to adjust for the fact that XCOL is not really ! -- at the image center XELE=XCOL-(XLON-ZSLON-180.0*IWEST)/(ZDLON*IWEST) ELSE XELE=XCOL-(XLON-ZSLON)/(ZDLON*IWEST) ENDIF NVXEAS=0 RETURN END FUNCTION NVXEAS INTEGER FUNCTION NVXOPT(IFUNC,XIN,XOUT) IMPLICIT NONE REAL:: XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST COMMON/RCTCOM/XROW,XCOL,ZSLAT,ZSLON,ZDLAT,ZDLON,ITYPE,IWEST REAL*4 XIN(*),XOUT(*) CHARACTER (LEN=4):: CFUNC INTEGER:: ifunc ! ! IFUNC= 'SPOS' SUBSATELLITE LAT/LON ! ! XIN - NOT USED ! XOUT - 1. STANDARD LATITUDE ! - 2. NORMAL LONGITUDE ! ! ! IFUNC= 'ORAD' OBLATE RADIUS ! ! XIN - LATITUDE ! XOUT - RADIUS IN KM ! CFUNC=CLIT(IFUNC) NVXOPT=0 IF(CFUNC.EQ.'SPOS') THEN !ccc XOUT(1)=ZSLAT !ccc XOUT(2)=ZSLON ! !--- NO SAT. POSITION FOR THIS PROJECTION NVXOPT = -1 ELSE IF(CFUNC.EQ.'ORAD') THEN CALL LLOBL(XIN,XOUT) ELSE NVXOPT=1 ENDIF RETURN END FUNCTION NVXOPT SUBROUTINE JDAYI(julday,iyear,imon,iday) ! This routine calculates the month (imon) and day (iday) ! from the Julian day (julday) and year (iyear). ! The appropriate correction is made for leap year. ! INTEGER:: ndmon(12),nsum(13) INTEGER:: iyear,imon,iday,julday INTEGER:: mxjul,iijdi ! ! Specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 ! ! Correct for leap year IF (MOD(iyear,4) .EQ. 0) ndmon(2)=29 ! ! Check for illegal input IF (MOD(iyear,4) .EQ. 0) THEN mxjul = 366 ELSE mxjul = 365 END IF IF (julday .LT. 1 .OR. julday .GT. mxjul) THEN imon = -1 iday = -1 return END IF ! Calculate the month and day nsum(1) = 0 DO iijdi=1,12 nsum(iijdi+1) = nsum(iijdi) + ndmon(iijdi) END DO DO iijdi=2,13 IF (julday .LE. nsum(iijdi)) THEN imon = iijdi-1 go to 1000 END IF END DO 1000 CONTINUE iday = julday - nsum(imon) RETURN END SUBROUTINE JDAYI DOUBLE PRECISION FUNCTION yearfrac(ifyear,ifjday,iftime) INTEGER:: ifyear,ifjday,iftime REAL:: fday INTEGER::NXJUL,hhh,mmm,sss IF (MOD(ifyear,4) .EQ. 0) THEN nxjul = 366 ELSE nxjul = 365 END IF hhh=(iftime/10000) mmm=(iftime-hhh*10000)/100 sss=iftime-hhh*10000-mmm*100 fday=hhh/24.+mmm/1400.+sss/86400. yearfrac=dble(ifyear)+(dble(ifjday)+fday)/dble(nxjul) RETURN END FUNCTION yearfrac SUBROUTINE bitcon(idata,idatau) ! converts bit to integer 0 64 ! from bit structure on modern hardware ! John Knaff, NOAA/NESDIS ! INTEGER (KIND=1) ::idata INTEGER::ikb,idatau,iii,idatat INTEGER :: ibytekk(8) idatat=idata+32 IF(idatat.gt.255) THEN idatat=idatat-255 END IF idatau=idatat END SUBROUTINE bitcon END MODULE onebyte_1band_rect_read