wocss.f

Go to the documentation of this file.

C
C********************************************************************
        PROGRAM  wocss
C********************************************************************
C
C       WINDS ON CRITICAL STREAMLINE SURFACES (WOCSS) PROGRAM FINDS 
C       TOPOGRAPHICALLY INDUCED WINDS BY MAKING THE ORIGINALLY ANALYZED 
C       WINDS NONDIVERGENT WITHIN FLOW SURFACES THAT ARE DEFINED USING 
C       A CONCEPT ANALAGOUS TO THE CRITICAL STREAMLINE WHEN THE 
C       ATMOSPHERE IS STABLY STRATIFIED. FLOW SURFACES CAN INTERSECT 
C       THETERRAIN RESULTING IN ZERO WINDS IN THESE 'UNDERGROUND CELLS.' 
C       ADJUSTMENT TOWARD NONDIVERGENCE IN BAL5 CAUSES FLOW AROUND THE 
C       OBSTACLES.  GEOSIG READS IN WIND SOUNDINGS AND SURFACE DATA.  
C       DIRECT VECTOR ALTERATIONS ARE USED IN SUBROUTINE BAL5. WIND 
C       COMPONENTS AT SPECIFIED POINTS CAN BE HELD CONSTANT OR ADJUSTED 
C       AT A FRACTION OF THE ADJUSTMENTS IN OTHER CELLS IN BAL5.
C
C       THE MODEL HAS BEEN DESCRIBED IN:
C                LUDWIG, F. L., J. M. LIVINGSTON, AND R. M. ENDLICH, 
C                1991:  'USE OF MASS CONSERVATION AND DIVIDING 
C                STREAMLINE CONCEPTS FOR EFFICIENT OBJECTIVE ANALYSIS 
C                OF WINDS IN COMPLEX TERRAIN,' 
C                  J. APPL. METEOROL., VOL. 30, PP. 1490-1499.
C
C       THIS VERSION HAS REMOVED SEVERAL OPTIONS FOUND IN SOME EARLIER 
C       VERSIONS.  ONLY ONE GRID IS USED, NO NESTING.   AT LEAST ONE 
C       SOUNDING IS REQUIRED; THE GEOSTROPHIC WIND PROVISION IS NO 
C       LONGER AVAILABLE.  THE LOW ALTITUDE GRID POINTS TO BE USED
C       FOR DEFINING CRITICAL STREAMLINE WINDS ARE NOW INPUT.  
C
C       NEW PROVISIONS INCLUDE A PARAMETER INPUT FILE SO THAT 
C       ADJUSTMENT FACTORS, FLOW SURFACE COMPRESSION AND OTHER FACTORS 
C       CAN BE CHANGED WITHOUT RECOMPILING.  THE NUMBER OF LABELED 
C       COMMONS HAS ALSO BEEN SUBSTATIALLY REDUCED.
C
C       JANUARY 2002   
C                    F. LUDWIG
C                    ENVIRONMENTAL FLUID MECHANICS LAB
C                    DEPT. OF CIVIL ENGINEERING
C                    STANFORD UNIVERSITY
C                    STANFORD, CA 94305-4020
C
C           BASED IN LARGE PART ON EARLIER WORK WITH:
C           R. ENDLICH, A. BECKER, D. SINTON, K. NITZ, 
C                     J. LIVINGSTON, B. MORLEY AND C. BHUMRALKAR
C
C*********************************************************************
C  ADJMAX           THE FRACTION OF THE USUAL ITERATIVE ADJUSTMENT 
C                   TOWARD NONDIVERGENCE THAT IS MADE AT GRID POINTS 
C                    NEAR OBSERVATIONS IN SUBROUTINE BAL5.
C  AVTHK            HEIGHT AGL OF TOP SURFACE OVER LOWEST TERRAIN.
C  CMPRES           MAXIMUM 'COMPRESSION' OF SURFACES -- 0 MEANS THAT  
C                   LOWER SFCS MUST PARALLEL THOSE ABOVE, 1 MEANS THAT  
C                   THE LOW SFC CAN TOUCH THE NEXT LEVEL -- DEFINES  
C                   INFLUENCE OF UPPER STABLE LAYERS ON LESS STABLE  
C                   LOWER ONES.
C  DBUG,DBUG2       LOGICAL FLAGS TO TRIGGER PRINTING OF STATEMENTS IN  
C                   PROGRAM
C  D2MIN            THE MINIMUM DIDTANCE ALLOWED IN THE INVERSE  
C                   WEIGHTING DENOMINATOR.
C  DPOTMIN      THE MINIMUM ALLOWED POTENT TEMP LAPSE RATE  
C                   (DEG/M) -- LIMITS INSTABILITY AND HENCE RISE OF  
C                   FLOW SFCS.
C  DS               GRID SIZE IN METERS 
C  DSCRS            GRID SIZE (KM)
C  DTWT             DISTANCE WEIGHT POWER --WEIGHT=1/(DISTANCE**DTWT)
C  DZMAX(JT,JZ)     MAXIMUM RISE FOR JZTH FLOW SFC AS DETERMINED FROM  
C                   T-SONDE JT 
C  GRDHI            HEIGHT MSL OF HIGHEST TERRAIN ON GRID.
C  MDATE            DAY OF MONTH
C  KGRIDX           X INDEX OF REFERENCE POINT           
C  KGRIDY           Y INDEX OF REFERENCE POINT
C  LOWIX(JL)        X INDICES OF LOWEST PTS 
C  LOWIY(JL)        Y INDICES OF LOWEST PTS 
C  NCOL             NUMBER OF COLUMNS (X INDEX) IN GRID             
C  NCOLM1           NCOL-1
C  NEND             NUMBER OF HOURS TO CALCULATE BEFORE QUITTING IN A  
C                   SEQUENCE OF DATA)
C  NLVL             NUMBER OF VERTICAL FLOW SURFACES
C  NROW             NUMBER OF ROWS (Y INDEX) IN GRID
C  NROWM1           NROW-1
C  NSITES           MAXIMUM NUMBER OF SITES (INCLUDING UPPER AIR) TO 
C                   BE USED.
C  NSNDHT           MAXIMUM NUMBER SOUNDING (WIND OR TEMP) LEVELS TO  
C                   BE USED.
C  NTSITE           MAXIMUM NUMBER T-SONDES TO BE USED.
C  NUMDOP           NUMBER OF UPPER WIND SITES
C  NUMNWS           NUMBER OF SURFACE WIND SITES
C  NUMTMP           NUMBER OF TEMEPRATURE SOUNDING SITES
C  NUMTOT           TOTAL NUMBER OF OBSERVING WIND SITES
C  NWSITE           MAXIMUM NUMBER WIND SOUNDINGS TO BE USED.
C  NXGRD            DIMENSION FOR GRIDS IN W-E DIRECTION.
C  NYGRD            DIMENSION FOR GRIDS IN S-N DIRECTION.
C  NZGRD            DIMENSION FOR FLOW SURFACES.
C  PWR          POWER TO WHICH SEPARATION IS RAISED FOR INVERSE  
C                   DISTANCE WTS.
C  RHS(JX,JY,JZ)    HEIGHT (M) ABOVE TERRAIN OF JZTH FLOW SFC ABOVE  
C                   PT JX,JY
C  RHSLO(JT,JZ)     HEIGHT (M) OF JZTH FLOW SFC ABOVE THE LOWEST  
C                   TERRAIN NEAR T-SONDE JT  
C  SFCHT(JX,JY)     TERRAIN HEIGHT (M MSL) AT PT JX,JY
C  SFCLOW           TERRAIN HEIGHT (M MSL) OF LOWEST PT IN THE DOMAIN
C  SIGMA(JZ)        FRACTION (OVER LOWEST TERRAIN) OF THE HEIGHT OF 
C                   THE TOP SURFACE FOR SFC JZ
C  SLFAC            CONTROLS THE DEGREE TO WHICH THE 1ST GUESS SURFACES  
C                   FOLLOW THE TERRAIN -- 0 GIVES FLAT SFCS & 1 GIVES  
C                   TERRAIN-FOLLOWING SFCS.
C  SPDCNV           CONVERSION FACTOR -- INPUT UNITS TO M/S
C  U(JX,JY,JZ)      WESTERLY COMPONENT (M/S) AT (JX,JY,JZ)
C  UCOMP(JOB)       OBSERVED WESTERLY COMPONENT (M/S) AT SITE JOB
C  USIG(JWS,JZ)     U COMPONENT AT FLOW SFC JZ ON WIND SOUNDING JWS
C  UTMAPX           UTM EASTING COORDINATE (KM) OF REFERENCE PT  
C                   (KGRIDX,KGRIDY).
C  UTMAPY           UTM NORTHING COORDINATE (KM) OF REFERENCE PT  
C                   (KGRIDX,KGRIDY).
C  V(JX,JY,JZ)      SOUTHERLY COMPONENT (M/S) AT (JX,JY,JZ)
C  VCOMP(JOB)        OBSERVED SOUTHERLY COMPONENT (M/S) AT SITE JOB
C  VSIG(JWS,JZ)     V COMPONENT AT FLOW SFC JZ ON WIND SOUNDING JWS
C  W(JX,JY,JZ)      UPWARD COMPONENT (M/S) AT (JX,JY,JZ)
C  XG(JSITE)        COORDINATE (GRID UNITS) OF OBSERVATION SITE JSITE.
C  XORIG            UTM EASTING COORDINATE (KM) OF THE GRID ORIGIN. 
C  YG(JSITE)        COORDINATE (GRID UNITS) OF OBSERVATION SITE JSITE.
C  YORIG            UTM NORTHING COORDINATE (KM) OF THE GRID ORIGIN.
C                       (NOTE: XORIG,YORIG ARE AT 1,1)
C  ZRISE            DIFFERENCE BETWEEN SFCHI -SFCLOW.
C  ZZERO            ROUGHNESS LENGTH (M) FOR LOG PROFILE ESTIMATES.    
C
    include 'NGRIDS.PAR'
