dopsig.f

Go to the documentation of this file.

C********************************************************************
        SUBROUTINE dopsig(ICALL)
C********************************************************************
C
C  ASSIGN WND STA WIND PROFILES TO FLOW SURFACES.MISSING WINDS ARE 
C  DENOTED BY -9999. IF SOUNDING IS NOT COMPLETE THE LAST REPORTED WIND 
C  IS USED AT THE HIGHEST ALTITUDES. AFTER FLOW SURFACES ARE REDEFINED 
C  (RESIG), DOPSIG IS RECALLED (ICALL >1) AND THE SOUNDINGS ARE 
C  REINTERPOLATED TO THE SURFACES. TOPWIND AND BETWIN ARE RECALLED TO 
C  PROVIDE THE WINDS THAT ARE TO BE BALANCED. THE ORIGINAL VERSION OF 
C  THIS SUBROUTINE WAS WRITTEN BY R.M. ENDLICH, SRI INTN'L, MENLO PARK 
C  CA 94025. 
C  IT WAS LARGELY REWRITTEN  BY LUDWIG IN NOV,1987. THIS VERSION 
C  INCLUDES CHANGES MADE IN APRIL, 1989.
C
C  FURTHER MODIFIED MAY 1989 TO PROVIDE FOR A SECOND CALL WHERE WINDS 
C  ARE REINTERPOLATED TO THE NEWLY DEFINED FLOW SURFACES; DATA READING
C  AND OTHER PARTS OF THE ROUTINE ARE SKIPPED -- F. LUDWIG
C
C
        include 'NGRIDS.PAR'
C
    include 'ANCHOR.CMM'
        include 'FLOWER.CMM'
        include 'LIMITS.CMM'
        include 'STALOC.CMM'
        include 'TSONDS.CMM'
C
    REAL DPHT(nsndht,nwsite), DPUC(nsndht,nwsite)
    REAL DPVC(nsndht,nwsite),RHS1(nzgrd),SYTHYT(nwsite)
    INTEGER LEVHI(nwsite)
C
        DATA uvo/0.0/
    SAVE
C
C  VARIABLES ARE:
C    DPUC=U COMPONENT OF WND STA WIND IN MPS
C    DPVC=V COMPONENT OF WND STA WIND IN MPS
C    NWHTS=NUMBER OF POINTS IN VERTICAL WIND PROFILE
C    NLVL=NUMBER OF FLOW LEVELS
C    RHS=HT OF FLOW SURFACES ABOVE TERRAIN (M)
C    XG,YG=STA. DIST IN X,Y IN GRID UNITS FROM 0,0 (SW CORNER)
C    Z0, UVO = ROUGHNESS HT AND ZERO-LEVEL WIND FOR INTERPOLATING
C    SYTHYT(IS)=SFC ELEVATION FOR SITE IS
C    LEVHI(IS)=HIGHEST FLOW SURFACE FOR WHICH WIND WAS OBSERVED AT SITE IS
C
C    ON 1ST CALL AT EACH TIME, READ IN WIND PROFILES.
C
    IF (icall .EQ. 1) THEN
       DO 45 it=1,numdop
C
C  ZERO VALUES ON FLOW SURFACES BEFORE CALCULATING
C
          DO 11 ilev=1,nlvl
             usig(it,ilev)=0.0
             vsig(it,ilev)=0.0
11        CONTINUE
C
C  SKIP BLANK LINE THEN READ COORDINATES OF SITE
C
              READ (12,*)
              READ (12,*)  dopx,dopy,dopz
          sythyt(it)=dopz
C
C  LOCATE GRID POINT NEAREST THE SOUNDING (NOTE ORIGIN AT 1,1).
C
          xdop(it)=1.0+(dopx-xorig)/dscrs
          ydop(it)=1.0+(dopy-yorig)/dscrs
          zdop(it)=dopz
              ix=max(1,nint(xdop(it)))
              jy=max(1,nint(ydop(it)))
              ix=min(ix,ncol)
              jy=min(jy,nrow)
          zfix=sfcht(ix,jy)-sythyt(it)
              READ (12,*) nwhts(it)
C
C  SKIP COLUMN HEADING LINE
C
              READ (12,*)
          DO 15 ll=1,nwhts(it)
                 READ (12,*) zzht,zzwd,zzws
C   
                 IF (ll .LT. nsndht) THEN
                    dpht(ll,it)=zzht-sythyt(it)
                    dpwd(ll)=zzwd
                    dpws(ll)=zzws
             ELSE
                    dpht(nsndht,it)=zzht-sythyt(it)
                    dpwd(nsndht)=zzwd
                    dpws(nsndht)=zzws
                 END IF
15        CONTINUE
          IF (nwhts(it) .GT. nsndht) nwhts(it)=nsndht
C
C  CONVERT WIND MEASUREMENT HEIGHTS IN METERS (MSL) TO METERS (AGL)
C  ADJUST FOR FACT THAT SOUNDING LOCATION NOT AT GRID PT.
C
          DO 25 ll=1,nlvl
                 rhs1(ll)=rhs(ix,jy,ll)
         zsigl(it,ll)=rhs1(ll)
25        CONTINUE
C
C
C  CHANGE DIRECTION AND SPEED (MPS) TO U AND V; CHECK FOR MISSING DATA 
C  OR BELOW GROUND HEIGHTS (AFTER CONVERSION FROM MSL) ON 1ST CALL.
C
          newll=0
              DO 40 ll=1,nwhts(it)
                 IF (dpwd(ll).GT.-9998.9) THEN
                    IF (dpht(ll,it) .GT. z0) THEN
                       newll=newll+1
                       dpuc(newll,it)=
     $                          -dpws(ll)*sin(dpwd(ll)/57.295)
                       dpvc(newll,it)=
     $                          -dpws(ll)*cos(dpwd(ll)/57.295)
                       dpht(newll,it)=dpht(ll,it)
                    END IF
                 END IF
