winterp.f

Go to the documentation of this file.

C*********************************************************************
        SUBROUTINE topwnd
C*********************************************************************
C
C  INTERPOLATES WINDS AT TOPMOST LEVEL BETWEEN OBSERVATION PTS--WHEN
C  ONLY ONE WIND AVAILABLE IT'S USED EVERYWHERE. LUDWIG, NOVEMBER,1987.
C
    include 'NGRIDS.PAR'
       
C
        include 'FLOWER.CMM'
        include 'LIMITS.CMM'
    include 'STALOC.CMM'
C
        dimension uu(nwsite),vv(nwsite)
C
C  USE THE ONE WIND AVAILABLE EVERYWHERE
C
        IF (numdop.EQ.1) THEN
          uu(1)=usig(1,nlvl)
          vv(1)=vsig(1,nlvl)
          DO 25 i = 1,ncol
            DO 25 j = 1,nrow
              u(i,j,nlvl)=uu(1)
              v(i,j,nlvl)=vv(1)
25        CONTINUE
        ELSE
          DO 30 i=1,numdop
            uu(i)=usig(i,nlvl)
            vv(i)=vsig(i,nlvl)
30        CONTINUE
C
          DO 65 j = 1,nrow
        xx=float(j)
        DO 60 i = 1,ncol
          yy=float(i)
          sumdop=0.0
              u(i,j,nlvl)=0.0
              v(i,j,nlvl)=0.0
              DO 50 ik=1,numdop
C
C  FOLLOWING STATEMENT CHANGED [ FROM IT=IDOP(JK) ] BY FLL 5/24/2000
C  BUG FOUND BY DOUG MILLER OF THE NAVAL POSTGRADUATE SCHOOL.  IT
C  APPEARS TO AFFECT ONLY RESULTS OBTAINED WHEN MORE THAN ONE
C  SOUNDING IS AVAILABLE.
C
                   it=ik
           dwate=wndwt(xx,yy,xdop(it),ydop(it))
           sumdop=sumdop+dwate
                   u(i,j,nlvl)=u(i,j,nlvl)+dwate*uu(it)
                   v(i,j,nlvl)=v(i,j,nlvl)+dwate*vv(it)
50      CONTINUE
        IF (sumdop.GT.0.0) THEN
                   u(i,j,nlvl)=u(i,j,nlvl)/sumdop
                   v(i,j,nlvl)=v(i,j,nlvl)/sumdop
        ELSE
           WRITE (*,*) 'NO TOP WINDS'
           stop
        END IF
60          CONTINUE
C
65    CONTINUE
        ENDIF
C
        RETURN
C
        END
C
C*********************************************************************
        SUBROUTINE tstwnd( CHMMDDTTTT)
C*********************************************************************
C
C  USES LEVWND TO INTERPOLATE MASS ADJUSTED FIELD TO ANEMOMETER 
C  HEIGHT (Z10, FOR INDEX=1) & NFLAT FLAT PLANES (FOR Z INDICES 2-6) 
C  SET AT THE HEIGHTS  ZCHOOZ ABOVE THE LOWEST TERRAIN GRID POINT IN 
C  THE COARSE GRID, AND TO CONVERT INTERPOLATED WINDS BACK TO 
C  METEOROLOGICAL SPEEDS AND ANGLES, THAT CAN BE PLOTTED IF DESIRED.  
C  
C   LUDWIG--FEB 2000
C
          include 'NGRIDS.PAR'
        
C
        parameter(rad2d=180./3.14159,zero=0.0,nmax=200)
    parameter(valmis=-9999.0)
C
    include 'ANCHOR.CMM'
    include 'FLOWER.CMM'
    include 'LIMITS.CMM'
        include 'STALOC.CMM'
C
    INTEGER NCALLS
    CHARACTER*8 CHMMDDTTTT
C
C  INTERPOLATES MASS ADJUSTED FIELD TO ANEMOMETER HEIGHT (Z10, FOR 
C  INDEX=1) THIS SUBROUTINE ALSO CONVERTS INTERPOLATED WINDS BACK TO 
C  METEOROLOGICAL SPEEDS AND ANGLES, THAT CAN BE PLOTTED IF DESIRED.  
C      LUDWIG--JANUARY 1988
C
C  REVISED FROM LEVWND MARCH 1997 TO JUST GET SFC LEVEL WINDS AT 
C  SPECIFIED GRID POINTS SO THAT THEY CAN BE COMPARED WITH THE 
C  OBSERVATIONS THAT WERE NOT USED IN CALCULATING THE WINDS.
C
C  MAKE SURE EVERYTHING IS STILL HERE ON SUBSEQUENT CALLS
C
    DATA ncalls /0/
C
    SAVE
C       
    CALL levwnd
C
C  FOR CHECKING, WE WILL COMPARE OBSERVATIONS WITH NEAREST GRID VALUES
C  FOR STATIONS ON GRID
C
    kase=0
C   IF (NCALLS .EQ.0) THEN
C      NCALLS=1
C      WRITE (33,6001)
C   END IF
    DO 75 it=1,numnws
C
C  GETTING NEAREST GRID POINT FOR EACH OBS
C
       klx=nint(xg(it))
           kly=nint(yg(it))
C
C  SKIP IF OUTSIDE CALCULATION GRID, OR OBSERVATION WAS MISSING
C
       IF (klx.GT.0 .AND. kly.GT.0 .AND.
     $      klx.LE.ncol .AND. kly.LE.nrow .AND.
     $                  nint(ucomp(it)) .NE. -9999) THEN
              kase=kase+1
          difdst=dscrs*sqrt((float(klx)-xg(it))**2+
     $                    (float(kly)-yg(it))**2)
          ztst=sfcht(klx,kly)
C
C CONVERTING COMPONENTS BACK TO SPD & DIRECTION FOR OBS 
C
              wsob=sp(ucomp(it),vcomp(it))
              wdob=dd(ucomp(it),vcomp(it))
C
C  GETTING CALCULATED WINDS AT Z10 HT  
C
          ulxly=0.01*float(iugraf(klx,kly,1))
          vlxly=0.01*float(ivgraf(klx,kly,1))
              wstst=sp(ulxly,vlxly)
              wdtst=dd(ulxly,vlxly)
C
C  IF OBSERVED WIND WAS CALM DIRECTION DIFFERENCE IS UNDEFINED
C
          IF (wsob .GT. 0.0) THEN
             deldir=degdif(wdob,wdtst)
          ELSE
             deldir=-999.
          END IF
C         WRITE(33,6002)IT,CHSTID(IT),WSOB,WSTST,WDOB,WDTST,
C     $            WSOB-WSTST,DELDIR,DIFDST,XG(IT),YG(IT),CHMMDDTTTT
       END IF
75      CONTINUE
C
6001    FORMAT ('  NUM  ID        SDPOB    SPDWOX  DIROB    DIRWOX ',
     $       '    SPDIF   DIRDIF   DISTDIF  XKM    YKM     DATE')
6002    FORMAT(1x,i4,1x,a5,2x,6f9.2,3f7.2,4x,a8)
C
        RETURN
C
        END