MODULE merc_4km ! ! This module reads the 4km mercator projection images contained in ! the cira tropical cyclone IR archive. ! ! Adapted from code that J. Kossin wrote while a post doc at CIRA (2002). ! ! 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 6/10/2009 ! ! 12/30/2008 - Take into consideration a calibration block (JFD) ! 06/10/2009 - Logical error fixed wrt calibration block (JAK) ! 06/10/2009 - Time from image restored (JAK) ! !********************************************************************* IMPLICIT NONE ! ! Global variables INTEGER :: dtgj ! julian date YYJJJ INTEGER :: dtgg ! gregorian date YYMMDD INTEGER :: time ! HHMMSS INTEGER :: nliny,nelex DOUBLE PRECISION :: fyear ! fractional year YYYY.(jday/365) REAL, ALLOCATABLE :: lat(:),lon(:) ! latitude (degree N) longitude (0-360E) REAL, ALLOCATABLE :: bt (:,:) ! brightness temperature CONTAINS LOGICAL FUNCTION irar_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 ! CHARACTER (LEN=1) , ALLOCATABLE :: imgblk(:) INTEGER (KIND=1) , ALLOCATABLE :: imgblk(:) ! BYTE, ALLOCATABLE :: imgblk(:) INTEGER, ALLOCATABLE :: navdir(:) ! JFD added 12/30/2008 INTEGER, ALLOCATABLE :: caldir(:) INTEGER :: icb,ios 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 LOGICAL :: lfn, swap ! Executable statements irar_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=51 INQUIRE (FILE=trim(cfn),exist=lfn) IF (lfn) THEN OPEN (UNIT=luarea,FILE=trim(cfn),ACCESS='direct',& STATUS='old',recl=256,err=900,IOSTAT=ios) ! 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(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(areadir(12).NE.4.OR.areadir(13).NE.4) THEN ! PRINT*, 'In file ',cfn ! PRINT*, '...x or y-resolution is not' ! PRINT*, 'equal to 4 km. 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 ! ! 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 mercator 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 ! ! 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=901) READ(luarea,rec=1)areadir,navdir,imgblk nav_beg=areadir(35)/4+1 - 64 nav_end=areadir(34)/4 - 64 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=902) READ(luarea,rec=1)areadir,navdir,caldir,imgblk nav_beg=areadir(35)/4+1 - 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 ! 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 ! print*,'imgblk =',imgblk(icount) CALL byttowrd(imgblk(icount),indwrd,ival) ! print*,iline,ielem,indwrd,imgblk(icount),ival ! print*,icount,imgblk(icount),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 ! print*,iline,ielem,indwrd,imgblk(icount),ival,bt(ielem,iline) END DO END DO ! PRINT*,'AAA_1= ',AAA/float(icount) ! 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 ! Assign lines and elements (GLOBAL) nliny = areadir( 9) nelex = areadir(10) ! 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) irar_get=.true. ELSE PRINT*, 'Failed to open ',TRIM (cfn) END IF RETURN 900 CONTINUE PRINT *, '900 Error while opening file named ',TRIM(cfn) PRINT *, 'IOSTAT = ',ios RETURN 901 CONTINUE PRINT *, '901 Error while opening file named ',TRIM(cfn) PRINT *, 'IOSTAT = ',ios RETURN 902 CONTINUE PRINT *, '902 Error while opening file named ',TRIM(cfn) PRINT *, 'IOSTAT = ',ios RETURN 910 CONTINUE RETURN END FUNCTION irar_get ! ! ***************S U B R O U T I N E S****************** ! SUBROUTINE swap4(in) ! ! swaps the byte order of a 4-byte integer (big endian<->little endian) ! INTEGER :: in,out INTEGER (KIND=1) :: bout(4),b ! BYTE :: bout(4),b 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. ! INTEGER:: i,j,idum INTEGER (KIND=1):: word(4) ! BYTE:: word(4) 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. ! 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 CHARACTER*4 FUNCTION CLIT(L) IMPLICIT NONE INTEGER L ! --- local variables CHARACTER*4 C INTEGER*4 L1 EQUIVALENCE(L1,C) L1=L CLIT=C RETURN END FUNCTION CLIT ! ! 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 ! 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 ! ! 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: nvxmerc.dlm,v 1.10 1999/04/28 17:50:57 rodr Exp $ *** ! SUBROUTINE NVXINI(IFUNC,IPARMS,nbeg,nend,fndum,indwrd) INTEGER:: IPARMS(*) REAL*8 DRAD,DECC CHARACTER*4 dumindw1,dumindw2 character*80 fndum REAL:: xrow,xcol,xlat1,xspace,xqlon,xblat,xblon,rad,r INTEGER:: itype,iwest,leftlon,ifunc,nbeg,nend,indwrd COMMON/MERCOM/XROW,XCOL,XLAT1,XSPACE,XQLON,XBLAT,XBLON, & ITYPE,IWEST,LEFTLON DATA RAD/.01745329/ ! if(indwrd.eq.1)dumindw='MERC' ! if(indwrd.eq.4)dumindw='CREM' dumindw1='MERC' dumindw2='CREM' IF (IFUNC.EQ.1) THEN IF (IPARMS(nbeg).NE. LIT(dumindw1).and. & IPARMS(nbeg) .NE. LIT(dumindw2)) THEN ! PRINT*,CLIT(IPARMS(nbeg)),dumindw1,dumindw2 PRINT*, '...is not a MERCATOR projection. Execution halted.' STOP ! RETURN ENDIF ITYPE=1 XROW=IPARMS(nbeg+1) XCOL=IPARMS(nbeg+2) XLAT1=FLALO(IPARMS(nbeg+3)) XSPACE=IPARMS(nbeg+4)/1000. XQLON=FLALO(IPARMS(nbeg+5)) DRAD=IPARMS(nbeg+6)/1000.D0 R=DRAD DECC=IPARMS(nbeg+7)/1.D6 IWEST=IPARMS(nbeg+9) IF(IWEST.GE.0) IWEST=1 CALL LLOPT(DRAD,DECC,IWEST,IPARMS(nbeg+8)) XBLAT=R*COS(XLAT1*RAD)/XSPACE XBLON=RAD*R/XSPACE ELSE IF (IFUNC.EQ.2) THEN IF(INDEX(CLIT(IPARMS(nbeg)),'XY').NE.0) ITYPE=1 IF(INDEX(CLIT(IPARMS(nbeg)),'LL').NE.0) ITYPE=2 ENDIF !--- Adding to make sure NVXSAE drops points outside of bounding area LEFTLON=XQLON-180*IWEST ! NVXINI=0 RETURN END SUBROUTINE NVXINI SUBROUTINE NVXSAE(XLIN,XELE,XDUM,XLAT,XLON,Z,ylat,ylon) REAL:: xrow,xcol,xlat1,xspace,xqlon,xblat,xblon,rad REAL:: xlin,xele,xdum,xlat,xlon,z,ylat,ylon,xldif REAL:: xedif,xrlon,xrlat INTEGER:: itype,iwest,leftlon COMMON/MERCOM/XROW,XCOL,XLAT1,XSPACE,XQLON,XBLAT,XBLON, & ITYPE,IWEST,LEFTLON DATA RAD/.01745329/ XLDIF=XROW-XLIN XEDIF=XCOL-XELE XRLON=IWEST*XEDIF/XBLON XLON=XRLON+XQLON XRLAT=ATAN(EXP(XLDIF/XBLAT)) XLAT=(XRLAT/RAD-45.)*2.+XLAT1 ! ! IF(XLON.GT.180.) XLON=XLON-360. ! IF(XLON.LT.-180.) XLON=XLON+360. IF(XLON.GT.(360+LEFTLON)) GO TO 20 IF(XLON.LT.(LEFTLON)) GO TO 20 IF(ITYPE.EQ.1) THEN YLAT=XLAT YLON=XLON CALL LLCART(YLAT,YLON,XLAT,XLON,Z) ENDIF RETURN 20 continue PRINT*, 'NVSSAE:: somethings funny, but i dont know what' RETURN END SUBROUTINE NVXSAE SUBROUTINE NVXEAS(ZLAT,ZLON,Z,XLIN,XELE,XDUM) REAL:: xrow,xcol,xlat1,xspace,xqlon,xblat,xblon,rad REAL:: zlat,zlon,z,xlin,xele,xdum,xlat,xlon,x,y,xrlon,xrlat INTEGER:: itype,iwest,leftlon COMMON/MERCOM/XROW,XCOL,XLAT1,XSPACE,XQLON,XBLAT,XBLON,& ITYPE,IWEST,LEFTLON DATA RAD/.01745329/ XLAT=ZLAT XLON=ZLON IF(ITYPE.EQ.1) THEN X=XLAT Y=XLON CALL CARTLL(X,Y,Z,XLAT,XLON) ENDIF XRLON=IWEST*(XLON-XQLON) IF(XRLON.GT.180.) XRLON=XRLON-360. IF(XRLON.LT.-180.) XRLON=XRLON+360. IF(XLAT.GE.90.) XLAT=89.99 IF(XLAT.LE.-90.) XLAT=-89.99 XRLAT=((XLAT-XLAT1)/2.+45.)*RAD IF (XRLAT.LE.0.0) THEN ! Need to fail the call here since ALOG(<=0) will be undefined! ! call edest('Failing with XRLAT<=0.0',0) ! NVXEAS=-1 print*,'xrlat<=0 in routine nvxeas. Halting.' stop ! RETURN ENDIF XLIN=XROW-XBLAT*ALOG(TAN(XRLAT)) XELE=XCOL-XRLON*XBLON ! NVXEAS=0 RETURN END SUBROUTINE NVXEAS FUNCTION NVXOPT(IFUNC,XIN,XOUT) REAL:: xrow,xcol,xlat1,xspace,xqlon,xblat,xblon INTEGER:: itype,iwest,leftlon,ifunc,nvxopt COMMON/MERCOM/XROW,XCOL,XLAT1,XSPACE,XQLON,XBLAT,XBLON, & ITYPE,IWEST,LEFTLON REAL*4:: XIN(*),XOUT(*) CHARACTER*4 :: CFUNC ! ! 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 XOUT(1)=XLAT1 XOUT(2)=XQLON 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) INTEGER (KIND=1) ::idata INTEGER::ikb,idatau,iii,idatat INTEGER :: ibytekk(8) ! idatau=INT(idata,kind=1) idatat=idata+32 IF(idatat.gt.255) THEN idatat=idatat-255 END IF idatau=idatat ! DO ikb=6,0,-1 ! ibytekk(ikb+2)=0 ! IF(idatat.GE.2**ikb)THEN ! idatat=idatat-2**ikb ! ibytekk(ikb+2)=1 ! END IF ! END DO ! idatau=0 ! DO ikb=1,8 ! IF(ibytekk(ikb).EQ.1) THEN ! print*,'ikb,ibytekk,iii =',ikb,ibytekk(ikb),iii ! idatau=idatau+2**(ikb-1) ! END IF ! END DO ! IF(ibytekk(1).EQ.1) THEN ! idatau=127+idatau ! ENDIF END SUBROUTINE bitcon END MODULE merc_4km