./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
http://www.angulartutorial.net/2017/09/calculate-epi-week-from-date-convert.html |
! 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 | Is this redundant opening of file? |
! 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 | INPUT FILE ILCD_divisions_ALL |
! 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 | CD_Daylight.csv |
! 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 | Is this needed to run model? |
! 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 | Needed for running model? |
! 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 | CD_Coef_2017.csv |