Fortran Model Code

program compute
    week
    implicit none
    ! compute: norms curr_week_diffs and MIR

    ! compute 30-yr (1981-2010) daily weather normal data from each Illinois CD

    ! for analysis year and previous year compute based on CDC EPI week
    ! in 2016 and 2017, begin day use = saturday to friday
    ! in 2018, change to sunday to Saturday

    ! CDC defines ending saturday: The first epi week of the year ends,
    !  by definition, on the first Saturday of January,
    !  as long as it falls at least four days into the month.
    !  Each epi week begins on a Sunday and ends on a Saturday.

    ! get daily data for analysis year.
    !   compute weekly normals and weekly current year observations.
    !   compute weekly difference between analysis year and normals.
    ! weekly temp avg ( C);
    ! weekly precip tot (cm);
    ! weekly degree day avg based on 22 C. (accumate over time from Jan 1 to dec 31)

    ! first saturday set up through 2020.
    ! 2017, change to first sunday
    ! VET-MED first epi week ends on first saturday of january, unless falls after
    ! January 3, then ends in corresponding december day.

    ! as per Karki:
    ! begin jan 1 all years and go thru 365 days.
    ! weeks: 49-52,1-9 = WInter
    !        10-22 = SPring (CDC spring ends week 22)
    !        23-35 = SUmmer
    !        36-48 = FAll

    ! bflg=0
    ! bflg=1

    character*1 c
    character*2 ST, CD, cstate
    character*4 ayear

    integer icd, iyear, year, mon, day, i, j, k, m, bflg
    integer mm(366), dd(366)
    integer icdi
    real t(366, 30), p(366, 30), d(366, 30), dum_d, dsum_d
    real wt(53, 30), wp(53, 30), wd(53, 30)

    real tmax, tmin, tmean, prcp
    real wknorm_t(53), wknorm_p(53), wknorm_d(53), wknorm_np(53)
    real wksum_t(53), wksum_p(53), wksum_d, wksum_np(53)

    real tmn(366), ppt(366), ddd(366)

    real wctmn(53), wcppt(53), wcddd(53), wcsum_d
    real df_wkt(53), df_wkp(53), df_wkd(53), df_wknp(53)

    !       integer begday(40),begmon(40)
    integer begyr(40), ndays(40)
    integer nbegday(40), nbegyr(40), nbegmon(40)
    integer bday, bmon, nday, byear, n, end, endk

    real WInorm_t, SPnorm_t, SUnorm_t, FAnorm_t
    real WIsum_t, SPsum_t, SUsum_t, FAsum_t
    real WInorm_d, SPnorm_d, SUnorm_d, FAnorm_d
    real WIsum_d, SPsum_d, SUsum_d, FAsum_d

    real WInorm_p, SPnorm_p, SUnorm_p, FAnorm_p
    real WIsum_p, SPsum_p, SUsum_p, FAsum_p

    real WInorm_np, SPnorm_np, SUnorm_np, FAnorm_np
    real WIsum_np, SPsum_np, SUsum_np, FAsum_np

    real SP_prev_yr_t, SP_prev_yr_d, SP_prev_yr_p, SP_prev_yr_np
    real SU_prev_yr_t, SU_prev_yr_d, SU_prev_yr_p, SU_prev_yr_np
    real FA_prev_yr_t, FA_prev_yr_d, FA_prev_yr_p, FA_prev_yr_np

    real SP_prev_sum_t, SP_prev_sum_d, SP_prev_sum_p, SP_prev_sum_np
    real SU_prev_sum_t, SU_prev_sum_d, SU_prev_sum_p, SU_prev_sum_np
    real FA_prev_sum_t, FA_prev_sum_d, FA_prev_sum_p, FA_prev_sum_np

    real WI_curr_yr_t, WI_curr_yr_d, WI_curr_yr_p, WI_curr_yr_np
    real SP_curr_yr_t, SP_curr_yr_d, SP_curr_yr_p, SP_curr_yr_np

    real WI_curr_sum_t, WI_curr_sum_d, WI_curr_sum_p, WI_curr_sum_np
    real SP_curr_sum_t, SP_curr_sum_d, SP_curr_sum_p, SP_curr_sum_np

    real df_WIc_t, df_WIc_p, df_WIc_d, df_WIc_np
    real df_SPc_t, df_SPc_p, df_SPc_d, df_SPc_np
    real df_SPp_t, df_SPp_p, df_SPp_d, df_SPp_np
    real df_SUp_t, df_SUp_p, df_SUp_d, df_SUp_np
    real df_FAp_t, df_FAp_p, df_FAp_d, df_FAp_np

    character *86 line1
    integer wk(53), cumwk(53)
    real cd_dylt(53, 9), cd_avg_mir(53, 13), cd_2012mir(53, 13)

    real dum1, dum2, dum3, dum4, dum5, dum6, dum7
    real mir(53), daylit(53), mir_cur(53)
    real mir2012(53, 13)
    real ACT_mir(52, 9), ACTm_mir(52, 9), ACTp_mir(52, 9)
    real cd_act_mir(53, 13)

    real intercept(13)
    real cDW_Lg1(13), cDW_Lg2(13), cDW_Lg3(13), cDW_Lg4(13)
    real cPr_Lg1(13), cPr_Lg2(13), cPr_Lg3(13), cPr_Lg4(13)
    real cDWLg1xPLg1(13), cDWLg1xPLg2(13)
    real cDWLg1xPLg3(13), cDWLg1xPLg4(13)
    real cDWLg2xPLg1(13), cDWLg2xPLg2(13)
    real cDWLg2xPLg3(13), cDWLg2xPLg4(13)
    real cDWLg3xPLg1(13), cDWLg3xPLg2(13)
    real cDWLg3xPLg3(13), cDWLg3xPLg4(13)
    real cDWLg4xPLg1(13), cDWLg4xPLg2(13)
    real cDWLg4xPLg3(13), cDWLg4xPLg4(13)
    real cSPc_temp(13), CSPc_prcp(13)
    real cWIc_temp(13), cSPp_temp(13), cSUp_temp(13), cFAp_temp(13)
    real cWIc_prcp(13), cSPp_prcp(13), cSUp_prcp(13), cFAp_prcp(13)
    real cDayLite(13), cDWLg1xDayLite_Lg1(13)
    real cDayLite_Lg1(13), cDayLite_Lg2(13)
    real cDayLite_Lg3(13), cDayLite_Lg4(13)

    integer ilcd, rundate
    character*10 valdate
    character*6 cprcp6
    character*5 lab1
    character*20 lab(40)
    real prcpd(10, 9), mxtemp(10, 9), mntemp(11, 9), prcp6(13, 9)
    integer n_mxT, n_mnT, n_ppt
    integer ncurr
    integer firstmon, firstday
    character*2 cfirstmon, cfirstday

    integer pday(53), pmon(53), pdd(366), pmm(366)
    character*2 cpdd, cpmm
    integer nwk(53)
    integer pm, pd
    integer wmonth(52), wday(52), wweek(52), wyear
    character*10 wdate
    character*3 wdow
    integer emon, eday

    ! used to comply with epi week beginning in December
    ! for previous year (iyear-1) and analysis year (iyear)
    ! for 2016 and 2017
    data nbegday/3, 2, 1, 31, 29, 28, 3, 2, 31, 30, &
            1, 31, 29, 28, 3, 2, 31, 30, 29, 28/
    data nbegmon/1, 1, 1, 12, 12, 12, 1, 1, 12, 12, &
            1, 12, 12, 12, 1, 1, 12, 12, 12, 12/
    data nbegyr/1981, 1982, 1983, 1983, 1984, 1985, 1987, 1988, 1988, 1989, &
            2011, 2011, 2012, 2013, 2015, 2016, 2016, 2017, 2018, 2019/


    ! for 2018
    ! defining years, months, days, epi weeks for previous year and analysis years
    !        data nbegday/4,3,2,1,30,29,4,3,1,31,
    !     !         30,29,3,2,1,31,30,4,3,2,
    !     !         31,30,29,4,2,1,31,30,4,3,
    !     !         2,1,30,29,4,3,1,31,30,29/
    !        data nbegmon/1,1,1,1,12,12,1,1,1,12,
    !     !         12,12,1,1,1,12,12,1,1,1,
    !     !         12,12,12,1,1,1,12,12,1,1,
    !     !         1,1,12,12,1,1,1,12,12,12/
    !        data nbegyr/1981,1982,1983,1984,1984,1985,1987,1988,1989,1989,
    !     !        1990,1991,1993,1994,1995,1995,1996,1998,1999,2000,
    !     !        2000,2001,2002,2004,2005,2006,2006,2007,2009,2010,
    !     !        2011,2012,2012,2013,2015,2016,2017,2017,2018,2019/

    ! USE REGARDLESS
    data begyr/1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, &
            2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020/
    data ndays/365, 365, 365, 366, 365, 365, 365, 366, 365, 365, &
            365, 366, 365, 365, 365, 366, 365, 365, 365, 366/
    call getarg (1, ST)
    call getarg (1, ST)
    call getarg (2, CD)
    call getarg (3, ayear)

    ! read from character to integer
    read(ayear, 25) iyear
    25    format(i4)
    read(CD, 26) icd
    26    format(i2)
    !        print *, 'CD = ',icd,' ST = ',ST,'iyear= ',iyear

    c = ','

    do j = 1, 30
        do i = 1, 366
            t(i, j) = 999.9
            d(i, j) = 0.0
            p(i, j) = 999.9
        enddo
    enddo

! input for days of week for current year
open(unit = 2, file = 'week_days_2016.csv', status = 'old', &
        form = 'formatted')
read (2, *) line1
do i = 1, 52
    read (2, *) wmonth(i), wday(i), wweek(i), wyear, wdate, wdow
enddo
close(2)
!MONTH,DAY,WEEK,YEAR,DATE,DOW
!1,2,1,2016,1/2/2016,SAT

! 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