PROGRAM STALST C---------------------------------------------------------------------- INCLUDE "parmeta" C---------------------------------------------------------------------- PARAMETER (IMJM=IM*JM-JM/2,IMT=IM*2-1,JMT=JM/2+1) PARAMETER (NSTA=2000) PARAMETER (JAM=6+2*(JM-10),LP1=LM+1) PARAMETER (D2R=1.745329252E-2,R2D=57.29577951) C---------------------------------------------------------------------- INCLUDE "LOOPS.com" C---------------------------------------------------------------------- INCLUDE "MAPOT.com" C---------------------------------------------------------------------- INCLUDE "MASKS.com" C---------------------------------------------------------------------- INCLUDE "INDX.com" C---------------------------------------------------------------------- C H A R A C T E R & AN*1,AW*1,CID*8,CIDDUM*8,CIDX*8 C---------------------------------------------------------------------- D I M E N S I O N & RLAT(NSTA),RLON(NSTA),TLAT(NSTA),TLON(NSTA) &,RLATD(NSTA),RLOND(NSTA) &,IHDUM(NSTA),JHDUM(NSTA),IVDUM(NSTA),JVDUM(NSTA) &,IHSTA(NSTA),JHSTA(NSTA),IVSTA(NSTA),JVSTA(NSTA) &,KCLAS(NSTA),KCLASD(NSTA),IDDUM(NSTA),ID(NSTA) &,CIDDUM(NSTA),CID(NSTA) C---------------------------------------------------------------------- R E A L & LAM0 C---------------------------------------------------------------------- DATA LNHB/12/,LUNLST/15/ C---------------------------------------------------------------------- C*** C*** WHEN CHANGING THE GRID, YOU MUST RESET VALUES IN THE BLOCK DATA C*** AT THE END OF THIS PROGRAM. C*** C---------------------------------------------------------------------- C*** C*** READ IN THE MASKS THAT ARE NEEDED C*** OPEN(UNIT=12,FILE='CNST.file',FORM='UNFORMATTED', + ACCESS='sequential') REWIND 12 READ(12) READ(12)LMH READ(12)LMV READ(12) READ(12) READ(12) READ(12)SM READ(12)SICE C*** C*** CALCULATE THE I-INDEX EAST-WEST INCREMENTS C*** DO J=1,JM IHE(J)=MOD(J+1,2) IHW(J)=IHE(J)-1 IVE(J)=MOD(J,2) IVW(J)=IVE(J)-1 ENDDO C write(6,*)' tph0d=',tph0d,' tlm0d=',tlm0d PHI0=TPH0D*D2R LAM0=-1.0*TLM0D*D2R DPHI=DPHD*D2R DLAM=DLMD*D2R C C SBDB AND WBDB ARE 2 ROWS INSIDE THE DOMAIN BOUNDARY C TO TRY TO EXCLUDE POINTS RIGHT ON THE BOUNDARY C SBDB=SBD+2.0*DPHD WBDB=WBD+2.0*DLMD write(6,*)' sbdb=',sbdb,' wbdb=',wbdb C*** C*** LOOP THROUGH ALL THE STATIONS TO COUNT HOW MANY INSIDE DOMAIN C*** NSUM=0 25 READ(LUNLST,50,END=100)IDX,RLATX,AN,RLONX,AW,CIDX,KCLS 50 FORMAT(I6,F6.2,A1,F7.2,A1,1X,A4,I3) IF (AN.EQ.'S') RLATX=-1.0*RLATX IF (AW.EQ.'E') RLONX=360.0-RLONX RLATX=RLATX*D2R RLONX=RLONX*D2R C*** C*** CONVERT GEODETIC TO TRANSFORMED COORDINATES OF THE STATION C*** X=COS(PHI0)*COS(RLATX)*COS(RLONX-LAM0) 1 +SIN(PHI0)*SIN(RLATX) Y=-COS(RLATX)*SIN(RLONX-LAM0) Z=COS(PHI0)*SIN(RLATX) 1 -SIN(PHI0)*COS(RLATX)*COS(RLONX-LAM0) TLATX=R2D*ATAN(Z/SQRT(X*X+Y*Y)) TLONX=R2D*ATAN(Y/X) IF(ABS(TLATX).GT.ABS(SBDB).OR.ABS(TLONX).GT.ABS(WBDB))GO TO 25 NSUM=NSUM+1 RLATD(NSUM)=RLATX RLOND(NSUM)=RLONX TLAT(NSUM)=TLATX TLON(NSUM)=TLONX IDDUM(NSUM)=IDX CIDDUM(NSUM)=CIDX KCLASD(NSUM)=KCLS GOTO 25 100 CONTINUE C C---------------------------------------------------------------------- DO 200 N2=1,NSUM C---------------------------------------------------------------------- C*** C*** FIND THE CLOSEST POINT REGARDLESS OF LAND/SEA MASK C*** ROW=TLAT(N2)/DPHD+JMT COL=TLON(N2)/DLMD+IM NROW=INT(ROW) NCOL=INT(COL) TLAT(N2)=TLAT(N2)*D2R TLON(N2)=TLON(N2)*D2R C*** C*** IDENTIFY THE ROW AND COLUMN OF THE LOWER LEFT CORNER OF THE BOX C*** THAT SURROUNDS EACH STATION INSIDE THE DOMAIN C*** C*** FIRST CONSIDER THE SITUATION WHERE THE STATION X IS AT C*** C*** H V C*** C*** C*** X C*** V H C*** IF(MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.1.OR. 1 MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.0)THEN TLAT1=(NROW-JMT+1)*DPHI TLAT2=TLAT1-DPHI TLON1=(NCOL-IM)*DLAM TLON2=TLON1+DLAM DLM1=TLON(N2)-TLON1 DLM2=TLON(N2)-TLON2 D1=ACOS(COS(TLAT(N2))*COS(TLAT1)*COS(DLM1) 1 +SIN(TLAT(N2))*SIN(TLAT1)) D2=ACOS(COS(TLAT(N2))*COS(TLAT2)*COS(DLM2) 1 +SIN(TLAT(N2))*SIN(TLAT2)) IF(D1.GT.D2)THEN NCOL=NCOL+1 ELSE NROW=NROW+1 ENDIF C*** C*** NOW CONSIDER THE SITUATION WHERE THE POINT X IS AT C*** C*** V H C*** C*** C*** X C*** H V C*** ELSE TLAT1=(NROW-JMT)*DPHI TLAT2=TLAT1+DPHI TLON1=(NCOL-IM)*DLAM TLON2=TLON1+DLAM DLM1=TLON(N2)-TLON1 DLM2=TLON(N2)-TLON2 D1=ACOS(COS(TLAT(N2))*COS(TLAT1)*COS(DLM1) 1 +SIN(TLAT(N2))*SIN(TLAT1)) D2=ACOS(COS(TLAT(N2))*COS(TLAT2)*COS(DLM2) 1 +SIN(TLAT(N2))*SIN(TLAT2)) IF(D1.GT.D2)THEN NROW=NROW+1 NCOL=NCOL+1 ENDIF ENDIF C*** C*** NOW WE CAN FIND THE K VALUE C*** KROWS=((NROW-1)/2)*IMT IF(MOD(NROW,2).EQ.1)THEN KMIN=KROWS+(NCOL+1)/2 ELSE KMIN=KROWS+IM+NCOL/2 ENDIF JX=(KMIN-1)/IMT+1 IX=KMIN-(JX-1)*IMT IF(IX.LE.IM)THEN IDUM=IX JDUM=2*JX-1 ELSE IDUM=IX-IM JDUM=2*JX ENDIF C IF (KCLASD(N2).GE.10.AND.KCLASD(N2).LT.20) THEN C C WE WANT THIS TO BE A LAND POINT, IF NOT, FIND THE CLOSEST LAND POINT C IF (SM(IDUM,JDUM)+SICE(IDUM,JDUM).GT.0.5) THEN DMIN=99999. DO J=1,JM TLAT1=(J-JMT)*DPHI DO I=1,IM IHV=I*2-MOD(J,2) TLON1=(IHV-IM)*DLAM DLM1=TLON(N2)-TLON1 IF(SM(I,J)+SICE(I,J).LT.0.5) THEN D1=ACOS(COS(TLAT(N2))*COS(TLAT1)*COS(DLM1)+ 1 SIN(TLAT(N2))*SIN(TLAT1)) IF(D1.LT.DMIN) THEN DMIN=D1 NCOL=IHV NROW=J ENDIF ENDIF ENDDO ENDDO ENDIF KROWS=((NROW-1)/2)*IMT IF(MOD(NROW,2).EQ.1)THEN KMIN=KROWS+(NCOL+1)/2 ELSE KMIN=KROWS+IM+NCOL/2 ENDIF ENDIF C IF (KCLASD(N2).GE.20.AND.KCLASD(N2).LT.30) THEN C C WE WANT THIS TO BE A SEA POINT, IF NOT, FIND THE CLOSEST SEA POINT C IF (SM(IDUM,JDUM)+SICE(IDUM,JDUM).LT.0.5) THEN DMIN=99999. DO J=1,JM TLAT1=(J-JMT)*DPHI DO I=1,IM IHV=I*2-MOD(J,2) TLON1=(IHV-IM)*DLAM DLM1=TLON(N2)-TLON1 IF(SM(I,J)+SICE(I,J).GT.0.5) THEN D1=ACOS(COS(TLAT(N2))*COS(TLAT1)*COS(DLM1)+ 1 SIN(TLAT(N2))*SIN(TLAT1)) IF(D1.LT.DMIN) THEN DMIN=D1 NCOL=IHV NROW=J ENDIF ENDIF ENDDO ENDDO ENDIF KROWS=((NROW-1)/2)*IMT IF(MOD(NROW,2).EQ.1)THEN KMIN=KROWS+(NCOL+1)/2 ELSE KMIN=KROWS+IM+NCOL/2 ENDIF ENDIF C*** C*** FIND THE GEODETIC LAT/LON C*** TLATS=(NROW-JMT)*DPHI TLONS=(NCOL-IM)*DLAM A=SIN(TLATS)*COS(PHI0)+COS(TLATS)*SIN(PHI0)*COS(TLONS) ALAT=ASIN(A) RLATD(N2)=ALAT B=COS(TLATS)*COS(TLONS)/(COS(ALAT)*COS(PHI0)) D=B-TAN(ALAT)*TAN(PHI0) IF (ABS(D).GT.1.0) D=ABS(D)/D C=ACOS(D) C=SIGN(C,-1.0*TLONS) RLOND(N2)=LAM0+C C*** C*** THIS IS THE H POINT, NOW FIND THE CLOSEST NEIGHBORING C*** V POINT THAT IS ON THE SAME STEP C*** C*** 2 KMIN+IM-1 C*** C*** 1 H 3 KMIN-1 H KMIN C*** C*** 4 KMIN-IM C*** C*** TLAT1=(NROW-JMT)*DPHI TLON1=(NCOL-IM-1)*DLAM TLAT2=TLAT1+DPHI TLON2=TLON1+DLAM TLAT3=TLAT1 TLON3=TLON2+DLAM TLAT4=TLAT1-DPHI TLON4=TLON2 DLM1=TLON(N2)-TLON1 DLM2=TLON(N2)-TLON2 DLM3=TLON(N2)-TLON3 DLM4=TLON(N2)-TLON4 D1=ACOS(COS(TLAT(N2))*COS(TLAT1)*COS(DLM1) 1 +SIN(TLAT(N2))*SIN(TLAT1)) D2=ACOS(COS(TLAT(N2))*COS(TLAT2)*COS(DLM2) 1 +SIN(TLAT(N2))*SIN(TLAT2)) D3=ACOS(COS(TLAT(N2))*COS(TLAT3)*COS(DLM3) 1 +SIN(TLAT(N2))*SIN(TLAT3)) D4=ACOS(COS(TLAT(N2))*COS(TLAT4)*COS(DLM4) 1 +SIN(TLAT(N2))*SIN(TLAT4)) DMIN=99999. C CALL K2IJH(KMIN,IMIN,JMIN) C IF(LMH(IMIN,JMIN).EQ.LMV(IMIN+IHW(JMIN),JMIN).AND. 1 D1.LT.DMIN)THEN DMIN=D1 KVMIN=KMIN-1 ENDIF C IF(LMH(IMIN,JMIN).EQ.LMV(IMIN,JMIN+1).AND.D2.LT.DMIN)THEN DMIN=D2 KVMIN=KMIN+IM-1 ENDIF C IF(LMH(IMIN,JMIN).EQ.LMV(IMIN+IHE(JMIN),JMIN).AND. 1 D3.LT.DMIN)THEN DMIN=D3 KVMIN=KMIN ENDIF C IF(LMH(IMIN,JMIN).EQ.LMV(IMIN,JMIN-1).AND.D4.LT.DMIN)THEN DMIN=D4 KVMIN=KMIN-IM ENDIF C IF(DMIN.EQ.99999.)THEN KVMIN=KMIN WRITE(6,180)IDDUM(N2),RLATD(N2)*R2D,RLOND(N2)*R2D,KMIN, 1 LMH(IMIN,JMIN),LMV(IMIN+IHW(JMIN),JMIN),LMV(IMIN,JMIN+1) 2, LMV(IMIN+IHE(JMIN),JMIN),LMV(IMIN,JMIN-1) ENDIF C KH=KMIN KV=KVMIN C JX=(KH-1)/IMT+1 IX=KH-(JX-1)*IMT IF(IX.LE.IM)THEN IHDUM(N2)=IX JHDUM(N2)=2*JX-1 ELSE IHDUM(N2)=IX-IM JHDUM(N2)=2*JX ENDIF C JX=(KV-1)/IMT+1 IX=KV-(JX-1)*IMT IF(IX.LE.IM-1)THEN IVDUM(N2)=IX JVDUM(N2)=2*JX-1 ELSE IVDUM(N2)=IX-IM+1 JVDUM(N2)=2*JX ENDIF C 180 FORMAT(I6,2F8.2,I8,5I4,' NO NEIGHBOR ON SAME STEP') 200 CONTINUE C*** NDUM=0 DO I=1,NSUM IF(IDDUM(I).NE.-99999)THEN NDUM=NDUM+1 IHSTA(NDUM)=IHDUM(I) JHSTA(NDUM)=JHDUM(I) IVSTA(NDUM)=IVDUM(I) JVSTA(NDUM)=JVDUM(I) ID(NDUM)=IDDUM(I) RLAT(NDUM)=RLATD(I) RLON(NDUM)=RLOND(I) CID(NDUM)=CIDDUM(I) KCLAS(NDUM)=MOD(KCLASD(I),10) ENDIF ENDDO C NSUM=NDUM C----------------------------------------------------------------------- WRITE(63)NSUM,ID,RLAT,RLON,IHSTA,JHSTA,IVSTA,JVSTA,CID DTR=1.74532925E-2 WRITE(6,20)NSUM open(unit=16,file='nsta.txt',form='formatted', + access='sequential') WRITE(16,643)NSUM 643 FORMAT(8x,'PARAMETER (NSTAT=',I4,')') 20 FORMAT('STALST: NUMBER OF PROFILE STATIONS ',I5) WRITE(6,29)'ID','LAT','LON','IHSTA','JHSTA','IVSTA','JVSTA','CID' WRITE(6,30)(ID(N),RLAT(N)/DTR,RLON(N)/DTR 1, IHSTA(N),JHSTA(N),IVSTA(N),JVSTA(N) 2, CID(N),N=1,NSUM) 29 FORMAT(7A8,A4) 30 FORMAT(2X,I6,2F8.2,4I8,1X,A8) WRITE(1001,31)(ID(N),CIDX,CID(N),N=1,NSUM) 31 FORMAT(I6,1X,A8,1X,A8) C----------------------------------------------------------------------- ! STOP END SUBROUTINE K2IJH(K,I,J) C---------------------------------------------------------- C*** C*** COMPUTE I,J VALUES OF H POINTS FROM A ONE DIMENSIONAL K C*** C---------------------------------------------------------- INCLUDE "parmeta" C---------------------------------------------------------- C*** C*** CONVERT FROM ONE DIMENSIONAL K TO I,J ON THE E-GRID C*** IMT=2*IM-1 JX=(K-1)/IMT+1 IX=K-(JX-1)*IMT IF(IX.LE.IM)THEN I=IX J=2*JX-1 ELSE I=IX-IM J=2*JX ENDIF C RETURN END C BLOCK DATA STN INCLUDE "parmeta" PARAMETER (LP1=LM+1) INCLUDE "MAPOT.com" C--------------------------FINAL THREE LINES ADDED by new_prep.sh job DATA WBD/-34.0384615430000/,SBD/-33.3846154180000/ DATA TPH0D/-19/,TLM0D/-50/ DATA DPHD/.26923076950000000000/,DLMD/.28846153850000000000/ END BLOCK DATA STN