C
    include 'ANCHOR.CMM'
    include 'FLOWER.CMM'
    include 'LIMITS.CMM'
    include 'STALOC.CMM'
C
        CHARACTER*1 CHGRID
        CHARACTER*6 CHDDHR
    CHARACTER*12 CHFNAM
    CHARACTER*8 CHMMDDTTTT
C
    INTEGER*4 NREAD,KEY
C
        DATA key /16/
C
        DATA rhslo /narray*0.0/
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C  MAKE SURE TO SET PATH DEFINITION APPROPRIATE TO YOUR JOB
C  STRING LENGTH FOR CHPATH IS DEFINED IN INCLUDE FILE, NGRIDS.PAR'
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
    chpath='slcin/'
C
C  GRID SIZE SELECTION OPTION DISABLED BY ASSIGNING A 1 KM MESH SIZE
C
    chmesh='1'
C
C       SET PATH & OPEN INPUT PARAMETER FILE
C
        CALL runrd (nit,kfyl1,kfyl2)
    nread=0
C
C ****************DISABLED******************************
C  OPEN FILES FOR KEEPING TRACK OF WHAT HAPPENED WHEN & TEST OUTPUT.
C
C   OPEN(33,FILE='CHPOUT//'TEST'//CHMESH//'KM',
C     $      FORM='FORMATTED', STATUS='UNKNOWN')
C ****************DISABLED******************************
C
    WRITE (*,*) 'STARTING'
