./cmp_n_nweek <state>
| run model | example |
|---|---|
| ./cmp_n_week <state> <climate division> <year> | ./cmp_n_week 11 09 2017 |
| Code | Comments |
|---|---|
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 FILE
|
! input for normals (1981-2010 daily)
open(unit = 3, file = ST // CD // 'w_' // ayear, &
status = 'unknown', form = 'formatted')
| OUTPUT FILE e.g., 1109w_2017 SPRING Prev_YR NORM DIF for T D and P 14.57, 13.22, 1.35, 0.00, 0.00, 0.00, 2.95, 2.52, 0.43 SUMMER Prev_YR NORM DIF for T D and P 25.48, 24.38, 1.09, 0.00, 15.34, 0.00, 3.47, 1.80, 1.67 FALL Prev_YR NORM DIF for T D and P 16.03, 14.05, 1.98, 0.00, 32.28, 0.00, 1.59, 2.03, -0.44 |
! 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'
| OUTPUT FILE |
! 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'
| OUTPUT FILE |
! 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 NORMALS |
! 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 WEEKLY DAILY AVERAGES d,t? |
! 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.
| COMPUTE NORMALS |
! 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
| GET PREVIOUS YEAR |
! 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
| |