40        CONTINUE
              nwhts(it)=newll
45     CONTINUE
    END IF
C
C  INTERPOLATE TO ORIGINAL, OR NEWLY DEFINED, FLOW SURFACE HEIGHTS
C  AT LEVELS WHERE DATA ARE AVAILABLE
C   
    DO 150 it=1,numdop
           ix=max(1,nint(xdop(it)))
           jy=max(1,nint(ydop(it)))
           ix=min(ix,ncol)
           jy=min(jy,nrow)
C
C  FIND HIGHEST FLOW SURFACE REACHED BY THIS SOUNDING. ALSO CHECK
C  FOR HIGHEST SFC REACHED BY ANY SNDING
C
       levhi(it)=0
           DO 100 kl=1,nlvl
              zlevkl=rhs(ix,jy,kl)
              rhs1(kl)=zlevkl
          zsigl(it,kl)=zlevkl
          DO 96 iz=1,nwhts(it)
             IF (dpht(iz,it).GT.zlevkl) levhi(it)=kl
96        CONTINUE
100    CONTINUE
C
C  ASSIGNING OBSERVED WINDS TO FLOW SURFACE HEIGHTS UP TO TOP 
C  OF SONDE
C
       DO 140 kl=1,levhi(it)
          zf=rhs1(kl)
C
          IF (zf .LE. z0) THEN
C
C  ZERO WIND WHEN FLOW SFC BELOW ROUGHNESS HEIGHT.
C
                 usig(it,kl)=0.0
                 vsig(it,kl)=0.0
          ELSE IF (zf.GT.z0 .AND. zf.LE. dpht(1,it)) THEN
C
C  FLOW SFC BELOW 1ST OBSERVATION HT -- LOG INTERP TO FLOW SFC.
C
             CALL lgntrp(usig(it,kl),z0,dpht(1,it),
     $                                  zf,uvo,dpuc(1,it))
                 CALL lgntrp(vsig(it,kl),z0,dpht(1,it),
     $                                  zf,uvo,dpvc(1,it))
          ELSE
C
C  FLOW SURFACE ABOVE LOWEST OB -- FIND WHICH OBS IT IS BETWEEN
C  AND LOG INTERPOLATE.
C
             DO 135 jht=1,nwhts(it)-1
                zw1=dpht(jht,it)
                    zw2=dpht(jht+1,it)
                    IF (zf.GT.zw1 .AND. zf.LE.zw2) THEN
C
C FLOW SFC BETWEEN OBSERVATION HTS
C
                       CALL lgntrp(usig(it,kl),zw1,zw2,zf,
     $                           dpuc(jht,it),dpuc(jht+1,it))
                       CALL lgntrp(vsig(it,kl),zw1,zw2,zf,
     $                             dpvc(jht,it),dpvc(jht+1,it))
                       GO TO 140
                    END IF
135      CONTINUE
C
          END IF
140    CONTINUE
150 CONTINUE
C
C  FIND HIGHEST LEVEL REACHED BY ANY SONDE
C
    levtop=0
    DO 160 jsond=1,numdop
       levtop=max(levtop,levhi(jsond))
160 CONTINUE
C
C  ALL OBSERVED DATA INTERPOLATED VERTICALLY TO FLOW SURFACES.  NOW,
C  INTERPOLATE HORIZONTALLY UP TO TOP OBSERVATION WHERE REQUIRED
C
    DO 200 jsond=1,numdop
           xx=xdop(jsond)
       yy=ydop(jsond)
       DO 180 il=levhi(jsond)+1,levtop
C
C  THIS ONE IS MISSING INTERPOLATE FROM OTHERS AVAILABLE AT THIS LEVEL
C
          sumu=0.0
          sumv=0.0
          sumwt=0.0
          DO 175 jd=1,numdop
C
C  SKIP WHEN LOOKING AT SAME SOUNDING BEING INTERPOLATED FOR.
C
             IF (jd.NE.jsond) THEN
            IF (levhi(jd).GE.il) THEN
               weight=wndwt(xx,yy,xdop(jd),ydop(jd))
               sumwt=sumwt+weight
               sumu=sumu+weight*usig(jd,il)
               sumv=sumv+weight*vsig(jd,il)
            END IF
         END IF
175       CONTINUE
C
          IF (sumwt.GT.0.0) THEN
             usig(jsond,il)=sumu/sumwt
             vsig(jsond,il)=sumv/sumwt
          ELSE
         WRITE (*,*) 'CHECK DOPSIG NEAR 180'
         stop
          END IF
180    CONTINUE
200     CONTINUE
C
C  IF NO SOUNDING GOES AS HIGH AS UPPERMOST LEVEL, EXTRAPOLATE
C  UPPERMOST VALUESTO TOP OF DOMAIN.
C
    IF (levtop.LT.nlvl) THEN
       DO 225 jsond=1,numdop
          DO 220 il=levtop+1,nlvl
             usig(jsond,il)=usig(jsond,levtop)
             vsig(jsond,il)=vsig(jsond,levtop)
220       CONTINUE
225    CONTINUE
    END IF
C
C  ALL MISSING UPPER WIND DATA HAVE BEEN ESTIMATED.
C
        RETURN
C
        END