* Cloudscope data reduce program.
* Written by Pat Arnott.
* Started 4/10/98.
* USE AND ABUSE AT YOUR OWN RISK.
* Written as one huge program for portability.
PROGRAM cloudscopereduce
IMPLICIT REAL (a-h,o-z)
LOGICAL slice1(640,480),slice2(640,480)
LOGICAL offedge
INTEGER sumarea,sumoverlap
REAL mass
REAL area(10000),frame(10000)
REAL xcoor(10000),ycoor(10000),perim(10000)
REAL major(10000),minor(10000),angle(10000),pixnum(10000)
* New particle descriptors.
REAL nuarea(10000),nuframe(10000),nupartnum(10000)
REAL nuxcoor(10000),nuycoor(10000),nuperim(10000)
REAL numajor(10000),numinor(10000),nuangle(10000)
REAL nupixnum(10000),svolcorr(10000)
* Size bin descriptors.
INTEGER binlo(500),binhi(500),number(500)
REAL scalemajor(10000),scaleminor(10000),scaleArea(10000)
REAL nperlium(500),areabin(500),binmass(500),meandiam
* Other variables.
CHARACTER*1 dummy,ans,setting,scope,html(40),fn(40)
CHARACTER*1 headchk(17)
CHARACTER*8 date,time,timedone
CHARACTER*17 headcheck
DOUBLE PRECISION usec,uhms
CHARACTER*40 filename,htmlfile
INTEGER hrs,mins,isecs
CHARACTER*2 hh,mm,ss
* Common statements.
COMMON /time/ hrs,mins,isecs
EQUIVALENCE (html,htmlfile)
EQUIVALENCE (fn,filename)
EQUIVALENCE (headchk,headcheck)
* External statements.
EXTERNAL mass
1 FORMAT(a1)
121 FORMAT(a17)
8 FORMAT(a8)
nummax = 10000
2 WRITE(*,*) 'OPEN THE FILE CONTAINING RESULTS'
OPEN(1,FILE='',STATUS='OLD')
READ(1,121,ERR=2) headcheck
IF (headchk(17).EQ.'a') THEN
READ(1,1,ERR=2) dummy
DO 10 i=1,100000
10 READ(1,*,END=20,ERR=2) dum,area(i),xcoor(i),ycoor(i),perim(i),
* major(i),minor(i),angle(i),frame(i),pixnum(i)
ELSE
BACKSPACE(1)
DO 11 i=1,100000
11 READ(1,*,END=20,ERR=2) area(i),xcoor(i),ycoor(i),perim(i),
* major(i),minor(i),angle(i),frame(i),pixnum(i)
END IF
20 CLOSE(1)
numdat=i-1
3 WRITE(*,*) 'OPEN THE NEW PARTICLE REPORT FILE'
OPEN(1,FILE='',STATUS='UNKNOWN',ERR=3)
WRITE(1,*) 'N Area X Y Length Major Minor Angle Frame PixNum'
4 WRITE(*,*) ''
WRITE(*,*) 'OPEN THE SIZE DISTRIBUTION FILE'
OPEN(3,FILE='',STATUS='UNKNOWN',ERR=4)
* Define the html output file and open it.
INQUIRE(3,NAME=filename)
DO 400 i=1,40
400 html(i)=' '
DO 410 i=40,1,-1
IF (fn(i).NE.' ') THEN
iend=i
GOTO 420
END IF
410 CONTINUE
420 DO 430 i=1,iend
430 html(i)=fn(i)
html(iend+1)='.'
html(iend+2)='h'
html(iend+3)='t'
html(iend+4)='m'
html(iend+5)='l'
OPEN(2,FILE=htmlfile)
* Set the time and date of the data.
WRITE(*,*) ''
WRITE(*,*) 'ENTER TIME DATA WAS OBTAINED as hh:mm:ss'
READ(*,8) time
WRITE(*,*) 'ENTER DATE DATA WAS OBTAINED as mm_dd_yy'
READ(*,8) date
* Set up the cloudscope being used.
12 WRITE(*,*) 'ENTER THE CLOUDSCOPE USED FROM THE FOLLOWING'
WRITE(*,*) ' a = Standard Aircraft Cloudscope (default)'
WRITE(*,*) ' b = Large Format Cloudscope'
WRITE(*,*) ' c = Lab Cloudscope'
WRITE(*,*) ' d = Medium Format Cloudscope'
READ(*,1) scope
IF ((scope.EQ.' ').OR.(scope.EQ.'a')) THEN
scope='a'
xmax = 420. ! resolution in microns for full screen
ELSE IF (scope.EQ.'b') THEN
xmax = 38000. ! Large cloudscope resolution
ELSE IF (scope.EQ.'c') THEN
xmax = 175. ! Lab cloudscope resolution.
ELSE IF (scope.EQ.'d') THEN
xmax = 19000. ! Lab cloudscope resolution.
ELSE
WRITE(*,*) 'Your choices are a b c d or enter, try again'
GOTO 12
END IF
* Set up the image size.
*13 WRITE(*,*) 'ENTER THE STACK SIZE FROM THE FOLLOWING'
* WRITE(*,*) ' a = 640 x 480 pixels (default)'
* WRITE(*,*) ' b = 320 X 240 pixels'
* READ(*,1) ans
* IF ((ans.EQ.'a').OR.(ans.EQ.' ')) THEN
* pixels=640.
* ELSE IF (ans.EQ.'b') THEN
* pixels=320.
* ELSE
* WRITE(*,*) 'Your choices are a or b, try again'
* GOTO 13
* END IF
pixels=640.
scalefac = xmax/pixels
* Set up the setting, aircraft or ground operation.
14 WRITE(*,*) 'ENTER THE SETTING FROM THE FOLLOWING'
WRITE(*,*) ' a = aircraft operation (default)'
WRITE(*,*) ' b = ground or lab operation'
READ(*,1) setting
IF (((setting.EQ.'a').OR.(setting.EQ.' ')).OR.
* (setting.EQ.'b')) THEN
IF (setting.EQ.' ') setting='a'
CALL getaircraft(time,setting,tas,press,altitude,
* xlatitude,xlongitude,airtemp,totaltemp)
ELSE
WRITE(*,*) ' Your choices are a or b try again'
GOTO 14
END IF
* Setup the number of frames per second
WRITE(*,*) ''
WRITE(*,*) 'ENTER THE FRAMES PER SECOND DURING DIGITIZATION'
READ(*,*) fps
WRITE(*,*) 'Enter the NUMBER OF FRAMES digitized in the stack'
READ(*,*) nframedig
dt=REAL(nframedig)/fps
* Calculate the sample volume or sample area, use samplevol for both.
samplevol = ((pixels*scalefac)**2)*0.75
WRITE(*,*) ' '
IF (setting.EQ.'a') THEN
samplevol = samplevol*tas*dt*1.e-9 ! In liters
WRITE(*,*) 'Sample volume in Liters',samplevol
ELSE
samplevol = samplevol*1.e-8 ! In square cm.
WRITE(*,*) 'Sample area in square cm',samplevol
END IF
WRITE(*,*) ' '
WRITE(*,*) '*****************************************'
WRITE(*,*) '**** Working .... '
WRITE(*,*) '*****************************************'
WRITE(*,*) ' '
* NOW COMPUTE THE NEW PARTICLES FROM THE VIDEO DATA.
* Initialize the slices.
DO 30 ifx=1,640
DO 40 ify=1,480
slice1(ifx,ify)=.FALSE.
40 slice2(ifx,ify)=.FALSE.
30 CONTINUE
np=0 ! Particle number counter.
nup=0 ! New particle number counter.
DO 60 i=1,nummax
IF (frame(i).EQ.1) THEN
np=np+1
CALL fillslice(slice2,slice1,xcoor(i),ycoor(i),angle(i),
* major(i),minor(i),sumarea,sumoverlap,offedge)
ELSE
GOTO 50
END IF
60 CONTINUE
50 numframe=2 ! Slice counter.
* Evaluate the slices for particles.
DO 100 ndum=1,nummax
np=np+1
IF (np.LE.numdat) THEN
IF (frame(np).NE.numframe) THEN
* Initialize the slices.
DO 110 ifx=1,640
DO 120 ify=1,480
slice1(ifx,ify)=slice2(ifx,ify)
120 slice2(ifx,ify)=.FALSE.
110 CONTINUE
numframe=frame(np)
END IF
CALL fillslice(slice1,slice2,xcoor(np),ycoor(np),angle(np),
* major(np),minor(np),sumarea,sumoverlap,offedge)
overlap=REAL(sumoverlap)/REAL(sumarea)
* Check for new particle.
IF ((overlap.EQ.0.0).AND.(offedge)) THEN ! Found new particle!
nup=nup+1
nupartnum(nup)=np
nuArea(nup)=Area(np)
nuxcoor(nup)=xcoor(np)
nuycoor(nup)=ycoor(np)
nuperim(nup)=perim(np)
numajor(nup)=major(np)
numinor(nup)=minor(np)
nuangle(nup)=angle(np)
nuframe(nup)=frame(np)
nupixnum(nup)=pixnum(np)
scalemajor(nup)=scalefac*numajor(nup)
scaleminor(nup)=scalefac*numinor(nup)
scaleArea(nup)=scalefac*scalefac*nuArea(nup)
del=(scalemajor(nup)+scaleminor(nup))/2.
svolcorr(nup)=1.-4.*del*(1.75*xmax-2.*del)/(3.*xmax*xmax)
* Cloud properties are now summed up.
END IF
ELSE
GOTO 200
END IF
100 CONTINUE
* Write out the nu particle report, and scale to microns.
200 DO 210 i=1,nup
WRITE(1,*) nupartnum(i),nuArea(i),nuxcoor(i),nuycoor(i),
* nuperim(i),numajor(i),numinor(i),nuangle(i),nuframe(i),
* nupixnum(i)
210 CONTINUE
* Now bin the data and produce size distributions.
* a = Standard Aircraft Cloudscope
* b = Large Format Cloudscope
* c = Lab Cloudscope
* d = Medium Format Cloudscope
CALL setbins(scope,nbins,binlo,binhi)
extcoef = 0.
meandiam = 0.
totalmass = 0.
totalperli = 0.
DO 300 i=1,nbins
number(i) = 0
nperlium(i) = 0.
areabin(i)= 0.
300 binmass(i)= 0.
DO 310 i=1,nup
meandiam=meandiam + scalemajor(i)/(samplevol*svolcorr(i))
extcoef =extcoef + 2.*scaleArea(i)/(samplevol*svolcorr(i))
totalperli=totalperli+1./(samplevol*svolcorr(i))
DO 320 j=1,nbins
IF ( (scalemajor(i).LE.FLOAT(binhi(j))).AND.
* (scalemajor(i).GE.FLOAT(binlo(j))) )THEN
number(j)=number(j) + 1
dD = FLOAT(binhi(j)-binlo(j))
nperlium(j)=nperlium(j) + 1./(samplevol*svolcorr(i)*dD)
areabin(j) =areabin(j) +
* (1.e-3)*scaleArea(i)/(samplevol*svolcorr(i)*dD)
GOTO 310
END IF
320 CONTINUE
310 CONTINUE
meandiam=meandiam/totalperli
extcoef=extcoef*1.e-6 ! 1/km units.
DO 330 i=1,nbins
binmass(i)=mass(binlo(i),binhi(i),nperlium(i))
330 totalmass=totalmass + binmass(i)*dD
* Now write out the data in regular and in html format.
* First decode the time and determine the end time. Then reencode it.
DECODE (8,370,time) ih,im,is
370 FORMAT(i2,1x,i2,1x,i2)
usec=DFLOAT(dt + FLOAT(ih)*3600. + FLOAT(im)*60. + FLOAT(is))
CALL gct(usec,uhms) ! hrs,mins,isecs are loaded by the common.
IF (hrs.LT.10) THEN
ENCODE(2,381,hh) hrs
381 FORMAT('0',i1)
ELSE
ENCODE(2,382,hh) hrs
382 FORMAT(i2)
END IF
IF (mins.LT.10) THEN
ENCODE(2,381,mm) mins
ELSE
ENCODE(2,382,mm) mins
END IF
IF (isecs.LT.10) THEN
ENCODE(2,381,ss) isecs
ELSE
ENCODE(2,382,ss) isecs
END IF
timedone=hh//':'//mm//':'//ss
WRITE(2,*) '<HTML>'
WRITE(2,*) '<HEAD>'
WRITE(2,340) date,time
340 FORMAT('<TITLE>','Cloudscope',1x,a8,1x,a8,' UTC',1x,'</TITLE>')
WRITE(2,*) '</HEAD>'
WRITE(2,*) '<BODY>'
WRITE(2,*) '<PRE>'
* Write out the data to file.
DO 390 j=2,3
IF (j.EQ.2) WRITE(2,350) date,time
IF (j.EQ.3) WRITE(3,351) date,time
350 FORMAT('<B><FONT SIZE=+1>','DRI CLOUDSCOPE DATA from',
* 1x,a8,1x,a8,' UTC','</FONT></B>')
351 FORMAT('DRI CLOUDSCOPE DATA from',1x,a8,1x,a8,' UTC')
WRITE(j,*) 'StartTime_UTC as hh:mm:ss ',time
WRITE(j,*) ' EndTime_UTC as hh:mm:ss ',timedone
WRITE(j,*) 'Seconds analyzed',dt
WRITE(j,*) 'frames per second digitization',fps
WRITE(j,*) 'True Air Speed in m/s for the aircraft',tas
WRITE(j,*) 'Altitude km AboveSeaLevel',altitude
WRITE(j,*) 'Air pressure in mb',press/100.
WRITE(j,*) 'Temperature C',airtemp-273.16
WRITE(j,*) 'Window Temperature C',totaltemp-273.16
WRITE(j,*) 'Latitude deg',xlatitude
WRITE(j,*) 'Longitude deg',xlongitude
WRITE(j,*) ' '
WRITE(j,*) 'Total Sample Volume in Liters',samplevol
WRITE(j,*) 'Minimum Concentration 1/Li',1./samplevol
WRITE(j,*) ' '
WRITE(j,*) 'Total Concentration as #/Li ',totalperli
WRITE(j,*) 'MeanDiameter um ',meandiam
WRITE(j,*) 'Extinction Coefficient as 1/km ',extcoef
WRITE(j,*) 'Total Mass as mg/m^3',totalmass
WRITE(j,*) ' '
WRITE(j,*) 'The bulk density of ice was assumed as 900 kg/m^3.'
WRITE(j,*) 'Crystals had mass estimated from the mass dimension'
WRITE(j,*) 'relation given as Eq. 3 in Arnott,et.al. JGR, 1994.'
WRITE(j,*) ''
WRITE(j,*) ' '
WRITE(j,*) ' BinUM Tot# #/Li/um mm^2/m^3/um mg/m^3/um '
IF (j.EQ.2) WRITE(2,*) '<FONT SIZE=-3>'
DO 360 i=1,nbins
WRITE(j,*) binLo(i),number(i),nperlium(i),areabin(i),binmass(i)
WRITE(j,*) binHi(i),number(i),nperlium(i),areabin(i),binmass(i)
360 CONTINUE
390 CONTINUE
WRITE(2,*) '</FONT>'
WRITE(2,*) '</PRE>'
WRITE(2,*) '</BODY>'
WRITE(2,*) '</HTML>'
CLOSE(1) ! particle report file.
CLOSE(2) ! html file with size distribution.
CLOSE(3) ! Size distribution file for plotting purposes.
WRITE(*,*) 'ENTER y TO EVALUATE ANOTHER FILE'
READ(*,1) ans
IF (ans.EQ.'y') GOTO 2
END
*********************************************************************
* Subroutine fillslice. Places logical true where particles are,
* according to treating the particles as ellipses.
* Original coordinates are 0,639 X 0,479, so xcoord and ycoord are
* having 1 added to them below.
*********************************************************************
SUBROUTINE fillslice(oldslice,slice,xcoor,ycoor,angle,max,min,
* sumarea,sumoverlap,offedge)
IMPLICIT REAL (a-h,o-z)
LOGICAL slice(640,480),oldslice(640,480)
LOGICAL offedge
INTEGER sumarea,sumoverlap
REAL max, min
a=1.00*max/2. ! enlarge the particle size by 10%
b=1.00*min/2.
b2=b*b
ia=INT(a+0.5)
ib=INT(b+0.5)
rad=3.141592654*angle/180.
cs=COS(rad)
sn=SIN(rad)
sumarea=0
sumoverlap=0
offedge=.TRUE.
DO 10 ix=-ia,ia
DO 20 iy=-ib,ib
y=FLOAT(iy)
IF (y.GT.b) y=b
x=a*SQRT(1.-y*y/b2)
xm=-x
xix=FLOAT(ix)
IF ((xix.GE.xm).AND.(xix.LE.x)) THEN
ifx=INT(xcoor + cs*xix - sn*y + 1.5)
ify=INT(ycoor + sn*xix + cs*y + 1.5)
IF (ifx.LE.1) THEN
ifx=1
offedge=.FALSE.
END IF
IF (ifx.GE.640) THEN
ifx=640
offedge=.FALSE.
END IF
IF (ify.LE.1) THEN
ify=1
offedge=.FALSE.
END IF
IF (ify.GE.480) THEN
ify=480
offedge=.FALSE.
END IF
slice(ifx,ify)=.TRUE.
sumarea=sumarea+1
IF ((oldslice(ifx,ify)).AND.(slice(ifx,ify))) THEN
sumoverlap=sumoverlap+1
END IF
END IF
20 CONTINUE
10 CONTINUE
RETURN
END
*********************************************************************
* Subroutine setbins, sets the bin sizes for the various cloudscopes.
* a = Standard Aircraft Cloudscope, same bins as standard replicator.
* b = Large Format Cloudscope.
* c = Lab Cloudscope.
* d = Medium Format Cloudscope.
* bin sizes are all in microns.
*********************************************************************
SUBROUTINE setbins(scope,nbins,binlo,binhi)
IMPLICIT REAL (a-h,o-z)
CHARACTER*1 scope
INTEGER binlo(500),binhi(500)
LOGICAL t
t=((scope.EQ.'a').OR.(scope.EQ.'b'))
IF (t.OR.(scope.EQ.'d')) THEN
nbins=19
binlo(1)=0
binlo(2)=11
binlo(3)=22
binlo(4)=33
binlo(5)=44
binlo(6)=55
binlo(7)=66
binlo(8)=77
binlo(9)=88
binlo(10)=99
binlo(11)=110
binlo(12)=131
binlo(13)=164
binlo(14)=197
binlo(15)=230
binlo(16)=263
binlo(17)=329
binlo(18)=395
binlo(19)=461
DO 10 i=2,19
10 binhi(i-1)=binlo(i)
binhi(19)=527
IF (scope.EQ.'b') THEN
DO 20 i=1,19
binlo(i)=binlo(i)*35
20 binhi(i)=binhi(i)*35
END IF
IF (scope.EQ.'d') THEN
DO 21 i=1,19
binlo(i)=binlo(i)*18
21 binhi(i)=binhi(i)*18
END IF
END IF
IF (scope.EQ.'c') THEN
nbins=18
DO 30 i=1,11
binlo(i)=(i-1)*5
30 binhi(i)=i*5
DO 40 i=12,18
binlo(i)=55 + (i-12)*10
40 binhi(i)=55 + (i-11)*10
END IF
RETURN
END
**************************************
* Function to compute the mass dimensional relationships
**************************************
FUNCTION mass(iDL,iDH,c)
IMPLICIT REAL (a-h,o-z)
INTEGER iDL,iDH
REAL mass,c
* Mass dimension relationship is assumed of the form
* m(D)=aD^b where D is maximum dimension in microns,
* a and b are constants.
* DL is the max dimension at the left end of the size bin, DH is the
* max dimension at the right end, both given in microns.
* c is the concentration in #/Li/um.
* mass will come out in units mg/m^3/um.
* mass dimension relations are from Arnott, 1995, JGR.
DH=FLOAT(iDH)
DL=FLOAT(iDL)
IF (DH.LE.200.) THEN
a=1.27e-10
b=3.
ELSE
a=6.4e-9
b=2.27
END IF
e=b+1.
mass=1000.*a*c*(DH**e - DL**e)/((b+1.)*(DH-DL))
RETURN
END
**************************************
* Subroutine to compute hhmmss.ss time from seconds after midnite.
**************************************
SUBROUTINE gct(usec,uhms)
IMPLICIT DOUBLE PRECISION (a-h,o-z)
INTEGER hrs,mins,isecs
COMMON /time/ hrs,mins,isecs
hrs=INT(usec/3600.)
mins=INT((usec - 3600.*DFLOAT(hrs))/60.)
secs=usec - 3600.*DFLOAT(hrs) - 60.*DFLOAT(mins)
isecs=INT(secs+0.5)
uhms=secs+DFLOAT(100*mins+10000*hrs)
RETURN
END
**************************************
* Subroutine to compute time in seconds from time in hhmmss.ss .
**************************************
SUBROUTINE timeconv(uhms,usec)
IMPLICIT REAL*4 (a-h,o-z)
INTEGER hrs,mins
hrs=INT(uhms/10000.)
mins=INT((uhms-10000.*FLOAT(hrs))/100.)
secs=uhms-10000.*FLOAT(hrs)-100.*FLOAT(mins)
usec=FLOAT(60*mins+3600*hrs)+secs
RETURN
END
*******************************************************************
* Subroutine getaircraft. Obtain the aircraft properties by direct
* input, or by reading them from a file.
*******************************************************************
SUBROUTINE getaircraft(time,setting,tas,press,altitude,
* xlatitude,xlongitude,airtemp,totaltemp)
IMPLICIT REAL (a-h,o-z)
CHARACTER*8 time
CHARACTER*1 tempset,setting,ans
REAL d1(51),d2(51)
1 FORMAT(a1)
IF (setting.EQ.'a') THEN
WRITE(*,*) 'TO READ AIRCRAFT INFO FROM FILE ENTER a (default)'
READ(*,1) ans
IF ((ans.EQ.'a').OR.(ans.EQ.' ')) THEN
DO 2 j=1,10
2 WRITE(*,*) ' '
WRITE(*,*) 'OPEN THE AIR DATA FILE FOR THIS FLIGHT'
OPEN(69,FILE='',STATUS='OLD')
DECODE (8,370,time) ih,im,is
370 FORMAT(i2,1x,i2,1x,i2)
usec=FLOAT(ih)*3600. + FLOAT(im)*60. + FLOAT(is)
READ(69,*) (d1(i),i=1,51)
DO 10 i=1,100000
READ(69,*,END=20) (d2(j),j=1,51)
IF ((d2(1).GE.usec).AND.(d1(1).LE.usec)) THEN
tas=d2(5)
press=100.*d2(3)
altitude=d2(2)/1000. ! in km.
airtemp=273.16+d2(9)
xlatitude=d2(44)
xlongitude=d2(45)
GOTO 20
ELSE
DO 30 ij=1,51
30 d1(ij)=d2(ij)
END IF
10 CONTINUE
20 CLOSE(69)
ELSE
WRITE(*,*) 'ENTER THE AIRCRAFT TRUE AIR SPEED in m/s'
READ(*,*) tas
WRITE(*,*)
WRITE(*,*) 'ENTER THE AIR PRESSURE in mb'
READ(*,*) press
press=press*100. ! convert to Pascal
WRITE(*,*) 'ENTER THE AIRCRAFT ALTITUDE IN km'
READ(*,*) altitude
WRITE(*,*) 'ENTER THE AIRCRAFT LATITUDE IN DEG xx.xxx'
READ(*,*) xlatitude
WRITE(*,*) 'ENTER THE AIRCRAFT LONGITUDE IN DEG xxx.xxx'
READ(*,*) xlongitude
WRITE(*,*) 'ENTER THE AIR TEMPERATURE in Celsius'
READ(*,*) airtemp
airtemp = airtemp + 273.16
END IF
15 WRITE(*,*) 'ENTER THE WINDOW TEMPERATURE FROM...'
WRITE(*,*) ' c = compute the total temperature (default)'
WRITE(*,*) ' d = known measured temperature '
READ(*,1) tempset
IF ((tempset.EQ.'c').OR.(tempset.EQ.' ')) THEN
tempset='c'
Rgas = 8.3143*1000./29. ! Joules / kg K for air.
cp = 7.*Rgas/2. ! for air, same units as Rgas.
totaltemp=tas*tas/cp + airtemp
ELSE IF (tempset.EQ.'d') THEN
WRITE(*,*) 'Enter the window temperature in Celsius'
READ(*,*) totaltemp
totaltemp=totaltemp + 273.16
ELSE
WRITE(*,*) ' Enter c or d, try again...'
GOTO 15
END IF
ELSE IF (setting.EQ.'b') THEN
WRITE(*,*) 'ENTER THE WINDOW TEMPERATURE IN CELSIUS'
READ(*,*) totaltemp
totaltemp=totaltemp + 273.16
END IF
RETURN
END