c 湿潤断熱線を引くためのデータを作るプログラム
c 2000/10/11
c 2000/11/02 等値線上数値出力追加
c 2000/12/27 上の改良

c Variable

      real*8 p,mintw,maxtw,minp,maxp,thetaw,thetae,e,w,po,l
      real*8 ii,jj,f,er,ab,t0,t1,cp,r,tc,step,et0,wt0
      real*8 steptw,stepp,pp,texst,texst1
      integer twdata,pdata,fl
      parameter(mintw=240.,maxtw=310.,minp=10.,maxp=1050.,r=287.)
      parameter(steptw=5.,stepp=10.,er=0.01)
      parameter(po=1000.,cp=1004.,l=2.5*10**6,texst1=30)
      parameter(twdata=(maxtw-mintw)/steptw+1,pdata=(maxp-minp)/stepp+1)
      dimension twt(pdata,twdata)

c Theta_w Calculation & Output

      open(10,file='thetaw.dat')
      open(20,file='thetawtext.dat')

      pp=750.
      texst=texst1

      if(mintw-int(mintw/10)*10.ne.0) then
	fl=1
	else
	  fl=0
      end if

      do j=1,twdata

         write(10,50)
 50      format('>')

         do i=1,pdata
            ii=i
            jj=j
            p=(ii-1)*stepp+minp
            thetaw=(jj-1)*steptw+mintw
            tc=thetaw-273.15
            e=6.112*exp((17.67*tc)/(tc+243.5))
            w=0.622*e/(po-e)
            thetae=thetaw*exp(l/cp*w/thetaw)

            t0=213.15
            step=100
 10         continue
            if(step.eq.100/2**10) go to 20
            et0=6.112*exp((17.67*(t0-273.15))/((t0-273.15)+243.5))
            wt0=0.622*et0/(p-et0)
            f=thetae-t0*(po/p)**(r/cp)*exp(l/cp*wt0/t0)
            ab=abs(f)
            if(ab.lt.er) go to 20
              t0=t0+sign(step,f)
              step=step/2
              go to 10
 20         continue
            twt(i,j)=t0-273.15

            write(10,200)twt(i,j),p
 200        format(f8.2,f8.2)

            if(fl.ne.1.and.p.eq.pp) then
               write(20,110)twt(i,j),p,8,0,0,'CM',int(thetaw)
 110           format(f8.2,f8.2,i8,i8,i8,A8,i8)
               pp=pp-texst
            endif
         enddo

	if(fl.ne.1) then
	  fl=fl+1
	  else
	    fl=fl-1
	endif

      enddo

      close(10)
      close(20)

      stop
      end








