program enclosure_convection
  implicit none
  integer :: etype,ftype,i
  double precision :: Lg,Hh,Ww,Di,Do,Thot,Tcold
  double precision :: rho,mu,kf,Pr,cp,beta
  double precision :: Tf,dT,nu_v,alp,Ra,Nu,h_eff,Q,keff,Lc,HL,g
  double precision :: dTs,Ras,Nus,Qs
  double precision, parameter :: pi=3.14159265358979d0
  g=9.80665d0
  read(*,*) etype; read(*,*) Lg; read(*,*) Hh; read(*,*) Ww
  read(*,*) Di; read(*,*) Do
  read(*,*) Thot; read(*,*) Tcold
  read(*,*) ftype; read(*,*) rho; read(*,*) mu
  read(*,*) kf; read(*,*) Pr; read(*,*) cp; read(*,*) beta
  Tf=(Thot+Tcold)/2d0; dT=abs(Thot-Tcold)
  if(ftype==1) then
    rho=1.177d0;mu=1.85d-5;kf=0.0263d0;Pr=0.71d0;cp=1007d0;beta=1d0/(Tf+273.15d0)
  elseif(ftype==2) then
    rho=997d0;mu=8.9d-4;kf=0.613d0;Pr=6.13d0;cp=4180d0;beta=2.1d-4
  elseif(ftype==3) then
    rho=870d0;mu=0.05d0;kf=0.14d0;Pr=500d0;cp=2000d0;beta=7d-4
  endif
  nu_v=mu/rho; alp=kf/(rho*cp)
  Ra=g*beta*dT*Lg**3/(nu_v*alp)
  Nu=1d0; Q=0d0; keff=kf; Lc=Lg; HL=1d0
  select case(etype)
  case(1)
    HL=Hh/max(Lg,1d-6)
    if(HL<=2d0) then
      Nu=0.18d0*(Pr*Ra/(0.2d0+Pr))**0.29d0
    elseif(HL<=10d0) then
      Nu=0.22d0*(Pr*Ra/(0.2d0+Pr))**0.28d0*(HL)**(-0.25d0)
    else
      Nu=0.42d0*Ra**0.25d0*Pr**0.012d0*(HL)**(-0.3d0)
    endif
    if(Nu<1d0) Nu=1d0
    h_eff=Nu*kf/Lg; Q=h_eff*Hh*Ww*dT
  case(2)
    Nu=0.069d0*Ra**(1d0/3d0)*Pr**0.074d0
    if(Nu<1d0) Nu=1d0
    h_eff=Nu*kf/Lg; Q=h_eff*Hh*Ww*dT
  case(3)
    Nu=1d0; h_eff=kf/Lg; Q=h_eff*Hh*Ww*dT
  case(4)
    Lc=(Do-Di)/2d0
    Ra=g*beta*dT*Lc**3/(nu_v*alp)
    keff=0.386d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ra**0.25d0
    if(keff<kf) keff=kf
    Q=2d0*pi*keff*dT/log(Do/max(Di,1d-6))
    Nu=keff/kf; h_eff=Nu*kf/Lc
  case default
    Lc=(Do-Di)/2d0
    Ra=g*beta*dT*Lc**3/(nu_v*alp)
    keff=0.74d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ra**0.25d0
    if(keff<kf) keff=kf
    Q=pi*Di*Do*keff*dT/Lc
    Nu=keff/kf; h_eff=Nu*kf/Lc
  end select
  write(*,'(A)') '============================================'
  write(*,'(A)') '  FREE CONVECTION IN ENCLOSURES'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,I4)') '  Enclosure Type          = ',etype
  write(*,'(A,F10.4,A)') '  Gap L                   = ',Lg,' m'
  write(*,'(A,F10.4,A)') '  Height H                = ',Hh,' m'
  write(*,'(A,F10.4,A)') '  Width W                 = ',Ww,' m'
  write(*,'(A,F10.2,A)') '  T_hot                   = ',Thot,' C'
  write(*,'(A,F10.2,A)') '  T_cold                  = ',Tcold,' C'
  write(*,'(A,F10.2,A)') '  Delta T                 = ',dT,' C'
  write(*,'(A)') ''
  write(*,'(A)') '--- FLUID ---'
  write(*,'(A,F10.4,A)') '  rho                     = ',rho,' kg/m3'
  write(*,'(A,ES12.4,A)') '  mu                      = ',mu,' Pa.s'
  write(*,'(A,F10.6,A)') '  k                       = ',kf,' W/mK'
  write(*,'(A,F10.4)') '  Pr                      = ',Pr
  write(*,'(A,ES12.4,A)') '  beta                    = ',beta,' 1/K'
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,ES14.4)') '  Rayleigh Number Ra      = ',Ra
  write(*,'(A,F12.4)') '  Nusselt Number Nu       = ',Nu
  write(*,'(A,F12.4,A)') '  Effective h             = ',h_eff,' W/m2K'
  write(*,'(A,F12.4,A)') '  k_effective             = ',keff,' W/mK'
  write(*,'(A,F12.4,A)') '  Heat Transfer Q         = ',Q,' W (or W/m)'
  if(etype<=3) write(*,'(A,F8.2)') '  Aspect H/L              = ',HL
  write(*,'(A)') ''
  write(*,'(A)') '--- DELTA-T SWEEP ---'
  write(*,'(A)') '  dT[C]      Ra           Nu        Q[W]'
  write(*,'(A)') '  ---------------------------------------------------'
  do i=1,25
    dTs=1d0+(dT*3d0-1d0)*dble(i-1)/24d0
    if(etype<=3) then; Ras=g*beta*dTs*Lg**3/(nu_v*alp)
    else; Ras=g*beta*dTs*Lc**3/(nu_v*alp); endif
    select case(etype)
    case(1)
      if(HL<=2d0) then; Nus=0.18d0*(Pr*Ras/(0.2d0+Pr))**0.29d0
      elseif(HL<=10d0) then; Nus=0.22d0*(Pr*Ras/(0.2d0+Pr))**0.28d0*HL**(-0.25d0)
      else; Nus=0.42d0*Ras**0.25d0*Pr**0.012d0*HL**(-0.3d0); endif
      if(Nus<1d0) Nus=1d0; Qs=Nus*kf/Lg*Hh*Ww*dTs
    case(2)
      Nus=0.069d0*Ras**(1d0/3d0)*Pr**0.074d0
      if(Nus<1d0) Nus=1d0; Qs=Nus*kf/Lg*Hh*Ww*dTs
    case(3); Nus=1d0; Qs=kf/Lg*Hh*Ww*dTs
    case(4)
      keff=0.386d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ras**0.25d0
      if(keff<kf) keff=kf; Nus=keff/kf
      Qs=2d0*pi*keff*dTs/log(Do/max(Di,1d-6))
    case default
      keff=0.74d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ras**0.25d0
      if(keff<kf) keff=kf; Nus=keff/kf
      Qs=pi*Di*Do*keff*dTs/Lc
    end select
    write(*,'(2X,F8.2,2X,ES10.3,2X,F9.3,2X,F12.4)') dTs,Ras,Nus,Qs
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Vertical rect: Catton (1978) correlations'
  write(*,'(A)') '  Horizontal hot bottom: Globe-Dropkin Nu=0.069*Ra^(1/3)*Pr^0.074'
  write(*,'(A)') '  Concentric cyl: k_eff/k=0.386(Pr/(0.861+Pr))^0.25*Ra^0.25'
  write(*,'(A)') '  Concentric sph: k_eff/k=0.74(Pr/(0.861+Pr))^0.25*Ra^0.25'
  write(*,'(A)') '  Ref: Incropera Ch.9 Sec.9.8, Kothandaraman Ch.10'
end program enclosure_convection