#include #include double tau1,delta,rcp,rad(),IPstrat(),IPtrop(); double solar ; double sig = 5.67e-8; main() { double Tg,Ttp; double delta1,delta2,rad1,radmid,rad2; int its; /*Get input parameters*/ /*printf("enter absorbed solar radiation (W/m**2)\n");*/ scanf("%lf",&solar); /*printf("enter tau inf\n");*/ scanf("%lf",&tau1); /*printf("enter R/Cp\n");*/ scanf("%lf",&rcp); if(rcp < 0.) rcp = 2./7.; html_header(); /*Write out input parameters*/ printf("

Absorbed Solar Flux (w/m**2)

  %lf

 Optical Thickness of Atmosphere

  %lf

 R/Cp

  %lf
",solar,tau1,rcp); /*Determine tropopause pressure (relative to surface pressure). delta is ratio (p(trop)/p(surf)). */ /* The following iteration loop solves for rad(delta)=1 by bisection. */ delta1 = 1.e-8;delta2=1. - 1.e-8; for(its=0;its<100;its++) { rad1 = rad(delta1);rad2=rad(delta2);radmid=rad(.5*(delta1+delta2)); if( (rad1-1.)*(radmid-1.) <= 0.) delta2 = .5*(delta1+delta2); if( (radmid-1.)*(rad2-1.) <= 0.) delta1 = .5*(delta1+delta2); if( fabs(delta1-delta2) < 1.e-7) break; } delta = .5*(delta1+delta2); /* Determine tropopause temperature from delta and solar input */ Ttp = pow(solar*(1.+tau1*delta) /(2.*sig), .25); /*Warning. The following estimate of surface temperature becomes innacurate if delta gets too large.*/ Tg = Ttp*pow(1./delta,rcp); /*Write out results*/ printf("

 Tropopause pressure

  %f times surface pressure

 Tropopause Temperature

 %f Kelvin

 Surface Temperature

 %f Kelvin

", delta, Ttp, Tg); printf("\n\n"); } double rad(delta) double delta; { return(IPstrat(delta)/IPtrop(delta));} double IPtrop(delta) double delta; /*This computes the upward IR radiation coming from the troposphere, assuming that the temperature is equal to the radiative equilibrium temperature at the tropopause. It is normalized to the solar flux.*/ { double sum,eta,deta,a; deta = .005; a = 1. - delta; /* a is the ratio of the optical depth of the tropopause to the optical depth of the whole atmosphere */ sum = 0.; for(eta=0.;eta<1.;eta += deta) sum += deta*( pow((1.-eta*a)/(1.-a),4.*rcp)*exp(-a*tau1*(1.-eta))); return( .5*(1+tau1*delta)*( a*tau1*sum + exp(-a*tau1)/pow(1.-a,4.*rcp) ) ); } double IPstrat(delta) double delta; { /*Computes upward IR flux based on the stratospheric radiative equilibrium value . Normalized to solar flux.*/ double a; a = 1. - delta; return( .5*(2+tau1*(1.-a)) ); } html_header() { printf("\n \n Greygas Results\n"); printf("\n \n"); printf("

Grey Gas Model Results\n"); }