program calc1dpbl parameter (n=100,ns=10) real ugi(ns),vgi(ns),tgi(ns),dsoil(ns),wsoil(ns),tsoil(ns) real zz(n),un(n),vn(n),ss(n),rh(n),pwfl(n),zsp(n),pblk(n), $ ug(n),vg(n) real sigguma (n,31), sigmaall(n,18,184,31) character*73 outfile,ascfile,fname,fext character*10 infile character*72 texti integer dcount, filecount, ifilecount c c program designed to output hourly files for 1D PBL, FSU v. 2.2 c July 2006 phr and Met 5920 Colloquium class c c query user for file name c c starting the count of files for statistics filecount=0 c the number of hour forecasted (include 0 in count) fcastmaxt=25 print *, ' Enter file name UPA Station from which to extract data' print *, ' ' read (5,*) infile OPEN (unit=12, form='formatted', file=infile, status='old') print*, infile 313 READ (12, 314, END=301) outfile print*, outfile c print *, ' Enter file name from which to extract data' c print *, ' ' c read (5,*) outfile lunb=2 OPEN (unit=LUNB,form='unformatted',file=outfile, $ status='old', ERR=313) filecount=filecount+1 print*, 'filecount is ',filecount c c input data from binary output file c print *, ' file ',outfile,'.bin is open.' print *, ' ' c c strip off any extra blanks on file name (if any added by compiler) c and create ascii version of .bin file (to be called .asc c lenf = len(outfile) kskip=0 kfb = 0 do k=1,lenf if (outfile(k:k).eq.']') kfb=k if (outfile(k:k).eq.' '.and.k.ne.1.and.kskip.eq.0) then kfe = k-kfb-1 kskip = 1 endif enddo lenf = kfe fname = outfile(kfb+1:lenf-4) lenf = lenf-4 fext(1:4) = '.asc' ascfile(1:lenf) = fname(1:lenf) ascfile(lenf+1:lenf+4) = fext(1:4) lenf = lenf+4 c c input control variables from binary file (same as in .dmp file) c read (lunb) ifhf,ifcri,ifsno,ifcld READ(LUNB) NOOFGR,DELTAT,TEND,RICR,PINK,KOOL READ(LUNB) TEXTI READ (LUNB)SLA,SLO,TZONE,Z0,Z0H,ZD0,ALBEDO READ (LUNB)MO,IDY,TIMEIS READ (LUNB)PSFC READ (LUNB)TREF READ (LUNB) CLC READ (LUNB) CMC READ (LUNB)PRST,PREND,PRCIP,ESD,TSNOW READ (LUNB)ISOIL,TWILT,SIGMAF,IFTC,TSO0,TSOREF,PC,RC print *,' ' print *,' Preliminaries written; ',sla,slo,timeis print *,texti open (unit=10,status='new',form='formatted',file=ascfile) write (10,100) texti write (10,*) ifhf,ifcri,ifsno,ifcld write (10,*) noofgr,deltat,tend,ricr,pink,kool write (10,*) sla,slo,tzone,zo,zoh,zdo,albedo write (10,*) mo,dy,timeis,psfc,tref,clc,cmc write (10,*) prst,prend,prcip,esd,tsnow write (10,*) isoil,twilt,sigmaf,iftc,tso0,tsoref,pc,rc write (10,*) c c Leaving print * statements in for debugging purposes in case changes c are needed later c READ (LUNB) NUG c print *,' # of geostrophic winds = ',nug do i=1,nug READ (LUNB) UGI(I),VGI(I),TGI(I) c print *,tgi(i),ugi(i),vgi(i) enddo READ (LUNB) NSOIL c print *,' ' c print *,' # of soil layers = ',nsoil do js=1,nsoil READ (LUNB) DSOIL(JS),WSOIL(JS),TSOIL(JS) c print *,dsoil(js),tsoil(js),wsoil(js) enddo READ (LUNB) mz1 print *,' # of input levels = ',mz1 write (10,101) nug,(tgi(i),ugi(i),vgi(i),i=1,nug) write (10,102) ! write soil init data now write (10,101) nsoil,(dsoil(i),tsoil(i),wsoil(i),i=1,nsoil) c READ (LUNB) MZ write (10,103) mz1,mz print *,' # of levels this step = ',mz c c End of preliminary header information/control variables; now get c hourly dump data and print out to ASCII file for convenient viewing c c c dcount=1 c c 10 DO jack=1,fcastmaxt read (lunb) mo,idy,ido,iho,min,sec c print *,' Date/time (mmddhhmm ss.ss) = ',mo,idy,iho,min, c $ sec,', day ',ido write (10,110) write (10,120) mo,idy,iho,min,sec,ido c c This loop gets all data for vertical profiles c do j=1,mz READ (LUNB) ZZ(J),UN(J),VN(J),TH,TC,QQ, $ P,SS(J),RH(J),PUFL,PVFL,PHFL,PWFL(J), $ Q1,Q2,TD,ZSP(J),PBLK(J) c print *,j c print *,zz(j),un(j),vn(j),tc,td,pblk(j) write (10,130) zz(j),un(j),vn(j),th,tc,qq,p*0.01,ss(j),rh(j), $ pufl,pvfl,phfl,pwfl(j),q1,q2,td,zsp(j),pblk(j) sigguma (j,dcount)=ss(j) enddo c DO jj=1,mz sigmaall(jj, 1, filecount, jack)=ss(jj) sigmaall(jj, 2, filecount, jack)=zz(jj) sigmaall(jj, 3, filecount, jack)=un(jj) sigmaall(jj, 4, filecount, jack)=vn(jj) sigmaall(jj, 5, filecount, jack)=th sigmaall(jj, 6, filecount, jack)=tc sigmaall(jj, 7, filecount, jack)=qq sigmaall(jj, 8, filecount, jack)=p*0.01 sigmaall(jj, 9, filecount, jack)=rh(jj) sigmaall(jj, 10, filecount, jack)=pufl sigmaall(jj, 11, filecount, jack)=pvfl sigmaall(jj, 12, filecount, jack)=phfl sigmaall(jj, 13, filecount, jack)=pwfl(jj) sigmaall(jj, 14, filecount, jack)=q1 sigmaall(jj, 15, filecount, jack)=q2 sigmaall(jj, 16, filecount, jack)=td sigmaall(jj, 17, filecount, jack)=zsp(jj) sigmaall(jj, 18, filecount, jack)=pblk(jj) c print*, sigguma(jj, dcount) ENDDO dcount=dcount+1 c c Now get all time series data - surface variables, energy, water, and rad c budgets and surface/boundary layer parameters c READ (LUNB) FDOWN,H,SSOIL,E,FUP,SNOFLX,SUM,FNET,EP1,THBAR, $ QBAR,ETOT c print *,fdown,h,ssoil,e,fup,snoflx,sum,fnet,ep1,thbar, c $ QBAR,ETOT READ (LUNB) TAOX,TAOY,USTAR,WSTAR,UG(1),VG(1),GEOS,SINA, $ anem,wsc,clc c print *,TAOX,TAOY,USTAR,WSTAR,UG(1),VG(1),GEOS,SINA, c $ anem,wsc,clc READ (LUNB) HPBL,XL,RIB,RIF,RIR,TSFC,TAIR,CM,CH,CG,THERM c print *,hpbl,xl,rib,rif,rir,tsfc,tair,cm,ch,cg,therm READ (LUNB) PRCP,DEW,ACCP,ACCD c print *,prcp,dew,accp,accd READ (LUNB) CMC,ESD,PC,RC,SOLARN c print *,cmc,esd,pc,rc,solarn write (10,140) fdown,h,ssoil,e,fup,snoflx,sum,fnet,ep1,thbar, $ qbar,etot,taox,taoy,ustar,wstar,ug(1),vg(1),geos,sina,anem, $ wsc,clc,hpbl,xl,rib,rif,rir,tsfc,tair,cm,ch,cg,therm,prcp, $ dew,accp,accd,cmc,esd,pc,solarn c c Finally, get soil output and print it; last record in binary file c is soil output, so use end of file trap to end data extraction c do js=1,nsoil read (lunb,end=99) dsoil(js),wsoil(js),tss write (10,150) dsoil(js),tss,wsoil(js) enddo c c go to 10 ! next hour dump c jack=jack+1 ENDDO c 99 print *,' ' print *,ascfile,' has been written; EOF on binary file reached' close (unit=lunb) go to 313 ! next file in the series 301 continue c this concludes reading in the files. We are left with the array of c sigmaall which is (levels, 18 parameters, # of times for file, hour of file forecast) c sigmaall has mz levels - currently 74 c sigmaall has 18 parameters - the sigma coordinate is the first one c sigmaall has filecount # of files for times (integer) c sigmaall has fcastmaxt for the maximum number of forecast hours (include 0 in count) DO i=1,100 print *, 'a', sigmaall(i,2,1,1) ENDDO close (unit=12) stop c 100 format ("1FSU PBL Model Version 2.2 July 2006",//,1x,a72,//, $ " Number of geostrophic wind levels:",/) 101 format (i5,/,10(i5,f12.4)) 102 format (" Number of soil levels: "/) 103 format (" Number of input levels = ",i5,/, $ " Number of model levels = ",i5) 110 format ("1FSU PBL v. 2.2",//) 120 format (1H1," Mo Dy Hr Mn Sc Day",//4i5,f8.1,i5,//, $ " Ht u v th", $ " T q p sigma RH uflx vflx ", $ "uvfx wflx Q1 Q2 Td zsp K",//) 130 format (f8.0,4f6.1,f7.2,f7.1,f6.2,f5.0,6f8.2,f6.1,2f8.1) 140 format ("0Radiation Budget: FDOWN SENS SOIL EVAP", $ " FUP SNOW SUM FNET",/,15x,8f10.2,//, $ " Surface Info: EP1 THBAR QBAR ETOT", $ " TAOX TAOY USTR WSTR",/, $ 15x,8f10.2,/,20x,"UGs VGs GEOS SINA ANEM ", $ " WSC CLC",/,15x,7f10.2,/," PBL Parameters: HPBL", $ 7x,"L RIB RIF RIR TSFC TAIR CM ", $ " CH CG THERM",/,15x,11f8.2,//," Water Bal:",/, $ 14x,"PRCP DEW ACCP ACCD CMC ESD ", $ "PC SOLARN",//,10X,8f8.2,//, $ " Soil Output Z T w",//) 150 format (15x,3f12.4) 314 format (A73) end