!-------------------- pdi_v2.f90 -------------------------------- PROGRAM PALMER ! ! This program computes the Palmer Drought Severity Index ! (PDSI), the Palmer Hydrologic Drought Index (PHDI), the ! real-time operational Palmer Drought Index (WPLM), and the Z ! Index. It is a modified version of the program in use at ! the National Climatic Data Center; the old version was written ! for mainframe computers with tape drives. Code taken from the ! old version is in upper case; new code is in lower case. This ! version is Fortran F77 compatible and is designed for disk rather ! than tape input/output/storage. As written, it accepts climatic ! division temperature and precipitation data obtainable from the ! NCDC web site under the ftp address: ! ! ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv ! ! This version of the program was written by Ned Guttman and ! modified per Neds suggestions by Richard Heim. ! ! It is a flexible program in that parameter settings control the ! input, output, years of analysis, etc. These parameter settings ! are set within the program and are not input from a read statement. ! Changing the settings therefore requires editing the parameter ! statement in the program and then recompiling. Note that the required ! input are temperature and precipitation data that are serially ! complete from the beginning to the end of the analysis period ! of record. The temperature and precipitation data must be in ! the same sort order. There is a simple option for filling in ! missing data with constant monthly values; the user can modify ! this routine as needed for serially completing the data record. ! Another required input is the soils constant file which contains ! under soil layer water capacity, Thornthwaite coefficients for ! computing evapotranspiration, and the negative tangent of the ! latitude of the location. ! ! A description of the parameter settings follows: ! ! nbcaly, necaly beginning, ending years of the period for ! which the calibration, or "climatologically ! appropriate" average values are computed ! nbegyr, nendyr beginning, ending years for the analysis ! nbcaly and necaly must be contained within ! nbegyr and nendyr ! nendmo last month in nendyr for which data are ! available ! tmsng, pmsng missing value code for temp, prec input files ! amsng missing value code for all output ! inunit basic input unit device number from which all ! input and output devices are determined ! output codes 0 = no, 1 = yes; all output is year/month array ! in the format identification, year, 12 monthly ! values ! output file (device) ! ioutp (iout21) precip used for analysis (should = prec input) ! ioutt (iout20) temp used for analysis (should = temp input) ! ioutsp (iout22) soil water content ! ioutpe (iout23) potential evapotranspiration ! ioutpl (iout24) potential loss ! ioutpr (iout25) potential recharge ! ioutr (iout26) recharge ! iouttl (iout27) total loss ! ioutet (iout28) evapotranspiration ! ioutro (iout29) runoff ! ioutsss (iout30) surface layer water content ! ioutssu (iout31) under layer water content ! iouteff (iout32) effective wetness or dryness in ending a ! dry/wet spell ! ioutpro (iout33) probability a spell has ended ! ioutx1 (iout34) index of incipient wetness (wet spell) ! ioutx2 (iout35) index of incipient dryness (dry spell) ! ioutx3 (iout36) index for established wet or dry spell ! ! iouthat (iout37) phat (CAFEC precipitation) ! ! ioutwp (iout11) operational Palmer (WPLM or PMDI) ! ioutpd (iout12) PDSI ! ioutz (iout13) Z Index ! ioutph (iout14) PHDI ! ioutcp (iout15) evap+recharge+runoff-loss ! ! ioutbyr,iouteyr beginning, ending year for output ! <<< WANT ENDING YEAR FOR OUTPUT TO BE ! NO EARLIER THAN ENDING YEAR FOR ANALYSIS >>> ! ! NOTE that the climatologically appropriate values, coefficents, ! etc. are output automatically to device iout10. If you do not want ! this output, comment out the appropriate write statement to device ! iout10. ! ! -------------------------------------------------------------- ! ! NOTE that the id of the soils constants file must match the id ! of the temp and precip id input. If there is not a match, the ! program will stop with an error message that says "no constants". ! This program is written for the climate division input data; ! if you want to run the program on site data, include the ! appropriate state-division code in the identification of the site. ! The present code calls the state-division code id for temp input, ! jd for precip input, jjd for soils constant input, and kd for ! keeping track of the identification. For site data, an approach is ! to keep the current state-division codes, but append a site ! identification such as station number with variable names istn for ! temp input and jstn for precip input. The input read statements and ! format as well as the output write statements and format will have ! to be changed accordingly. The file of soil constants for climatic ! divisions used in this program are available from NCDC. The user can ! provide his own constants if desired with appropriate modifications ! to the read/format statements. ! ! The device numbers XX (ioutXX, inXX) are coded inunit+XX ! ! ! National Climatic Data Center ! Asheville, NC ! December 1997 !--------------------------------------------------------------------- ! INTEGER :: ISTAT(12,4) real :: prec(200,14),temp(200,14), prec1(200,14),temp1(200,14) integer :: indexj(40),indexm(40) real :: pdat(200,12),spdat(200,12),pedat(200,12),pldat(200,12), & prdat(200,12),rdat(200,12),tldat(200,12),etdat(200,12), & rodat(200,12),tdat(200,12),sssdat(200,12),ssudat(200,12), & tsum(12),psum(12),spsum(12),pesum(12),plsum(12), & prsum(12),rsum(12),tlsum(12),etsum(12),rosum(12), & pdsi(200,13),phdi(200,13),wplm(200,13), & phi(12),days(12),alp(12),bet(12), & gam(12),del(12),trat(12),sabsd(12),akhat(12),ak(12), & cp(200,13),z(200,13),ppr(200,12),px1(200,12), & px2(200,12),px3(200,12),eff(200,12), & sx(96),sx1(40),sx2(40),sx3(40),x(200,12) real :: phatsav(12) real :: amsng, b, cd, cet, cl, cr, cro, d, dbar, dk, dum real :: et, h integer :: iass, idv, iendmo, iendyr, imult, inunit integer :: iout10, iout11, iout12, iout13, iout14, iout15, iout20, iout21 integer :: iout22, iout23, iout24, iout25, iout26, iout27, iout28, iout29 integer :: iout30, iout31, iout32, iout33, iout34, iout35, iout36, iout37 integer :: ioutbyr, ioutcp, iouteff, ioutet, iouteyr, iouthat, ioutp, ioutpd, ioutpe integer :: ioutph, ioutpl, ioutpr, ioutpro, ioutr, ioutro, ioutsp, ioutsss, ioutssu integer :: ioutt, iouttl, ioutwp, ioutx1, ioutx2, ioutx3, ioutz integer :: j, k integer :: ist, iy, iyrpr integer :: jbcaly, jbegyr, jj, joutbyr, k8, k8max, lmnpr, lyear, m, mo integer :: nbcaly, nbegyr, necaly, nendmo, nendyr, numyr real :: p, pe, phat, pl, pr, pro, pv, q, r, ro, rs, ru, sl, sp, ss, sss real :: ssu, su, swtd, t, tl, tla, ud, ul, uw, v, wcbot, wctop, wctot real :: x1, x2, x3, ze integer :: iok character :: inputpath*100, metadir*100, pdireg*5, cfmt118*21, cfmt4321*21, soilfmt*13 DATA DAYS /31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./ DATA PHI /-.3865982,-.2316132,-.0378180, .1715539, .3458803, & .4308320, .3916645, .2452467, .0535511,-.15583436, & -.3340551,-.4310691 / !************************************************************ !************************************************************ !---GET STATUS INFO FOR PUERTO RICO & OTHER TERRITORIES ! OPEN (4,FILE='prolog.status',STATUS='OLD') call getenv("SRC",inputpath) call getenv("metadir",metadir) call getenv("pdireg",pdireg) inputpath = trim(inputpath)//'/input/' metadir = trim(metadir)//'/' IYRPR=1998 LMNPR=6 !401 READ (4,404,END=499) IY,((ISTAT(M,I),I=1,4),M=1,12) !404 FORMAT (1X,I4,12(1X,4I1)) ! DO 410 M=1,12 ! IF (ISTAT(M,2).EQ.0 .AND. ISTAT(M,3).EQ.0 .AND. ! * ISTAT(M,4).EQ.0) GO TO 499 ! IYRPR=IY ! LMNPR=M !410 CONTINUE ! GO TO 401 !499 CLOSE (4) !************************************************************ !************************************************************ ! set input/output device assignments ! read the latest month to process from ../metadata/control.status CALL COMPIM (IMULT,nendyr,nendmo) iendyr=nendyr iendmo=nendmo ! The pdinew.parameter file is the same for each run. call load_params (nbcaly, necaly, nbegyr, ioutbyr, iouteyr, amsng, & inunit, ioutp, ioutt, ioutsp, ioutpe, ioutpl, ioutpr, ioutr, iouttl, & ioutet, ioutro, ioutsss, ioutssu, ioutwp, ioutpd, ioutz, ioutph, & ioutcp, iouteff, ioutpro, ioutx1,ioutx2, ioutx3, iouthat) jbcaly=nbcaly jbegyr=nbegyr joutbyr=ioutbyr if (iouteyr < nendyr) iouteyr=nendyr open (9, file='processing.log',status='replace') ! open (7, file=trim(inputpath)//'pdinew.parameter',status='old') write (9,*) ' running pdi_v2.exe ' write (6,*) ' running pdi_v2.exe ' write (6,*) ' pdireg = '//trim(pdireg) call open_iouts (inunit, iout10, iout11, iout12, iout13, iout14, iout15, & iout20, iout21, iout22, iout23, iout24, iout25, iout26, iout27, iout28, & iout29, iout30, iout31, iout32, iout33, iout34, iout35, iout36, iout37) ! Assign different formats for counties vs. state/div if ( pdireg == 'cty' ) then soilfmt = "(i2,i3,4f8.4)" cfmt118 = "(i2.2,i3.3,i4,13f7.2)" cfmt4321 = "(i2.2,i3.3,i5,12f7.2)" else soilfmt = "(i3,i2,4f8.4)" cfmt118 = "(i3.3,i2.2,i4,13f7.2)" cfmt4321 = "(i3.3,i2.2,i5,12f7.2)" endif ! The pdinew.soilconst file changes for state, divisional, and ! river basin options. Control of this is handled by the palmer.csh ! script. ! The contents of pdinew.soilconst determines what areas are processed. open (14,file=trim(metadir)//'pdinew.soilconst',status='old') ! beginning of iteration over all regions (states, divisions, counties) do read(14,soilfmt,END=16)IST,IDV,wcbot,b,h,tla !---STATES ALASKA TO THE PACIFIC ISLANDS DO NOT HAVE DATA PRIOR !---TO 1931. SO, NEED TO SET START YEARS TO 1931 FOR THOSE !---STATE CODES. if (IST > 48 .and. IST < 100) then if (nbcaly < 1931) nbcaly=1931 if (nbegyr < 1931) nbegyr=1931 if (ioutbyr < 1931) ioutbyr=1931 nendyr=IYRPR nendmo=LMNPR else nbcaly=jbcaly nbegyr=jbegyr ioutbyr=joutbyr nendyr=iendyr nendmo=iendmo endif IF (IST.EQ.1 .AND. IDV.EQ.1 .and. trim(pdireg) == "div" ) THEN write (9,*) ' processing divisions' write (6,*) ' processing divisions' else IF (IST.EQ.1 .AND. IDV.EQ.1 .and. trim(pdireg) == "cty" ) THEN write (9,*) ' processing counties' write (6,*) ' processing counties' ELSE IF (IST.EQ.1 .AND. IDV.EQ.0) THEN write (9,*) ' processing states & T.K. regions' write (6,*) ' processing states & T.K. regions' ELSE IF (IST.EQ.111 .AND. IDV.EQ.0) THEN write (9,*) ' processing river basins and special regions' write (6,*) ' processing river basins and special regions' ELSE IF (IST.EQ.121 .AND. IDV.EQ.0) THEN write (9,*) ' processing NWS regions' write (6,*) ' processing NWS regions' ENDIF ! initialize variables tsum(:) = 0. ; psum(:) = 0. ; spsum(:) = 0. ; pesum(:) = 0. ; plsum(:) = 0. prsum(:) = 0. ; rsum(:) = 0. ; tlsum(:) = 0. ; etsum(:) = 0. ; rosum(:) = 0. pdsi(:,:)=amsng ; phdi(:,:)=amsng ; z(:,:)=amsng ; wplm(:,:)=amsng cp(:,:)=amsng ; pdat(:,:)=amsng ; spdat(:,:)=amsng ; pedat(:,:)=amsng pldat(:,:)=amsng ; prdat(:,:)=amsng ; rdat(:,:)=amsng ; tldat(:,:)=amsng etdat(:,:)=amsng ; rodat(:,:)=amsng ; tdat(:,:)=amsng ; sssdat(:,:)=amsng ssudat(:,:)=amsng eff(:,:)=0.0 ; PPR(:,:)=0.0 ; PX1(:,:)=0.0 ; PX2(:,:)=0.0 ; PX3(:,:)=0.0 call GETPOR(IST,IDV,1,prec1,iok) if (iok .eq. 0) then write (9,35) IST,IDV write (*,35) IST,IDV 35 format(3X,'***ERROR GETPOR(precip.): IST=',I3,' DIV=',I2) call exit(1) endif call GETPOR(IST,IDV,2,temp1,iok) if (iok .eq. 0) then write (9,836) IST,IDV write (*,836) IST,IDV 836 format(3X,'***ERROR GETPOR(temp.): IST=',I3,' DIV=',I2) call exit(1) endif ! sort through the input, keeping only the years that are within range j = 0 do k=1,200 iy=nint(prec1(k,14)) if (iy < nbegyr) cycle j = j + 1 do m=1,14 prec(j,m)=prec1(k,m) temp(j,m)=temp1(k,m) enddo enddo WCTOP = 1.0 SS = WCTOP SU = WCBOT ! comes from the soil constants file WCTOT = WCBOT + WCTOP ! loop on years (j) and months (m) do 70 j=1,nendyr-nbegyr+1 iy=nINT(prec(j,14)) ! check for leap year LYEAR = (IY/4*4)/IY ! fails on 1800, 1900, 2100, etc. if (mod(iy,100) .eq. 0 .and. mod(iy,400) .ne. 0) lyear=0 do 60 m=1,12 if ((j.eq.nendyr-nbegyr+1).and.(m.gt.nendmo)) goto 60 t=temp(j,m) p=prec(j,m) !----------------------------------------------------------------------- ! HERE START THE WATER BALANCE CALCULATIONS !----------------------------------------------------------------------- SP = SS + SU PR = WCBOT + WCTOP - SP !----------------------------------------------------------------------- ! 1 - CALCULATE PE (POTENTIAL EVAPOTRANSPIRATION) !----------------------------------------------------------------------- IF (T.LE.32.) THEN PE = 0.0 ELSE DUM = PHI(M)*TLA DK = ATAN(SQRT(1. - DUM*DUM)/DUM) IF (DK.LT.0.) DK = 3.141593 + DK DK = (DK + .0157)/1.57 IF (T.GE.80.) THEN PE = (SIN(T/57.3 - .166) - .76)*DK ELSE DUM=log(t-32) PE = (EXP(-3.863233 + B*1.715598 -B*log(h) + B*DUM))*DK ENDIF ENDIF !----------------------------------------------------------------------- ! CONVERT DAILY TO MONTHLY !----------------------------------------------------------------------- IF (M.EQ.2.AND.LYEAR.EQ.1) THEN PE = PE * 29. ELSE PE = PE * DAYS(M) ENDIF !----------------------------------------------------------------------- ! 2 - PL POTENTIAL LOSS !----------------------------------------------------------------------- IF (SS.GE.PE) THEN PL = PE ELSE PL = ((PE - SS) * SU) / (WCBOT + WCTOP) + SS PL = AMIN1(PL,SP) ENDIF !----------------------------------------------------------------------- ! 3 - CALCULATE RECHARGE, RUNOFF, RESIDUAL MOISTURE, LOSS TO BOTH ! SURFACE AND UNDER LAYERS, DEPENDING ON STARTING MOISTURE ! CONTENT AND VALUES OF PRECIPITATION AND EVAPORATION. !----------------------------------------------------------------------- IF (P.GE.PE) THEN ! ----------------- PRECIP EXCEEDS POTENTIAL EVAPORATION ET = PE TL = 0.0 IF ((P - PE).GT.(WCTOP - SS)) THEN ! ------------------------------ EXCESS PRECIP RECHARGES ! UNDER LAYER AS WELL AS UPPER RS = WCTOP - SS SSS = WCTOP IF ((P - PE -RS).LT.(WCBOT - SU)) THEN ! ---------------------------------- BOTH LAYERS CAN TAKE ! THE ENTIRE EXCESS RU = P - PE - RS RO = 0.0 ELSE ! ---------------------------------- SOME RUNOFF OCCURS RU = WCBOT - SU RO = P - PE - RS - RU ENDIF SSU = SU + RU R = RS + RU ELSE ! ------------------------------ ONLY TOP LAYER RECHARGED R = P - PE SSS = SS + P - PE SSU = SU RO = 0.0 ENDIF ELSE ! ----------------- EVAPORATION EXCEEDS PRECIPITATION R = 0.0 IF (SS.GE.(PE - P)) THEN ! ----------------------- EVAP FROM SURFACE LAYER ONLY SL = PE - P SSS = SS - SL UL = 0.0 SSU = SU ELSE ! ----------------------- EVAP FROM BOTH LAYERS SL = SS SSS = 0.0 UL = (PE - P - SL) * SU / (WCTOT) UL = AMIN1(UL,SU) SSU = SU - UL ENDIF TL = SL + UL RO = 0.0 ET = P + SL + UL ENDIF ! -------------------------------------------------------------- ! FOR CALIBRATION YEARS, SUM VALUES NEEDED TO CALCULATE THE ! NORMAL CLIMATE COEFFICIENTS (ALPHA, BETA, ETC.) ! -------------------------------------------------------------- IF(IY.GE.NBCALY.AND.IY.LE.NECALY) THEN TSUM (M) = TSUM (M) + T PSUM (M) = PSUM (M) + P SPSUM(M) = SPSUM(M) + SP PESUM(M) = PESUM(M) + PE PLSUM(M) = PLSUM(M) + PL PRSUM(M) = PRSUM(M) + PR RSUM (M) = RSUM (M) + R TLSUM(M) = TLSUM(M) + TL ETSUM(M) = ETSUM(M) + ET ROSUM(M) = ROSUM(M) + RO ENDIF pdat(j,m)=p spdat(j,m)=sp pedat(j,m)=pe pldat(j,m)=pl prdat(j,m)=pr rdat(j,m)=r tldat(j,m)=tl etdat(j,m)=et rodat(j,m)=ro tdat(j,m)=t sssdat(j,m)=sss ssudat(j,m)=ssu SS = SSS SU = SSU 60 CONTINUE ! -------------------------- END OF MONTH LOOP ! output the intermediate variables if(iy.ge.ioutbyr.and.iy.le.iouteyr)then if(ioutt.eq.1)write(iout20,cfmt4321)IST,IDV,iy,(tdat(j,m),m=1,12) if(ioutp.eq.1)write(iout21,cfmt4321)IST,IDV,iy,(pdat(j,m),m=1,12) if(ioutsp.eq.1)write(iout22,cfmt4321)IST,IDV,iy,(spdat(j,m),m=1,12) if(ioutpe.eq.1)write(iout23,cfmt4321)IST,IDV,iy,(pedat(j,m),m=1,12) if(ioutpl.eq.1)write(iout24,cfmt4321)IST,IDV,iy,(pldat(j,m),m=1,12) if(ioutpr.eq.1)write(iout25,cfmt4321)IST,IDV,iy,(prdat(j,m),m=1,12) if(ioutr.eq.1)write(iout26,cfmt4321)IST,IDV,iy,(rdat(j,m),m=1,12) if(iouttl.eq.1)write(iout27,cfmt4321)IST,IDV,iy,(tldat(j,m),m=1,12) if(ioutet.eq.1)write(iout28,cfmt4321)IST,IDV,iy,(etdat(j,m),m=1,12) if(ioutro.eq.1)write(iout29,cfmt4321)IST,IDV,iy,(rodat(j,m),m=1,12) if(ioutsss.eq.1)write(iout30,cfmt4321)IST,IDV,iy,(sssdat(j,m),m=1,12) if(ioutssu.eq.1)write(iout31,cfmt4321)IST,IDV,iy,(ssudat(j,m),m=1,12) endif write (60,cfmt118) IST,IDV,iy,(pedat(j,m), m=1,12),amsng write (61,cfmt118) IST,IDV,iy,(pldat(j,m), m=1,12),amsng write (62,cfmt118) IST,IDV,iy,(prdat(j,m), m=1,12),amsng write (63,cfmt118) IST,IDV,iy,(etdat(j,m), m=1,12),amsng write (64,cfmt118) IST,IDV,iy,(tldat(j,m), m=1,12),amsng write (65,cfmt118) IST,IDV,iy,(rdat(j,m), m=1,12),amsng write (66,cfmt118) IST,IDV,iy,(rodat(j,m), m=1,12),amsng write (67,cfmt118) IST,IDV,iy,(sssdat(j,m),m=1,12),amsng write (68,cfmt118) IST,IDV,iy,(ssudat(j,m),m=1,12),amsng write (69,cfmt118) IST,IDV,iy,(spdat(j,m), m=1,12),amsng 70 CONTINUE ! ------------------------- END OF YEAR LOOP !----------------------------------------------------------------------- ! LOOP 90 IS BY MONTH: CAFEC COEFFICIENTS CALCULATIONS AND ! 'T' RATIO FOR EACH MONTH !----------------------------------------------------------------------- DO 90 M = 1, 12 ! ! ALPHA CALCULATION ! ----------------- IF(PESUM(M).NE.0.) THEN ALP(M) = ETSUM(M) / PESUM(M) ELSE IF (ETSUM(M).EQ.0.) THEN ALP(M) = 1.0 ELSE ALP(M) = 0.0 ENDIF ENDIF ! ! BETA CALCULATION ! ---------------- IF (PRSUM(M).NE.0.) THEN BET(M) = RSUM(M)/PRSUM(M) ELSE IF (RSUM(M).EQ.0.) THEN BET(M) = 1.0 ELSE BET(M) = 0.0 ENDIF ENDIF ! ! GAMMA CALCULATION ! ----------------- IF (SPSUM(M).NE.0.) THEN GAM(M) = ROSUM(M) / SPSUM(M) ELSE IF (ROSUM(M).EQ.0.) THEN GAM(M) = 1. ELSE GAM(M) = 0.0 ENDIF ENDIF ! ! DELTA CALCULATION ! ----------------- IF (PLSUM(M).NE.0.) THEN DEL(M) = TLSUM(M) / PLSUM(M) ELSE DEL(M) = 0.0 ENDIF !----------------------------------------------------------------------- ! CALCULATION OF Z-INDEX WEIGHTING FACTORS (VARIABLE AK) ! TRAT(M) IS THE 'T' RATIO OF AVERAGE MOISTURE DEMAND TO AVERAGE ! MOISTURE SUPPLY IN MONTH M. !----------------------------------------------------------------------- TRAT(M) = (PESUM(M)+RSUM(M)+ROSUM(M)) / (PSUM(M)+TLSUM(M)) !----------------------------------------------------------------------- ! 'SUM' VARIABLES NOW BECOME AVERAGES OVER THE CALIBRATION PERIOD !----------------------------------------------------------------------- numyr=necaly-nbcaly+1 TSUM (M) = TSUM(M) / NUMYR PSUM (M) = PSUM (M) / NUMYR SPSUM(M) = SPSUM(M) / NUMYR PESUM(M) = PESUM(M) / NUMYR PLSUM(M) = PLSUM(M) / NUMYR PRSUM(M) = PRSUM(M) / NUMYR RSUM (M) = RSUM (M) / NUMYR TLSUM(M) = TLSUM(M) / NUMYR ETSUM(M) = ETSUM(M) / NUMYR ROSUM(M) = ROSUM(M) / NUMYR 90 CONTINUE ! ----------------------END OF MONTHLY COEFFICIENT LOOP SABSD(:) = 0.0 !----------------------------------------------------------------------- ! REREAD MONTHLY PARAMETERS FOR CALCULATION OF ! THE 'K' MONTHLY WEIGHTING FACTORS USED IN Z-INDEX CALCULATION !----------------------------------------------------------------------- do j=nbcaly-nbegyr+1,necaly-nbegyr+1 iy=j+nbegyr-1 do m=1,12 PHAT = ALP(M)*PEdat(j,m) + BET(M)*PRdat(j,m) + & GAM(M)*SPdat(j,m) - DEL(M)*PLdat(j,m) phatsav(m) =PHAT D = Pdat(j,m) - PHAT SABSD(M)= SABSD(M) + ABS(D) enddo ! next month !---save the PHAT values write (iout37,cfmt4321) IST,IDV,iy,(phatsav(m),m=1,12) enddo ! next year SWTD = 0.0 DO MO = 1,12 DBAR = SABSD(MO) / NUMYR AKHAT(MO) = 1.5 * log10((trat(mo)+2.8)/dbar) + .5 SWTD = SWTD + DBAR*AKHAT(MO) enddo !----------------------------------------------------------------------- ! PRINT OUT COEFFICIENTS AND OTHER STUFF for calibration years !----------------------------------------------------------------------- DO MO = 1,12 AK(MO) = 17.67 * AKHAT(MO) / SWTD write(iout10,6060) IST,IDV,MO, ALP(MO), BET(MO), GAM(MO), DEL(MO), & AK(MO),PSUM (MO), SPSUM(MO), PRSUM(MO), PESUM(MO), & PLSUM(MO), ETSUM(MO), RSUM (MO), TLSUM(MO), & ROSUM(MO) enddo 6060 format(i3.3,i2.2,i3,4f9.4,f9.3,9f9.2) CONTINUE ! ?? ok V = 0.0 PRO = 0.0 X1 = 0.0 X2 = 0.0 X3 = 0.0 K8 = 1 k8max = 1 !----------------------------------------------------------------------- ! LOOP FROM 160 TO 230 REREADS data FOR CALCULATION OF ! THE Z-INDEX (MOISTURE ANOMALY) AND PDSI (VARIABLE X). ! THE FINAL OUTPUTS ARE THE VARIABLES PX3, X, AND Z WRITTEN ! TO FILE 11. !----------------------------------------------------------------------- do 231 j=1,nendyr-nbegyr+1 iy=j+nbegyr-1 do 230 m=1,12 if ((j.eq.nendyr-nbegyr+1).and.(m.gt.nendmo)) goto 2311 ! reset indentation indexj(k8)=j indexm(k8)=m ZE = 0.0 UD = 0.0 UW = 0.0 CET = ALP(M) * PEdat(j,m) CR = BET(M) * PRdat(j,m) CRO = GAM(M) * SPdat(j,m) CL = DEL(M) * PLdat(j,m) CP(j,m) = CET + CR + CRO - CL CD = Pdat(j,m)- CP(j,m) Z(j,m) = AK(M) * CD IF (PRO.EQ.100. .OR. PRO.EQ.0.) THEN ! ------------------------------------ NO ABATEMENT UNDERWAY ! WET OR DROUGHT WILL END IF ! -.5 =< X3 =< .5 IF (X3.GE.-.5 .AND. X3.LE..5) THEN ! ---------------------------------- END OF DROUGHT OR WET PV = 0.0 PPR(j,m) = 0.0 PX3(j,m) = 0.0 ! ------------ BUT CHECK FOR NEW WET OR DROUGHT START FIRST GO TO 200 ELSE IF (X3.GT..5) THEN ! ----------------------- WE ARE IN A WET SPELL IF (Z(j,m).GE..15) THEN ! ------------------ THE WET SPELL INTENSIFIES GO TO 210 ELSE ! ------------------ THE WET STARTS TO ABATE (AND MAY END) GO TO 170 ENDIF ELSE IF (X3.LT.-.5) THEN ! ------------------------- WE ARE IN A DROUGHT IF (Z(j,m).LE.-.15) THEN ! -------------------- THE DROUGHT INTENSIFIES GO TO 210 ELSE ! ------------------ THE DROUGHT STARTS TO ABATE (AND MAY END) GO TO 180 ENDIF ENDIF ELSE ! ------------------------------------------ABATEMENT IS UNDERWAY IF (X3.GT.0.) GO TO 170 ! ----------------------- WE ARE IN A WET SPELL IF (X3.LE.0.) GO TO 180 ! ----------------------- WE ARE IN A DROUGHT ENDIF !----------------------------------------------------------------------- ! WET SPELL ABATEMENT IS POSSIBLE !----------------------------------------------------------------------- 170 UD = Z(j,m) - .15 PV = UD + AMIN1(V,0.) IF (PV.GE.0) GO TO 210 ! ---------------------- DURING A WET SPELL, PV => 0 IMPLIES ! PROB(END) HAS RETURNED TO 0 ZE = -2.691 * X3 + 1.5 GO TO 190 !----------------------------------------------------------------------- ! DROUGHT ABATEMENT IS POSSIBLE !----------------------------------------------------------------------- 180 UW = Z(j,m) + .15 PV = UW + AMAX1(V,0.) IF (PV.LE.0) GO TO 210 ! ---------------------- DURING A DROUGHT, PV =< 0 IMPLIES ! PROB(END) HAS RETURNED TO 0 ZE = -2.691 * X3 - 1.5 !----------------------------------------------------------------------- ! PROB(END) = 100 * (V/Q) WHERE: ! V = SUM OF MOISTURE EXCESS OR DEFICIT (UD OR UW) ! DURING CURRENT ABATEMENT PERIOD ! Q = TOTAL MOISTURE ANOMALY REQUIRED TO END THE ! CURRENT DROUGHT OR WET SPELL !----------------------------------------------------------------------- 190 IF (PRO.EQ.100.) THEN ! --------------------- DROUGHT OR WET CONTINUES, CALCULATE ! PROB(END) - VARIABLE ZE Q = ZE ELSE Q = ZE + V ENDIF PPR(j,m) = (PV / Q) * 100. IF (PPR(j,m).GE.100.) THEN PPR(j,m) = 100. PX3(j,m) = 0. ELSE PX3(j,m) = .897 * X3 + Z(j,m)/3. ENDIF !----------------------------------------------------------------------- ! CONTINUE X1 AND X2 CALCULATIONS. ! IF EITHER INDICATES THE START OF A NEW WET OR DROUGHT, ! AND IF THE LAST WET OR DROUGHT HAS ENDED, USE X1 OR X2 ! AS THE NEW X3. !----------------------------------------------------------------------- 200 PX1(j,m) = .897 * X1 + Z(j,m)/3. PX1(j,m) = AMAX1(PX1(j,m),0.) IF (PX1(j,m).GE.1.) THEN IF (PX3(j,m).EQ.0.) THEN ! ------------------- IF NO EXISTING WET SPELL OR DROUGHT ! X1 BECOMES THE NEW X3 X(j,m) = PX1(j,m) PX3(j,m) = PX1(j,m) PX1(j,m) = 0. iass=1 call assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,& j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) GO TO 220 ENDIF ENDIF PX2(j,m) = .897 * X2 + Z(j,m)/3. PX2(j,m) = AMIN1(PX2(j,m),0.) IF (PX2(j,m).LE.-1.) THEN IF (PX3(j,m).EQ.0.) THEN ! ------------------- IF NO EXISTING WET SPELL OR DROUGHT ! X2 BECOMES THE NEW X3 X(j,m) = PX2(j,m) PX3(j,m) = PX2(j,m) PX2(j,m) = 0. iass=2 call assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,& j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) GO TO 220 ENDIF ENDIF IF (PX3(j,m).EQ.0.) THEN ! -------------------- NO ESTABLISHED DROUGHT (WET SPELL), BUT X3=0 ! SO EITHER (NONZERO) X1 OR X2 MUST BE USED AS X3 IF (PX1(j,m).EQ.0.) THEN X(j,m) = PX2(j,m) iass=2 call assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,& j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) GO TO 220 ENDIF IF (PX2(j,m).EQ.0) THEN X(j,m) = PX1(j,m) iass=1 call assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm,& j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) GO TO 220 ENDIF ENDIF !----------------------------------------------------------------------- ! AT THIS POINT THERE IS NO DETERMINED VALUE TO ASSIGN TO X, ! ALL VALUES OF X1, X2, AND X3 ARE SAVED IN FILE 8. AT A LATER ! TIME X3 WILL REACH A VALUE WHERE IT IS THE VALUE OF X (PDSI). ! AT THAT TIME, THE ASSIGN SUBROUTINE BACKTRACKS THROUGH FILE ! 8 CHOOSING THE APPROPRIATE X1 OR X2 TO BE THAT MONTHS X. !----------------------------------------------------------------------- IF (K8.GT.40) then write (9,*) 'pdi_v2.f90 STOP X STORE ARRAYS FULL' write (6,*) 'pdi_v2.f90 STOP X STORE ARRAYS FULL' call exit(1) endif SX1(K8) = PX1(j,m) SX2(K8) = PX2(j,m) SX3(K8) = PX3(j,m) X(j,m) = PX3(j,m) K8 = K8 + 1 k8max=k8 GO TO 220 !----------------------------------------------------------------------- ! PROB(END) RETURNS TO 0. A POSSIBLE ABATEMENT HAS FIZZLED OUT, ! SO WE ACCEPT ALL STORED VALUES OF X3 !----------------------------------------------------------------------- 210 PV = 0.0 PX1(j,m) = 0.0 PX2(j,m) = 0.0 PPR(j,m) = 0.0 PX3(j,m) = .897 * X3 + Z(j,m)/3. X(j,m) = PX3(j,m) IF (K8.EQ.1) THEN PDSI(j,m) = X(j,m) PHDI(j,m) = PX3(j,m) IF ( PX3(j,m).EQ. 0.0 ) PHDI(j,m)=X(j,m) call casepdi (PPR(j,m),PX1(j,m),PX2(j,m),PX3(j,m),WPLM(j,m)) ELSE iass=3 call assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm, & j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) ENDIF !----------------------------------------------------------------------- ! SAVE THIS MONTHS CALCULATED VARIABLES (V,PRO,X1,X2,X3) FOR ! USE WITH NEXT MONTHS DATA !----------------------------------------------------------------------- 220 V = PV eff(j,m)=PV PRO = PPR(j,m) X1 = PX1(j,m) X2 = PX2(j,m) X3 = PX3(j,m) 230 continue ! next month (m) !---store these intermediate variables 2311 continue if (iy.eq.nendyr) then call setmis(j,nendmo,eff,amsng) call setmis(j,nendmo,PPR,amsng) call setmis(j,nendmo,PX1,amsng) call setmis(j,nendmo,PX2,amsng) call setmis(j,nendmo,PX3,amsng) endif if(iy.ge.ioutbyr.and.iy.le.iouteyr)then if(iouteff.eq.1)write(iout32,cfmt4321)IST,IDV,iy,(eff(j,m),m=1,12) if(ioutpro.eq.1)write(iout33,cfmt4321)IST,IDV,iy,(PPR(j,m),m=1,12) if(ioutx1.eq.1)write(iout34,cfmt4321)IST,IDV,iy,(PX1(j,m),m=1,12) if(ioutx2.eq.1)write(iout35,cfmt4321)IST,IDV,iy,(PX2(j,m),m=1,12) if(ioutx3.eq.1)write(iout36,cfmt4321)IST,IDV,iy,(PX3(j,m),m=1,12) endif write (79,cfmt118) IST,IDV,iy,(eff(j,m), m=1,12),amsng write (80,cfmt118) IST,IDV,iy,(PPR(j,m), m=1,12),amsng write (81,cfmt118) IST,IDV,iy,(PX1(j,m), m=1,12),amsng write (82,cfmt118) IST,IDV,iy,(PX2(j,m), m=1,12),amsng write (83,cfmt118) IST,IDV,iy,(PX3(j,m), m=1,12),amsng !------------------------- READ THE NEXT MONTHS DATA 231 continue ! next year (j) !----------------------------------------------------------------------- ! FINISH UP BY CLEARING OUT FILE 8 AND WRITING IT TO FILE 11 !----------------------------------------------------------------------- do k8=1,k8max-1 PDSI(indexj(k8),indexm(k8)) = X(indexj(k8),indexm(k8)) PHDI(indexj(k8),indexm(k8)) = PX3(indexj(k8),indexm(k8)) IF (PX3(indexj(k8),indexm(k8)).EQ. 0.0 ) & PHDI(indexj(K8),indexm(k8))=X(indexj(k8),indexm(k8)) call casepdi(PPR(indexj(k8),indexm(k8)),PX1(indexj(k8),indexm(k8)), & PX2(indexj(k8),indexm(k8)),PX3(indexj(k8),indexm(k8)), & WPLM(indexj(k8),indexm(k8))) enddo ! ! 118 format (I3.3,I2.2,I4,13F7.2) ! JJ=nendyr-nbegyr+1 do j=1,JJ iy=j+nbegyr-1 write (75,cfmt118) IST,IDV,iy,(PDSI(j,k),k=1,12),amsng write (76,cfmt118) IST,IDV,iy,(PHDI(j,k),k=1,12),amsng write (77,cfmt118) IST,IDV,iy,(Z(j,k), k=1,12),amsng write (78,cfmt118) IST,IDV,iy,(WPLM(j,k),k=1,12),amsng write (70,cfmt118) IST,IDV,iy,(cp(j,k), k=1,12),amsng enddo !--------------------------------------------------------------- ! output the indices do j=1,nendyr-nbegyr+1 iy=j+nbegyr-1 if(iy.ge.ioutbyr.and.iy.le.iouteyr.and.ioutwp.eq.1) & write (iout11,cfmt4321) IST,IDV,IY,(WPLM(j,m),m=1,12) if(iy.ge.ioutbyr.and.iy.le.iouteyr.and.ioutpd.eq.1)& write (iout12,cfmt4321) IST,IDV,IY,(pdsi(j,m),m=1,12) if(iy.ge.ioutbyr.and.iy.le.iouteyr.and.ioutz.eq.1)& write (iout13,cfmt4321) IST,IDV,IY,(z(j,m),m=1,12) if(iy.ge.ioutbyr.and.iy.le.iouteyr.and.ioutph.eq.1)& write (iout14,cfmt4321) IST,IDV,IY,(phdi(j,m),m=1,12) if(iy.ge.ioutbyr.and.iy.le.iouteyr.and.ioutcp.eq.1)& write (iout15,cfmt4321) IST,IDV,IY,(cp(j,m),m=1,12) enddo !-------------------------------------------------------------- enddo ! next row in soilconst. file 16 close(14) call close_iouts (inunit) call final_output(pdireg) write (9,*)' Successful completion of pdi_v2.exe' close (9) write (6,*)' Successful completion of pdi_v2.exe' stop end program palmer !----------------------------------------------------------------------- subroutine compim (imult,lyr,lmon) ! COMPUTE THE IMULT MULTIPLIER, WHICH IDENTIFIES IF A NEW YEAR IS ! TO BE ADDED TO THE DATA BASE FILE, AND DETERMINE THE LATEST YEAR ! AND MONTH TO BE UPDATED. ! In other words, imult = -1 if lmon is Jan and imult = 1 otherwise. ! Actually gets this info from control.status ! WRITTEN 2-22-1999 BY RICHARD HEIM JR. ! imult no longer has any use here, but we still need to determine lyr and lmon character :: metadir*100 integer, intent(out) :: imult, lyr, lmon call getenv ("metadir", metadir) !---READ control.status TO DETERMINE LYR AND LMON imult=1 open (2,file=trim(metadir)//'/control.status',status='old') read (2,'(i4,x,i2)') lyr,lmon if(lmon.eq.1) imult=-1 close (2) return end subroutine compim !--------------------------------------------------------------- subroutine load_params (nbcaly, necaly, nbegyr, ioutbyr, iouteyr, amsng, & inunit, ioutp, ioutt, ioutsp, ioutpe, ioutpl, ioutpr, ioutr, iouttl, & ioutet, ioutro, ioutsss, ioutssu, ioutwp, ioutpd, ioutz, ioutph, & ioutcp, iouteff, ioutpro, ioutx1,ioutx2, ioutx3, iouthat) integer, intent(out)::nbcaly, necaly, nbegyr, ioutbyr, iouteyr, & inunit, ioutp, ioutt, ioutsp, ioutpe, ioutpl, ioutpr, ioutr, iouttl, & ioutet, ioutro, ioutsss, ioutssu, ioutwp, ioutpd, ioutz, ioutph, & ioutcp, iouteff, ioutpro, ioutx1,ioutx2, ioutx3, iouthat real, intent(out) :: amsng character :: inputpath*100 call getenv("SRC",inputpath) inputpath = trim(inputpath)//'/input/' open (7, file=trim(inputpath)//'pdinew.parameter',status='old') read (7,'(i4)') nbcaly read (7,'(i4)') necaly read (7,'(i4)') nbegyr read (7,'(i4)') ioutbyr read (7,'(i4)') iouteyr read (7,'(f6.2)') amsng read (7,'(i2)') inunit ! value is usually 10 read (7,'(i1)') ioutp read (7,'(i1)') ioutt read (7,'(i1)') ioutsp read (7,'(i1)') ioutpe read (7,'(i1)') ioutpl read (7,'(i1)') ioutpr read (7,'(i1)') ioutr read (7,'(i1)') iouttl read (7,'(i1)') ioutet read (7,'(i1)') ioutro read (7,'(i1)') ioutsss read (7,'(i1)') ioutssu read (7,'(i1)') ioutwp read (7,'(i1)') ioutpd read (7,'(i1)') ioutz read (7,'(i1)') ioutph read (7,'(i1)') ioutcp read (7,'(i1)') iouteff read (7,'(i1)') ioutpro read (7,'(i1)') ioutx1 read (7,'(i1)') ioutx2 read (7,'(i1)') ioutx3 read (7,'(i1)') iouthat close (7) return end subroutine load_params !----------------------------------------------------------------------- subroutine open_iouts (inunit, iout10, iout11, iout12, iout13, iout14, iout15, & iout20, iout21, iout22, iout23, iout24, iout25, iout26, iout27, iout28, & iout29, iout30, iout31, iout32, iout33, iout34, iout35, iout36, iout37) integer, intent(in) :: inunit integer, intent(out):: iout10, iout11, iout12, iout13, iout14, iout15, & iout20, iout21, iout22, iout23, iout24, iout25, iout26, iout27, iout28, & iout29, iout30, iout31, iout32, iout33, iout34, iout35, iout36, iout37 CHARACTER(len=30) :: name10,name11,name12,name13,name14,name15,& name20,name21,name22,name23,name24,name25,& name26,name27,name28,name29,name30,name31,& name32,name33,name34,name35,name36,name37 integer :: k character :: metadir*100 call getenv("metadir",metadir) metadir = trim(metadir)//'/' ! Read all the files names ! These file names include the relative path ../metadata/ which will ! be removed by putting 12x in the format below. open(7,file=trim(metadir)//'pdinew.outputfiles',status='old') read (7,'(12x,a30)') name10 read (7,'(12x,a30)') name11 read (7,'(12x,a30)') name12 read (7,'(12x,a30)') name13 read (7,'(12x,a30)') name14 read (7,'(12x,a30)') name15 read (7,'(12x,a30)') name20 read (7,'(12x,a30)') name21 read (7,'(12x,a30)') name22 read (7,'(12x,a30)') name23 read (7,'(12x,a30)') name24 read (7,'(12x,a30)') name25 read (7,'(12x,a30)') name26 read (7,'(12x,a30)') name27 read (7,'(12x,a30)') name28 read (7,'(12x,a30)') name29 read (7,'(12x,a30)') name30 read (7,'(12x,a30)') name31 read (7,'(12x,a30)') name32 read (7,'(12x,a30)') name33 read (7,'(12x,a30)') name34 read (7,'(12x,a30)') name35 read (7,'(12x,a30)') name36 read (7,'(12x,a30)') name37 close (7) ! define the output unit numbers iout10=inunit+10 iout11=inunit+11 iout12=inunit+12 iout13=inunit+13 iout14=inunit+14 iout15=inunit+15 iout20=inunit+20 iout21=inunit+21 iout22=inunit+22 iout23=inunit+23 iout24=inunit+24 iout25=inunit+25 iout26=inunit+26 iout27=inunit+27 iout28=inunit+28 iout29=inunit+29 iout30=inunit+30 iout31=inunit+31 iout32=inunit+32 iout33=inunit+33 iout34=inunit+34 iout35=inunit+35 iout36=inunit+36 iout37=inunit+37 ! set file assignments for ioutNN files and open open(unit=iout10,file=trim(metadir)//trim(name10),status='replace') open(unit=iout11,file=trim(metadir)//trim(name11),status='replace') open(unit=iout12,file=trim(metadir)//trim(name12),status='replace') open(unit=iout13,file=trim(metadir)//trim(name13),status='replace') open(unit=iout14,file=trim(metadir)//trim(name14),status='replace') open(unit=iout15,file=trim(metadir)//trim(name15),status='replace') open(unit=iout20,file=trim(metadir)//trim(name20),status='replace') open(unit=iout21,file=trim(metadir)//trim(name21),status='replace') open(unit=iout22,file=trim(metadir)//trim(name22),status='replace') open(unit=iout23,file=trim(metadir)//trim(name23),status='replace') open(unit=iout24,file=trim(metadir)//trim(name24),status='replace') open(unit=iout25,file=trim(metadir)//trim(name25),status='replace') open(unit=iout26,file=trim(metadir)//trim(name26),status='replace') open(unit=iout27,file=trim(metadir)//trim(name27),status='replace') open(unit=iout28,file=trim(metadir)//trim(name28),status='replace') open(unit=iout29,file=trim(metadir)//trim(name29),status='replace') open(unit=iout30,file=trim(metadir)//trim(name30),status='replace') open(unit=iout31,file=trim(metadir)//trim(name31),status='replace') open(unit=iout32,file=trim(metadir)//trim(name32),status='replace') open(unit=iout33,file=trim(metadir)//trim(name33),status='replace') open(unit=iout34,file=trim(metadir)//trim(name34),status='replace') open(unit=iout35,file=trim(metadir)//trim(name35),status='replace') open(unit=iout36,file=trim(metadir)//trim(name36),status='replace') open(unit=iout37,file=trim(metadir)//trim(name37),status='replace') ! Open scratch unit files. ! This should not interfere with all the other open files since inunit=10 do k=60,70 open (k,status='scratch') enddo do k=75,83 open (k,status='scratch') enddo return end subroutine open_iouts !----------------------------------------------------------------- ! Subroutine getpor returns the data for one DRD element for a ST or DV ! matrix passed must be (200,14) - 200 years x 14 values ! Values are monthly (1-12), annual(13) and year(14) the latter as a real. subroutine getpor(nst,ndv,nel,dara,iok) integer, intent(in) :: nst, ndv, nel integer, intent(out) :: iok real, intent(out) :: dara(200,14) integer :: eltnum(14), iyear, miss, i, j, ios character :: stabv*2 character :: ename*4 character :: filen*9 character :: grepfn*43 character :: datadir*100, pdireg*5 character :: cdum*3, cfmt*17 real :: amiss(14) data eltnum / 1, 2, 3, 4, 5, 6, & 7, 8, 9, 11, 12, 24, & 25, 26 / data amiss / -9.99, -99.9, -9999., -9999., -99.99, -99.99, & -99.99, -99.99, -99.9, -99.9, -9.999, -9999., & -9999., -9999. / ! grepfn = 'grep stnumdv /raid3/drd964x/elemnrst.' // call getenv("datadir",datadir) call getenv("pdireg",pdireg) iok = 1 ! grepfn = 'grep stnumdv data/elemnrst. > yearmodahrmmss.elemnr_0' grepfn = 'grep stnumdv data/elemnnxx. > getpor.tmp' ! ,--------------------------------------------------------------------, ! | Get State abbreviation from state number | ! `--------------------------------------------------------------------' if ( trim(pdireg) == "div" .or. trim(pdireg) == "state" ) then call regid(nst,stabv,nel,iok) if (iok .eq. 0) then write (9,*) ' error in regid called from getpor ' write (6,*) ' error in regid called from getpor ' return endif else stabv = 'CC' endif ! ,--------------------------------------------------------------------, ! | Get abbreviated element name | ! `--------------------------------------------------------------------' call elem_abv(nel,ename,iok) if (iok .eq. 0) then write (9,*) ' error in elem_abv called from getpor ' write (6,*) ' error in elem_abv called from getpor ' return endif ! ,--------------------------------------------------------------------, ! | Initialize matrix to missing | ! `--------------------------------------------------------------------' miss = -1 select case (nel) case (55:77) ! spi & palmer model miss=5 case (13:21) ! cir miss=1 case (81:87) ! prob miss=11 case (91:97) ! pacc miss=1 case default do i = 1, 14 if (nel == eltnum(i)) then miss = i exit endif enddo end select if ( miss == -1 ) then write (9,'(a,i3)') ' getpor unable to look-up missing value for ',nel write (6,'(a,i3)') ' getpor unable to look-up missing value for ',nel iok = 0 return endif dara(1:200,1:13) = amiss(miss) dara(1:200,14) = amiss(3) ! ,--------------------------------------------------------------------, ! | Construct input file name | ! `--------------------------------------------------------------------' ! we would need to adjust this for counties select case (trim(pdireg)) case ("div") write (filen(1:9),'(a4,i2.2,a3)') ename, nel, 'dv.' case ("state") write (filen(1:9),'(a4,i2.2,a3)') ename, nel, 'st.' case ("cty") write (filen(1:9),'(a4,i2.2,a3)') ename, nel, 'cy.' case default write (6,*) ' unknown pdireg value:'//pdireg//'.' iok = 0 return end select ! if (ndv == 0) then ! write (filen(1:9),'(a4,i2.2,a3)') ename, nel, 'st.' ! else ! write (filen(1:9),'(a4,i2.2,a3)') ename, nel, 'dv.' ! endif !! The temporary file is just getpor.tmp ! ! ,--------------------------------------------------------------------, ! | Create grep call to get data requested | ! `--------------------------------------------------------------------' ! worksheet for figuring out characters in grepfn ! 111111111122222222223333333333444444444455555555556666666666777 !23456789012345678901234567890123456789012345678901234567890123456789012 !rep stnumdv data/elemnrst. > getpor.tmp' ! we would need to adjust this for counties if ( trim(pdireg) == "div" .or. trim(pdireg) == "state" ) then write (grepfn(6:12),'(a2,i3.3,i2.2)') stabv, nst, ndv else write (grepfn(6:12),'(a2,i2.2,i3.3)') stabv, nst, ndv endif write (grepfn(22:30),'(a9)') filen ! ,--------------------------------------------------------------------, ! | issue grep call then read data | ! `--------------------------------------------------------------------' ! write (6,*) ' getpor issuing system command:' ! write (6,'(a)') grepfn(1:13)//trim(datadir)//grepfn(21:43) call system (grepfn(1:13)//trim(datadir)//grepfn(21:43)) ! Changed unit number from 102 to 87. unit > 100 can be problematic. open (87,file='getpor.tmp',iostat=ios) if ( ios /= 0 ) then write (9,'(a,i3,a)') 'iostat ',ios,' for file getpor.tmp' write (6,'(a,i3,a)') 'iostat ',ios,' for file getpor.tmp' iok = 0 return endif ! determine the actual format of the input data ! (previous format (t10,f4.0,13f7.0) is incorrect but apparently worked) read (87,'(17x,a3)') cdum if (cdum(1:1) == '.') cfmt = '(t10,i4,13f7.2)' if (cdum(2:2) == '.') cfmt = '(t10,i4,13f7.1)' if (cdum(3:3) == '.') cfmt = '(t10,i4,13f7.0)' rewind (87) do i = 1, 200 ! read(87,16,end=18) dara(i,14),(dara(i,j),j=1,13) read(87,cfmt,end=18) iyear, (dara(i,j),j=1,13) dara(i,14) = real(iyear) enddo ! 16 format(t10,f4.0,13f7.0) ! This format is incorrect. It could be f7.2 or f7.1 also. 18 close(87) ! Throw error flag if a low number of rows are read in. if ( i < 60 ) then write (6,'(a,i2,a)') 'Only ',i,' values read from getpor.tmp' iok = 0 return endif call system('\rm -f getpor.tmp') return end subroutine getpor !----------------------------------------------------------------------- !---- subroutine formerly known as assign ------------------------------ subroutine assignpdi (iass,k8,ppr,px1,px2,px3,x,pdsi,phdi,wplm, & j,m,nendyr,nbegyr,sx1,sx2,sx3,sx,indexj,indexm) ! intent ??? real :: SX1(40), SX2(40), SX3(40), SX(96), & pdsi(200,13),phdi(200,13),wplm(200,13),ppr(200,12), & px1(200,12),px2(200,12),px3(200,12),x(200,12) integer :: indexj(40),indexm(40) integer :: iass, isave, k8, j, m, nendyr, nbegyr integer :: k8max, mm, n ! nbegyr and nendyr are not used in this subroutine which yields a ! warning message. Do a trivial operation with them to silence the warning. mm = nendyr + nbegyr ! !----------------------------------------------------------------------- ! FIRST FINISH OFF FILE 8 WITH LATEST VALUES OF PX3, Z,X ! X=PX1 FOR I=1, PX2 FOR I=2, PX3, FOR I=3 !----------------------------------------------------------------------- SX(K8) = X(j,m) ISAVE = iass IF (K8.EQ.1) THEN PDSI(j,m) = X(j,m) PHDI(j,m) = PX3(j,m) IF ( PX3(j,m) .EQ. 0.0 ) PHDI(j,m)=X(j,m) call casepdi (PPR(j,m),PX1(j,m),PX2(j,m),PX3(j,m),WPLM(j,m)) RETURN ELSE ENDIF IF (Iass.EQ.3) THEN ! ---------------- USE ALL X3 VALUES DO 10 Mm = 1, K8-1 10 SX(Mm) = SX3(Mm) GO TO 70 ELSE ! -------------- BACKTRACK THRU ARRAYS, STORING ASSIGNED X1 (OR X2) ! IN SX UNTIL IT IS ZERO, THEN SWITCHING TO THE OTHER ! UNTIL IT IS ZERO, ETC. DO 60 Mm = K8-1, 1, -1 GO TO (20,40) ISAVE 20 IF (SX1(Mm).EQ.0) GO TO 50 30 ISAVE = 1 SX(Mm) = SX1(Mm) GO TO 60 40 IF (SX2(Mm).EQ.0) GO TO 30 50 ISAVE = 2 SX(Mm) = SX2(Mm) 60 CONTINUE ENDIF !----------------------------------------------------------------------- ! PROPER ASSIGNMENTS TO ARRAY SX HAVE BEEN MADE, ! OUTPUT THE MESS !----------------------------------------------------------------------- 70 DO 80 N = 1,K8 PDSI(indexj(n),indexm(n)) = SX(N) PHDI(indexj(n),indexm(n)) = PX3(indexj(n),indexm(n)) IF(PX3(indexj(n),indexm(n)).EQ. 0.0)PHDI(indexj(n),indexm(n))=SX(N) call casepdi (PPR(indexj(n),indexm(n)),PX1(indexj(n),indexm(n)),& PX2(indexj(n),indexm(n)),PX3(indexj(n),indexm(n)),& WPLM(indexj(n),indexm(n))) 80 CONTINUE K8 = 1 k8max=k8 RETURN end subroutine assignpdi !----------------------------------------------------------------------- !----- subroutine formerly known as case -------------------------------- subroutine casepdi (PROB,X1,X2,X3,PALM) real, intent(in) :: prob, x1, x2, x3 real, intent(out):: palm real :: pro ! ! THIS SUBROUTINE SELECTS THE PRELIMINARY (OR NEAR-REAL TIME) ! PALMER DROUGHT SEVERITY INDEX (PDSI) FROM THE GIVEN X VALUES ! DEFINED BELOW AND THE PROBABILITY (PROB) OF ENDING EITHER A ! DROUGHT OR WET SPELL. ! ! X1 - INDEX FOR INCIPIENT WET SPELLS (ALWAYS POSITIVE) ! X2 - INDEX FOR INCIPIENT DRY SPELLS (ALWAYS NEGATIVE) ! X3 - SEVERITY INDEX FOR AN ESTABLISHED WET SPELL (POSITIVE) ! OR DROUGHT (NEGATIVE) ! PALM - THE SELECTED PDSI (EITHER PRELIMINARY OR FINAL) ! ! This subroutine written and provided by CPC (Tom Heddinghaus ! & Paul Sabol). ! only the finest spaghetti code ! IF (X3.NE.0.) GO TO 10 ! IF X3=0 THE INDEX IS NEAR NORMAL AND EITHER A DRY OR WET SPELL ! EXISTS. CHOOSE THE LARGEST ABSOLUTE VALUE OF X1 OR X2. PALM = X1 IF ( ABS(X2) .GT. ABS(X1) ) PALM=X2 GO TO 100 ! 10 CONTINUE IF ( PROB.GT.0. .AND. PROB.LT.100. ) GO TO 20 ! A WEATHER SPELL IS ESTABLISHED AND PALM=X3 AND IS FINAL PALM = X3 GO TO 100 ! 20 CONTINUE PRO = PROB/100. IF (X3 .GT. 0.) GO TO 30 ! TAKE THE WEIGHTED SUM OF X3 AND X1 PALM = (1.-PRO)*X3 + PRO*X1 GO TO 100 ! 30 CONTINUE ! TAKE THE WEIGHTED SUM OF X3 AND X2 PALM = (1.-PRO)*X3 + PRO*X2 100 RETURN END subroutine casepdi !--------------------------------------------------------------- subroutine setmis(j,nendmo,d,amsng) ! fill in the remainder of the current year with missing values integer, intent(in) :: j, nendmo real, intent(in) :: amsng real, intent(inout) :: d(200,12) integer :: m if (nendmo.ge.12) return do m=nendmo+1,12 d(j,m)=amsng enddo return end subroutine setmis !--------------------------------------------------------------- subroutine close_iouts (inunit) integer, intent(in):: inunit integer :: k logical :: itsopen do k = 10, 37 inquire(unit=inunit+k, opened=itsopen) if ( itsopen ) then close(inunit+k) ! default status for non-scratch files is "keep" end if enddo ! Closing these scratch units will make them disappear do k= 60, 70 close (k) enddo do k= 75, 83 close (k) enddo return end subroutine close_iouts !--------------------------------------------------------------- ! It is assumed that a standard set of "intermediate" files has been ! created in the metadir, and now we just need to do things like add ! "state" abbreviations, element number, and an annual or "missing" ! value in the 13th data column. subroutine final_output(pdireg) character, intent(in) :: pdireg*5 character :: metadir*100, datadir*100 character :: inputstub(20)*10, outputstub(20)*6 character :: extin*6, extout*3, cfmt*16, dfmt*11, creg*5, stid*2, efmt*29 integer :: iyr, ist, idv, k, ielem, iok real :: mval(12), rsum data inputstub / 'pmdi.index','pdsi.index','z.index','phdi.index', & 'cp.index','spdat','pedat','pldat', & 'prdat','rdat','tldat','etdat','rodat','sssdat', & 'ssudat','effdat','probdat','x1dat','x2dat','x3dat'/ ! the element number will be extracted from the outputstub data outputstub /'pmdi08','pdsi05','zndx07','phdi06', & 'pmcp70','pmsw69','pmpe60','pmpl61', & 'pmpr62','pmre65','pmtl64','pmet63','pmro66','pmsl67', & 'pmul68','pmef55','pmpb56','pmx157','pmx258','pmx359'/ call getenv("metadir",metadir) call getenv("datadir",datadir) metadir = trim(metadir)//'/' datadir = trim(datadir)//'/' ! establish file names and formats depending on pdireg select case (pdireg) case ("div") extin = ".div" extout = "dv." dfmt = "(i3.3,i2.2)" efmt = "(a2,i3.3,i2.2,i2.2,i4,13f7.2)" case ("state") extin = ".state" extout = "st." dfmt = "(i3.3,i2.2)" efmt = "(a2,i3.3,i2.2,i2.2,i4,13f7.2)" case ("cty") extin = ".cty" extout = "cy." dfmt = "(i2.2,i3.3)" efmt = "(a2,i2.2,i3.3,i2.2,i4,13f7.2)" end select cfmt = "(a5,x,i4,12f7.2)" do k = 1, 20 open (14,file=trim(metadir)//trim(inputstub(k))//trim(extin), & status='old',form='formatted') open (24,file=trim(datadir)//trim(outputstub(k))//trim(extout), & status='replace',form='formatted') if (k==1) write (6,*) ' creating files like '//trim(outputstub(k))//trim(extout) read (outputstub(k),'(4x,i2)') ielem do read (14,cfmt,end=99) creg, iyr, mval read (creg,dfmt) ist, idv if ( trim(pdireg) == "cty" ) then stid = "CC" else call regid (ist, stid, ielem, iok) if (iok /= 1) call exit(1) endif ! Most variables get annual value = -99.99, but a few get the sum of the ! montly values rsum = -99.99 select case (ielem) case (60:61,63:64,66) ! If the December value is NOT missing, that indicates no missing values ! for the year, so make the output the sum. Otherwise -99.99. if ( abs(mval(12)+99.99) > 0.01 ) rsum = sum(mval) end select write (24,efmt) stid, ist, idv, ielem, iyr, mval, rsum enddo 99 close (14) close (24) enddo ! next file (k) end subroutine final_output !--------------------------------------------------------------- ! Return a 2 character abbreviation for given numeric state/region code ! ! Author: Roger Winchell ! Date: July 16, 1999 ! ! This routine replaces stateid with simpler code. ! stateid is complicated because the codes are not sequential (there are ! gaps) and because there are two sets of abbreviations for codes ! 101-110, depending on element: HDDU(3), CDDU(4), HDDC(25), CDDC(26) ! use the codes 151-160 below, so the subroutine adds 50 to the ! input number for these elements. ! ! If the input number is invalid, a subroutine returns a 1 error code. ! ! Main reason for this subroutine is to get data code out of main program. ! ! Have allowed 500 abbreviations for future growth. ! subroutine regid(nreg,abb,nel,iok) integer, intent(in) :: nreg, nel character, intent(out) :: abb*2 integer, intent(out) :: iok integer :: n character(len=2) :: regab1(100),regab2(100), & regab3(100),regab4(100),regab5(100) data regab1 / & 'AL','AZ','AR','CA','CO','CT','DE','FL','GA','ID', & 'IL','IN','IA','KS','KY','LA','ME','MD','MA','MI', & 'MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY', & 'NC','ND','OH','OK','OR','PA','RI','SC','SD','TN', & 'TX','UT','VT','VA','WA','WV','WI','WY','**','AK', & 'HI','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','PR','VI','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & 'PI','**','**','**','**','**','**','**','**','**' / data regab2 / & 'A1','A2','A3','A4','A5','A6','A7','A8','A9','AS', & 'AG','**','AN','AM','AU','**','**','**','**','AW', & 'ER','SR','CR','WR','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & 'P1','P2','P3','P4','P5','P6','P7','P8','P9','PS', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**' / data regab3 / & 'Q1','Q2','Q3','Q4','Q5','Q6','Q7','Q8','Q9','QA', & 'QB','QC','QD','QE','QF','QG','QH','QI','**','QJ', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','B1', & '**','**','**','**','B2','B4','**','**','**','B3', & 'B5','B6','**','**','B7','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**' / data regab4 / & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','E1', & '**','**','**','**','**','E4','**','**','**','**', & 'E5','E6','**','**','E7','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**' / data regab5 / & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','H1', & '**','**','**','**','**','H4','**','**','**','**', & 'H5','H6','**','**','H7','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**', & '**','**','**','**','**','**','**','**','**','**' / iok = 1 if ( nreg < 1 .OR. nreg > 500 ) then write (6,*)'regid says nreg is out of range (1-500): ', nreg iok = 0 return endif if ( nreg >= 151 .and. nreg <= 160 ) then write (6,*)'regid says nreg is in invalid range (151-160): ', nreg iok=0 return endif ! Fix the index for HDD/CDD elements n = nreg if ( ( nel == 3 .or. nel == 4 .or. nel == 25 .OR. nel == 26 ) & .and. ( nreg >= 101 .and. nreg <= 110 ) ) n = nreg + 50 select case (n) case (1:100) abb = regab1(n) case (101:200) abb = regab2(n-100) case (201:300) abb = regab3(n-200) case (301:400) abb = regab4(n-300) case (401:500) abb = regab5(n-400) end select if (abb == '**') then write (6,*) 'regid says abb = ** and that is a problem :(' iok=0 return endif return end subroutine regid !-------------------------------------- subroutine elem_abv(nel,aname,iok) ! ,--------------------------------------------------------------------, ! | Return an alpha name abbreviation for a given element code | ! `--------------------------------------------------------------------' integer, intent(in) :: nel character, intent(out) :: aname*4 integer, intent(out) :: iok character(len=1) :: cir(9) character(len=2) :: spi(7) character(len=4) :: elem(31),palelm(16) integer :: i, nel_code(31) data cir / 'a','b','c','d','e','f','h','l','r' / data elem / & 'pcpn','tmpc','hddu','cddu','pdsi','phdi','zndx','pmdi', & 'tmpu','tzcc','pzcc','torn','hddc','cddc','patw','patc', & 'papw','papd','paed','pasd','pamd','pald','paid','pann', & 'paiw','pasw','pamw','pavw','paew','tmax','tmin'/ data nel_code / & 1, 2, 3, 4, 5, 6, 7, 8, & 9, 11, 12, 24, 25, 26, 31, 32, & 33, 34, 41, 42, 43, 44, 45, 46, & 47, 48, 49, 50, 51, 27, 28 / data palelm /'pmef','pmpb','pmx1','pmx2','pmx3', & 'pmpe','pmpl','pmpr','pmet','pmtl','pmre','pmro','pmsl', & 'pmul','pmsw','pmcp'/ data spi / '01','02','03','06','09','12','24'/ iok = 1 aname = 'xxxx' select case (nel) case (55:70) aname = palelm(nel-54) case (13:21) aname = 'cir'//cir(nel-12) case (71:77) aname = 'sp'//spi(nel-70) case (81:87) aname = 'pr'//spi(nel-80) case (91:97) aname = 'ac'//spi(nel-90) end select if ( aname /= 'xxxx' ) return do i = 1, 31 if ( nel == nel_code(i) ) then aname = elem(i) return endif end do ! if this point is reached, an nel value was passed in that did not ! have a matching rule to return an aname write (6,*) ' no match in elem_abv for nel value: ', nel iok = 0 return end subroutine elem_abv !-------------------------