program compute week
implicit none
c compute: norms curr_week_diffs and MIR
c compute 30-yr (1981-2010) daily weather normal data from each Illinois CD
c for analysis year and previous year compute based on CDC EPI week
c in 2016 and 2017, begin day use = saturday to friday
c in 2018, change to sunday to Saturday
c CDC defines ending saturday: The first epi week of the year ends,
c by definition, on the first Saturday of January,
c as long as it falls at least four days into the month.
c Each epi week begins on a Sunday and ends on a Saturday.
c get daily data for analysis year.
c compute weekly normals and weekly current year observations.
c compute weekly difference between analysis year and normals.
c weekly temp avg ( C);
c weekly precip tot (cm);
c weekly degree day avg based on 22 C. (accumate over time from Jan 1 to dec 31)
c first saturday set up through 2020.
c 2017, change to first sunday
c VET-MED first epi week ends on first saturday of january, unless falls after
c January 3, then ends in corresponding december day.
c as per Karki:
c begin jan 1 all years and go thru 365 days.
c weeks: 49-52,1-9 = WInter
c 10-22 = SPring (CDC spring ends week 22)
c 23-35 = SUmmer
c 36-48 = FAll
c bflg=0
c 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)
c 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
c used to comply with epi week beginning in December
c for previous year (iyear-1) and analysis year (iyear)
c for 2016 and 2017
data nbegday/3,2,1,31,29,28,3,2,31,30,
! 29,28,2,1,31,30,29,3,2,1,
! 30,29,28,3,1,31,30,29,3,2,
! 1,31,29,28,3,2,31,30,29,28/
data nbegmon/1,1,1,12,12,12,1,1,12,12,
! 12,12,1,1,12,12,12,1,1,1,
! 12,12,12,1,1,12,12,12,1,1,
! 1,12,12,12,1,1,12,12,12,12/
data nbegyr/1981,1982,1983,1983,1984,1985,1987,1988,1988,1989,
! 1990,1991,1993,1994,1994,1995,1996,1998,1999,2000,
! 2000,2001,2002,2004,2005,2005,2006,2007,2009,2010,
! 2011,2011,2012,2013,2015,2016,2016,2017,2018,2019/
c for 2018
c defining years, months, days, epi weeks for previous year and analysis years
c data nbegday/4,3,2,1,30,29,4,3,1,31,
c ! 30,29,3,2,1,31,30,4,3,2,
c ! 31,30,29,4,2,1,31,30,4,3,
c ! 2,1,30,29,4,3,1,31,30,29/
c data nbegmon/1,1,1,1,12,12,1,1,1,12,
c ! 12,12,1,1,1,12,12,1,1,1,
c ! 12,12,12,1,1,1,12,12,1,1,
c ! 1,1,12,12,1,1,1,12,12,12/
c data nbegyr/1981,1982,1983,1984,1984,1985,1987,1988,1989,1989,
c ! 1990,1991,1993,1994,1995,1995,1996,1998,1999,2000,
c ! 2000,2001,2002,2004,2005,2006,2006,2007,2009,2010,
c ! 2011,2012,2012,2013,2015,2016,2017,2017,2018,2019/
c USE REGARDLESS
data begyr/1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,
! 1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,
! 2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,
! 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,
! 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)
c read from character to integer
read(ayear,25) iyear
25 format(i4)
read(CD,26) icd
26 format(i2)
c 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
c 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)
cMONTH,DAY,WEEK,YEAR,DATE,DOW
c1,2,1,2016,1/2/2016,SAT
c input for normals (1981-2010 daily)
open(unit=3,file=ST//CD//'w_'//ayear,
! status='unknown',form='formatted')
c 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_
!mir day_lt df_wkd df_wkp aSPc_p aWIc_p aFAp_p aSUp_p aSPp_p
! aSPc_t aWIc_t aFAp_t aSUp_t aSPp_t'
c 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'
c NORMALS computations
c put data in daily,yearly array to compute normalso
c same normals for each analysis year
c week 1, 8 days, dec 29-Jan 5;
c all other weeks=7 days (skip day 365 on new year)
c read into array
print *, '1 CD = ',CD
print *, ST//CD//'n'
open(unit=10,file=ST//CD//'n',status='old', form='formatted')
c for normals
bmon = 12
bday = 29
emon = 12
eday = 28
nday=365
c becasue skip 2/29
do j=1,30
bflg=0
c 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
c 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
c write (*,*) i,year,mon,day,tmax,tmin,prcp,byear,bmon,bday
c compute daily mean
c convert to celcius from farenheit
c compute degree days, base 22 C
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
c go to next year at nday=365
998 enddo
c close(10)
enddo
999 close(10)
c 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
c 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
c 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))
c print *, i,j,k,t(i,j),t((i+1),j),t((i+2),j),t((i+3),j),t((i+4),j),
c ! t((i+5),j),t((i+6),j)
c print *,'weekly sum and avg ppt_cm',i,j,k,wp(k,j),(wp(k,j)/7.0)
c print *,i,j,k,p(i,j),p((i+1),j),p((i+2),j),p((i+3),j),p((i+4),j),
c ! p((i+5),j),p((i+6),j),wp(k,j),(wp(k,j)/7.0)
enddo
enddo
c 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
c 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)
c wksum_d(k)=wksum_d(k)+wd(k,j)
enddo
wknorm_t(k)=wksum_t(k)/30.
wknorm_p(k)=wksum_p(k)/30.
c wknorm_d(k)=wksum_d(k)/30.
c 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)
c wknorm_d(53)=wknorm_d(1)
c compute weekly dw 30-yr normals
do k=1,53
wknorm_d(k)=0.0
enddo
do k=1,53
c 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
c compute seasonal weekly 30-year normals
do k=1,9
WIsum_t =WIsum_t +wknorm_t(k)
WIsum_p =WIsum_p +wknorm_p(k)
c 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)
c 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)
c 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)
c 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)
c FAsum_np=FAsum_np+wknorm_np(k)
FAsum_d =FAsum_d +wknorm_d(k)
enddo
c 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.
print *, 'seasonal Weekly Norm Temp deg_day prcp'
print *, 'WInter norm',' ',WInorm_t,' ',WInorm_d,' ',WInorm_p
print *, 'SPring norm',' ',SPnorm_t,' ',SPnorm_d,' ',SPnorm_p
print *, 'SUmmer norm',' ',SUnorm_t,' ',SUnorm_d,' ',SUnorm_p
print *, 'FAll norm ',' ',FAnorm_t,' ',FAnorm_d,' ',FAnorm_p
c print *, 'WInter norm',' ',WInorm_p,' ',WInorm_np
c print *, 'SPring norm',' ',SPnorm_p,' ',SPnorm_np
c print *, 'SUmmer norm',' ',SUnorm_p,' ',SUnorm_np
c print *, 'FAll norm',' ',FAnorm_p,' ',FAnorm_np
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 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
c normals are done
c previouw year = iyear-1
c 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
c18 write (*,*) i,iyear,year,mon,day,tmax,tmin,prcp,bflg
18 n=n+1
if (tmax .eq. 999.9) go to 19
c mean temperature
c convert to centigrade
c 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
c write (*,*) i,iyear,year,mon,day,tmn(i),ppt(i)
enddo
19 close(1)
print *,'iyear,n,nday = ',iyear,n,nday
c 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
c write (*,*) k,i,i+1,i+2,i+3,i+4,i+5,i+6
c print *, k, tmn(i),tmn(i+1),tmn(i+2),tmn(i+3),tmn(i+4),
c ! tmn(i+5),tmn(i+6)
c print *, k, ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
c ! 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))
c print *, 'prev_ yr weekly sum avg ppt ',wcppt(k),(wcppt(k)/7.0)
c write (*,*) ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
c ! ppt(i+5),ppt(i+6),wcppt(k),(wcppt(k)/7.0)
c 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
c 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
c 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
c 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
c 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
c 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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 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
c 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
c write (*,*) cstate,icd,year,mon,day,tmax,tmin,prcp,'begin',
c ! byear,bday,bmon
if ((tmax .eq. -99 .or. tmin .eq. -99 .or. prcp .lt. 0.0)
! .and. n.gt.140) then
c for forecast
firstday=day
firstmon=mon
write(cfirstmon,11) firstmon
write(cfirstday,11) firstday
11 format(i2)
c 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
c 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
c mean temperature
c convert to centigrade
c compute degree days
c 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
c changed to accomodate temperature good, but precipitation missing!
c5-19-2017
if (ppt(n) .lt. 0.0) ppt(n)=0.0
c if (ppt(n) .gt. 1000) ppt(n)=0.0
c if (ppt(n) .lt. 0) ppt(n)=999.9
c write (*,*) i,iyear,year,mon,day,tmn(i),ppt(i)
enddo
c missing data or end of year, so close
88 close(1)
ncurr=n
print *,'iyear,n,ncurr,nday = ',iyear,n,ncurr,nday
c ADD FORECAST DATA
c input - get 10-day max and min temp forecasts, and 3 day precip forecasts
c for this year 2017:
if (iyear .eq. 2017) then
c 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
c set later days of precip to 0.0
do i=4,10
prcpd(i,j)=0.0
enddo
enddo
c see /home/nan/WNV_WX/Not_Smooth/run.get.forecast
c open(unit=7,file='ILCD_divisions_ALL',
c ! 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
cDIV,CYCLE,ELEMENT,VALIDTIME,MEAN-VALUE
c1,2014041412,MINT,2014041512,23.115
c 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
c 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'
c 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
c go to next cliamte division
enddo
c 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
c end of forecast data for 2017
c compute weekly averages for model analysis year
do k=1,53
wctmn(k)=999.9
wcppt(k)=999.9
enddo
k=1
c do i=1,n-1,7
do i=1,ncurr+8,7
c only go out 8 days (should match dupage model, will have 3 days with precip forecast
pday(k)=pdd(i)
pmon(k)=pmm(i)
c pday(k)=pdd(i+4)
c 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+2) .ge. 999.9 .or. tmn(i+3) .ge. 999.9 .or.
! tmn(i+4) .ge. 999.9 .or. tmn(i+5) .ge. 999.9 .or.
! tmn(i+6) .ge. 999.9) then
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)
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
c wcppt(k) = (ppt(i)+ppt(i+1)+ppt(i+2)+ppt(i+3)+ppt(i+4)
c ! +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)
c write (*,*) k,i,i+1,i+2,i+3,i+4,i+5,i+6
c print *,'cur_yr weekly avg sum ppt ',k,wcppt(k),(wcppt(k)/7.0)
c print *, k, tmn(i),tmn(i+1),tmn(i+2),tmn(i+3),tmn(i+4),
c ! tmn(i+5),tmn(i+6),wctmn(k),pmon(k),pday(k)
c print *, k, ddd(i),ddd(i+1),ddd(i+2),ddd(i+3),ddd(i+4),
c ! ddd(i+5),ddd(i+6),wcddd(k)
c print *, k, ppt(i),ppt(i+1),ppt(i+2),ppt(i+3),ppt(i+4),
c ! ppt(i+5),ppt(i+6),wcppt(k)
k=k+1
79 continue
endk=k-1
print *,'endk = ',endk
enddo
c 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
c summed december already from previous year
c WI_curr_sum_t=0.0
c WI_curr_sum_d=0.0
c 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
c 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
c WI_curr_yr_p=WI_curr_sum_p
endif
c print *, 'DecFeb WI_curr_sum_p = ',WI_curr_sum_p
c 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
c 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))
c print *,'k,df_wkp ',k,df_wkp(k)
c print *,'year,endk,k,wcppt(k),wknorm_p(k),df_wkp(k) ',
c ! year,endk,k,wcppt(k),wknorm_p(k),df_wkp(k)
c print *,'year,endk,k,wcddd(k),wknorm_d(k) ',
c ! year,endk,k,wcddd(k),wknorm_d(k)
enddo
c what is this!!
c if (iyear .eq. 2015)then
c df_wkd(endk)=0.0
c df_wkp(endk)=0.0
cc print *,'endk,df_wkp ',k,df_wkp(endk)
c 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)
c 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)
c 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
c 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,
! wcddd(k),c,wknorm_d(k),c,df_wkd(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,
! wcddd(k),c,wknorm_d(k),c,df_wkd(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
c 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'
c51 format (1x,3(f6.2,a1),4x,3(f6.2,a1),4x,3(f6.2,a1),4x,3(f6.2,a1))
c write (3,*) 'SP_curr_t SPnrm_t df_SPc_t SP_curr_d SPnrm_d df_q
c !2c_d SP_curr_p SPnrm_p df_SPc_p'
c write (3,51) SP_curr_yr_t,c,SPnorm_t,c,df_SPc_t,c, SP_curr_yr_d,
c ! c,SPnorm_d,c,df_SPc_d,c, SP_curr_yr_p,c,SPnorm_p,c,df_SPc_p
c ! c,SP_curr_yr_np,c,SPnorm_np,c,df_SPc_np,c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input - climate division daylight data
c 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)
c print *, wk(i),(cd_dylt(i,j),j=1,9)
enddo
close(14)
cWEEK,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
c1,9.228,9.228,9.367,9.367,9.350,9.467,9.489,9.600,9.600
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input - MIR 9-year average by 9 CDs and (4) combined areas
c 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)
c print *, 'each year',wk(i),(cd_avg_mir(i,j),j=1,13)
enddo
close(15)
cWeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08
c1,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 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)
c print *, '2012 mir1',wk(i),(cd_2012mir(i,j),j=1,13)
enddo
print *,'act_mir year and iyear = ',year,' ',iyear
c input - 2005-2013 MIR by CD
c do i=1,53
c do i=1,13
c cd_act_mir=0.0
c enddo
c enddo
c
c if (iyear .lt. 2005) go to 228
c if (iyear .gt. 2013) go to 228
c open(unit=26,file='CD_ACT_mir_2005_13.csv',status='old',
c ! form='formatted')
c read (26,*) line1
c i=1
c226 read (26,*,end=228) year,wk(i),cumwk(i),(cd_act_mir(i,j),j=1,9)
c print *, i,year,wk(i),(cd_act_mir(i,j),j=1,9)
c if (year .lt. iyear) go to 226
c if (year .gt. iyear) go to 227
c print *, year,' ',iyear,'cd =',j,'act_mir = ',cd_act_mir(i,j)
c i=i+1
c go to 226
c228 close(26)
cWeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08
c1,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
c input - current year (2017) ACTUAL MIR by CD
c 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
cWeekNum,CD_02,CD_02,CD_05,CD_05,CD_05,Cd_08,Cd_08,Cd_08,Cd_08
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 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)
c if (icd .eq. 1 .or. icd .eq. 2) icdi=10
c if (icd .eq. 3 .or. icd .eq. 4 .or. icd .eq. 5) icdi=11
c 2016
c if (icd .eq. 1) icdi=1
c if (icd .eq. 2) icdi=2
c if (icd .eq. 3) icdi=11
c if (icd .eq. 4) icdi=11
c if (icd .eq. 5) icdi=11
c if (icd .eq. 6 .or. icd .eq. 7) icdi=12
c if (icd .eq. 8 .or. icd .eq. 9) icdi=13
c 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),
! cDW_Lg1(icdi),cDW_Lg2(icdi),cDW_Lg3(icdi),cDW_Lg4(icdi),
! cSPc_temp(icdi),cWIc_temp(icdi),
! cFAp_temp(icdi),cSUp_temp(icdi),cSPp_temp(icdi),
! cPr_Lg1(icdi),cPr_Lg2(icdi),cPr_Lg3(icdi),cPr_Lg4(icdi),
! cSPc_prcp(icdi),cWIc_prcp(icdi),
! cFAp_prcp(icdi),cSUp_prcp(icdi),cSPp_prcp(icdi),
! cDWLg1xDayLite_Lg1(icdi),
! cDayLite_Lg1(icdi),
! cDWLg1xPLg1(icdi),cDWLg1xPLg2(icdi),
! cDWLg1xPLg3(icdi),cDWLg1xPLg4(icdi),
! cDWLg2xPLg1(icdi),cDWLg2xPLg2(icdi),
! cDWLg2xPLg3(icdi),cDWLg2xPLg4(icdi),
! cDWLg3xPLg1(icdi),cDWLg3xPLg2(icdi),
! cDWLg3xPLg3(icdi),cDWLg3xPLg4(icdi),
! cDWLg4xPLg1(icdi),cDWLg4xPLg2(icdi),
! cDWLg4xPLg3(icdi),cDWLg4xPLg4(icdi)
c compute mir
c future week = Lag0 (i+1)
c Lg1 = this past week, the most recent week (i)
c Lg2 = 2 weeks back (i-1)
c Lg3 = 3 weeks back (i-2)
c Lg4 = 4 weeks back (i-3)
c Lg0 = next week
c MIR computed for next week
c 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 +
!act_mir +day_lt dw_wkd df_wkp 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'
dum1=0.0
dum2=0.0
dum3=0.0
dum4=-99.99
dum5=-99.99
dum6=-99.99
dum7=-99.99
c 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,
! 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),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)'
c lag1 (i)
c lag2 (i-1)
c lag3 (i-2)
c lag4 (i=3)
do i=17,endk
mir(i+1)= ( intercept(icdi) +
! cDW_Lg1(icdi)*df_wkd(i) + cDW_Lg2(icdi)*df_wkd(i-1) +
! cDW_Lg3(icdi)*df_wkd(i-2) + cDW_Lg4(icdi)*df_wkd(i-3) +
! cSPc_temp(icdi)*df_SPc_t +
! cWIc_temp(icdi)*df_WIc_t + cFAp_temp(icdi)*df_FAp_t +
! cSUp_temp(icdi)*df_SUp_t + cSPp_temp(icdi)*df_SPp_t +
! cPr_Lg1(icdi)*df_wkp(i) + cPr_Lg2(icdi)*df_wkp(i-1) +
! cPr_Lg3(icdi)*df_wkp(i-2) + cPr_Lg4(icdi)*df_wkp(i-3) +
! cSPc_prcp(icdi)*df_SPc_p +
! cWIc_prcp(icdi)*df_WIc_p + cFAp_prcp(icdi)*df_FAp_p +
! cSUp_prcp(icdi)*df_SUp_p + cSPp_prcp(icdi)*df_SPp_p +
! cDWLg1xDayLite_Lg1(icdi)*df_wkd(i)*cd_dylt(i,icd) +
! cDayLite_Lg1(icdi)*cd_dylt(i,icd) +
! cDWLg1xPLg1(icdi)*df_wkd(i)*df_wkp(i) +
! cDWLg1xPLg2(icdi)*df_wkd(i)*df_wkp(i-1) +
! cDWLg1xPLg3(icdi)*df_wkd(i)*df_wkp(i-2) +
! cDWLg1xPLg4(icdi)*df_wkd(i)*df_wkp(i-3) +
! cDWLg2xPLg1(icdi)*df_wkd(i-1)*df_wkp(i) +
! cDWLg2xPLg2(icdi)*df_wkd(i-1)*df_wkp(i-1) +
! cDWLg2xPLg3(icdi)*df_wkd(i-1)*df_wkp(i-2) +
! cDWLg2xPLg4(icdi)*df_wkd(i-1)*df_wkp(i-3) +
! cDWLg3xPLg1(icdi)*df_wkd(i-2)*df_wkp(i) +
! cDWLg3xPLg2(icdi)*df_wkd(i-2)*df_wkp(i-1) +
! cDWLg3xPLg3(icdi)*df_wkd(i-2)*df_wkp(i-2) +
! cDWLg3xPLg4(icdi)*df_wkd(i-2)*df_wkp(i-3) +
! cDWLg4xPLg1(icdi)*df_wkd(i-3)*df_wkp(i) +
! cDWLg4xPLg2(icdi)*df_wkd(i-3)*df_wkp(i-1) +
! cDWLg4xPLg3(icdi)*df_wkd(i-3)*df_wkp(i-2) +
! cDWLg4xPLg4(icdi)*df_wkd(i-3)*df_wkp(i-3) )
c 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
c weekly weather output
c open(unit=3,file=ST//CD//'w_'//ayear
write (3,60)i,i+1,icd,iyear,year, wk(i+1),
! mir(i+1), cd_avg_mir(i+1,icdi), mir_cur(i+1),
! cd_2012mir(i+1,icd), cd_dylt(i+1,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
c weekly MIR file
c open(unit=4,file=ST//CD//'_'//ayear//'_mir'
write (4,60)i,i+1,icd,iyear,year, wk(i+1),
! mir(i+1), cd_avg_mir(i+1,icdi), mir_cur(i+1),
! cd_2012mir(i+1,icd), cd_dylt(i+1,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
60 format (2(1x,i2),4(1x,i4),7(1x,f6.2,1x),1x,8(1x,f6.2))
if (iyear .eq. 2017) then
c MIR plotting file
c open(unit=8,file=ST//CD//'_'//ayear//'_plot'
write (8,83)icd,year,wk(i+1),wmonth(i+1),wday(i+1),
! cd_2012mir(i+1,icd),cd_avg_mir(i+1,icdi),mir_cur(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
c if (iyear .eq. 2017) then
c 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
c 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))
c write (4,60)i,i,icd,iyear,year,wk(i),real(wmonth(i)),
c ! real(wday(i))
c 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)),
! real(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
c write (4,60)i,i,icd,iyear,year,wk(i),wmonth(i),wday(i),
c ! dum3,dum4,cd_dylt(i,icd),df_wkd(i),df_wkp(i),
c ! df_SPc_p,df_WIc_p,df_FAp_p,df_SUp_p,df_SPp_p,
c ! 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
Overview
Content Tools