character*14 input character*18 output external f call getarg(1,input) open(unit=10,file=input,status='unknown') write(output,1000) input 1000 format(a14,'_adj') open(unit=20,file=output,status='unknown') c******************************************************** c Obtain buoy height c******************************************************** if(input.eq.'mdrm1h1997.txt') z=22.6 if(input.eq.'mism1h1997.txt') z=16.5 if(input.eq.'44007h1997.txt') z=5.0 if(input.eq.'iosn3h1997.txt') z=19.2 if(input.eq.'44005h1997.txt') z=5.0 if(input.eq.'44013h1997.txt') z=5.0 if(input.eq.'44028h1997.txt') z=13.8 if(input.eq.'44011h1997.txt') z=5.0 if(input.eq.'buzm3h1997.txt') z=24.8 if(input.eq.'alsn6h1997.txt') z=49.1 if(input.eq.'44025h1997.txt') z=5.0 if(input.eq.'44008h1997.txt') z=5.0 if(input.eq.'tplm2h1997.txt') z=18.0 if(input.eq.'44009h1997.txt') z=5.0 if(input.eq.'44004h1997.txt') z=5.0 if(input.eq.'chlv2h1997.txt') z=43.3 if(input.eq.'44014h1997.txt') z=5.0 if(input.eq.'ducn7h1997.txt') z=20.4 if(input.eq.'dsln7h1997.txt') z=46.6 if(input.eq.'clkn7h1997.txt') z=9.8 if(input.eq.'41001h1997.txt') z=5.0 if(input.eq.'fpsn7h1997.txt') z=44.2 if(input.eq.'41004h1997.txt') z=5.0 if(input.eq.'fbis1h1997.txt') z=9.8 if(input.eq.'41002h1997.txt') z=5.0 if(input.eq.'41008h1997.txt') z=5.0 if(input.eq.'sauf1h1997.txt') z=16.5 if(input.eq.'41010h1997.txt') z=5.0 if(input.eq.'41009h1997.txt') z=5.0 if(input.eq.'lkwf1h1997.txt') z=13.7 if(input.eq.'spgf1h1997.txt') z=9.8 if(input.eq.'fwyf1h1997.txt') z=43.9 if(input.eq.'mlrf1h1997.txt') z=15.8 c******************************************************** c Define constants c******************************************************** cp=1004 rho=1.275 cH=0.003 von=0.4 g=9.8 pi=3.14159 zten=10.0 alpha=0.015 write(20,1500) 1500 format('YY MM DD hh WD WSPD GST WVHT DPD APD MWD BAR AT &MP WTMP DEWP VIS') c print*,'z=',z c******************************************************* c Read in buoy data c******************************************************* read(10,*) do 5000 i=1,10000 read(10,*,end=500) yy,mm,dd,hh,wd,wspd,gst,wvht,dpd, & apd,mwd,bar,atmp,wtmp,dewp,vis c******************************************************** c If data is not in March, April, May, or June, skip c******************************************************** if(mm.lt.3.or.mm.gt.6) goto 5000 c******************************************************** c If wind not available, skip whole normalization process c and go to line 50 c******************************************************** if(wspd.gt.98.0) goto 50 c******************************************************** c If wind is less than 0.5 m/s, don't convert to 10 meters c******************************************************** if(wspd.lt.0.5) adjwspd=wspd if(wspd.lt.0.5) goto 35 c******************************************************** c Iterate for ustar c******************************************************** xleft=0.000001 xright=5.0 call bisect(f,xleft,xright,ustar,z,wspd) c******************************************************** c Compute roughness length zo using Charnock's equation c******************************************************** zo=alpha*ustar*ustar/g c******************************************************* c Adjust wind speed to 10 meters height c******************************************************* adjwspd=ustar/von*(log(zten/zo)) 35 continue c print*,'for nonneutral old=',wspd,' m/s and new=',adjwspd c******************************************************* c Increase 8-minute averaged wind by 9% for moored buoys c******************************************************* if(input.eq.'44007h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44005h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44013h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44028h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44011h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44025h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44008h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44009h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44004h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'44014h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41001h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41004h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41002h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41008h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41010h1997.txt') adjwspd=adjwspd*1.09 if(input.eq.'41009h1997.txt') adjwspd=adjwspd*1.09 50 write(20,200) int(yy),mm,int(dd),int(hh),int(wd),adjwspd, & gst,wvht,dpd,apd,mwd,bar,atmp,wtmp,dewp,vis 200 format(i2,1x,i2,1x,i2,1x,i2,1x,i3,1x,f4.1,1x,f4.1,1x,f5.2, & 1x,f5.2,1x,f5.2,1x,i3,1x,f6.1,1x,f5.1,1x,f5.1,1x, & f5.1,1x,f4.1) 5000 continue 500 end subroutine bisect(f,xleft,xright,xmid,a,b) c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c Bisection program c c Written by Pat Fitzpatrick, Jackson State University c c Originally written February 14, 1996 c Updated May 14, 1998 c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c This program calculates the root of a function by solving f(x)=0 c with an error tolerance of 0.001 (this can be modified by changing c eps). For example, if the user wants to find the cube root of 10, c f(x) = x**3 - 10 c c User must define the function f(x) and the interval for within which c the root exists (xleft & xright) in a function subprogram. c c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c Variables used: c c f(x) --- function to be solved for f(x)=0 c eps --- error tolerance c xleft --- left side of interval c xright --- right side of interval c xmid --- midpoint of interval c c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c********************************************************************* c Set error tolerance c********************************************************************* eps=0.001 c********************************************************************* c Check to see if interval properly chosen. If c f(xright,a,b)*f(xleft,a,b) > 0 , c then interval improperly chosen because either xleft and xright are c both left of the root or both right of the root c c If interval incorrect, program stops, and a message prints to c screen notifying the user of the problem. c********************************************************************* if( f(xright,a,b)*f(xleft,a,b).gt.0.0) then print*,f(xright,a,b),f(xleft,a,b) print*,a,b print*,'Improper interval chosen; program stops' stop endif c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c Main bisection routine exists here. There are two separate c sections, since either f(xleft,a,b) < 0 and f(xright,a,b) > 0 , c or f(xleft,a,b)) > 0 and f(xright,a,b)) > 0 . c If after 40 iterations no c answer is found, then either eps has been set too small or c more loops are required, and the program stops with a warning c message. c c The steps are as follows: c c 1) The midpoint value of the interval "xmid" is computed. c c 2) f(xmid,a,b) is then computed. c c 3) If f(xmid,a,b) < absolute value of error tolerance, c calculation is completed. Goto line 20 and print the c answer. Otherwise, continue to 4) c c 4) Insert xmid into the interval limit which has the same function c sign as f(xmid,a,b) . Goto part 1. c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& c This section is used if f(xleft,a,b)<0 and f(xright,a,b)>0 c********************************************************************* if(f(xleft,a,b).lt.0.0) then do i=1,40 xmid=(xright+xleft)/2.0 if(f(xmid,a,b).lt.eps.and.f(xmid,a,b).gt. & (-1.0*eps))goto 20 if(f(xmid,a,b).lt.0.0) then xleft=xmid else xright=xmid endif enddo else c********************************************************************* c Otherwise this section is used since f(xleft,p,thetae)>0 and c f(xright,p,thetae)<0 c********************************************************************* do i=1,40 xmid=(xright+xleft)/2.0 if(f(xmid,a,b).lt.eps.and.f(xmid,a,b).gt. & (-1.0*eps)) goto 20 if(f(xmid,a,b).gt.0.0) then xleft=xmid else xright=xmid endif enddo endif print*,'Iterated solution not found in subroutine bisect' stop 20 return end function f(ustar,z,wspd) von=0.4 g=9.8 alpha=0.015 dum=z*g/(alpha*ustar**2) f=(von*wspd)**2/(log(dum))**2-ustar**2 return end