strat.f

Go to the documentation of this file.

C**********************************************************************
         SUBROUTINE strat
C**********************************************************************
C
C  READS TEMPERATURE  SOUNDING INFORMATION,CALCULATES LAPSE RATES & 
C  OTHER PARAMETERS FOR VARIOUS FLOW LEVELS. 
C   MODIFIED BY LUDWIG, JAN 2002.
C
          include 'NGRIDS.PAR'
C
          include 'ANCHOR.CMM'
          include 'FLOWER.CMM'
          include 'LIMITS.CMM'
          include 'STALOC.CMM'
          include 'TSONDS.CMM'
C
    REAL RHS1(ntsite,nzgrd),SYTHYT(ntsite),PMIDS(nsndht,ntsite)
    INTEGER LOWZEE(ntsite,nzgrd),NXTZEE(ntsite,nzgrd)
    INTEGER LEVHI(ntsite),LEVTOP,IXRASS(2),IYRASS(2)
C
C  INTEGER COORDINATES OF RASS SITES -- NOT USED FOR UPPER PRESSURE
C
    DATA ixrass,iyrass /423,432,4497,4499/
C
        IF(numtmp.LE.0) THEN
       WRITE (*,*) 'NO TEMP SONDE DATA -- CANNOT RUN'
       stop
    END IF
C
    levtop=0
C
    DO 100 jtsond=1,numtmp
C
C  FILL ARRAYS WITH MSG DATA IDENTIFIER
C
       DO 8 i=1,nsndht
              zsnd(i,jtsond)=-9999.
              tsnd(i,jtsond)=-9999.
              psnd(i,jtsond)=-9999.
          altmtr(i,jtsond)=-9999.
              zmids(i,jtsond)=-9999.
              pmids(i,jtsond)=-9999.
              dptdzs(i,jtsond)=-9999.
              potemp(i,jtsond)=-9999.
8      CONTINUE
           ii=jtsond
C
C  READ  STATION NO., NO. OF HEIGHTS (NHITES) IN SOUNDING & DATA INPUT 
C  TYPE (ITYP) DATA ARE READ AS FOLLOWS FOR ITYP= :
C       1--HEIGHT(M--MSL),TEMPERATURE(C), POTENTIAL TEMPERATURE
C          LAPSE RATE(K/M) -- THIS OPTION NOT AVAILABLE HERE
C       2--HEIGHT(M--MSL), TEMP(C) & PRESSURE(MB)
C       3--HEIGHT(IN FT. MSL), TEMP(C), PRESSURE (MB)
C
C  SKIP BLANK LINE
C
           READ(12,*) 
           READ(12,*) xs,ys,sythyt(jtsond)
       IF ((nint(xs) .EQ. ixrass(1) .AND. nint(ys).EQ.iyrass(1))
     $          .OR. (nint(xs).EQ.ixrass(2)  
     $                            .AND. nint(ys).EQ.iyrass(2)) ) THEN
          goodt(jtsond)=.false.
       ELSE
          goodt(jtsond)=.true.
       END IF
       xtmp(jtsond)=1.0+(xs-xorig)/dscrs
           ytmp(jtsond)=1.0+(ys-yorig)/dscrs
           READ(12,*) nthts(jtsond),ityp
C
C  SKIP COLUMN HEADERS
C
           READ(12,*) 
C
           nhites=nthts(jtsond)
       IF (nhites .GT.0) THEN
C
          i=0
C
          DO 10 j=1,nhites
                 READ(12,*,end=186) z,t,p
                 IF(nint(z) .LE.-9999 .OR. nint(t).LE.-9999 
     $                .OR. nint(p) .LE.-9999 .OR.   
     $                    z-sythyt(jtsond).LT. -1.0) GO TO 10    
                 i=i+1
             IF (ityp.GE.2) THEN
C
C  CONVERT FEET TO METERS & T TO POT TMP, IF REQUIRED (ITYP=3)
C  CONVERT TO HEIGHT ABOVE SFC-- 
C
                IF (ityp .EQ. 3 ) z=0.3048*z
                IF (i .LE. nsndht) THEN
                       zsnd(i,jtsond)=z-sythyt(jtsond)
                       potemp(i,jtsond)=tpot(p,t)
                       tsnd(i,jtsond)=t+273.13
                       psnd(i,jtsond)=p
               altmtr(i,jtsond)=altset(z,p)
               IF (altmtr(i,jtsond).LT. 28.) WRITE(*,*)
     $                    i,' CHECK T-SONDE ',jtsond,z,p,altset(z,p)
                ELSE
                       zsnd(nsndht,jtsond)=z-sythyt(jtsond)
                       potemp(nsndht,jtsond)=tpot(p,t)
                       tsnd(nsndht,jtsond)=t+273.13
                       psnd(nsndht,jtsond)=p
               altmtr(i,jtsond)=altset(z,p)
               IF (altmtr(i,jtsond).LT. 28.) WRITE(*,*)
     $                    i,' CHECK T-SONDE ',jtsond,z,p,altset(z,p)
                END IF
             ELSE IF (ityp.EQ.1) THEN
