C Last change: AU 29 Feb 2000 11:22 am c cfeb_gen_6.for c c ################################################################## c ################################################################## c ###### ###### c ###### ###### c ###### PROGRAM FOR TEMPERATURE AND SOLAR ###### c ###### RADIATION PARAMETEIZATION OF USCLIMATE MODEL ###### c ###### ###### c ###### ARS - 114 ###### c ###### (c) Jan. 1995 ###### c ###### ###### c ###### UNITED STATES DEPARTMENT OF AGRICULTUE ###### c ###### AGRICULTURAL RESEARCH SERVICE. ###### c ###### All rights reserved. ###### c ###### ###### c ################################################################## c ################################################################## PROGRAM GENPAR6 c c####################################################################### c c PURPOSE: c c This program generates the parameters w(i) in "*.TEMP" file c as an input for GEM6 climate model. w(i) consists of the annual c means, amplitudes of mean value of Tmax, Tmin and Rad, and their c coefficient of variation for wet/dry days. c Note: The mean and amplitude of standard deviation is not needed c for USCLIMATE model as input, they are evaluateD here anyway. c Users can modify this program's output very easily to use the std. c dev. value for their further study. c c####################################################################### c c AUTHOR: YUNYUN LU (based on old program on VAX from c CLARENCE RICHARDSON - USDA_ARS, TEMPLE, TX (1983) c Modify history: c Yunyun Lu (Jan,1995) c transfered program from machine VAX to HP c reduced the unnecessary procedure c added the full documentation ! very minimal documentation c c Ward Ballard (Sep 1997) c DOS version using Lahey F77L3 compiler c without graphics part c Just changed file names for DOS compatibility in low_gen_m_6.f c Also change sind,cosd,tand (vax trig functions for degrees c to sin,cos,tan, and change arguments to radians c c c Will Frymire (Jan 2000) c c !!! 3 !!! indicated areas that need work c CCC indicates code that can probably be removed c c Added test "MEANS" routine for checking and verifying i/o data. c c Modified code to start from Jan instead of March. Pervious c GenPar_6 output parameters did not match WGENPar & GenPar (3x3) c output. This version of 6x6 is consistant with the original c WGENPAR and GENPAR used with USCLIMATE. c c Verified the output - checked it using WGENPAR and GENPAR (3x3). c c Added comments for variables: c !!! PROBLEM area : some variable names may be inconstistant c c !!! COMMON BLOCKs need to be reworked as Includes c !!! OUTPUT comments need to be checked and corrected c !!! Temperatures: change inputs to Kelvin or Rankin then get rid of c the adding and subtracting of 100. c c !!!Comments and documentation needs to be verified and corrected!!!! c c Eventually we will be c adding a control file, default or on command line c ie: genpar_6 (control file ) or look for default control file c named genpar.inp for genpar input. c contents of genpar input will consist of c 1. Name of input weather file to read c Note that the weather file contains: c Year, month, day, daily precipitation, daily max temp, daily c Min temperature, dewpoint/relative humidity,total solar c radiation, and daily mean windrun c 2. Number of years in that file, user responsible that this c is accurate c 3. Latitude associated with the site c 4. output file names for means and such c 5. input units for precip,temperature,radiation,dew/rh,solar c radiation, and wind c 6. Output units for precip,temperature,radiation,dew/rh,solar c radiation, and wind. Choices should be the common units c of measure: inches/mm, F/C, Dew F/C: RH %/decimal radiation c Ly/day or w/m2, mph / m/sec c c####################################################################### c c Variable Declarations. c c####################################################################### c c OUTPUT: c c xdata(i), parameter array for 18 snotel stations c temperatures in F, radiation in Langleys c 1. Annual mean of tmax for dry days; c 2. Amplitude of tmax for wet or dry days; c 3. Annual mean of the coefficient of variation of tmax c for wet or dry days; c 4. Amplitude of the coefficient of variation of tmax c for wet or dry days; c 5. Annual mean of tmax for wet days; c 6. Annual mean of tmax for wet or dry days; c 7. Amplitude of tmin for wet or dry days; c 8. Annual mean of the coefficient of variation of tmin c for wet or dry days; c 9. Amplitude of the coefficient of variation of tmin c for wet or dry days; c 10. Annual mean of radiation for dry days; c 11. Amplitude of radiation for dry days; c 12. Annual mean of radiation for wet days; c 13. Amplitude of radiation for wet days; c 14. Annual mean of the coefficient of variation of c radiation for dry days; c 15. Amplitude of the coefficient of variation of c radiation for dry days; c !!! 16. Annual mean of the coefficient of variation of c radiation for dry days; c !!! 17. Amplitude of the coefficient of variation of c radiation for dry days; c c INPUT: (comment used to say M years instead of N years) c c tmax maximum temperature (daily) c tmin minimum temperature (daily) c rad radiation (daily) c rain precipitation (daily) c dew dew point temperature (daily) c wind wind speed (daily) c c####################################################################### c c Variable Declarations. for 30 years (was 56) c c####################################################################### c parameter(n=30,m=365,ni=6) ! 30 years data ?? If through Feb?? common /jj/jct common /const/pi common /mm1/rr0(ni,ni),rr1(ni,ni) common /mm2/am(ni,ni),bm(ni,ni) common /mm3/nyrs integer jct,nyrs integer mon(13) ccc End of MONTH 03,04,05, 06, 07, 08, 09, 10, 11, 12, 01, 02 March--Feb ccc data mon/0,31,61,92,122,153,184,214,245,275,306,337,365/ !start March ccc data mon is julian style but starts from March 1st c FOR: normal start in January data mon/0,31,59,90,120,151,181,212,243,273,304,334,365/ !Jan to Jan real pi,angle dimension tmax(n,m),tmin(n,m),rad (n,m),rain(n,m) ! rain=ppt+snow dimension dew(n,m),wind(n,m) dimension rr3(ni,ni),rr4(ni,ni) real rr0,rr1,am,bm dimension xdata(40),rc(m) character*40 header ! used to read header for boise61a.txt c character name*7,filename*60 c angle is a temporary variable to convert degrees to radians character*25 filename,temp,max,infile ! to build output filenames c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c pi=4.*atan(1.) ! give a constant value jct=0 nyrs=n C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C write (6,*)'Enter Location file name (NO extension-no blanks) ' read(5,*) filename ilen = index(filename,' ') ! string position of 1st ' ' c infile=filename(1:ilen-1)// '.txt' ! example boise.txt temp=filename(1:ilen-1)// '.temp' ! annual means boise.temp c open(unit=12,file='boise61a.txt', status='old') C open(unit=8,file=infile, status='unknown') !example: boise.txt C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C open(unit=8,file='boise61a.txt',status='old') ! Boise air port c alat=43.+43./60. !boise bbb.out (boisewx9.dat - english units) c alat=43.+4./60. !Reynolds mountain (176.out - mixed english/metric) c alat=43.+7./60. !lower sheep (127.out - mixed english/metric) c Temperatures at Reynolods Creek watershed are sent back as c degrees C. C======================================================================= c The 1st line of the record should contain the # of year in the record c and the latitude as a decimal valure. c c The 2nd line contains text information for each field. c====================================================================== c BOISE data set boise61a.txt: c VARIABLE MEANING UNITS c v1 - !tmax deg F c v2 - !tmin deg F c v3 - !dewpoint deg F c v4 - !solar langleys (watts/m^2 may range to 1500) c v5 - !wind meter/second c v6 - !precipitation inches (example 0.01) c======================================================================= read(8,*) numyears,alat ! decimal value for alat c read(8,*) header ! header for boise61a.txt !skip the header c ! for data input write(6,*) header c======================================================================== c####################################################################### c Notice: following loop reads 2 month data in order to let proc. start c from Mar. 1 c If your data set already start from Mar. 1, jump to next step. c This simplifys leap year problems - the 366th day is skipped c which should be Feb 29. c####################################################################### c CCCCC DO 20 iuu=1,59 CCCCC 21 read(8,*,end=440)iyr,imo,ida,v6,v1,v2,v3,v4,v5 ! for boise61a.txt CCC 21 read(8,*,end=440)iyr,imo,ida,v6,v1,v2,v3,v4,v5 !bbb.out CCC v6=v6/100. !bbb.out uses 100th of inches CCCCC IF(imo.eq.2.and.ida.eq.29) GOTO 21 ! If Feb 29 skip it too. CCCCC 20 CONTINUE c####################################################################### c Following loop gives a constrain for extreme value of radiation c for each day. c REF "Physical Climatology, 65-24983, 1965 c William D. Sellers, The University of Chicago Press, p 16 eq 3.7 c summary: page 232 c Where do the following constants used below come from?: c 80.25 ! 6:00 A.M. March 21 ( equinox) Julian date c 0.4102 ! c 88.2 c .0172 -- .017204 is 2pi/365.25 ie # of radian earth moves in a day c 1 c 0.0335 c 889.2305 (I think this is a constant maximum solar radiation -WARD) c (1440/pi)*1.94 && 1.94 is the solar constant in langleys. c 0.90 c####################################################################### DO 40 i=1,m ! m = 365 days / year angle= (float(i)-80.25) * (pi/180.0) ! degrees to radians !! WGENPar had 2pi/365.25 sd=0.4102*sin(angle) ! sin in radians ? .4102? angle= alat * (pi/180.0) ! convert degrees to radians ch=-tan(angle)*tan(sd) IF(ch.gt.1.) THEN h=0. goto 30 ELSE IF(ch.lt.-1.) THEN h=pi goto 30 ELSE h=acos(ch) !acos is Arccosine - old arcos(x) END IF END IF angle= (float(i) + 88.2) * (pi/180.0) ! wgenpar 2pi/365.25 30 dd=1.+0.0335*sin(angle) ! ? .0335 ? angle= alat * (pi/180.0) ! convert degrees to radians rc(i)=889.2305*dd*((h*sin(angle)*sin(sd)) : +(cos(angle)*cos(sd)*sin(h))) ! ? 889.2305 ? rc(i)=rc(i)*0.90 ! ? .90 ? Rich. 84 has .80 40 CONTINUE c####################################################################### c start to read data c####################################################################### DO 70 i=1,nyrs ! was march to march (nyrs-1) DO 77 j=1,m ! Checks for missing of out of range data v1=-999. !tmax ! needs to be changed from abs(-999) to v2=-999. !tmin ! checking the -999 itself. Reasoning: v3=-999. !dewpoint ! rad in watts/m^2 may range to 1500. v4=-999. !solar v5=-999. !wind v6=-999. !precipitation tmax(i,j)=-999. tmin(i,j)=-999. rain(i,j)=-999. rad(i,j) =-999. ! max value (watts/meter^2) about 1500 dew(i,j) =-999. wind(i,j)=-999. c 60 read(8,*,end=440)imo,ida,iyr,v6,v1,v2,v3,v4,v5 ! boise61a.txt c CCC 60 read(8,*,end=440)iyr,imo,ida,v6,v1,v2,v3,v4,v5 !bbb.out CCC v6=v6/100. !bbb.out uses 100th of inches ccccc write(6,*)'',iyr,imo,ida,v1,v2,v3,v4,v5,v6 c ! read in data day by day IF(imo.eq.2.and.ida.eq.29) GOTO 60 ! skip record for leap year dew(i,j)=v3 IF(abs(v3).lt.999..and.abs(v1).lt. @ 999..and.abs(v2).lt.999.) THEN ! if missing data (-999 marker) dew(i,j)=dew(i,j)+100. ! unit change and add 100 F if data is celcus (reynolds sites) c dew(i,j)=(dew(i,j)*9./5.)+132. ! unit change and add 100 F ?? c END IF c wind(i,j)=v5 c tmax(i,j)=v1 IF(abs(v1).lt.999.) THEN ccc tmax(i,j)=(tmax(i,j)*9./5.)+132. ! unit change and add 100 F (CELCUS???) tmax(i,j)=tmax(i,j)+100. ! add 100 F END IF c tmin(i,j)=v2 IF(abs(v2).lt.999.) THEN ccc tmin(i,j)=(tmin(i,j)*9./5.)+132. ! unit change and add 100 F tmin(i,j)=tmin(i,j)+100. ! add 100 F END IF c rain(i,j)=v6 c c below line converted ly/day to w/m2. ! max (solar constant) in watts/meter^2 1368 ?? c if(abs(v4).lt.999.) v4=v4*0.4847 ! max langely/day = 2822 rad(i,j)=min(v4,rc(j)) c write(6,'(3i5,4f7.1)'),iyr,imo,ida,tmax(i,j),tmin(i,j),dew(i,j) c !,rain(i,j) c 77 CONTINUE c write(6,'(10f8.3)')(rain(i,j),j=1,m) ! NEEDS to BE WRITTEN TO A FILE 70 CONTINUE close(8) ! CLOSE INPUT FILE DO 50 i=1,40 ! initilize xdata array to 0 xdata(i)=0. 50 CONTINUE open(unit=21,file='s1',status='replace') ! for creating MO matrix ? open(unit=22,file='s2',status='replace') ! for creating M1 matrix ? open(unit=23,file='s3',status='replace') ! for creating A matrix ? open(unit=24,file='s4',status='replace') ! for creating B matrix ? write(21,'(a10)')'M0 matrix:' write(22,'(a10)')'M1 matrix:' write(23,'(a10)')'A matrix:' write(24,'(a10)')'B matrix:' do i=1,6 ! initillize rr3 & rr4 do j=1,6 rr3(i,j)=0. rr4(i,j)=0. end do end do do icc=1,12 c CCCCC igg=icc+2 ! offset to start in March CCCCC if(icc.gt.10) igg=icc-10 CCCCC March offset removed to allow use of all years of data. CCCCC The March offset does not use the last 10 months of input data. CCCCC THUS the user must arrange the A & B matrixs for march start for CCCCC for use in GEM6. (*.max) c igg=icc write(21,'(a10,i5)')'month: ',igg write(22,'(a10,i5)')'month: ',igg write(23,'(a10,i5)')'month: ',igg write(24,'(a10,i5)')'month: ',igg nj1=mon(icc)+1 ! from March julian style date: first day of month icc nj2=mon(icc+1) ! from March julian style date: last day of month icc c nj1=1 c nj2=365 CALL M01_MATRIX(tmax,tmin,dew,rad,wind,nj1,nj2) write(21,*) write(21,'(6f7.3)') ((rr0(i,j),j=1,6),i=1,6) write(22,*) write(22,'(6f7.3)') ((rr1(i,j),j=1,6),i=1,6) write(23,*) write(23,'(6f7.3)') ((am(i,j),j=1,6),i=1,6) write(24,*) write(24,'(6f7.3)') ((bm(i,j),j=1,6),i=1,6) c do i=1,6 do j=1,6 rr3(i,j)=rr3(i,j)+rr0(i,j) rr4(i,j)=rr4(i,j)+rr1(i,j) end do end do end do write(21,'(a20)')'annual mean: ' write(22,'(a20)')'annual mean: ' write(23,'(a20)')'annual mean: ' write(24,'(a20)')'annual mean: ' do i=1,6 do j=1,6 rr0(i,j)=rr3(i,j)/12. rr1(i,j)=rr4(i,j)/12. end do end do CALL ABMATRIX write(21,*) write(21,'(6f7.3)') ((rr0(i,j),j=1,6),i=1,6) write(22,*) write(22,'(6f7.3)') ((rr1(i,j),j=1,6),i=1,6) write(23,*) write(23,'(6f7.3)') ((am(i,j),j=1,6),i=1,6) write(24,*) write(24,'(6f7.3)') ((bm(i,j),j=1,6),i=1,6) CALL MSD(tmax,rain,xdata) ! for statistic cal. c CALL MSD1(tmax,rain,xdata) ! for statistic cal. CALL MSD(tmin,rain,xdata) CALL MSD(dew,rain,xdata) CALL MSD(rad,rain,xdata) CALL MSD(wind,rain,xdata) c####################################################################### c For final output phase angles shifted by 180 degrees and c signs changed on amplitudes c####################################################################### c####################################################################### c Temperatures subtracted by 100 F. c####################################################################### do ic=1,10 io=(ic-1)*4+1 if(ic.le.6) then xdata(io)=xdata(io)-100. end if xdata(io+1)=-xdata(io+1) if(xdata(io+3).gt.0.)xdata(io+3)=-xdata(io+3) end do cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(6,*) 'ANNUAL MEANS at end of program ' call MEANS (tmax,'FOR TMAX ',1) call MEANS (tmin,'FOR TMIN ',1) call MEANS (rain,'For RAIN ',-1) call MEANS (dew,'For dew ',1) call MEANS (rad,'For RAD ',3) call MEANS (wind,'For Wind ',-1) write (6,*)' ' ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c print*,(xdata(ip),ip=1,40) C write(6,'(8f10.5)')(xdata(ip),ip=1,40) OPEN(UNIT=19,FILE=temp,STATUS='replace') write(19,'(8f10.5)')(xdata(ip),ip=1,40) c write(lunw,1000)(xdata(ic),ic=1,5), c !(xdata(ic),ic=9,14),(xdata(ic),ic=17,18) 1000 format(2i4,2f7.3,3i4,2f7.3,4i4) GOTO 450 440 write(6,*)'eof' 450 CONTINUE c close (21) close (22) close (23) close (24) close (16) c c Write the GEM6 input matrix file starting with march instead of c January. call GEMatrix(filename) c STOP END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE M01_MATRIX ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE M01_MATRIX(tmax,tmin,dew,rad,wind,nj1,nj2) parameter(n=30,m=365,ni=6,nj=36) !n was 56 years dimension tmax(n,m),tmin(n,m),rad(n,m) dimension dew(n,m),wind(n,m) common /mm1/rr0(ni,ni),rr1(ni,ni) common /mm2/am(ni,ni),bm(ni,ni) common /mm3/nyrs do ijk=1,ni rr0(ijk,ijk)=1 end do CALL PEARSN1(tmax,tmax,r,r1,r2,r3,r4,nj1,nj2) ! tmax tmax ??? rr0(1,2)=r rr0(2,1)=r rr1(1,1)=r1 rr1(2,2)=r2 rr1(1,2)=r3 rr1(2,1)=r4 CALL PEARSN1(tmax,tmin,r,r1,r2,r3,r4,nj1,nj2) rr0(1,3)=r rr0(3,1)=r rr1(3,3)=r2 rr1(1,3)=r3 rr1(3,1)=r4 CALL PEARSN1(tmax,rad,r,r1,r2,r3,r4,nj1,nj2) rr0(1,4)=r rr0(4,1)=r rr1(4,4)=r2 rr1(1,4)=r3 rr1(4,1)=r4 CALL PEARSN1(tmax,dew,r,r1,r2,r3,r4,nj1,nj2) rr0(1,5)=r rr0(5,1)=r rr1(5,5)=r2 rr1(1,5)=r3 rr1(5,1)=r4 CALL PEARSN1(tmax,wind,r,r1,r2,r3,r4,nj1,nj2) rr0(1,6)=r rr0(6,1)=r rr1(6,6)=r2 rr1(1,6)=r3 rr1(6,1)=r4 CALL PEARSN(tmax,tmin,r,r1,r2,r3,r4,nj1,nj2) rr0(2,3)=r rr0(3,2)=r rr1(2,3)=r3 rr1(3,2)=r4 CALL PEARSN(tmax,rad,r,r1,r2,r3,r4,nj1,nj2) rr0(2,4)=r rr0(4,2)=r rr1(2,4)=r3 rr1(4,2)=r4 CALL PEARSN(tmax,dew,r,r1,r2,r3,r4,nj1,nj2) rr0(2,5)=r rr0(5,2)=r rr1(2,5)=r3 rr1(5,2)=r4 CALL PEARSN(tmax,wind,r,r1,r2,r3,r4,nj1,nj2) rr0(2,6)=r rr0(6,2)=r rr1(2,6)=r3 rr1(6,2)=r4 CALL PEARSN(tmin,rad,r,r1,r2,r3,r4,nj1,nj2) rr0(3,4)=r rr0(4,3)=r rr1(3,4)=r3 rr1(4,3)=r4 CALL PEARSN(tmin,dew,r,r1,r2,r3,r4,nj1,nj2) rr0(3,5)=r rr0(5,3)=r rr1(3,5)=r3 rr1(5,3)=r4 CALL PEARSN(tmin,wind,r,r1,r2,r3,r4,nj1,nj2) rr0(3,6)=r rr0(6,3)=r rr1(3,6)=r3 rr1(6,3)=r4 CALL PEARSN(rad,dew,r,r1,r2,r3,r4,nj1,nj2) rr0(4,5)=r rr0(5,4)=r rr1(4,5)=r3 rr1(5,4)=r4 CALL PEARSN(rad,wind,r,r1,r2,r3,r4,nj1,nj2) rr0(4,6)=r rr0(6,4)=r rr1(4,6)=r3 rr1(6,4)=r4 CALL PEARSN(dew,wind,r,r1,r2,r3,r4,nj1,nj2) rr0(5,6)=r rr0(6,5)=r rr1(5,6)=r3 rr1(6,5)=r4 CALL ABMATRIX RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE PEARSN ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE PEARSN(xx,yy,r,r1,r2,r3,r4,nj1,nj2) common /mm3/nyrs dimension xx(30,365),yy(30,365) ! 30 was 56 years CCCCC real x(10585),y(10585) real x(10950),y(10950) n=0 do i=1,nyrs c do j=1,365 do j=nj1,nj2 jj=j+1 ii=i if(jj.gt.365) then ii=i+1 jj=1 end if if(ii.le.nyrs) then if(abs(xx(i,j)).lt.999..and.abs(yy(i,j)). ! lt.999..and.abs(xx(ii,jj)).lt.999.. ! and.abs(yy(ii,jj)).lt.999.) then n=n+1 x(n)=xx(i,j) ! In Pearsn1 x(n)=xx(ii,jj)-xx(i,j) y(n)=yy(i,j) end if end if end do end do ! at end of loop n = nyrs(nj2-nj1) r=0. r1=0. r2=0. r3=0. r4=0. ax=0. ax1=0. ay=0. ay1=0. ! NOTE n = nyrs(nj2-nj1) do i=2,n ! Get the sums for computing the means ax=ax+x(i) ! WHY SKIP x(1) ??? etc. ay=ay+y(i) ax1=ax1+x(i-1) ay1=ay1+y(i-1) end do ax=ax/float(n-1) ! determine the means ax1=ax1/float(n-1) ay=ay/float(n-1) ay1=ay1/float(n-1) ! end of the means determination sxx=0. syy=0. sxy=0. sx1x1=0. sy1y1=0. sxx1=0. syy1=0. sxy1=0. syx1=0. ! NOTE n = nyrs(nj2-nj1) do i=2,n ! Compute the correlation coefficient xt=x(i)-ax yt=y(i)-ay xt1=x(i-1)-ax1 yt1=y(i-1)-ay1 sxx=sxx+xt**2 syy=syy+yt**2 sxy=sxy+xt*yt sx1x1=sx1x1+xt1**2 sy1y1=sy1y1+yt1**2 sxx1=sxx1+xt*xt1 syy1=syy1+yt*yt1 sxy1=sxy1+xt*yt1 syx1=syx1+yt*xt1 end do r=sxy/sqrt(sxx*syy) ! r values to be returned to calling routine r1=sxx1/sqrt(sxx*sx1x1) r2=syy1/sqrt(syy*sy1y1) r3=sxy1/sqrt(sxx*sy1y1) r4=syx1/sqrt(syy*sx1x1) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE PEARSN1 ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## c c Only one line is different from routine PEARSN (xx(ii,jj)-(xx(i,j))) c SUBROUTINE PEARSN1(xx,yy,r,r1,r2,r3,r4,nj1,nj2) common /mm3/nyrs dimension xx(30,365),yy(30,365) ! 30 was 56 years CCCCC real x(10585),y(10585) real x(10950),y(10950) c n=0 do i=1,nyrs c do j=1,365 do j=nj1,nj2 jj=j+1 ii=i if(jj.gt.365) then ii=i+1 jj=1 end if if(ii.le.nyrs) then if(abs(xx(i,j)).lt.999..and.abs(yy(i,j)). ! lt.999..and.abs(xx(ii,jj)).lt.999.. ! and.abs(yy(ii,jj)).lt.999.) then n=n+1 x(n)=xx(ii,jj)-xx(i,j) ! PEARSN does not have differences ! ?? goodness of fit ??? y(n)=yy(i,j) end if end if end do end do r =0. r1 =0. r2 =0. r3 =0. r4 =0. ax =0. ax1=0. ay =0. ay1=0. do i=2,n ax =ax+x(i) ay =ay+y(i) ax1=ax1+x(i-1) ay1=ay1+y(i-1) end do ax =ax/float(n-1) ax1=ax1/float(n-1) ay =ay/float(n-1) ay1=ay1/float(n-1) sxx=0. syy=0. sxy=0. sx1x1=0. sy1y1=0. sxx1=0. syy1=0. sxy1=0. syx1=0. do i=2,n xt=x(i)-ax yt=y(i)-ay xt1=x(i-1)-ax1 yt1=y(i-1)-ay1 sxx=sxx+xt**2 syy=syy+yt**2 sxy=sxy+xt*yt sx1x1=sx1x1+xt1**2 sy1y1=sy1y1+yt1**2 sxx1=sxx1+xt*xt1 syy1=syy1+yt*yt1 sxy1=sxy1+xt*yt1 syx1=syx1+yt*xt1 end do r=sxy/sqrt(sxx*syy) r1=sxx1/sqrt(sxx*sx1x1) r2=syy1/sqrt(syy*sy1y1) r3=sxy1/sqrt(sxx*sy1y1) r4=syx1/sqrt(syy*sx1x1) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE ABMATRIX ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE ABMATRIX c####################################################################### c c PURPOSE: c c From M0 and M1 matrix to got A and B matrix. c c####################################################################### c c AUTHOR: Yunyun Lu c c####################################################################### c c Variable Declarations. c c####################################################################### c integer nmax,n PARAMETER (NMAX=60, n=6) dimension ipiv(nmax),indxr(nmax),indxc(nmax) common /mm1/rr0(n,n),rr1(n,n) common /mm2/am(n,n),bm(n,n) dimension a(n,n),c(n,n),d(n,n) c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% do 750 i1=1,n ipiv(i1)=0 do 750 i2=1,n a(i1,i2)=rr0(i1,i2) 750 continue DO 40 i=1,n big=0. DO 50 j=1,n IF(ipiv(j).ne.1) THEN DO 60 k=1,n IF (ipiv(k).eq.0) THEN IF (abs(a(j,k)).ge.big) THEN big=abs(a(j,k)) irow=j icol=k ENDIF ELSE IF (ipiv(k).gt.1) THEN pause 'Singular matrix' END IF 60 CONTINUE END IF 50 CONTINUE ipiv(icol)=ipiv(icol)+1 IF(irow.ne.icol) THEN DO 70 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 70 CONTINUE END IF indxr(i)=irow indxc(i)=icol IF (a(icol,icol).eq.0.) pause 'Singular matrix.' pivinv=1./a(icol,icol) a(icol,icol)=1. DO 80 l=1,n a(icol,l)=a(icol,l)*pivinv 80 CONTINUE DO 90 ll=1,n IF(ll.ne.icol) THEN dum=a(ll,icol) a(ll,icol)=0. DO 100 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 100 CONTINUE END IF 90 CONTINUE 40 CONTINUE DO 110 l=n,1,-1 IF(indxr(l).ne.indxc(l)) THEN DO 120 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 120 CONTINUE END IF 110 CONTINUE DO 130 i=1,n DO 130 j=1,n am(i,j)=0. bm(i,j)=0. c(j,i)=0. DO 140 l=1,n am(i,j)=am(i,j)+rr1(i,l)*a(l,j) 140 CONTINUE c(j,i)=am(i,j) 130 CONTINUE DO 150 i=1,n DO 150 j=1,n d(i,j)=0. DO 160 l=1,n d(i,j)=d(i,j)+rr1(i,l)*c(l,j) 160 CONTINUE d(i,j)=rr0(i,j)-d(i,j) 150 CONTINUE bm(1,1)=sqrt(d(1,1)) bm(1,2)=0. bm(1,3)=0. bm(1,4)=0. bm(1,5)=0. bm(1,6)=0. bm(2,1)=d(1,2)/bm(1,1) bm(2,2)=sqrt(d(2,2)-bm(2,1)*bm(2,1)) bm(2,3)=0. bm(2,4)=0. bm(2,5)=0. bm(2,6)=0. bm(3,1)=d(1,3)/bm(1,1) bm(3,2)=(d(3,2)-bm(3,1)*bm(2,1))/bm(2,2) bm(3,3)=sqrt(d(3,3)-bm(3,2)*bm(3,2) ! -bm(3,1)*bm(3,1)) bm(3,4)=0. bm(3,5)=0. bm(3,6)=0. bm(4,1)=d(1,4)/bm(1,1) bm(4,2)=(d(4,2)-bm(4,1)*bm(2,1))/bm(2,2) bm(4,3)=(d(4,3)-bm(4,1)*bm(3,1) ! -bm(4,2)*bm(3,2) ! )/bm(3,3) bm(4,4)=sqrt(d(4,4)-bm(4,3)*bm(4,3) ! -bm(4,2)*bm(4,2) ! -bm(4,1)*bm(4,1)) bm(4,5)=0. bm(4,6)=0. bm(5,1)=d(1,5)/bm(1,1) bm(5,2)=(d(5,2)-bm(5,1)*bm(2,1))/bm(2,2) bm(5,3)=(d(5,3)-bm(5,1)*bm(3,1) ! -bm(5,2)*bm(3,2) ! )/bm(3,3) bm(5,4)=(d(5,4)-bm(5,1)*bm(4,1) ! -bm(5,2)*bm(4,2) ! -bm(5,3)*bm(4,3) ! )/bm(4,4) bm(5,5)=sqrt(d(5,5)-bm(5,4)*bm(5,4) ! -bm(5,3)*bm(5,3) ! -bm(5,2)*bm(5,2) ! -bm(5,1)*bm(5,1)) bm(5,6)=0. bm(6,1)=d(1,6)/bm(1,1) bm(6,2)=(d(6,2)-bm(6,1)*bm(2,1))/bm(2,2) bm(6,3)=(d(6,3)-bm(6,1)*bm(3,1) ! -bm(6,2)*bm(3,2) ! )/bm(3,3) bm(6,4)=(d(6,4)-bm(6,1)*bm(4,1) ! -bm(6,2)*bm(4,2) ! -bm(6,3)*bm(4,3) ! )/bm(4,4) bm(6,5)=(d(6,5)-bm(6,1)*bm(5,1) ! -bm(6,2)*bm(5,2) ! -bm(6,3)*bm(5,3) ! -bm(6,4)*bm(5,4) ! )/bm(5,5) bm(6,6)=sqrt(d(6,6)-bm(6,5)*bm(6,5) ! -bm(6,4)*bm(6,4) ! -bm(6,3)*bm(6,3) ! -bm(6,2)*bm(6,2) ! -bm(6,1)*bm(6,1)) RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE FOUR ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE FOUR(xm,sd,cv,xdata) c c####################################################################### c c PURPOSE: c c CALCULATE the mean and amplitude for each inputed variable. c c####################################################################### c c AUTHOR: CLARENCE RICHARDSON - USDA_ARS, TEMPLE, TX (1983) c c Modify history: c Yunyun Lu (Jan,1995) c transfered program from machine VAX to HP c reduced the unnecessary procedure c added the full documentation c c####################################################################### c c INPUT: c c xm mean value for 26 periods c sd standard deviation for 26 period c cv coefficient of variation c c OUTPUT: c c xdata parameters c c####################################################################### c c Variable Declarations. c c####################################################################### c common /jj/jct common /const/pi parameter (m=26) dimension xm(m),sd(m),cv(m),xdata(40) c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c####################################################################### c Calculate the mean of mean (s), std.dev. (s1) and coefficient c of variation (s2). c####################################################################### s=0. s1=0. s2=0. DO 10 i=1,m s=s+xm(i) s1=s1+sd(i) s2=s2+cv(i) 10 CONTINUE s=s/float(m) s1=s1/float(m) s2=s2/float(m) c####################################################################### c Calculate the amplitude of mean (c), std.dev. (c1) and coefficient c of variation (c2). c####################################################################### sa=0. sa1=0. sa2=0. sb=0. sb1=0. sb2=0. DO 20 k=1,m ! fourier analysis sa = sa +(xm(k)-s)*cos(2*pi*float(K)/float(m)) sa1 = sa1 +(sd(k)-s1)*cos(2*pi*float(K)/float(m)) sa2 = sa2 +(cv(k)-s2)*cos(2*pi*float(K)/float(m)) sb = sb +(xm(k)-s)*sin(2*pi*float(K)/float(m)) sb1 = sb1 +(sd(k)-s1)*sin(2*pi*float(K)/float(m)) sb2 = sb2 +(cv(k)-s2)*sin(2*pi*float(K)/float(m)) 20 CONTINUE a=sa*2/float(m) ! coeff. of cos term a1=sa1*2/float(m) a2=sa2*2/float(m) b=sb*2/float(m) ! coeff. of sin term b1=sb1*2/float(m) b2=sb2*2/float(m) t=atan(-b/a) ! phase angle t1=atan(-b1/a1) t2=atan(-b2/a2) c=a/cos(t) ! amplitude c1=a1/cos(t1) c2=a2/cos(t2) jct=jct+1 xdata(jct)=s jct=jct+1 xdata(jct)=c jct=jct+1 xdata(jct)=s2 jct=jct+1 xdata(jct)=c2 RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE MSD ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE MSD(w,rain,xdata) c c####################################################################### c c PURPOSE: c c CALCULATE the mean, standard deviation and coefficient variation c value for each 26 days period from 30 years daily data. c Separate the wet and dry days for Tmax and Rad calculation. c c####################################################################### c c AUTHOR: CLARENCE RICHARDSON - USDA_ARS, TEMPLE, TX (1983) c c Modify history: c Yunyun Lu (Jan,1995) c transfered program from machine VAX to HP c reduced the unnecessary procedure c added the full documentation c cc####################################################################### c c INPUT: c c w 30 years daily data. c rain 30 years daily precipitation. c c OUTPUT: c c xm mean value for 26 periods c sd starnard deviation for 26 period c cv coefficient of variation for 26 periods c c####################################################################### c c Variable Declarations. c c####################################################################### c common /jj/jct common /const/ pi common /mm3/nyrs parameter (m=26,n=30) ! 30 was 56 years dimension w(n,365), rain(n,365) dimension xm(m),xm1(m),sd(m),sd1(m),cx(m),cx1(m) dimension xdata(40) c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sumxx=0 sumx1=0 xcount=0 x2count=0 DO 60 i=1,26 nf=i*14 ni=nf-13 xn=0. xn1=0. sum=0. sum1=0. ss=0. ss1=0. DO 30 jd=ni,nf DO 30 jy=1,nyrs IF(abs(rain(jy,jd)).ge.999.) GOTO 30 ! not count for missing prec. IF(rain(jy,jd))10,10,20 10 CONTINUE ! for Tmax and Rad count dry day. ! for Tmin count for all. IF(abs(w(jy,jd)).ge.999.) GOTO 30 ! not count for missing data xn=xn+1. sum=sum+w(jy,jd) ss=ss+(w(jy,jd)*w(jy,jd)) GOTO 30 20 CONTINUE ! for Tmax and Rad count wet day. IF(abs(w(jy,jd)).ge.999.) GOTO 30 ! not count for missing data xn1=xn1+1. sum1=sum1+w(jy,jd) ss1=ss1+(w(jy,jd)*w(jy,jd)) 30 CONTINUE xm(i)=0. sd(i)=0. cx(i)=0. xm(i)=sum/xn ! dry day's or {Tmin (all)'s} sd(i)=sqrt((ss-sum*sum/xn)/(xn-1.)) ! mean, std. dev. and c.v. cx(i)=sd(i)/xm(i) 40 CONTINUE xm1(i)=0. sd1(i)=0. cx1(i)=0. xm1(i)=sum1/xn1 sd1(i)=sqrt((ss1-sum1*sum1/xn1)/(xn1-1.)) cx1(i)=sd1(i)/xm1(i) 60 CONTINUE c c c sumxx = sumxx + sum sumx1 = sumx1 + sum1 xcount = xcount + xn x2count=x2count + xn1 c c xmm=0 xmm1=0 do kk=1,26 xmm=xmm+xm(kk) xmm1=xmm1+xm1(kk) end do c CCC write(6,*) 'xn=',xn,'sums XN1=',XN1,'xm=',xmm,'XM1=',xmm1 xdsum=0 xwsum=0 ! running sum do jj=1,26 xdsum=xdsum+xm(jj) xwsum=xwsum+xm1(jj) xdry=xdsum/jj xwet=xwsum/jj end do ccc write(6,*) 'I=',jj,' XM=',xm(jj),' XM1=',xm1(jj),xdry,' ',xwet ccc write(6,*)'sum=',sumxx,'sum1=',SUMx1,'xc=',xcount,'xc2=',x2count c c c####################################################################### c Call subroutine FOUR to calculate the mean and amplitude of c 1) xm, sd, cv for dry day of Tmax and Rad, dry/wet day for Tmin; c 2) xm1, sd1, cv1 for wet day of Tmax, Rad. c####################################################################### CALL FOUR(xm,sd,cx,xdata) CALL FOUR(xm1,sd1,cx1,xdata) 70 CONTINUE RETURN END c ################################################################## c ################################################################## c ###### ###### c ###### SUBROUTINE MSD1 ###### c ###### ###### c ###### (c) 1995 ###### c ###### ###### c ################################################################## c ################################################################## SUBROUTINE MSD1(w,rain,xdata) c c####################################################################### c c PURPOSE: c c CALCULATE the mean, standard deviation and coefficient variation c value for each 26 days period from 30 years daily data. c Separate the wet and dry days for Tmax and Rad calculation. c c####################################################################### c c AUTHOR: CLARENCE RICHARDSON - USDA_ARS, TEMPLE, TX (1983) c c Modify history: c Yunyun Lu (Jan,1995) c transfered program from machine VAX to HP c reduced the unnecessary procedure c added the full documentation c cc####################################################################### c c INPUT: c c w 30 years daily data. c rain 30 years daily precipitation. c c OUTPUT: c c xm mean value for 26 periods c sd starnard deviation for 26 period c cv coefficient of variation for 26 periods c c####################################################################### c c Variable Declarations. c c####################################################################### c common /jj/jct common /const/ pi common /mm3/nyrs parameter (m=26,n=30) ! 30 was 56 years dimension w(n,365), rain(n,365) dimension xm(m),xm1(m),sd(m),sd1(m),cx(m),cx1(m) dimension xdata(40) c####################################################################### c c C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Beginning of executable code... C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO 60 i=1,26 nf=i*14 ni=nf-13 xn=0. xn1=0. sum=0. sum1=0. ss=0. ss1=0. DO 30 jd=ni,nf DO 30 jy=1,nyrs IF(abs(rain(jy,jd)).ge.999.) GOTO 30 ! not count for missing prec. IF(rain(jy,jd))10,10,20 10 CONTINUE ! for Tmax and Rad count dry day. ! for Tmin count for all. IF(abs(w(jy,jd)).ge.999.) GOTO 30 ! not count for missing data IF(abs(w(jy,jd+1)).ge.999.) GOTO 30 ! not count for missing data xn=xn+1. sum=sum+(w(jy,jd+1)-w(jy,jd)) ss=ss+((w(jy,jd+1)-w(jy,jd))*(w(jy,jd+1)-w(jy,jd))) GOTO 30 20 CONTINUE ! for Tmax and Rad count wet day. IF(w(jy,jd).eq.-999.) GOTO 30 xn1=xn1+1. sum1=sum1+(w(jy,jd+1)-w(jy,jd)) ss1=ss1+((w(jy,jd+1)-w(jy,jd))*(w(jy,jd+1)-w(jy,jd))) 30 CONTINUE xm(i)=0. sd(i)=0. cx(i)=0. xm(i)=sum/xn ! dry day's or {Tmin (all)'s} sd(i)=sqrt((ss-sum*sum/xn)/(xn-1.)) ! mean, std. dev. and c.v. cx(i)=sd(i)/xm(i) 40 CONTINUE xm1(i)=0. sd1(i)=0. cx1(i)=0. xm1(i)=sum1/xn1 sd1(i)=sqrt((ss1-sum1*sum1/xn1)/(xn1-1.)) cx1(i)=sd1(i)/xm1(i) 60 CONTINUE c####################################################################### c Call subroutine FOUR to calculate the mean and amplitude of c 1) xm, sd, cv for dry day of Tmax and Rad, dry/wet day for Tmin; c 2) xm1, sd1, cv1 for wet day of Tmax, Rad. c####################################################################### CALL FOUR(xm,sd,cx,xdata) CALL FOUR(xm1,sd1,cx1,xdata) 70 CONTINUE RETURN END c********************************************************************** cwlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf c********************************************************************** c TEST routines : This area contains routines for data testing and c verification c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%% %%%% c%%% Subroutine: MEANS %%%% c%%% %%%% c%%% Jan 200 %%%% c%%% %%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c c PURPOSE: c This routine tests tmax array elements for the mean values passed c in each array. To test at various points in the program copy this c block to the location you wish to test, Then remove the "c" comment in c front of the statements. c c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%c ccccSTART of BLOCKccccccccccccccccccccccccccccccccccccccccccccccccccccc c write(6,*) ' Location comment ' c call MEANS (tmax,'FOR TMAX ',1) ! 1 flag for temperature c call MEANS (tmin,'FOR TMIN ',1) c call MEANS (rain,'For RAIN ',-1) ! -1 flag for precipitation c call MEANS (dew,'For dew ',1) c call MEANS (rad,'For RAD ',3) ! 3 flag for radiation c call MEANS (wind,'For Wind ',-1) cccccEND of BLOCKcccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c C AUTHOR: W.L. Frymire: Jan 2000 c********************************************************************** SUBROUTINE MEANS(xarray,String,iflag) DIMENSION xarray(30,365),yarray(30,365) INTEGER iflag CHARACTER*10 String c open(unit=16,file='MEANS_Test.txt',status='replace') open(unit=36,file='MEANS_Array.txt',status='replace') kkount = 0 total =0 do 10 k = 1,30 !years do 20 kk= 1,364 !365 c if (iflag.eq.1) then yarray(k,kk)=xarray(k,kk)-100 kkount=kkount+1 total = yarray(k,kk)+ total endif c if (iflag.eq.-1) then yarray(k,kk)=xarray(k,kk) kkount=kkount+1 total = Yarray(k,kk)+ total endif c if (iflag.eq.3) then yarray(k,kk)=xarray(k,kk) kkount=kkount+1 total = yarray(k,kk)+ total end if c if (iflag.eq.3) goto 20 c if(xarray(k,kk).gt.211.or.xarray(k,kk).lt.-25) then write(16,*)String,' k=',k,'kk=',kk,'xarray=',xarray(k,kk) end if c 20 continue 10 continue xmean=total/kkount write(6,*)String,'Count=',kkount,'Tot=',total,'Mn=',xmean,iflag write(36,*)String write(36,37) xarray 37 format(16f7.2) c return END c******************************************************************** c********************************************************************** cwlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf_wlf c********************************************************************** c******************************************************************** c C C ################################################################## C ################################################################## C ###### ###### C ###### G E M 6 T O O L S ###### C ###### ###### C ###### USE with GENPAR6 ###### C ###### ###### C ###### PROGRAM FOR DAILY WEATHER SIMULATION ###### C ###### IN THE CONTIGOUS UNITED STATES ###### C ###### ###### C ###### ARS - 114 ###### c ###### Feburay 2000 ###### C ###### ###### C ###### UNITED STATES DEPARTMENT OF AGRICULTUE ###### C ###### AGRICULTURAL RESEARCH SERVICE. ###### C ###### All rights reserved. ###### C ###### ###### C ################################################################## C ################################################################## SUBROUTINE GEMatrix(filename) c c###################################################################### c PURPOSE: c c This program takes the S3 and S4 files from GENPAR6 output c to create the *.max (A & B) matrix file for input into GEM6. c c##################################################################### c N O T E: c c GEM6 uses data starting from March 1st. c GENPAR6 uses all available data thus outputs from Jan thru Dec. c c#################################################################### c CHARACTER*25 filename,outfile CHARACTER*43 Amatrix(106) CHARACTER*43 Bmatrix(106) CHARACTER*43 GEMatrx (194) c open(unit=33,file='s3',status='unknown') ! A matrix open(unit=34,file='s4',status='unknown') ! B matrix c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c B U I L D O U T P U T F I L E N A N E c ilen = index(filename,' ') ! string position of 1st ' ' c outfile=filename(1:ilen-1)// '.max' ! example IDboise.max c open (UNIT=31,FILE=outfile,STATUS='replace') c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c do j=1,97 read (33,10)Amatrix(j) 10 format (a42) c write (6,11)Amatrix(j),'Amatrix ' 11 format (a42,1x,a8) read (34,10)Bmatrix(j) c write (6,11)Bmatrix(j),'Bmatrix ' end do c c A matrix header ( get the 1st line of the s3 file) c GEMatrx(1) = Amatrix(1) c c month: 3 thru 12 (set months 3-12 in front of 1 & 2) c do i=2,81 GEMatrx(i)=Amatrix(i+16) end do c c month 1 - 2 (place Jan & Feb after the rest of the months c do k=82,97 GEMatrx(k)=Amatrix(k-80) end do c c B matrix (concat. B matrix into same file) c GEMatrx(98) = Bmatrix(1) c c c month: 3 thru 12 c do i=99,178 GEMatrx(i)=Bmatrix(i-81) end do c c month 1 - 2 c do k=179,194 GEMatrx(k)=Bmatrix(k-177) end do c c c write the GEMartix to a *.max file do k=1,194 write (31,*)GEMatrx(k) end do c close (33) close (34) close (31) c return END