	program alpa2010
	character fname*10
	print*,'################################'
	print*,'#    Alpa2002 - versao 2010    #'
	print*,'#    Experimento  PRE-CHUVA    #' 
	print*,'################################'
	print*
	print*,'Bora, macho! Entre com o prefixo'
	print*,'dos arquivos a processar! =)'
	print*
	print*,'Formato: AAAAMMDD-N'
	print*,'(A=ano, M=mes, D=dia, N=numero de ordem)'
	read*,fname
	call alpa(fname)
	call readfssp(fname)
	call alpa1hz(fname)
	stop
	end

	subroutine alpa(fname)
	implicit real*8(a-h,o-z)
	character fname*10,fname1*14,fname2*14,
     1          fname3*12,exten1*4,exten2*4,
     2          exten3*4
	PARAMETER(CMB2P=100.0,D027=1.176,CK2MS=0.514444)
	PARAMETER(GAMMA=1.399721448468,RD=287.05)
	PARAMETER(Cv=718.0,Cp=1005.0,rf=0.892)
	PARAMETER(P0=101325.0,Aa0=-6.5E-3)
	PARAMETER(G0=9.80665,T0=288.15,RE=6.3755E6)
	dimension var(300000,20),i200x(300000,16),
     1            i200y(300000,16),ivar(300000,10),
     2		    a0(4),a1(4),vargps(8,20000),
     3            spp(100,50000),s200x(15),d200x(15),
     4            s200y(15),d200y(15),c200x(15),c200y(15),
     5            dfssp(30),cfssp(30),ss0(15),v0(15)
	data exten1/'.GPS'/exten2/'.dat'/exten3/'.spp'/
	data a0/-56.7177775,-49.406,-1.34,0.141146/
	data a1/21.6192383,20.0676,108.54,17.2166/
	data s200x/0.561,4.234,12.188,22.355,32.413,31.11,29.28,
     1     27.45,25.62,23.79,21.96,20.13,18.3,16.47,14.64/
	data d200x/36.15,65.85,95.55,125.25,154.95,184.5,213.9,243.45,
     1     274.35,304.05,332.25,361.65,391.2,420.75,450.3/
	data s200y/1485.0,1657.0,1578.0,1499.0,1420.0,1341.0,1262.0,
     1     1183.0,1105.0,1026.0,947.0,868.0,789.0,710.0,631.0/
	data d200y/361.5,661.5,958.5,1254.0,1551.0,1845.0,2140.5,2437.5,
     1     2731.5,3028.5,3322.5,3607.5,3900.0,4200.0,4500.0/
	data ss0/.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,1.2,1.4,1.6,1.8,2./
	data v0/.35,.51,.63,.72,.81,.89,.96,1.02,1.08,1.14,1.25,1.35,
     1        1.43,1.51,1.59/
	data alfa/4.313/beta/-0.161/

	dfssp(1)=2.75
	do ifssp=2,30
	   dfssp(ifssp)=dfssp(ifssp-1)+1.5
	enddo
	pi=4.*atan(1.)
	pi6=pi/6.
	umter=1./3.

c	Correções introduzidas:
c	Nova calibracao do CCNC: alfa=3.359 / beta=-0.443
	print*

	fname1=fname//exten1
	fname2=fname//exten2
	fname3=fname(3:10)//exten3
	open(1,file=fname2,status='old')
	i=1
199	read(1,*,end=299)(ivar(i,j),j=1,7),(var(i,j),j=1,16),
     1           (i200y(i,j),j=1,16),(i200x(i,j),j=1,16)
	i=i+1
	go to 199
299	continue
	close(1)
	write(*,*)'dados analogicos e das sondas opticas lidos',i
	open(1,file=fname1,status='old')
	igps=0
399	read(1,*,end=499)idata,ihora,imin,isec,
     1     (vargps(ii,igps+1),ii=3,8)
	igps=igps+1
	vargps(1,igps)=float(idata)
	vargps(2,igps)=float(ihora)*10000+float(imin)*100+
     1   float(isec)
c	write(*,*)igps,(vargps(i,igps),i=1,8)
	go to 399
499	close(1)
	write(*,*)'dados do gps lidos',igps