C
C  REMOVED LAPSE RATE INPUT OPTION FOR SLC VTMX APPLICATION
C    FLUDWIG 2/02
C*********************************************************
            WRITE (*,*) 'TEMP SOUNDINGS MUST PROVIDE HEIGHT ',
     $                  '(M, TYP2 OR FT, TYP3), TEMP (DEG C) & ',
     $                  'PRESS (MB/HP)'
                    stop
                 END IF
C
10        CONTINUE
C
C  EMPTY SOUNDING
C      
       ELSE
          WRITE (*,*) 'TEMPERATURE SONDE ', jtsond, ' HAS NO OBS.'
          stop
       END IF
186    nhites=i
       nthts(jtsond)=i
C
C  SOUNDING READ; FIND HEIGHTS OF FLOW SURFACES
C  FIRST, FIND NEAREST LOW REFERENCE POINT 
C
       idis2=(nint(xtmp(jtsond))-lowix(1))**2+
     $                    (nint(ytmp(jtsond))-lowiy(1))**2
       jnear(jtsond)=1
       DO 17 jd=2,5
          id2here=(nint(xtmp(jtsond))-lowix(jd))**2+
     $                     (nint(ytmp(jtsond))-lowiy(jd))**2
          IF (id2here .LT. idis2) THEN
                 idis2=id2here
                 jnear(jtsond)=jd
              END IF
17     CONTINUE
C
       lllx=lowix(jnear(jtsond))
       llly=lowiy(jnear(jtsond))
C
C  CHECK TO SEE IF TOP OBSERVATION IS ABOVE THIS SURFACE.  RECORD 
C  HIGHEST LEVELS REACHED FOR INDIVIDUAL SONDES AND ALL SONDES.
C 
       stoplev=zsnd(nthts(jtsond),jtsond)
       DO 25 kl=1,nlvl
              zlevkl=rhs(lllx,llly,kl)
          rhs1(jtsond,kl)=zlevkl
C
C  FINDING HIGHEST FLOW SURFACE REACHED BY THIS SONDE
C         
          IF (stoplev .GE. zlevkl) THEN
             levhi(jtsond)=kl
         IF (kl .GT. levtop) levtop=kl
          END IF
25     CONTINUE
C      
C  GETTING INDICES OF OBS THAT ARE NEAREST MIDPOINTS ABOVE & BELOW
C  THE FLOW SURFACES.
C
       lowzee(jtsond,1)=1
       DO 45 kl=1,levhi(jtsond)
          zlevkl=rhs1(jtsond,kl)
          IF (kl.LT.nlvl) THEN
                 zbar=0.5*(zlevkl+rhs1(jtsond,kl+1))
          ELSE IF (kl.EQ.nlvl) THEN
                 zbar=zlevkl+0.5*(zlevkl-rhs1(jtsond,kl-1))
          END IF
          close1=2000000.
          DO 40 iz=1,nthts(jtsond)-1
         IF (abs(zbar-zsnd(iz,jtsond)).LT.close1) THEN
            close1=abs(zsnd(iz,jtsond)-zbar)
            nxtzee(jtsond,kl)=iz
            IF (nxtzee(jtsond,kl).EQ.
     $                  lowzee(jtsond,kl)) nxtzee(jtsond,kl)=iz+1
C
C  LOWER POINT FOR NEXT SURFACE WILL BE SAME AS UPPER POINT FOR THIS 
C  ONE
C
            IF (kl .LT. nlvl) THEN
               lowzee(jtsond,kl+1)=iz
            END IF
             END IF
40        CONTINUE
45     CONTINUE
C
C  CHECKING TO SEE THAT FLOW SFCS ARE BETWEEN BOUNDS WHEN POSSIBLE.
C
       DO 70 kl=2,levhi(jtsond)
          IF (zsnd(lowzee(jtsond,kl),jtsond) .GE.  
     $            rhs1(jtsond,kl) .AND. lowzee(jtsond,kl) .GT. 1) THEN
C
C  PICK 1ST PLACE YOU COME TO BELOW THE PROPER MIDPOINT
C
             DO 60 iz=lowzee(jtsond,kl),1,-1
                zbar=0.5*(rhs1(jtsond,kl-1)+rhs1(jtsond,kl))
                IF (zsnd(iz,jtsond) .LT. zbar) THEN
               lowzee(jtsond,kl)=iz
               GO TO 70
            END IF
