【急】麻烦帮我看一下这段C语言程序有什么错误 多谢

倔强的葱头

程序如下，多谢多谢！！！！

#include <stdio.h>
#include <math.h>
#include "types.h"
#include "date.h"
#include "numeric.h"
#include "units.h"

#define a 107.0
#define b 0.2
#define qoo 1367.0        /* solar constant (1367 W/m2) */
#define beta 0.17        /* Albedo */
#define c 0.25
#define d 0.5
#define k 13750.98708        /* conversion factor from solar angular units to
                           seconds (12/pi*3600) */

/*
*       Function petpar
*
*       Calculates mid-month daily photosynthetically active radiation flux,
*       daylength and daily potential evapotranspiration given temperature,
*       sunshine percentage and latitude
*
*       CALCULATION OF NET DOWNWARD SHORT-WAVE RADIATION FLUX
*       (also known as insolation or incident solar radiation)
*       Refs: Prentice et al 1993, Monteith & Unsworth 1990,
*             Henderson-Sellers & Robinson 1886
*
*        (1) rs = (c + d*ni) * (1 - beta) * Qo * cos Z * k
*              (Eqn 7, Prentice et al 1993)
*        (2) Qo = Qoo * ( 1 + 2*0.01675 * cos ( 2*pi*i / 365) )
*              (Eqn 8, Prentice et al 1993; angle in radians)
*        (3) cos Z = sin(lat) * sin(delta) + cos(lat) * cos(delta) * cos h
*              (Eqn 9, Prentice et al 1993)
*        (4) delta = -23.4 * pi / 180 * cos ( 2*pi*(i+10) / 365 )
*              (Eqn 10, Prentice et al 1993, angle in radians)
*        (5) h = 2 * pi * t / 24 = pi * t / 12
*
*            where rs    = instantaneous net downward shortwave radiation
*                          flux (W/m2 = J/m2/s)
*                  c, d  = empirical constants (c+d = clear sky
*                          transmissivity)
*                  ni    = proportion of bright sunshine
*                  beta  = average 'global' value for shortwave albedo
*                         (not associated with any particular vegetation)
*                  i     = julian date, (1-365, 1=1 Jan)
*                  Qoo   = solar constant, 1360 W/m2
*                 Z     = solar zenith angle (angular distance between the
*                          sun's rays and the local vertical)
*                  k     = conversion factor from solar angular units to
*                          seconds, 12 / pi * 3600
*                  lat   = latitude (+=N, -=S, in radians)
*                  delta = solar declination (angle between the orbital
*                         plane and the Earth's equatorial plane) varying
*                          between +23.4 degrees in northern hemisphere
*                          midsummer and -23.4 degrees in N hemisphere
*                          midwinter
*                  h     = hour angle, the fraction of 2*pi (radians) which
*                          the earth has turned since the local solar noon
*                  t     = local time in hours from solar noon
*
*       From (1) and (3), shortwave radiation flux at any hour during the
*       day, any day of the year and any latitude given by
*        (6) rs = (c + d*ni) * (1 - beta) * Qo * ( sin(lat) * sin(delta) +
*                 cos(lat) * cos(delta) * cos h ) * k
*       Solar zenith angle equal to -pi/2 (radians) at sunrise and pi/2 at
*       sunset.  For Z=pi/2 or Z=-pi/2,
*        (7) cos Z = 0
*       From (3) and (7),
*        (8)  cos hh = - sin(lat) * sin(delta) / ( cos(lat) * cos(delta) )
*             where hh = half-day length in angular units
*       Define
*        (9) u = sin(lat) * sin(delta)
*       (10) v = cos(lat) * cos(delta)
*       Thus
*       (11) hh = acos (-u/v)
*       To obtain the net DAILY downward short-wave radiation flux, integrate
*       equation (6) from -hh to hh with respect to h,
*       (12) rs_day = 2 * (c + d*ni) * (1 - beta) * Qo *
*                     ( u*hh + v*sin(hh) ) * k
*       Define
*       (13) w = (c + d*ni) * (1 - beta) * Qo
*       From (12) & (13),
*         (14) rs_day = 2 * w * ( u*hh + v*sin(hh) ) * k
*/

#define gamma(temp)  (65.05+temp*0.064)
#define lambda(temp) (2.495e6-temp*2380)