c	open(1,file=fname//'.spp')
c	ispp=1
c	do j=1,34
c	   read(1,*)
c	enddo
c599	read(1,*,end=699)ihspp,imspp,(spp(ii,ispp),ii=1,9),nn
c	write(*,*)ihspp,imspp,(spp(ii,ispp),ii=1,9),nn
c	read(1,*)(spp(ii+9,ispp),ii=1,nn)
c	spp(100,ispp)=float(ihspp)*3600.+float(imspp)*60.+spp(1,ispp)
c	ispp=ispp+1
c	go to 599
c699	close(1)
	write(*,*)'dados da fssp lidos',ispp
	open(2,file=fname//'.cld')
	open(12,file=fname//'.200x')
c	open(22,file=fname//'.200y')
c	open(32,file=fname//'.fssp')
	write(2,9)'dia','mes','ano','hora','tempo','lon','lat','ps','qc',
     1  'ttf','tdp','temp','theta','rm','q','LWCk','ias','tas','gsgps',
     2  'h','altbar','altgps','Nfssp','N200x','N200y','dbarf','dbarx',
     3  'dbary','LWCf','LWCx','LWCy','ssat','ccnc'
	write(12,*)'dia mes ano hora tempo lon lat altbar altgps '//
     1  'N200x dbarx LWCx dmassx deffx 200x01 200x02 200x03 200x04 '//
     2  '200x05 200x06 200x07 200x08 200x09 200x10 200x11 200x12 '//
     3  '200x13 200x14 200x15'
c	write(22,*)'dia mes ano hora tempo lon lat altbar altgps '//
c     1  'N200y dbary LWCy dmassy deffy 200y01 200y02 200y03 200y04 '//
c     2  '200y05 200y06 200y07 200y08 200y09 200y10 200y11 200y12 '//
c	3  '200y13 200y14 200y15'
c
c	write(32,*)'dia mes ano hora tempo lon lat altbar altgps '//
c     1  'Nfssp dbarf LWCf dmassf defff fssp01 fssp02 fssp03 fssp04 '//
c	2  'fssp05 fssp06 fssp07 fssp08 fssp09 fssp10 fssp11 fssp12 '//
c	3  'fssp13 fssp14 fssp15 fssp16 fssp17 fssp18 fssp19 fssp20 '//
c	4  'fssp21 fssp22 fssp23 fssp24 fssp25 fssp26 fssp27 fssp28 '//
c	5  'fssp29 fssp30'
	do k=1,i-1
c	write(*,*)k
	   time=float(ivar(k,4))*3600.+float(ivar(k,5))*60.+
     1          float(ivar(k,6))+float(ivar(k,7))/1000.
c		write(*,*)time
	   do jgps=1,igps-1
	       call timegps(vargps(2,jgps),tgps1)
c	write(*,*)tgps1
c	pause
	       call timegps(vargps(2,jgps+1),tgps2)
		   if (tgps1.le.time.and.tgps2.gt.time) then
	          dt=time-tgps1
	          dt0=tgps2-tgps1
	          dv=vargps(3,jgps+1)-vargps(3,jgps)
	          alat=vargps(3,jgps)+dv*dt/dt0
		      dv=vargps(4,jgps+1)-vargps(4,jgps)
	          alon=vargps(4,jgps)+dv*dt/dt0
	          dv=vargps(5,jgps+1)-vargps(5,jgps)
	          altgps=vargps(5,jgps)+dv*dt/dt0
		      dv=vargps(6,jgps+1)-vargps(6,jgps)
		      gspeed=vargps(6,jgps)+dv*dt/dt0
			  gspeed=gspeed*ck2ms
	          go to 799
	       endif
799	       continue
	   enddo
	   itime=ivar(k,6)+100*ivar(k,5)+10000*ivar(k,4)
	   ps=a0(3)+a1(3)*var(k,1)
	   qc=a0(4)+a1(4)*var(k,2)
	   ttf=a0(1)+a1(1)*var(k,3)
	   tdp=a0(2)+a1(2)*var(k,4)
	   alwc=var(k,5)
! No procedimento combinado, AIAS=75 m/s
! qc=0.5*75*75/100.
	   aIAS=2.0*QC*CMB2P
	   aias=sqrt(max(aias,0.))
	   aias=75.
	   qc=0.5*75*75/100.
	   aMACH=2.0*Cv/RD * ((QC/PS + 1.0)**(RD/Cp) - 1.0)
	   amach=sqrt(max(amach,0.))
	   ttfk=ttf+273.15
	   TAMBxK=TTFK/( 1.0 +aMACH*aMACH*RF*(GAMMA-1.0)/2.0 )
	   tambx=tambxk-273.15
	   TAS=SQRT(GAMMA*RD*tambxk) * aMACH
	   TSL = TAMBXK *  ( (PS*CMB2P)/P0 )**(Aa0*RD/G0)
	   H = (TAMBXK-TSL)/Aa0 + (TSL-T0)/Aa0 * LOG(1.0+(TAMBXK-TSL)/T0)
	   ALT = (RE*H) / (RE-H)
	   pumpccn=var(k,15)
	   tempccn=var(k,16)
	   if (pumpccn.lt.1.) detccn=var(k,14)
	   if (pumpccn.gt.3.) detccn0=var(k,14)
	   voltccn=0.
	   if (pumpccn.lt.1.) voltccn=max(0.,detccn-detccn0)
	   do iss=1,14
	      if (tempccn.gt.v0(iss).and.tempccn.lt.v0(iss+1)) then
	         dv=tempccn-v0(iss)
	         dv0=v0(iss+1)-v0(iss)
	         ssat=ss0(iss)+dv*(ss0(iss+1)-ss0(iss))/dv0
	      endif
	   enddo
c	   ssat=0.2
c	   if (tempccn.gt.0.55) ssat=0.3
c	   if (tempccn.gt.0.7) ssat=0.5
c	   if (tempccn.gt.0.85) ssat=0.7
c	   ssat=tempccn*tempccn/2.5
	   do jspp=1,ispp-1
	      if (spp(100,jspp).le.time.and.spp(100,jspp+1).gt.time) then
c	         write(*,*)jspp,spp(100,jspp),time,spp(100,jspp+1)
	         dt=time-spp(100,jspp)
		     dt0=spp(100,jspp+1)-spp(100,jspp)
	         dv=spp(4,jspp+1)-spp(4,jspp)
			 fssp=(spp(4,jspp)+dv*dt/dt0)*tas/80.
		     dv=spp(5,jspp+1)-spp(5,jspp)
	         alwcf=(spp(5,jspp)+dv*dt/dt0)*tas/80.
	         dv=spp(7,jspp+1)-spp(7,jspp)
	         corrf=spp(7,jspp)+dv*dt/dt0
	         do ifssp=1,30
	            dv=spp(ifssp+9,jspp+1)-spp(ifssp+9,jspp)
	            cfssp(ifssp)=(spp(ifssp+9,jspp)+dv*dt/dt0)*tas/80.
	         enddo
	         go to 899
	      endif
899	      continue
	   enddo	      
	   ccn=1000*voltccn*(alfa+beta/ssat)
	   if (ssat.eq.0) ccn=0.
	   concx=0.
	   concy=0.
	   alwcx=0.
	   alwcy=0.
	   dbarx=0.
	   dmassx=0.
	   deffx=0.
	   dbary=0.
	   dmassy=0.
	   deffy=0.
	   dbarf=0.
	   dmassf=0.
	   defff=0.
	   do i200=1,15
c	write(*,*)concx,concy,alwcx,alwcy,dbarx,dbary
	      svol=max(tas*0.05*s200x(i200)*1.e-3,1.e-15)
	      c200x(i200)=float(i200x(k,i200+1))/svol
c	write(*,*)i200,k,i200x(k,i200+1),svol,c200x(i200)
	      concx=concx+c200x(i200)
		  alwcx=alwcx+c200x(i200)*pi6*(d200x(i200)**3)/1.e9
	      dbarx=dbarx+c200x(i200)*d200x(i200)/max(concx,1.e-10)
	      deffx=deffx+c200x(i200)*d200x(i200)*d200x(i200)/
     1            max(concx,1.e-10)
	      dmassx=dmassx+c200x(i200)*d200x(i200)*d200x(i200)*
     1            d200x(i200)/max(concx,1.e-10)
	      svol=max(tas*0.05*s200y(i200)*1.e-3,1.e-15)
	      c200y(i200)=float(i200y(k,i200+1))/svol
	      concy=concy+c200y(i200)
	      alwcy=alwcy+c200y(i200)*pi6*(d200y(i200)**3)/1.e9
	      dbary=dbary+c200y(i200)*d200y(i200)/max(concy,1.e-10)
	      deffy=deffy+c200y(i200)*d200y(i200)*d200y(i200)/
     1            max(concy,1.e-10)
	      dmassy=dmassy+c200y(i200)*d200y(i200)*d200y(i200)*
     1            d200y(i200)/max(concy,1.e-10)
	      radar=radar+1.e-15*(c200x(i200)*(d200x(i200)**6.)+
     1            c200y(i200)*(d200y(i200)**6.))
	   enddo
	   do ifssp=1,30
		  dbarf=dbarf+cfssp(ifssp)*dfssp(ifssp)/max(fssp,1.e-10)
	      defff=defff+cfssp(ifssp)*dfssp(ifssp)*dfssp(ifssp)
     1            /max(fssp,1.e-10)
	      dmassf=dmassf+cfssp(ifssp)*dfssp(ifssp)*dfssp(ifssp)*
     1            dfssp(ifssp)/max(fssp,1.e-10)
	   enddo
	   deffx=sqrt(deffx)
	   deffy=sqrt(deffy)
	   defff=sqrt(defff)
	   dmassx=dmassx**umter
	   dmassy=dmassy**umter
	   dmassf=dmassf**umter
c	write(*,*)radar
c	   radar=10.*log(radar+1.e-10)
c	write(*,*)tsl,h,alt
	   call termo(Tambx,Tambxk,tdp,ps,ew,rm,q,U,theta,
     1        thetae,tv,thetav)
c	   if (fssp.gt.20.) 
	write(2,10)(ivar(k,j),j=1,3),itime,time,
     1        alon,alat,ps,qc,ttf,tdp,tambx,theta,rm,q,
     2        alwc,aias,tas,gspeed,h,alt,altgps,fssp,
     3        concx,concy,dbarf,dbarx,dbary,alwcf,alwcx,alwcy,
     4        ssat,ccn
	   write(12,10)(ivar(k,j),j=1,3),itime,time,alon,alat,alt,altgps,
     1        concx,dbarx,alwcx,dmassx,deffx,(c200x(i200),i200=1,15)
c	   write(22,10)(ivar(k,j),j=1,3),itime,time,alon,alat,alt,altgps,
c     1        concy,dbary,alwcy,dmassy,deffy,(c200y(i200),i200=1,15)
c	   write(32,10)(ivar(k,j),j=1,3),itime,time,alon,alat,alt,altgps,
c     1        fssp,dbarf,alwcf,dmassf,defff,(cfssp(ifssp),ifssp=1,30)
9	format(3a6,a7,41a12)
10	format(3i6,i7,f12.3,40f12.4)
c	write(*,*)k
	enddo
	close(2)
	close(12)
	close(22)
	return
	end

	subroutine timegps(thms,t)
	implicit real*8(a-h,o-z)
	itgps=int(thms)
	ih=itgps/10000
	im=itgps/100-100*ih
	is=itgps-100*im-10000*ih
c	write(*,*)ih,im,is
	t=3600.*float(ih)+60.*float(im)+is
	return
	end

c******************************************************************************************
c
c Esta subrotina é  parte do programa ALPA-2002 e calcula parametros termodinamicos:
c pressão do vapor d'água ew(mb), razão de mistura rm(g/Kg), umidade específica q(g/Kg),
c umidade absoluta U(g/m3), temperatura potencial theta(K), temperatura potencial 
c equivalente thetae(K), temperatura virtual tv(C), temperatura potencial virtual 
c thetav (K), Umidade absoluta U(g/m3), Umidade especifica(g/kg), razao de mistura (g/kg).
c
c                          By Emerson Mariano da Silva - Fortaleza/2002.
c*****************************************************************************************
 

      subroutine termo(T,Tk,tdp,ps,ew,rm,q,U,theta,thetae,tv,thetav)
	implicit real*8(a-h,o-z)

c####################################################################
c
c Esta função calcula a pressão de vapor d'água (ew) em mb, obtido pelo
c metodo de Paul R. Lowe(1977). O erro nos calculos é menor que 1% para 
c faixa de temperatura do ambiente (T) entre -50 e 50 Celsius.
c
c Fonte: http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html 
c###################################################################

      
       if(T.lt.-50.0) then
       

        ews=(4.4685+T*(0.27347+T*(6.8381e-7+T*(8.7094e-5+T*(5.63513e-7
     &         +T*1.47796e-9))))) 
        
        end if 

        if (T.ge.-50) then
	 
	 ews=(6.107799961+T*(0.4436518521+T*(0.01428945805+
     &      T*(2.650648471e-4
     &      +T*(3.031240396e-6+T*(2.034080948e-8+T*6.136820929e-11))))))

	 end if
       if(Tdp.lt.-50.0) then
       

        ew=(4.4685+Tdp*(0.27347+Tdp*(6.8381e-7+Tdp*(8.7094e-5+
     1     Tdp*(5.63513e-7+Tdp*1.47796e-9))))) 
        
        end if 

        if (Tdp.ge.-50) then
	 
	 ew =(6.107799961+Tdp*(0.4436518521+Tdp*(0.01428945805+
     1     Tdp*(2.650648471e-4+Tdp*(3.031240396e-6+Tdp*
     2	 (2.034080948e-8+Tdp*6.136820929e-11))))))

	end if

c####################################################################
c
c Esta função calcula a razão de mistura (g/Kg). Onde ew(mb) é a 
c pressão do vapor d'água e ps(mb) e a pressão estática.
c
c Fonte: http//raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c###################################################################

       rm=622.0*(ew/(ps-ew))     

c####################################################################
c
c Esta função calcula a umidade específica q (g/Kg). Onde ew(mb) é a 
c pressão do vapor d'água e ps(mb) e a pressão estática.
c
c Formula obtida em http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c###################################################################  

          q=622.0*(ew/(ps-0.378*ew))

c####################################################################
c
c Esta função calcula a umidade absoluta (g/m3). Onde ew(mb) é a 
c pressão do vapor d'água e Tk a temperatura do ambiente(K).
c
c Fonte: http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c################################################################### 
     
       U=216.68*(ew/Tk)

c####################################################################
c
c Esta função calcula a temperatura potencial(theta) em Kelvin(K).
c Onde Tk é a temperatura do ambinte em Kelvin (K).
c
c  theta=T(1000/ps)**R/cp ; R=287,05 e cp=1005 ; R/cp=0,286 
c###################################################################

       theta=Tk*(1000.0/ps)**0.286       

c####################################################################
c
c Esta função calcula a temperatura potencial equivalente (thetae) em kelvin.  
c derivada pelo método de Bolton (1980). Onde rm(g/kg), tlcl é 
c a temperatura no nível de condensação livre em Kelvin (K), T(k),
c e ew(mb).
c   Fonte: http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c###################################################################
       
       tlcl=((2840.0/(3.51*log(Tk)-log(ew)-4.805))+55.0)
       thetae=theta*(((3.376/Tlcl)-0.00254)*(rm*(1.0+0.00081*rm)))       

c####################################################################
c
c Esta função calcula a temperatura virtual (tv) em celsius (C).  
c Smithsonian Meteorologycal Tables (1958).
c Onde q é a umidade específica em g/kg.
c
c   Fonte: http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c###################################################################

        tv=(Tk*(1.0+1.6078*q)/(1.0+q))-273.16      

c####################################################################
c
c Esta função calcula a temperatura potencial virtual (tv) em celsius (C).  
c ps (mb), R/cp=0.286.
c
c   Fonte: http://raf.atd.ucar.edu/Bulletins/b9appdx_B.html
c###################################################################

      thetav=(tv+273.16)*(1000/ps)**0.286
       
	return
        end

	subroutine readfssp(fname)
	dimension ncount(20),fdist(20),gdist(20),diam(20),dd(20),amass(20)
	character*80 aline
	character*8 date
	character*6 hour
	character fname*10,fname1*14,fname2*14,
     1          fname3*12,exten1*4,exten2*4,
     2          exten3*4

	data exten1/'.txt'/exten2/'.spp'/

	tas=7500. !true air speed (cm/s)
	beam=0.02 !beam diameter (cm)
	dof=0.25  !depth of field (cm)
	fr=1.     !frequency (Hz) 
	pi=4.*atan(1.)
	pi6=pi/6. !pi/6
	row=1.e3  !liquid water density (kg/m3)

	diam(1)=5.
	dd(1)=2.
	do i=2,15
	   diam(i)=diam(i-1)+2.
	   dd(i)=2.
	enddo
	do i=16,20
	   diam(i)=diam(i-1)+3.
	   dd(i)=3.
	enddo

	do i=1,20
	   dmeter=diam(i)*1.e-6
	   amass(i)=pi6*row*(dmeter**3) !droplet mass in kg
	enddo

!	open FSSP raw data
	open(1,file=fname//exten1,status='old')
!	create file for FSSP processed data
	open(2,file=fname//exten2)
	write(2,*)
     1 'date    hour    conc    lwc    number and mass distributions'

!	reads FSSP counts
19	continue
	read(1,'(a25)')aline
	write(*,*)aline,aline(2:5)
	if (aline(2:5).eq.'DATA') go to 29
	go to 19
29	continue
!	write(*,*)'reading counts'
	read(1,'(a25)',end=39)aline
!	write(*,*)aline(5:8)
!	write(*,*)aline(10:11)
!	write(*,*)aline(13:14)
!	write(iyear,'(i4.4)')aline(5:8)
!	write(*,*)iyear
!	write(imon,'(i2.2)')aline(10:11)
!	write(idate,'(i2.2)')aline(13:14)
!	write(ihour,'(i2.2)')aline(16:17)
!	write(imin,'(i2.2)')aline(19:20)
!	write(isec,'(i2.2)')aline(22:23)
!	write(*,*)iyear,imon,idate,ihour,imin,isec
	read(1,*)
	read(1,*)
	read(1,*)(ncount(i),i=1,10)
	read(1,*)(ncount(i),i=11,20)
!	write(*,*)aline(5:24),(ncount(i),i=1,20)

!	Sample volume
	sv=beam*dof*tas/fr

!	Total concentration
	tcount=0.
	do i=1,20
	  tcount=tcount+float(ncount(i))
	enddo
	conc=tcount/sv !concentration in cm-3

!	DSD in cm-3.um-1
	conc1=0.
	do i=1,20
	   fdist(i)=float(ncount(i))/sv/dd(i)
	   conc1=conc1+fdist(i)*dd(i)
	enddo

!	DSD in g.m-3.um-1
	do i=1,20
	   gdist(i)=fdist(i)*amass(i)*1.e9
	enddo

!	LWC
	alwc=0.
	do i=1,20
	   alwc=alwc+gdist(i)*dd(i)
	enddo
c	write(*,*)aline(16:17),aline(19:20),aline(22:23),conc,alwc

	date=aline(4:7)//aline(9:10)//aline(12:13)
	hour=aline(15:16)//aline(18:19)//aline(21:22)

	write(2,'(a8,1x,a6,1x,f7.2,f8.4,40e10.2)')
     1 date,hour,conc,alwc,(fdist(i),i=1,20),(gdist(i),i=1,20)

	go to 29
39	continue
	close(1)
	close(2)
	return
	end

	subroutine alpa1hz(fname)
	dimension ivar(4),var(29),vart(29)
	dimension jvar(4),varx(29),vartx(29)
	dimension hgt(500),press(500),temp(500)
	character fname*10,fname1*14,fname2*14,
     1          fname3*12,exten1*4,exten2*4,
     2          exten3*4
	PARAMETER(CMB2P=100.0,D027=1.176,CK2MS=0.514444)
	PARAMETER(GAM=1.399721448468,RD=287.05)
	PARAMETER(Cv=718.0,Cp=1005.0,rf=0.892)
	PARAMETER(P0=101325.0,Aa0=-6.5E-3)
	PARAMETER(G0=9.80665,T0=288.15,RE=6.3755E6)
	data exten1/'.cld'/exten2/'.out'/

	open(1,file=fname//exten1,status='old')
	open(2,file=fname//exten2)
	open(9,file=fname//'.snd')
!	open(3,file=fname//'.200x',status='old')
!	open(4,file=fname//'.oap')
	read(1,*)
!	read(3,*)
	write(2,9)'hora','tempo','lon','lat','psnd','tsnd','pair','qc',
     1  'ttf','tdp','temp','theta','rm','q','LWCk','ias','tas','gsgps',
     2  'h','altbar','altgps','Nfssp','N200x','N200y',
     3  'LWCf','LWCx','LWCy','ssat','ccnc'

	do i=1,10
	read(9,*)
	enddo

	do i=1,500
	read(9,*)itime,x1,x2,hgt(i),press(i),temp(i),rh,td,dir,vel
c	write(*,*)itime,hgt(i),press(i),temp(i)
	enddo	

19	continue
	do k=1,29
	   vart(k)=0.
	enddo
	do ii=1,20
	read(1,*,end=99)(ivar(j),j=1,4),(var(k),k=1,29)
c	time,
c     1        alon,alat,ps,qc,ttf,tdp,tambx,theta,rm,q,
c     2        alwc,aias,tas,gspeed,h,alt,altgps,fssp,
c     3        concx,concy,dbarf,dbarx,dbary,alwcf,alwcx,alwcy,
c     4        ssat,ccn
!	read(3,*,end=99)(jvar(j),j=1,4),(varx(k),k=1,25)

	do k=1,29
	   vart(k)=vart(k)+0.05*var(k)
	enddo

!	do k=1,25
!	   vartx(k)=vartx(k)+0.05*varx(k)
!	enddo

	enddo

!	conc=0.	
!	do i=1,15
!	   if (vartx(i+10).ne.0) write(*,*)ivar(4),vartx(i+10)
!	   conc=conc+vartx(i+10)*30.
!	enddo
!	if (conc.gt.0) write(*,*)'200x',conc,var(20),var(26)

c	vart(3)=1008.43-0.10743*(990.009-vart(3))/0.082248
	do ii=1,499
	 if (vart(18).gt.hgt(ii).and.vart(18).le.hgt(ii+1)) then
	   write(*,*)hgt(ii),vart(18),hgt(ii+1)
	   d1=vart(18)-hgt(ii)
	   d2=hgt(ii+1)-vart(18)
	   pr=press(ii)+(press(ii+1)-press(ii))*d1/(d1+d2)
	   tc=temp(ii)+(temp(ii+1)-temp(ii))*d1/(d1+d2)
	   go to 49
	 endif
	enddo
	   aias=75.
	   qc=0.5*75*75/100.
	   aMACH=2.0*Cv/RD * ((QC/PS + 1.0)**(RD/Cp) - 1.0)
	   amach=sqrt(max(amach,0.))
	   ttfk=tc+273.15
	   TAMBxK=TTFK/( 1.0 +aMACH*aMACH*RF*(GAM-1.0)/2.0 )
	   tambx=tambxk-273.15
	   TAS=SQRT(GAM*RD*tambxk) * aMACH
c     1        alon,alat,ps,qc,ttf,tdp,tambx,theta,rm,q,
c     2        alwc,aias,tas,gspeed,h,alt,altgps,fssp,
c     3        concx,concy,dbarf,dbarx,dbary,alwcf,alwcx,alwcy,
c     4        ssat,ccn
	vart(5)=qc
	vart(13)=aias
	vart(14)=tas
49	continue

	write(2,10)ivar(4),(vart(k),k=1,3),pr,tc,
     1  (vart(k),k=4,21),(vart(k),k=25,29)
!	write(4,11)ivar(4),conc,alwc,(vartx(k),k=11,25)
	
	go to 19
99	continue
	close(1)
	close(2)
	
 9	format(a7,40a12)
10	format(i7,f12.3,39f12.4)  
11	format(i7,17f10.4)     
	return
	end
	