60           CONTINUE
          END IF
70     CONTINUE
C
C  FINALLY, CHECK TO SEE IF TOP PT IS HIGH ENOUGH (ABOVE MIDPT TO
C  NEXT SURFACE) TO EXTRAPOLATE TO ONE MORE LEVEL
C
       IF (levhi(jtsond).LT.nlvl) THEN
          zbar=0.5*(rhs1(jtsond,levhi(jtsond))+
     $                        rhs1(jtsond,levhi(jtsond)+1))
          IF (zsnd(nthts(jtsond),jtsond) .GE. zbar) THEN
             levhi(jtsond)=1+levhi(jtsond)
         IF (levhi(jtsond).GT.levtop) levtop=levhi(jtsond)
             lowzee(jtsond,levhi(jtsond))=
     $                   (lowzee(jtsond,levhi(jtsond)-1) + 
     $                    nxtzee(jtsond,levhi(jtsond)-1))/2 
             nxtzee(jtsond,levhi(jtsond))=nthts(jtsond)
          END IF
       END IF
100 CONTINUE
C
C  ALL GOOD DATA READ; CONVERT WHERE NECESSARY & FIND HEIGHTS NEAREST 
C  THE POINTS TO USE FOR ESTIMATING GRADIENTSAT FLOW SURFACE HEIGHTS 
C  HAVE BEEN FOUND.  NOW, GET POTENTIAL TEMP LAPSE RATES TO TOP OF 
C  EACH SOUNDING AND ESTIMATE LAPSE RATES ON THE FLOW SFCS.
C
    DO 170 jtsond=1,numtmp
       DO 160 kl=levhi(jtsond),1,-1
C
C  USE OBS UP TO THE HIGHEST OB. INTERPOLATE STD ATMOS HTS 
C  CORRESPONDING TO OBSERVED PRESSURES -- THEN GET PRESSURE 
C  CORRSPONDING TO THE INTERPOLATED HT.  OTHER VARIABLES ARE OBTAINED 
C  BY SIMPLE INTERP, OR FINITE DIFF ESIMATE.
C
          IF (kl.GT.1) THEN
                 zz0=zsnd(lowzee(jtsond,kl),jtsond)+sythyt(jtsond)
                 z1=zsnd(nxtzee(jtsond,kl),jtsond)+sythyt(jtsond)
             zlvl=rhs1(jtsond,kl)+sythyt(jtsond)
         aaa1=altmtr(nxtzee(jtsond,kl),jtsond)
         aaa0=altmtr(lowzee(jtsond,kl),jtsond)
             tt0=tsnd(lowzee(jtsond,kl),jtsond)
             tt1=tsnd(nxtzee(jtsond,kl),jtsond)
             pt0=potemp(lowzee(jtsond,kl),jtsond)
             pt1=potemp(nxtzee(jtsond,kl),jtsond)
             ptlaps(jtsond,kl)=(pt1-pt0)/(z1-zz0)
         altsig(jtsond,kl)=xlintr(aaa0,aaa1,zlvl,zz0,z1)
         psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
             ptsigl(jtsond,kl)=xlintr(pt0,pt1,zlvl,zz0,z1)
             tsigl(jtsond,kl)=
     $                     pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
C
          ELSE IF (kl.EQ.1) THEN
C
C  EXTRAPOLATE TO 10 M LEVEL
C
             z1=zsnd(2,jtsond)+sythyt(jtsond)
             zz0=zsnd(1,jtsond)+sythyt(jtsond)
             zlvl=sythyt(jtsond)
             ptlaps(jtsond,kl)=
     $               (potemp(2,jtsond)-potemp(1,jtsond))/(z1-zz0)
         IF (abs(zlvl-zz0).LT.0.1) THEN
            altsig(jtsond,kl)=altmtr(1,jtsond)
                ptsigl(jtsond,kl)=potemp(1,jtsond)
            tsigl(jtsond,kl)=tsnd(1,jtsond)
            psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
         ELSE IF (abs(zlvl-z1) .LT. 0.1) THEN
            altsig(jtsond,kl)=altmtr(2,jtsond)
                ptsigl(jtsond,kl)=potemp(2,jtsond)
            psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
            tsigl(jtsond,kl)=
     $                     pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
         ELSE
                    pt1=potemp(2,jtsond)  
                    pt0=potemp(1,jtsond)
                    tt1=tsnd(2,jtsond)
                    tt0=tsnd(1,jtsond)
            aaa1=altmtr(2,jtsond)
            aaa0=altmtr(1,jtsond)
            altsig(jtsond,kl)=xlintr(aaa0,aaa1,zlvl,zz0,z1)
            ptsigl(jtsond,kl)=xlintr(pt0,pt1,zlvl,zz0,z1)
            psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
                tsigl(jtsond,kl)=
     $                   pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
         END IF
              END IF