C
C    OPEN FILE 11 (TOPOGRAPHY)
C    ASSUMED SQUARE GRID, DEFAULT 80 x 80 FOR VTMX CASE
C
    lgrd=nint(10.0*dscrs)
    IF (lgrd.EQ.10 .OR. lgrd.EQ.20 .OR. lgrd.EQ.50) THEN
       chgrid=char(nint(48+dscrs))
    ELSE
       chgrid='H'
    END IF
C
    WRITE (*,*)'OPENING '//chpath//'SLC1KM.DAT'
C
        OPEN (11,file=chpath//'SLC1KM.DAT',status='OLD',
     $         form='FORMATTED')
C
C  READ TERRAIN HEIGHTS FOR ALL GRIDS USING TOPO
C
    CALL topo
        ds=dscrs*1.0e3
        ncolm1=ncol-1
        nrowm1=nrow-1
C
    CLOSE (11)
C
        OPEN(11,file=chpath//'SLCFILES',status='OLD',
     $         form='FORMATTED')
C
    nread=0
C
44  READ (11,6003,end=8586) chfnam
    chmmddtttt=chfnam(1:8)
C
    WRITE (*,*) 'OPENING '//chpath//chfnam
C
    OPEN (12,file=chpath//chfnam,status='OLD',
     $         form='FORMATTED')
C
C* READ AND ANALYZE WIND DATA USING WXANAL
C* MAKE INITIAL WIND ANALYSIS ON MESH
C
    READ (12,6001,end=8686) ihour,imin,imo,mdate,iyear
C
    nread=nread+1
    IF (nread .GT. nend)  GO TO 8586
C
    khr=ihour
    kmon=imo
    kday=mdate
    lastmo=imo
    minutz=julmin(iyear,imo,mdate,ihour,imin)
C
C
C  GEOSIG READS METEOROLOGICAL INPUTS & CALLS ROUTINES FOR
C  INTERPOLATION AND DEFINING SURFACE SHAPES.
C
    CALL geosig
C
C CALL SUBROUTINE TO MAKE WINDS NONDIVERGENT
C
    CALL bal5(nit)
C
C  IF NEEDED, GET PRESSURES AND TEMPERATURES OON FLOW SURFACES.
C
    IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) CALL betpar
C
C  GET PARAMETER VALUES ON FLAT SURFACES BY INTERPOLATIO IN Z.
C
        CALL tstwnd(chmmddtttt)
C
C  IF REQUESTED, CALCULATE BULK RICHARDSON NUMBER AND BRUNT-VAISALA 
C  PERIOD ON PLOTTING SURFACES.
C
    IF (dobri .OR. dobvpd) CALL getrich
C
    chddhr=chmmddtttt(3:8)
C
C  WRITE OUTPUT FILES
C
    CALL putout (chddhr,chfnam,minutz,nread)
C
    GO TO 44
C
8586    CONTINUE
C
8686    CONTINUE
C
    IF (nread .GT. nend)  GO TO 8786
C
8786    IF (nread .GT. nend) nread=nend
    IF (dowvsz) CALL putout (chddhr,chfnam,minutz,-nread)
    WRITE (*,*) 'FINISHED '
C
6001    FORMAT (4i2,1x,i4)  
6003    FORMAT (a12)
C
    stop
C
    END