void petpar(Real *daylength, /* daylength (h) */
            Real *par,       /* photosynthetic active radiation flux */
            Real *pet,       /* potential evapotranspiration */
            Real lat,        /* latitude (deg) */
            Real day,        /* day (1..365) */
            Real temp,       /* temperature (deg C) */
            Real alt,          /*altitude (m) */
            Real temp_max,
            Real temp_min,
            Real wind,
            Real rhumid,   /* relative humid (%) */
            Real sun)       /* sunshine (%) */      
{
        temp_max = 30.0;
        temp_min = 10.0;
        wind = 2.0;
        rhumid = 0.3;       
  Real delta,u,v,hh,w,s,es,ea,ra,xxx,rnl,rns,rn,decli,psy,tav,etx,etn,ws,dr,slope,z;
  delta=deg2rad(-23.4*cos(2*M_PI*(day+10.0)/NDAYYEAR));
  u=sin(deg2rad(lat))*sin(delta);
  v=cos(deg2rad(lat))*cos(delta);
  dr=1+0.033*cos(0.0172*day);
  w=(c+d*sun)*(1-beta)*qoo*dr;
  decli= 0.409*sin(2*M_PI*day/NDAYYEAR-1.39);
  tav=(temp_max+temp_min)/2; //%=temp
  slope=(2504*exp((17.27*tav)/(tav+237.3)))/((tav+237.3)*(tav+237.3)); //^2);  %kPa  =s
  etx=0.6108*exp((17.27*temp_max)/(temp_max+237.3)); //% kPa
  etn=0.6108*exp((17.27*temp_min)/(temp_min+237.3));
  es=(etx+etn)/2.0;
  ea=es*30/100;
  z=pow( 293-0.0065*500/293, 5.26 );
  psy=0.000665*101.3*z; //% kPa
  sun*=M_PI/24/hh;
  ws=acos(-tan(lat)*tan(delta));   // %ws=hh
  xxx=ws*sin(lat)*sin(decli)+cos(lat)*cos(decli)*sin(ws); //% xxx=xx [x]
  ra=24*60/M_PI*0.082*dr*xxx;     //%unit MJ/m2/day                       
  rns=(1-0.23)*(0.2+0.79*sun)*(0.75+2e-5*500)*ra;
  rnl=(2.45e-3/24/3600)*(0.1+0.9*sun)*(0.56-0.25*sqrt(ea))*(temp_max*temp_max*temp_max*temp_max+temp_min*temp_min*temp_min*temp_min);
  rn=rns-rnl;
  if(u>=v)
  {
    *daylength=24;
    *par=rn/2;
  }
  else if(u<=-v)
   {
            *daylength=*par=0;
  }
  else
  {
    hh=acos(-u/v);   
    *par=rn/2;
    *daylength=24*hh*M_1_PI;
  }
  u=w*u-(b+(1-b)*sun)*(a-temp);
  v*=w;
  if(u<=-v)                /*polar night*/
  {
            *pet=0;
    }
  else
  {
    s=2.503e6*exp(17.269*temp/(237.3+temp))/
            ((237.3+temp)*(237.3+temp));
    if(u>=v)
     {
         *pet=((0.408*slope*rn+900*psy*3*0.75*(es-ea)/(temp+273))/(slope+psy*(1+0.34*3*0.75)))*28.4;  //% mm/day to W/m2
    }
    else
    {
      hh=acos(-u/v);
      *pet=((0.408*slope*rn+900*psy*3*0.75*(es-ea)/(temp+273))/(slope+psy*(1+0.34*3*0.75)))*28.4;  //% mm/day to W/m2
    }
  }
} /* of 'petpar' */

空谷幽竹

编译运行不就知道哪出问题了么

才仔

完全不懂额

sarcr

这是部分的程序，你具体指的是这里面的运算错误？编译错误？还是程序是做什么用的？
只贴出部分程序，调试起来会比较麻烦的。像Real这样的还要typedef，deg2rad这样的还要自己写然后嵌进去。
表达清楚大概是哪部分出了问题，纠错起来会比较方便~

倔强的葱头

空谷幽竹 发表于 2013-6-24 10:47

                               
登录/注册后可看大图

编译运行不就知道哪出问题了么

编译起来是没有问题的，但是整体运行起来就出现了错误，显示“n a n”

倔强的葱头

sarcr 发表于 2013-6-24 13:48

                               
登录/注册后可看大图

这是部分的程序，你具体指的是这里面的运算错误？编译错误？还是程序是做什么用的？
只贴出部分程序，调试 ...

您说的很对，我在做模型，这是指模型中的一个文件，我对这个文件进行修改后，编译可以通过，但是运行起来就是老是显示“n a n”，我也不知道哪里错了，会不会是逻辑错误？

空谷幽竹

nan 是not a number 的意思 比如除以零的时候容易出现 这个需要仔细找啊

liutaoz21

前来围观~~~~~~

倔强的葱头

空谷幽竹 发表于 2013-6-24 22:22

                               
登录/注册后可看大图

nan 是not a number 的意思 比如除以零的时候容易出现 这个需要仔细找啊

多谢啊~！！！！！