C
160    CONTINUE
170 CONTINUE
C
C  NOW GET THE VALUES ABOVE THOSE OBSERVED AT EACH SITE, BUT BELOW
C  THE HIGHEST LEVEL OBSERVED ANYWHERE.
C
    DO 250 jtsond=1,numtmp
       IF (levhi(jtsond)+1.LE.levtop) THEN
          DO 240 kl=levhi(jtsond)+1,levtop
C
C  INVERSE DISTANCE HORIZ INTERP FROM THE HIGHEST OB HERE TO 
C  HIGHEST OB ANYWHERE
C
                 xx=xtmp(jtsond)
             yy=ytmp(jtsond)
         sumpt=0.0
             sumdpt=0.0
             sumwt=0.0
         sumalt=0.0
             DO 210 jd=1,numtmp
C
C  SKIP WHEN LOOKING AT SAME SOUNDING BEING INTERPOLATED FOR.
C
                IF (jd.NE.jtsond .AND. levhi(jd).GE.kl) THEN
               weight=wndwt(xx,yy,xtmp(jd),ytmp(jd))
               sumwt=sumwt+weight
               sumdpt=sumdpt+weight*ptlaps(jd,kl)
               sumpt=sumt+weight*ptsigl(jd,kl)
C
C  USE INTERPOLATED ALTIMETER SETTING TO GET ESTIMATE TO BE USED
C  FOR PRESSURE.
C
               sumalt=sumalt+weight*altsig(jd,kl)
            END IF
210          CONTINUE
C
             IF (sumwt.GT.0.0) THEN
                ptlaps(jtsond,kl)=sumdpt/sumwt
            altsig(jtsond,kl)=sumalt/sumwt
                ptsigl(jtsond,kl)=sumpt/sumwt
            psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),
     $                             rhs1(jtsond,kl)+sythyt(jtsond))
            tsigl(jtsond,kl)=pt2t(psigl(jtsond,kl),
     $                                        ptsigl(jtsond,kl))
C
C  CHECK FOR SUPERADIABATIC LAPSE RATE
C
                IF (ptsigl(jtsond,kl).LT.ptsigl(jtsond,kl-1)) THEN
               ptsigl(jtsond,kl)=ptsigl(jtsond,kl-1)
               tsigl(jtsond,kl)=
     $                      pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
            END IF
C                   
             ELSE
            WRITE (*,*) 'CHECK STRAT NEAR 210'
            stop
             END IF
240       CONTINUE
       END IF
C
250 CONTINUE
C
C  FOR LEVELS ABOVE THE HIGHEST OB, USE LAST VALUE OF LAPSE 
C  BASED ON OBS AND EXTRPOLATE TEMPERATURE.
C
    IF (levtop.LT.nlvl) THEN
       DO 275 jtsond=1,numtmp
C
C  IF NO SONDE REACHES LOWEST LEVEL, USE A MODERATELY STABLE LAPSE
C  RATE OF 0.005 DEGK/M AND SFC ALTIMETER SETTING
C
              IF (levtop.LT.2) THEN
             ptlaps(jtsond,1)=0.005
         levtop=1
          ELSE
             padjust=altsig(jtsond,levtop)
          END IF
C
          DO 265 kl=levtop+1,nlvl
         altsig(jtsond,kl)=altsig(jtsond,levtop)
             psigl(jtsond,kl)=cvt2p(altsig(jtsond,levtop),
     $                           rhs1(jtsond,kl)+sythyt(jtsond))
             ptlaps(jtsond,kl)=ptlaps(jtsond,levtop)
C
C  CHECK FOR SUPERADIABATIC LAPSE RATE
C
         IF (ptlaps(jtsond,kl).LT.0.0) ptlaps(jtsond,kl)=0.0
C
C  ESTIMATE PRESSURE AT LAST LEVEL FROM STAND ATMOS HEIGHT
C  USE IT TO ESTIMATE PTENTIAL TEMP AT THIS LEVEL BY EXTRAPOL
C  FROM PREVIOUS LEVEL.
C
         ptsigl(jtsond,kl)=ptsigl(jtsond,levtop) +       
     $              ((rhs1(jtsond,kl)-rhs1(jtsond,levtop))*
     $                    ptlaps(jtsond,levtop))
C
C NOW CONVERT POTENT TEMP BACK TO TEMP
C
                 tsigl(jtsond,kl)=
     $                 pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
265       CONTINUE
275    CONTINUE
    END IF
C
        RETURN
C
        END