Page History

program compute
   

! input for normals (1981-2010 daily)
open(unit = 3, file = ST // CD // 'w_' // ayear, &
        status = 'unknown', form = 'formatted')

! ananlyis file
open(unit = 4, file = ST // CD // '_' // ayear // '_mir', &
        status = 'unknown', form = 'formatted')
write (4, *) '  i  i+1  CD iyr yr  week  MIR avg_mir mir_cur act_&
        aSPc_t aWIc_t aFAp_t aSUp_t aSPp_t'

! MIR file
open(unit = 8, file = ST // CD // '_' // ayear // '_plot', &
        status = 'unknown', form = 'formatted')
write (8, *) ' CD yr week month day 2012_mir avg_mir curr_mir&
        ACT_mir ACTm_mir ACTp_mir'

! NORMALS computations
! put data in daily,yearly array to compute normalso
! same normals for each analysis year
! week 1, 8 days, dec 29-Jan 5;
! all other weeks=7 days (skip day 365 on new year)

! read into array
print *, '1 CD = ', CD
print *, ST // CD // 'n'
open(unit = 10, file = ST // CD // 'n', status = 'old', form = 'formatted')
! for normals
bmon = 12
bday = 29
emon = 12
eday = 28
nday = 365
! becasue skip 2/29

do j = 1, 30
    bflg = 0

    !        print *, 'byear ',byear,'bmon ',bmon,'bday ',bday,' nday ',nday

    do i = 1, nday
        1        read (10, 5, end = 999) cstate, icd, year, mon, day, tmax, tmin, prcp
        5        format (a2, 1x, i2, 1x, i4, 1x, i2, 1x, i2, 1x, f5.1, 1x, f5.1, 1x, f8.2)
        if (mon .eq. 2 .and. day .eq. 29) go to 1

        if (bflg .eq. 1) go to 3
        if (bflg .eq. 0 .and.&
                year.lt.byear .or. mon.lt.bmon .or. day.lt.bday) then
            !           print *,'passed test ',year,byear,mon,bmon,day,bday
            go to 1
        endif

        3         bflg = 1
        dd(i) = day
        mm(i) = mon

        if (year.gt.2010) then
            go to 999
        endif
        !          write (*,*) i,year,mon,day,tmax,tmin,prcp,byear,bmon,bday

        ! compute daily mean
        ! convert to celcius from farenheit
        ! compute degree days, base 22 C
        ! convert to cm from inches
        t(i, j) = (5.0 / 9.0) * (((tmax + tmin) / 2.0) - 32.0)
        p(i, j) = prcp * 2.54

        if (day.eq.eday .and. mon.eq.emon) then
            print *, i, ' ', j, ' ', year, ' ', mon, ' ', day
            go to 998
        endif
        ! go to next year at nday=365
    998     enddo

    !        close(10)
enddo

999    close(10)

! compute weekly daily averages for d and t (and weekly sum for precip) for each 30 years of data (1981-2010)

do j = 1, 30

    ! compute average T and sum P for 7 day weeks
    do i = 1, 8
        k = 1
        wt(1, j) = (t(i, j) + t((i + 1), j) + t((i + 2), j) + t((i + 3), j) + t((i + 4), j)&
                + t((i + 5), j) + t((i + 6), j) + t((i + 7), j)) / 8.0
        wp(1, j) = (p(i, j) + p((i + 1), j) + p((i + 2), j) + p((i + 3), j) + p((i + 4), j)&
                + p((i + 5), j) + p((i + 6), j) + p((i + 7), j))

        print *, j, k, t(i, j), t((i + 1), j), t((i + 2), j), t((i + 3), j), t((i + 4), j), &
                t((i + 5), j), t((i + 6), j)
    enddo

    do i = 9, 365, 7
        k = k + 1
        ! compute average T and sum P for 7 day weeks

        wt(k, j) = (t(i, j) + t((i + 1), j) + t((i + 2), j) + t((i + 3), j) + t((i + 4), j)&
                + t((i + 5), j) + t((i + 6), j)) / 7.0
        wp(k, j) = (p(i, j) + p((i + 1), j) + p((i + 2), j) + p((i + 3), j) + p((i + 4), j)&
                + p((i + 5), j) + p((i + 6), j))

        !        print *, i,j,k,t(i,j),t((i+1),j),t((i+2),j),t((i+3),j),t((i+4),j),
        !     !  t((i+5),j),t((i+6),j)
        !       print *,'weekly sum and avg ppt_cm',i,j,k,wp(k,j),(wp(k,j)/7.0)
        !       print *,i,j,k,p(i,j),p((i+1),j),p((i+2),j),p((i+3),j),p((i+4),j),
        !     !  p((i+5),j),p((i+6),j),wp(k,j),(wp(k,j)/7.0)

    enddo
enddo

! compute 30-year weekly normals
do k = 1, 53
    wknorm_t(k) = 0.0
    wknorm_p(k) = 0.0
    wknorm_d(k) = 0.0
    wksum_t(k) = 0.0
    wksum_p(k) = 0.0
enddo

! compute weekly 30-year normals
do k = 1, 52
    do j = 1, 30
        wksum_t(k) = wksum_t(k) + wt(k, j)
        wksum_p(k) = wksum_p(k) + wp(k, j)
        !            wksum_d(k)=wksum_d(k)+wd(k,j)
    enddo
    wknorm_t(k) = wksum_t(k) / 30.
    wknorm_p(k) = wksum_p(k) / 30.
    !            wknorm_d(k)=wksum_d(k)/30.
    !        print *,'wk,week_norms_p,weeksum_p',k,wknorm_p(k),wksum_p(k)
enddo
wknorm_t(53) = wknorm_t(1)
wknorm_p(53) = wknorm_p(1)
!           wknorm_d(53)=wknorm_d(1)

! compute weekly dw 30-yr normals
do k = 1, 53
    wknorm_d(k) = 0.0
enddo

do k = 1, 53
    ! if weekly average temperature is greater than 22:
    if(wknorm_t(k) .gt. 22.0) then
        wksum_d = wksum_d + (wknorm_t(k) - 22.0)
    elseif (wknorm_t(k) .le. 22.0) then
        wksum_d = wksum_d + 0.0
    endif
    wknorm_d(k) = wksum_d
enddo

! compute seasonal weekly 30-year normals
do k = 1, 9
    WIsum_t = WIsum_t + wknorm_t(k)
    WIsum_p = WIsum_p + wknorm_p(k)
    !            WIsum_np=WIsum_np+wknorm_np(k)
    WIsum_d = WIsum_d + wknorm_d(k)
enddo
do k = 49, 52
    WIsum_t = WIsum_t + wknorm_t(k)
    WIsum_p = WIsum_p + wknorm_p(k)
    !            WIsum_np=WIsum_np+wknorm_np(k)
    WIsum_d = WIsum_d + wknorm_d(k)
enddo

do k = 10, 22
    SPsum_t = SPsum_t + wknorm_t(k)
    SPsum_p = SPsum_p + wknorm_p(k)
    !            SPsum_np=SPsum_np+wknorm_np(k)
    SPsum_d = SPsum_d + wknorm_d(k)
enddo
do k = 23, 35
    SUsum_t = SUsum_t + wknorm_t(k)
    SUsum_p = SUsum_p + wknorm_p(k)
    !            SUsum_np=SUsum_np+wknorm_np(k)
    SUsum_d = SUsum_d + wknorm_d(k)
enddo
do k = 36, 48
    FAsum_t = FAsum_t + wknorm_t(k)
    FAsum_p = FAsum_p + wknorm_p(k)
    !            FAsum_np=FAsum_np+wknorm_np(k)
    FAsum_d = FAsum_d + wknorm_d(k)
enddo

! seasonal weekly normals
WInorm_t = WIsum_t / 13.
SPnorm_t = SPsum_t / 13.
SUnorm_t = SUsum_t / 13.
FAnorm_t = FAsum_t / 13.

WInorm_d = WIsum_d / 13.
SPnorm_d = SPsum_d / 13.
SUnorm_d = SUsum_d / 13.
FAnorm_d = FAsum_d / 13.

WInorm_p = WIsum_p / 13.
SPnorm_p = SPsum_p / 13.
SUnorm_p = SUsum_p / 13.
FAnorm_p = FAsum_p / 13.

! input - get previous year quarterly temperature and precipitation diffeences

print *, '2 CD = ', CD
open(unit = 1, file = ST // CD, status = 'old', form = 'formatted')
do i = 1, 366
    tmn(i) = 999.9
    ppt(i) = 999.9
enddo

! normals are done

! previouw year = iyear-1
! good for 1981 to 2020 (i=1 to 40)
do i = 1, 40
    if (begyr(i) .eq. iyear - 1) then
        byear = nbegyr(i)
        nday = ndays(i)
        bday = nbegday(i)
        bmon = nbegmon(i)
    endif
enddo

n = 0
bflg = 0
dsum_d = 0
do i = 1, nday
    9       read (1, 5, end = 19) cstate, icd, year, mon, day, tmax, tmin, prcp
    if (bflg.eq. 1) go to 18
    if (bflg.eq. 0 .and.&
            year.lt.byear .or. mon.lt. bmon .or. day.lt.bday) then
        dsum_d = 0
        go to 9
    endif
    bflg = 1
    !18         write (*,*) i,iyear,year,mon,day,tmax,tmin,prcp,bflg
    18          n = n + 1

    if (tmax .eq. 999.9) go to 19
    ! mean temperature
    ! convert to centigrade
    ! convert to cm
    write (*, 5) cstate, icd, year, mon, day, tmax, tmin, prcp
    tmn(i) = (5.0 / 9.0) * (((tmax + tmin) / 2.0) - 32.0)
    ppt(i) = prcp * 2.54
    if (ppt(i) .lt. 0) ppt(i) = 999.9
    !        write (*,*) i,iyear,year,mon,day,tmn(i),ppt(i)
enddo

19      close(1)
print *, 'iyear,n,nday = ', iyear, n, nday

! compute weekly averages for previous year
do k = 1, 53
    wctmn(k) = 999.9
    wcppt(k) = 999.9
enddo

k = 0
do i = 1, 366, 7
    k = k + 1

    !         write (*,*) k,i,i+1,i+2,i+3,i+4,i+5,i+6
    !         print *, k, tmn(i),tmn(i+1),tmn(i+2),tmn(i+3),tmn(i+4),
    !     ! tmn(i+5),tmn(i+6)
    !         print *, k, ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
    !     ! ppt(i+5),ppt(i+6)

    wctmn(k) = (tmn(i) + tmn(i + 1) + tmn(i + 2) + tmn(i + 3) + tmn(i + 4)&
            + tmn(i + 5) + tmn(i + 6)) / 7.0
    wcppt(k) = (ppt(i) + ppt(i + 1) + ppt(i + 2) + ppt(i + 3) + ppt(i + 4)&
            + ppt(i + 5) + ppt(i + 6))

    !        print *, 'prev_ yr weekly sum avg ppt ',wcppt(k),(wcppt(k)/7.0)
    !        write (*,*) ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
    !     ! ppt(i+5),ppt(i+6),wcppt(k),(wcppt(k)/7.0)
    !         print *, ' k,wcppt(k) = ',k,wcppt(k)
enddo

SP_prev_sum_t = 0.0
SP_prev_sum_p = 0.0
SP_prev_yr_t = 0.0
SP_prev_yr_p = 0.0

SU_prev_sum_t = 0.0
SU_prev_sum_p = 0.0
SU_prev_yr_t = 0.0
SU_prev_yr_p = 0.0

FA_prev_sum_t = 0.0
FA_prev_sum_p = 0.0
FA_prev_yr_t = 0.0
FA_prev_yr_p = 0.0

! get december from previous year
WI_curr_sum_t = 0.0
WI_curr_sum_p = 0.0
WI_curr_yr_t = 0.0
WI_curr_yr_p = 0.0

do k = 10, 22
    SP_prev_sum_t = SP_prev_sum_t + wctmn(k)
    SP_prev_sum_p = SP_prev_sum_p + wcppt(k)
enddo
SP_prev_yr_t = SP_prev_sum_t / 13.0
SP_prev_yr_p = SP_prev_sum_p / 13.0
!        SP_prev_yr_p=SP_prev_sum_p

do k = 23, 35
    SU_prev_sum_t = SU_prev_sum_t + wctmn(k)
    SU_prev_sum_p = SU_prev_sum_p + wcppt(k)
enddo
SU_prev_yr_t = SU_prev_sum_t / 13.0
SU_prev_yr_p = SU_prev_sum_p / 13.0
!        SU_prev_yr_p=SU_prev_sum_p

do k = 36, 48
    FA_prev_sum_t = FA_prev_sum_t + wctmn(k)
    FA_prev_sum_p = FA_prev_sum_p + wcppt(k)
enddo
FA_prev_yr_t = FA_prev_sum_t / 13.0
FA_prev_yr_p = FA_prev_sum_p / 13.0
!        FA_prev_yr_p=FA_prev_sum_p

do k = 49, 52
    WI_curr_sum_t = WI_curr_sum_t + wctmn(k)
    WI_curr_sum_p = WI_curr_sum_p + wcppt(k)
enddo

! compute quarterly previous year differences

df_SPp_t = SP_prev_yr_t - SPnorm_t
df_SPp_p = (SP_prev_yr_p - SPnorm_p)

df_SUp_t = SU_prev_yr_t - SUnorm_t
df_SUp_p = (SU_prev_yr_p - SUnorm_p)

df_FAp_t = FA_prev_yr_t - FAnorm_t
df_FAp_p = (FA_prev_yr_p - FAnorm_p)

write (3, *)  'SPRING Prev_YR NORM DIF for T D and P'
write (3, 51) SP_prev_yr_t, c, SPnorm_t, c, df_SPp_t, c, SP_prev_yr_d, &
        c, SPnorm_d, c, df_SPp_d, c, SP_prev_yr_p, c, SPnorm_p, c, df_SPp_p
write (3, *)  'SUMMER Prev_YR NORM DIF for T D and P'
write (3, 51) SU_prev_yr_t, c, SUnorm_t, c, df_SUp_t, c, SU_prev_yr_d, &
        c, SUnorm_d, c, df_SUp_d, c, SU_prev_yr_p, c, SUnorm_p, c, df_SUp_p
write (3, *)  'FALL Prev_YR NORM DIF for T D and P'
write (3, 51) FA_prev_yr_t, c, FAnorm_t, c, df_FAp_t, c, FA_prev_yr_d, &
        c, FAnorm_d, c, df_FAp_d, c, FA_prev_yr_p, c, FAnorm_p, c, df_FAp_p

write (*, *)  'SPRING Prev_YR NORM DIF for T D and P'
write (*, 51) SP_prev_yr_t, c, SPnorm_t, c, df_SPp_t, c, SP_prev_yr_d, &
        c, SPnorm_d, c, df_SPp_d, c, SP_prev_yr_p, c, SPnorm_p, c, df_SPp_p
write (*, *)  'SUMMER Prev_YR NORM DIF for T D and P'
write (*, 51) SU_prev_yr_t, c, SUnorm_t, c, df_SUp_t, c, SU_prev_yr_d, &
        c, SUnorm_d, c, df_SUp_d, c, SU_prev_yr_p, c, SUnorm_p, c, df_SUp_p
write (*, *)  'FALL Prev_YR NORM DIF for T D and P'
write (*, 51) FA_prev_yr_t, c, FAnorm_t, c, df_FAp_t, c, FA_prev_yr_d, &
        c, FAnorm_d, c, df_FAp_d, c, FA_prev_yr_p, c, FAnorm_p, c, df_FAp_p

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! input - get model analysis year

print *, '3 CD = ', CD
open(unit = 1, file = ST // CD, status = 'old', form = 'formatted')
do i = 1, 366
    tmn(i) = 999.9
    ppt(i) = 999.9
    ddd(i) = 0.0
enddo

do i = 1, 53
    df_wkd(k) = 0.0
enddo

! set current analysis year = iyear
do i = 1, 40
    if (begyr(i) .eq. iyear) then
        byear = nbegyr(i)
        nday = ndays(i)
        bday = nbegday(i)
        bmon = nbegmon(i)
    endif
enddo
print *, 'byear nday bday bmon ', byear, ' ', nday, ' ', bday, ' ', bmon

n = 0
bflg = 0
dsum_d = 0
do i = 1, nday
    cfirstmon = '01'
    cfirstday = '01'
    10       read (1, 5, end = 88) cstate, icd, year, mon, day, tmax, tmin, prcp
    if (year.lt.byear .and. day.lt.bday .and. mon.lt.bmon) go to 10
    !         write (*,*) cstate,icd,year,mon,day,tmax,tmin,prcp,'begin',
    !     !  byear,bday,bmon

    if ((tmax .eq. -99 .or. tmin .eq. -99 .or. prcp .lt. 0.0)&
            .and. n.gt.140) then
        ! for forecast
        firstday = day
        firstmon = mon
        write(cfirstmon, 11) firstmon
        write(cfirstday, 11) firstday
        11       format(i2)
        ! to compensate for blanks ...
        if (cfirstmon .eq. ' 1') cfirstmon = '01'
        if (cfirstmon .eq. ' 2') cfirstmon = '02'
        if (cfirstmon .eq. ' 3') cfirstmon = '03'
        if (cfirstmon .eq. ' 4') cfirstmon = '04'
        if (cfirstmon .eq. ' 5') cfirstmon = '05'
        if (cfirstmon .eq. ' 6') cfirstmon = '06'
        if (cfirstmon .eq. ' 7') cfirstmon = '07'
        if (cfirstmon .eq. ' 8') cfirstmon = '08'
        if (cfirstmon .eq. ' 9') cfirstmon = '09'
        if (cfirstday .eq. ' 1') cfirstday = '01'
        if (cfirstday .eq. ' 2') cfirstday = '02'
        if (cfirstday .eq. ' 3') cfirstday = '03'
        if (cfirstday .eq. ' 4') cfirstday = '04'
        if (cfirstday .eq. ' 5') cfirstday = '05'
        if (cfirstday .eq. ' 6') cfirstday = '06'
        if (cfirstday .eq. ' 7') cfirstday = '07'
        if (cfirstday .eq. ' 8') cfirstday = '08'
        if (cfirstday .eq. ' 9') cfirstday = '09'
        go to 88
    endif

    if (bflg.eq. 1) go to 8
    if (bflg.eq. 0 .and.&
            year.lt.byear .or. mon.lt. bmon .or. day.lt.bday) then
        dsum_d = 0
        !          print *, 'got iyear',iyear,'mon ',mon,'day ',day,'eday ',eday
        go to 10
    endif

    bflg = 1
    write (*, *)'first bflg=1,  i year mon day  tmax tmin prcp'
    write (*, *) i, year, mon, day, ' ', tmax, tmin, prcp
    8         write (*, *) i, iyear, year, mon, day, tmax, tmin, prcp, bflg
    n = n + 1

    ! mean temperature
    ! convert to centigrade
    ! compute degree days
    ! convert to cm
    pdd(n) = day
    pmm(n) = mon
    tmn(n) = (5.0 / 9.0) * (((tmax + tmin) / 2.0) - 32.0)
    dum_d = 0
    if (tmn(n) .ge. 22.0) dum_d = (tmn(n) - 22.0)

    ppt(n) = prcp * 2.54

    ! changed to accomodate temperature good, but precipitation missing!
    !5-19-2017

    if (ppt(n) .lt. 0.0) ppt(n) = 0.0
    !         if (ppt(n) .gt. 1000) ppt(n)=0.0
    !         if (ppt(n) .lt. 0) ppt(n)=999.9

    !        write (*,*) i,iyear,year,mon,day,tmn(i),ppt(i)
enddo

! missing data or end of year, so close
88      close(1)

ncurr = n
print *, 'iyear,n,ncurr,nday = ', iyear, n, ncurr, nday

! ADD FORECAST DATA

! input - get 10-day max and min temp forecasts, and 3 day precip forecasts

! for this year 2017:

if (iyear .eq. 2017) then
    ! inintalize forecast data
    do j = 1, 9
        do i = 1, 11
            mntemp(i, j) = 999.9
        enddo
        do i = 1, 10
            mxtemp(i, j) = 999.9
        enddo
        do i = 1, 3
            prcpd(i, j) = 999.9
        enddo
        ! set later days of precip to 0.0
        do i = 4, 10
            prcpd(i, j) = 0.0
        enddo
    enddo

    !  see /home/nan/WNV_WX/Not_Smooth/run.get.forecast
    !       open(unit=7,file='ILCD_divisions_ALL',
    !    !  status='old',form='formatted')

    print *, 'ILCD_' // ayear // cfirstmon // cfirstday

    open(unit = 7, file = 'ILCD_' // ayear // cfirstmon // cfirstday, &
            status = 'old', form = 'formatted')
    print *, '-----------------------------'
    print *, 'first day with missing data = '
    print *, 'ILCD_' // ayear // cfirstmon // cfirstday
    read (7, *) line1
    !DIV,CYCLE,ELEMENT,VALIDTIME,MEAN-VALUE
    !1,2014041412,MINT,2014041512,23.115

    ! should have 11 min temps, 10 max temps, 13 precips, otherwise bombs
    n_mnT = 0
    n_mxT = 0
    n_ppt = 0
    do j = 1, 9
        do i = 1, 11
            read (7, *, end = 6) ilcd, rundate, lab1, valdate, mntemp(i, j)
            cpdd = valdate(7 : 8)
            cpmm = valdate(5 : 6)
            read(cpdd, 27) pdd(n + i)
            read(cpmm, 27) pmm(n + i)
            27      format(i2)

            print *, '----------------------'
            print *, ' n= ', n, 'ncurr = ', ncurr
            write (*, *) 'fcst date', ilcd, rundate, lab1, valdate, mntemp(i, j), &
                    pmm(n + i), pdd(n + i), i, (n + i), n_mnT
            n_mnT = n_mnT + 1
        enddo

        6        do i = 1, 10
            read (7, *, end = 66) ilcd, rundate, lab1, valdate, mxtemp(i, j)
            write (*, *) ilcd, rundate, lab1, valdate, mxtemp(i, j), &
                    pmm(n + i), pdd(n + i)
            n_mxT = n_mxT + 1
        enddo

        66       do i = 1, 13
            read (7, *, end = 666) ilcd, rundate, lab1, valdate, cprcp6
            write (*, *) ilcd, rundate, lab1, valdate, cprcp6
            n_ppt = n_ppt + 1
            ! read from character to integer
            if(cprcp6(1 : 3).eq.'NaN') cprcp6 = '999.90'
            read(cprcp6, 28) prcp6(i, j)
            28      format(f6.2)
            if (prcp6(i, j) .ge. 999.0) prcp6(i, j) = 0.0
        enddo

        666        if (n_mnT.lt.11 .or. n_mxT.lt.10 .or. n_ppt.lt.13)&
                print *, 'too little data: n_mnT,n_mxT,n_ppt ', n_mnT, n_mxT, n_ppt
        print*, 'should have 11 n_mnT, 10 n_mxT, 13 n_ppt'
        prcpd(1, j) = prcp6(2, j) + prcp6(3, j) + prcp6(4, j) + prcp6(5, j)
        prcpd(2, j) = prcp6(6, j) + prcp6(7, j) + prcp6(8, j) + prcp6(9, j)
        prcpd(3, j) = prcp6(10, j) + prcp6(11, j) + prcp6(12, j) + prcp6(13, j)
        prcpd(4, j) = 0.0
        prcpd(5, j) = 0.0
        prcpd(6, j) = 0.0
        prcpd(7, j) = 0.0
        prcpd(8, j) = 0.0
        prcpd(9, j) = 0.0
        prcpd(10, j) = 0.0

        print *, 'read temp and precp forecasts'

        ! tack on to endof ytd data
        6666     print *, 'too little data? n_mnT,n_mxT,n_ppt ', n_mnT, n_mxT, n_ppt
        if (ilcd .eq. icd) then
            print *, 'ilcd = ', ilcd, icd
            dum_d = 0
            k = 1
            do i = ncurr + 1, ncurr + 8
                print *, 'ncurr+1 to ncurr+10', i, k, j, tmn(i), mntemp(k, j)
                tmn(i) = (5.0 / 9.0) * (((mxtemp(k, j) + mntemp(k, j)) / 2.0) - 32.0)
                dum_d = 0
                if (tmn(i) .ge. 22.0) dum_d = (tmn(i) - 22.0)
                print *, 'dum_d', dum_d
                print *, 'ncurr+1 to ncurr+10', i, k, j, tmn(i), mntemp(k, j)
                k = k + 1
            enddo

            k = 1
            do i = ncurr + 1, ncurr + 8
                ppt(i) = prcpd(k, j) * 2.54
                print *, 'ncurr+1 to ncurr+10', i, k, j, ppt(i), prcpd(k, j)
                k = k + 1
            enddo

        endif
        ! go to next cliamte division
    enddo

    ! write out after forecast
    do i = ncurr + 1, ncurr + 8
        write (*, *) 'FORECAST DAYS ', pmm(i), pdd(i)
        write (*, 15) year, i, ncurr, pday(i), pmm(i), tmn(i), ppt(i)
        15       format (i4, 1x, i4, 1x, i4, 1x, i2, 1x, i2, 1x, f7.3, 1x, f7.4)
    enddo
endif
! end of forecast data for 2017

! compute weekly averages for model analysis year
do k = 1, 53
    wctmn(k) = 999.9
    wcppt(k) = 999.9
enddo

k = 1
!        do i=1,n-1,7
do i = 1, ncurr + 8, 7
    ! only go out 8 days (should match dupage model, will have 3 days with precip forecast

    pday(k) = pdd(i)
    pmon(k) = pmm(i)
    !         pday(k)=pdd(i+4)
    !         pmon(k)=pmm(i+4)
    print *, 'meanT @7day ', i, pmm(i), pdd(i), &
            tmn(i), tmn(i + 1), tmn(i + 2), tmn(i + 3), tmn(i + 4), tmn(i + 5), tmn(i + 6)
    if (tmn(i) .ge. 999.9 .or. tmn(i + 1) .ge. 999.9 .or.&
            tmn(i + 6) .ge. 999.9) then
        print *, 'meanT @7day ', i, pmm(i), pdd(i), &
                tmn(i + 4), tmn(i + 5), tmn(i + 6)
        print *, '999', pmm(i + 4), pdd(i + 4)
        endk = k
        go to 79
    endif

    wctmn(k) = (tmn(i) + tmn(i + 1) + tmn(i + 2) + tmn(i + 3) + tmn(i + 4)&
            + tmn(i + 5) + tmn(i + 6)) / 7.0
    wcddd(k) = (ddd(i) + ddd(i + 1) + ddd(i + 2) + ddd(i + 3) + ddd(i + 4)&
            + ddd(i + 5) + ddd(i + 6)) / 7.0
    !         wcppt(k) = (ppt(i)+ppt(i+1)+ppt(i+2)+ppt(i+3)+ppt(i+4)
    !     ! +ppt(i+5)+ppt(i+6))/7.0
    wcppt(k) = (ppt(i) + ppt(i + 1) + ppt(i + 2) + ppt(i + 3) + ppt(i + 4)&
            + ppt(i + 5) + ppt(i + 6))

    print *, 'k,pmon(k),pday(k),i,pmm(i),pdd(i) = ', &
            k, pmon(k), pday(k), ' ', i, pmm(i), pdd(i)

    !         write (*,*) k,i,i+1,i+2,i+3,i+4,i+5,i+6
    !         print *,'cur_yr weekly avg sum ppt ',k,wcppt(k),(wcppt(k)/7.0)

    !         print *, k, tmn(i),tmn(i+1),tmn(i+2),tmn(i+3),tmn(i+4),
    !     ! tmn(i+5),tmn(i+6),wctmn(k),pmon(k),pday(k)
    !         print *, k, ddd(i),ddd(i+1),ddd(i+2),ddd(i+3),ddd(i+4),
    !     ! ddd(i+5),ddd(i+6),wcddd(k)
    !         print *, k, ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
    !     ! ppt(i+5),ppt(i+6),wcppt(k)

    k = k + 1
    79       continue
    endk = k - 1
    print *, 'endk = ', endk
enddo

! compute current weekly dw
do k = 1, endk
    wcddd(k) = 0
enddo
wksum_d = 0

do k = 1, endk
    if(wctmn(k) .gt. 22.0) then
        wcsum_d = wcsum_d + (wctmn(k) - 22.0)
    elseif (wctmn(k) .le. 22.0) then
        wcsum_d = wcsum_d + 0.0
    endif
    wcddd(k) = wcsum_d
enddo

! summed december already from previous year
!        WI_curr_sum_t=0.0
!        WI_curr_sum_d=0.0
!        WI_curr_sum_p=0.0
SP_curr_sum_t = 0.0
SP_curr_sum_d = 0.0
SP_curr_sum_p = 0.0

WI_curr_yr_t = 0.0
WI_curr_yr_d = 0.0
WI_curr_yr_p = 0.0
SP_curr_yr_t = 0.0
SP_curr_yr_d = 0.0
SP_curr_yr_p = 0.0

if (endk .ge. 9) then
    ! k=week, weeks 1-9 to previously computed WI_curr
    print *, 'Dec WI_curr_sum_p = ', WI_curr_sum_p
    do k = 1, 9
        WI_curr_sum_t = WI_curr_sum_t + wctmn(k)
        WI_curr_sum_d = WI_curr_sum_d + wcddd(k)
        WI_curr_sum_p = WI_curr_sum_p + wcppt(k)
    enddo
    WI_curr_yr_t = WI_curr_sum_t / 13.0
    WI_curr_yr_d = WI_curr_sum_d / 13.0
    WI_curr_yr_p = WI_curr_sum_p / 13.0
    !        WI_curr_yr_p=WI_curr_sum_p
endif
!        print *, 'DecFeb WI_curr_sum_p = ',WI_curr_sum_p
!        print *, 'DecFeb WI_curr_yr_p =  ',WI_curr_yr_p

if (endk .ge. 22) then
    do k = 10, 22
        SP_curr_sum_t = SP_curr_sum_t + wctmn(k)
        SP_curr_sum_d = SP_curr_sum_d + wcddd(k)
        SP_curr_sum_p = SP_curr_sum_p + wcppt(k)
    enddo
    SP_curr_yr_t = SP_curr_sum_t / 13.0
    SP_curr_yr_d = SP_curr_sum_d / 13.0
    SP_curr_yr_p = SP_curr_sum_p / 13.0
    !        SP_curr_yr_p=SP_curr_sum_p
endif

do k = 1, endk
    df_wkt(k) = wctmn(k) - wknorm_t(k)
    df_wkd(k) = wcddd(k) - wknorm_d(k)
    df_wkp(k) = (wcppt(k) - wknorm_p(k))

    !            print *,'k,df_wkp ',k,df_wkp(k)
    !            print *,'year,endk,k,wcppt(k),wknorm_p(k),df_wkp(k) ',
    !     !               year,endk,k,wcppt(k),wknorm_p(k),df_wkp(k)
    !            print *,'year,endk,k,wcddd(k),wknorm_d(k) ',
    !     !               year,endk,k,wcddd(k),wknorm_d(k)
enddo
! what is this!!
!        if (iyear .eq. 2015)then
!            df_wkd(endk)=0.0
!            df_wkp(endk)=0.0
!c            print *,'endk,df_wkp ',k,df_wkp(endk)
!        endif

df_WIc_t = WI_curr_yr_t - WInorm_t
df_WIc_d = WI_curr_yr_d - WInorm_d
df_WIc_p = (WI_curr_yr_p - WInorm_p)
!        df_WIc_np=(WI_curr_yr_p-WInorm_np)

df_SPc_t = SP_curr_yr_t - SPnorm_t
df_SPc_d = SP_curr_yr_d - SPnorm_d
df_SPc_p = (SP_curr_yr_p - SPnorm_p)
!        df_SPc_np=(SP_curr_yr_p-SPnorm_np)

if (endk .lt. 13) then
    df_WIc_t = 0.0
    df_WIc_d = 0.0
    df_WIc_p = 0.0
endif
if (endk .lt. 26) then
    df_SPc_t = 0.0
    df_SPc_d = 0.0
    df_SPc_p = 0.0
endif

! write weekly normals prevous and current quarts data
write (3, *) 'year week wctmn wknrm_t df_wkt wcddd wknrm_d df_wkd&
        wcppt wknrm_p df_wkp'
write (*, *) 'year week wctmn wknrm_t df_wkt wcddd wknrm_d df_wkd&
        wcppt wknrm_p df_wkp'

do k = 1, endk
    write (3, 50) iyear, c, k, c, wctmn(k), c, wknorm_t(k), c, df_wkt(k), c, &
            wcppt(k), c, wknorm_p(k), c, df_wkp(k), c
    write (*, 50) iyear, c, k, c, wctmn(k), c, wknorm_t(k), c, df_wkt(k), c, &
            wcppt(k), c, wknorm_p(k), c, df_wkp(k), c
    50      format (i4, a1, i4, 9(a1, f6.2, a1, f6.2, a1, f6.2))
enddo

! write quartly normals prevous and current quarts data
write (3, *)  'WI_curr_t WInrm_t df_WIc_t  WI_curr_d WInrm_d df_q&
        1c_d WI_curr_p WInrm_p df_WIc_p'
write (*, *)  'WI_curr_t WInrm_t df_WIc_t  WI_curr_d WInrm_d df_q&
        1c_d WI_curr_p WInrm_p df_WIc_p'

write (3, 51) WI_curr_yr_t, c, WInorm_t, c, df_WIc_t, c, WI_curr_yr_d, &
        c, WInorm_d, c, df_WIc_d, c, WI_curr_yr_p, c, WInorm_p, c, df_WIc_p
write (*, 51) WI_curr_yr_t, c, WInorm_t, c, df_WIc_t, c, WI_curr_yr_d, &
        c, WInorm_d, c, df_WIc_d, c, WI_curr_yr_p, c, WInorm_p, c, df_WIc_p
51      format (1x, 3(f6.2, a1), 4x, 3(f6.2, a1), 4x, 3(f6.2, a1))
print *, 'end of quart info'

!51      format (1x,3(f6.2,a1),4x,3(f6.2,a1),4x,3(f6.2,a1),4x,3(f6.2,a1))

!        write (3,*)  'SP_curr_t SPnrm_t df_SPc_t  SP_curr_d SPnrm_d df_q
!     !2c_d SP_curr_p SPnrm_p df_SPc_p'
!        write (3,51) SP_curr_yr_t,c,SPnorm_t,c,df_SPc_t,c, SP_curr_yr_d,
!     !      c,SPnorm_d,c,df_SPc_d,c, SP_curr_yr_p,c,SPnorm_p,c,df_SPc_p
!     !      c,SP_curr_yr_np,c,SPnorm_np,c,df_SPc_np,c

! input - climate division daylight data
! for cd's 1-9, cd = icd

open(unit = 14, file = 'CD_Daylight.csv'&
        , status = 'old', form = 'formatted')
read (14, *) line1
do i = 1, 53
    read (14, *) wk(i), (cd_dylt(i, j), j = 1, 9)
    !         print *, wk(i),(cd_dylt(i,j),j=1,9)
enddo
close(14)

!WEEK,HRS_DYLT_DIV1,HRS_DYLT_DIV2,HRS_DYLT_DIV3,HRS_DYLT_DIV4,HRS_DYLT_DIV5,HRS_DYLT_DIV6,HRS_DYLT_DIV7,HRS_DYLT_DIV8,HRS_DYLT_DIV9
!1,9.228,9.228,9.367,9.367,9.350,9.467,9.489,9.600,9.600

! input - MIR 9-year average by 9 CDs and (4) combined areas
! for cd's 1-9, cd = icd

open(unit = 15, file = 'Weekly_MIRavg_05_12.csv', &
        status = 'old', form = 'formatted')
read (15, *) line1
do i = 1, 52
    read (15, *) wk(i), (cd_avg_mir(i, j), j = 1, 13)
    !         print *, 'each year',wk(i),(cd_avg_mir(i,j),j=1,13)
enddo
close(15)

!WeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08
!1,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00

! input - 2012 MIR by CD

open(unit = 16, file = 'Weekly_MIR_2012.csv', status = 'old', &
        form = 'formatted')
read (16, *) line1
do i = 15, 45
    read (16, *) wk(i), (cd_2012mir(i, j), j = 1, 13)
    !          print *, '2012 mir1',wk(i),(cd_2012mir(i,j),j=1,13)
enddo
print *, 'act_mir year and iyear = ', year, ' ', iyear

! input - 2005-2013 MIR by CD
!        do i=1,53
!        do i=1,13
!        cd_act_mir=0.0
!        enddo
!        enddo
!
!        if (iyear .lt. 2005) go to 228
!        if (iyear .gt. 2013) go to 228
!        open(unit=26,file='CD_ACT_mir_2005_13.csv',status='old',
!     ! form='formatted')
!        read (26,*) line1
!        i=1
!226      read (26,*,end=228) year,wk(i),cumwk(i),(cd_act_mir(i,j),j=1,9)
!         print *, i,year,wk(i),(cd_act_mir(i,j),j=1,9)
!         if (year .lt. iyear) go to 226
!         if (year .gt. iyear) go to 227
!        print *, year,' ',iyear,'cd =',j,'act_mir = ',cd_act_mir(i,j)
!         i=i+1
!         go to 226
!228     close(26)

!WeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08
!1,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00

! input - current year (2017) ACTUAL MIR by CD

!        open(unit=18,file='CD_ACT_MIR_'//ayear//'.prn',status='old',
open(unit = 18, file = 'CD_ACT_MIR_2017.prn', status = 'old', &
        form = 'formatted')
read (18, *) line1
do i = 1, 52
    read(18, *, end = 118)wk(i), &
            (ACT_mir(i, j), ACTm_mir(i, j), ACTp_mir(i, j), j = 1, 9)
enddo
118     continue

do i = 1, 52
    write(*, *)&
            wk(i), icd, ACT_mir(i, icd), ACTm_mir(i, icd), ACTp_mir(i, icd)
enddo

!WeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08

    ! read coeffientes for 9 CD and 4 combined areas
    open(unit = 17, file = 'CD_Coef_2017.csv', &
            status = 'old', form = 'formatted')
    read (17, *) line1

    read (17, *) lab(1), (intercept(j), j = 1, 9)
    print *, lab(1), (intercept(icd))
    read (17, *) lab(2), (cDW_Lg1(j), j = 1, 9)
    print *, lab(2), (cDW_Lg1(icd))
    read (17, *) lab(3), (cDW_Lg2(j), j = 1, 9)
    print *, lab(3), (cDW_Lg2(icd))
    read (17, *) lab(4), (cDW_Lg3(j), j = 1, 9)
    print *, lab(4), (cDW_Lg3(icd))
    read (17, *) lab(5), (cDW_Lg4(j), j = 1, 9)
    print *, lab(5), (cDW_Lg4(icd))

    read (17, *)  lab(6), (cSPc_temp(j), j = 1, 9)
    print *, lab(6), (cSPc_temp(icd))
    read (17, *)  lab(7), (cWIc_temp(j), j = 1, 9)
    print *, lab(7), (cWIc_temp(icd))
    read (17, *)  lab(9), (cSUp_temp(j), j = 1, 9)
    print *, lab(9), (cSUp_temp(icd))
    read (17, *)  lab(8), (cFAp_temp(j), j = 1, 9)
    print *, lab(8), (cFAp_temp(icd))
    read (17, *) lab(10), (cSPp_temp(j), j = 1, 9)
    print *, lab(10), (cSPp_temp(icd))

    read (17, *) lab(11), (cPr_Lg1(j), j = 1, 9)
    print *, lab(11), (cPr_Lg1(icd))
    read (17, *) lab(12), (cPr_Lg2(j), j = 1, 9)
    print *, lab(12), (cPr_Lg2(icd))
    read (17, *) lab(13), (cPr_Lg3(j), j = 1, 9)
    print *, lab(13), (cPr_Lg3(icd))
    read (17, *) lab(14), (cPr_Lg4(j), j = 1, 9)
    print *, lab(14), (cPr_Lg4(icd))

    read (17, *) lab(15), (cSPc_prcp(j), j = 1, 9)
    print *, lab(15), (cSPc_prcp(icd))
    read (17, *) lab(16), (cWIc_prcp(j), j = 1, 9)
    print *, lab(16), (cWIc_prcp(icd))
    read (17, *) lab(18), (cSUp_prcp(j), j = 1, 9)
    print *, lab(18), (cSUp_prcp(icd))
    read (17, *) lab(17), (cFAp_prcp(j), j = 1, 9)
    print *, lab(17), (cFAp_prcp(icd))
    read (17, *) lab(19), (cSPp_prcp(j), j = 1, 9)
    print *, lab(19), (cSPp_prcp(icd))

    read (17, *) lab(36), (cDWLg1xDayLite_Lg1(j), j = 1, 9)
    read (17, *) lab(37), (cDayLite_Lg1(j), j = 1, 9)

    read (17, *) lab(20), (cDWLg1xPLg1(j), j = 1, 9)
    read (17, *) lab(21), (cDWLg1xPLg2(j), j = 1, 9)
    read (17, *) lab(22), (cDWLg1xPLg3(j), j = 1, 9)

    print *, "lab(23).........."
    print *, lab(23)

    read (17, *) lab(23), (cDWLg1xPLg4(j), j = 1, 9)

    print *, "lab(23).........."
    print *, lab(23)
    print *, "lab(24).........."
    print *, lab(24)

    do i = 1, 9
        print *, cDWLg1xPLg4(i), cDWLg2xPLg1(i)
    enddo

    read (17, *) lab(24), (cDWLg2xPLg1(j), j = 1, 9)

    read (17, *) lab(25), (cDWLg2xPLg2(j), j = 1, 9)
    read (17, *) lab(26), (cDWLg2xPLg3(j), j = 1, 9)
    read (17, *) lab(27), (cDWLg2xPLg4(j), j = 1, 9)

    read (17, *) lab(28), (cDWLg3xPLg1(j), j = 1, 9)
    read (17, *) lab(29), (cDWLg3xPLg2(j), j = 1, 9)
    read (17, *) lab(30), (cDWLg3xPLg3(j), j = 1, 9)
    read (17, *) lab(31), (cDWLg3xPLg4(j), j = 1, 9)

    read (17, *) lab(32), (cDWLg4xPLg1(j), j = 1, 9)
    read (17, *) lab(33), (cDWLg4xPLg2(j), j = 1, 9)
    read (17, *) lab(34), (cDWLg4xPLg3(j), j = 1, 9)
    read (17, *) lab(35), (cDWLg4xPLg4(j), j = 1, 9)
    close(17)

    !        if (icd .eq. 1 .or. icd .eq. 2) icdi=10
    !        if (icd .eq. 3 .or. icd .eq. 4 .or. icd .eq. 5) icdi=11
    ! 2016
    !        if (icd .eq. 1) icdi=1
    !        if (icd .eq. 2) icdi=2
    !        if (icd .eq. 3) icdi=11
    !        if (icd .eq. 4) icdi=11
    !        if (icd .eq. 5) icdi=11
    !        if (icd .eq. 6 .or. icd .eq. 7) icdi=12
    !        if (icd .eq. 8 .or. icd .eq. 9) icdi=13
    ! 2017
    if (icd .eq. 1) icdi = 1
    if (icd .eq. 2) icdi = 2
    if (icd .eq. 3) icdi = 3
    if (icd .eq. 4) icdi = 4
    if (icd .eq. 5) icdi = 5
    if (icd .eq. 6) icdi = 6
    if (icd .eq. 7) icdi = 7
    if (icd .eq. 8) icdi = 8
    if (icd .eq. 9) icdi = 9

    print *, 'intercept(icdi),icdi = ', intercept(icdi), ' ', icdi
    print *, intercept(icdi), &
            cDWLg4xPLg3(icdi), cDWLg4xPLg4(icdi)

    ! compute mir
    ! future week = Lag0 (i+1)
    ! Lg1 = this past week, the most recent week (i)
    ! Lg2 = 2 weeks back (i-1)
    ! Lg3 = 3 weeks back (i-2)
    ! Lg4 = 4 weeks back (i-3)
    ! Lg0 = next week

    ! MIR computed for next week
    ! CD 2: Ind Model weeks 18-38

    if (endk .gt. 44) endk = 44
    year = iyear

    do i = 1, 53
        mir(i) = 0.0
        mir_cur(i) = -99.99
    enddo

    write (3, *) '  i  i+1  CD iyr yr  week  +MIR +avg_mir +mir_cur +&
            df_SPc_t df_WIc_t df_FAp_t df_SUp_t df_SPp_t'
    dum1 = 0.0
    dum2 = 0.0
    dum3 = 0.0
    dum4 = -99.99
    dum5 = -99.99
    dum6 = -99.99
    dum7 = -99.99
    ! print week 1-17, current week weather and daylight and zero MIR
    do i = 1, 17
        write (4, 60)i, i, icd, iyear, year, wk(i), dum1, dum2, dum3, dum4, &
                df_SUp_p, df_SPp_p, df_SPc_t, df_WIc_t, df_FAp_t, df_SUp_t, df_SPp_t

        if (iyear .eq. 2017) then
            write (8, 83)icd, year, wk(i), pmon(i), pday(i), dum1, dum2, dum3, &
                    dum5, dum6, dum7
        endif
    enddo

    print *, 'computing mir for next week with this weeks weather'
    print *, 'and with next weeks daylight'

    print *, 'predicted (next week (i+1), start at 18 (17+1) to 44+1'
    print *, 'printing next week mir(i+1) wx(i) and daylight(i+1)'
    print *, ' and this weeks weather(i)'

    ! lag1 (i)
    ! lag2 (i-1)
    ! lag3 (i-2)
    ! lag4 (i=3)

    do i = 17, endk
        mir(i + 1) = (intercept(icdi) + &
                cDWLg4xPLg4(icdi) * df_wkd(i - 3) * df_wkp(i - 3))

        ! adjust next weeks mir(i+1) with next weeks avg_mir

        mir_cur(i + 1) = mir(i + 1) + cd_avg_mir(i + 1, icdi)
        print *, 'i+1, mir_cur, mir, avg_mir = ', i + 1, &
                mir_cur(i + 1), mir(i + 1), cd_avg_mir(i + 1, icdi)

        if(mir_cur(i + 1) .lt. 0.0) mir_cur(i + 1) = 0.0

        ! weekly weather output
        ! open(unit=3,file=ST//CD//'w_'//ayear
        write (3, 60)i, i + 1, icd, iyear, year, wk(i + 1), &
                df_spc_t, df_WIc_t, df_FAp_t, df_SUp_t, df_SPp_t

        ! weekly MIR file
        ! open(unit=4,file=ST//CD//'_'//ayear//'_mir'
        write (4, 60)i, i + 1, icd, iyear, year, wk(i + 1), &
                df_SPc_t, df_WIc_t, df_FAp_t, df_SUp_t, df_SPp_t
        60      format (2(1x, i2), 4(1x, i4), 7(1x, f6.2, 1x), 1x, 8(1x, f6.2))

        if (iyear .eq. 2017) then
            ! MIR plotting file
            ! open(unit=8,file=ST//CD//'_'//ayear//'_plot'
            write (8, 83)icd, year, wk(i + 1), wmonth(i + 1), wday(i + 1), &
                    ACT_mir(i + 1, icd), ACTm_mir(i + 1, icd), ACTp_mir(i + 1, icd)
            print *, icd, year, wk(i + 1), ACT_mir(i + 1, icd), ACTm_mir(i + 1, icd), &
                    ACTp_mir(i + 1, icd)
            83      format (i4, 4(1x, i4), 6(1x, f6.2))
        endif
    enddo

    !        if (iyear .eq. 2017) then
    ! post weeks
    dum1 = -99.99
    dum5 = -99.99
    dum6 = -99.99
    dum7 = -99.99
    do i = endk + 2, 45
        write (8, 83)icd, year, wk(i), wmonth(i), wday(i), &
                cd_2012mir(i, icd), cd_avg_mir(i, icdi), dum1, dum5, dum6, dum7
    enddo
    !        endif

    dum1 = 0.0
    dum2 = 0.0
    dum3 = 0.0
    dum4 = -99.99
    do i = 46, 52
        print *, "here now"
        print *, i, i, icd, iyear, year, wk(i), wmonth(i), wday(i), real(wday(i))


        !        write (4,60)i,i,icd,iyear,year,wk(i),real(wmonth(i)),
        !      ! real(wday(i))

        !        write (4,60)i,i,icd,iyear,year,wk(i),real(wmonth(i))

        write (4, 60)i, i, icd, iyear, year, wk(i), real(wmonth(i)), &
                df_SPc_t, df_WIc_t, df_FAp_t, df_SUp_t, df_SPp_t

        !        write (4,60)i,i,icd,iyear,year,wk(i),wmonth(i),wday(i),
        !     ! dum3,dum4,cd_dylt(i,icd),df_wkd(i),df_wkp(i),
        !     ! df_SPc_p,df_WIc_p,df_FAp_p,df_SUp_p,df_SPp_p,
        !     ! df_SPc_t,df_WIc_t,df_FAp_t,df_SUp_t,df_SPp_t

        if (iyear .eq. 2017) then
            write (8, 83)icd, year, wk(i), wmonth(i), wday(i), dum1, dum2, dum3, &
                    dum5, dum6, dum7
        endif
    enddo

    stop
end