ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   本程序用于2001年气象数据空间插值
c            2006-9-28
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
  parameter(n_sta=35,n=365,nn=1,n_pix=584915)
c  n－－－一年的总天数  nn－－－需要插值的天数
  dimension elev_sta(n_sta),temp_sta(n,n_sta)
  real mxtemp_sta(n,n_sta),mntemp_sta(n,n_sta)
  dimension vpre_sta(n,n_sta),suntime_sta(n,n_sta)
  dimension wind_sta(n,n_sta),p_sta(n,n_sta)
   character*50  No_sta(n_sta),Name_sta(n_sta),Name_prov,temps

c  n_sta----------------------气象站站点数
c  NO_sta---------------------气象站号
c  elev_sta-------------------气象站高程
c  temp_sta-------------------气象站气温
c  mxtemp_sta,mntemp_sta------气象站最高、最低气温
c  vpre_sta-------------------气象站水气压
c  suntime_sta----------------气象站日照时数
c  wind_sta-------------------气象站风速
c      P_sta----------------------气象站降雨
c  Name_sta,Name_prov---------气象站站名和所在省名

  dimension cor_x(n_sta),cor_y(n_sta)
  integer year(n,n_sta),month(n,n_sta),day(n,n_sta)
c  cor_x,cor_y ---------气象站经纬度坐标
c     
  dimension x(n_pix),y(n_pix),dem(n_pix)
c  x，y---------网格投影坐标（单位：米）
c  dem----------网格高程
  integer ye(nn),mo(nn),da(nn)
c  插值日期
  dimension x0(n_sta),y0(n_sta)
c  x0，y0---------气象站投影坐标(Alberes,WGS84)

  real temperature(n_sta),mxtemp(n_sta),mntemp(n_sta)
  dimension vpre(n_sta),suntime(n_sta),wind(n_sta),p(n_sta)
  dimension Vpre_out(n_pix),wind_out(n_pix),suntime_out(n_pix)
  dimension  temperature_out(n_pix),P_out(n_pix)
  real lanbuda,Ning

c  读入气象站点各项数据
       open(1,file='2001.txt',status='old')
       read(1,*)
  do i=1,n_sta
    do j=1,n
       read(1,*)NO_sta(i),cor_x(i),cor_y(i),Name_sta(i),Name_prov,temp1,
     $ elev_sta(i),temp2,temps,temp3,temp4,temp5,temp6,temp7,temp8,
     $ year(j,i),month(j,i),day(j,i),temp9,temp_sta(j,i),
     $ mxtemp_sta(j,i),mntemp_sta(j,i),
     $ vpre_sta(j,i),temp10,temp11,wind_sta(j,i),p_sta(j,i),
     $ temp12,suntime_sta(j,i)
    enddo
  enddo
  close(1)
c do i=1,n_sta
c write(*,*)cor_x(i),cor_y(i)
c enddo

c      读入栅格坐标
  open(2,file='coordinate.txt',status='old')
  do i=1,n_pix
   read(2,*) x(i),y(i),dem(i)
  enddo
  close(2)

c  输入日期
  open(3,file='input_data.txt',status='old')
  do i=1,nn
  read(3,*)ye(i),mo(i),da(i)
  enddo
  close(3)
 
c open(211,file='wind.txt')
c open(333,file='relation.txt')
c  找出该天的所有气象站点的数据并进行插值
  do ii=1,1  !nn
   
    !挑出当天各站的气象数据
    do i=1,n_sta
      do j=1,n
    if(year(j,i).eq.ye(ii).and.month(j,i).eq.mo(ii).
     $   and.day(j,i).eq.da(ii))then
       temperature(i)=temp_sta(j,i)*0.1      !0.1--scale factor
       mxtemp(i)=mxtemp_sta(j,i)*0.1
    mntemp(i)=mntemp_sta(j,i)*0.1
       vpre(i)=vpre_sta(j,i)*0.1
    suntime(i)=suntime_sta(j,i)*0.1
    wind(i)=wind_sta(j,i)*0.1
c       降雨对前五天累计进行插值,要保证j>5，否则不行
    p(i)=0.0
    do kk=1,5
     if(P_sta(j-kk,i).gt.3000) P_sta(j-kk,i)=0.0
     p(i)=p(i)+P_sta(j-kk,i)*0.1
    enddo
         
       endif
      enddo

 write(*,*)temperature(i),vpre(i),suntime(i),P(i)

c write(333,*)elev_sta(i),suntime(i)*0.1

  ! 找出各气象站点的风速，不插值
      lanbuda=cor_x(i)
    fi=cor_y(i)
 write(*,*)lanbuda,fi,wind(i)
   call proj(lanbuda,fi,Eing,Ning)
      x0(i)=Eing
   y0(i)=Ning
 
c     write(211,212)lanbuda*180/3.141593,fi*180/3.141593,x0(i)
c     $,y0(i),wind(i),name_sta(i)
c212 format(2(1x,f10.4),2(1x,f10.1),3x,f4.1,2x,a20)
  enddo
c close(211)

 do i=1,n_pix
   wind_out(i)=-10.0
   do j=1,n_sta
  if(abs(x(i)-x0(j)).le.500.and.abs(y(i)-y0(j)).le.500) then
    wind_out(i)=wind(j)
   write(*,*) 'find',j,wind_out(i)
  endif
   enddo
 enddo

  !------------------------------------------------------------------
  ! 对网格点插值
     ! 水汽压插值，对高程有明显的趋势，随高程指数下降
  do i=1,n_sta
       vpre(i)=vpre(i)*10**(elev_sta(i)/5000.0)
  enddo

  call DDWA(x0,y0,vpre,n_sta,x,y,n_pix,vpre_out)

  do i=1,n_pix
    if(dem(i).eq.-32768)then
   vpre_out(i)=-999
       else 
      vpre_out(i)=vpre_out(i)*10**(-dem(i)/5000.0)
    endif
  enddo

  !------------------------------------------------------------------
  ! 风速插值，对比高程和风速没有明显的趋势
c  call DDWA(x0,y0,wind,n_sta,x,y,n_pix,wind_out)

  !-------------------------------------------------------------------
  ! 太阳日照时数插值, 对高程没有趋势
  call DDWA(x0,y0,suntime,n_sta,x,y,n_pix,suntime_out)

  !-------------------------------------------------------------------
  ! 气温插值，需要将气象站高度订正到同一平面再插值，然后网格再算回去
     do i=1,n_sta
       offset=elev_sta(i)/100.0*0.6   !海拔每上升100米，气温平均下降0.6度
    temperature(i)=temperature(i)+offset
     enddo

  call DDWA(x0,y0,temperature,n_sta,x,y,n_pix,temperature_out)
  
  do i=1,n_pix
    if(dem(i).eq.-32768)then
   temperature_out(i)=-999
       else 
   offset=dem(i)/100*0.6
      temperature_out(i)=temperature_out(i)-offset
    endif
  enddo
c
  !------------------------------------------------------------------------
   ! 降雨插值
  call DDWA(x0,y0,P,n_sta,x,y,n_pix,P_out)

  do i=1,n_pix
    if(dem(i).eq.-32768)then
   p_out(i)=-999
       else 
      p_out(i)=p_out(i)
    endif
  enddo
  !------------------------------------------------------------------------

     open(4,file='Vpre.txt')
     do i=1,n_pix
          write(4,111)x(i),y(i),Vpre_out(i)
111  format(1x,f10.1,1x,f10.1,1x,f10.3)
     enddo
     close(4)

     open(5,file='Wind.txt')
     do i=1,n_pix
          write(5,112)x(i),y(i),wind_out(i)
112  format(1x,f10.1,1x,f10.1,1x,f10.3)
     enddo
     close(5)

  open(6,file='suntime.txt')
     do i=1,n_pix
      if(dem(i).eq.-32768) then
    suntime_out(i)=-999
   else
    suntime_out(i)=suntime_out(i)
   endif
          write(6,113)x(i),y(i),suntime_out(i)
113  format(1x,f10.1,1x,f10.1,1x,f10.3)
     enddo
     close(6)

  open(7,file='Temperature.txt')
     do i=1,n_pix
          write(7,114)x(i),y(i),Temperature_out(i)
114  format(1x,f10.1,1x,f10.1,1x,f10.3)
     enddo
     close(7)

  open(8,file='precipitation.txt')
     do i=1,n_pix
          write(8,115)x(i),y(i),P_out(i)
115  format(1x,f10.1,1x,f10.1,1x,f10.3)
     enddo
     close(8)

  enddo

 end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  该子程序用于将点的经纬度坐标转化为Albers的投影坐标,用于将气象站的经纬度转化为米
c      经纬度为度的形式，非度分秒形式
c  取中国区的参数为：fi1-------第一标准纬度25
c                    fi2-------第二标准纬度47
c                    lanbuda0--中央经线105
c                        fi0-------中央纬线0
c                    EF--------东偏0,单位：米
c                        NF--------北偏0,单位：米
c      WGS84椭球体参数       长半轴a=6378137
c                            短半轴b=6356752.3142
c      Karasovsky椭球体参数  长半轴a=6378245
c                            短半轴b=6356863.0188
c  偏心率e=SQRT(2f-f^2),1/f=a/(a-b)                    
cc     这里采用WGS84

  subroutine proj(lanbuda,fi,Eing,Ning)
c   输入（地理经纬度）lanbuda------精度
c                         fi-----------纬度
c       输出（投影坐标m） Eing---------东的坐标
c                         Ning---------北的坐标
  real lanbuda,lanbuda0,NF,m1,m2,n,Ning 
  pi=4*atan(1.0)
       fi0=0*pi/180
  lanbuda0=105*pi/180
  fi1=25*pi/180
  fi2=47*pi/180
  EF=0
  NF=0

  fi=fi*pi/180
  lanbuda=lanbuda*pi/180
 
  a=6378137
  b=6356752.3142
       e1=SQRT(1-b/a*b/a)
  e2=1-e1*e1

  alfa=e2*(sin(fi)/(1-e1*e1*sin(fi)*sin(fi))-0.5/e1
     $     *log((1-e1*sin(fi))/(1+e1*sin(fi))))
       alfa0=e2*(sin(fi0)/(1-e1*e1*sin(fi0)*sin(fi0))-0.5/e1
     $      *log((1-e1*sin(fi0))/(1+e1*sin(fi0))))
       alfa1=e2*(sin(fi1)/(1-e1*e1*sin(fi1)*sin(fi1))
     $      -0.5/e1*log((1-e1*sin(fi1))/(1+e1*sin(fi1))))
       alfa2=e2*(sin(fi2)/(1-e1*e1*sin(fi2)*sin(fi2))
     $      -0.5/e1*log((1-e1*sin(fi2))/(1+e1*sin(fi2))))
       m1=cos(fi1)/sqrt(1-e1*e1*sin(fi1)*sin(fi1))
       m2=cos(fi2)/sqrt(1-e1*e1*sin(fi2)*sin(fi2))
       n=(m1*m1-m2*m2)/(alfa2-alfa1)
       C=m1*m1+n*alfa1
       rou0=a*sqrt(C-n*alfa0)/n
       rou=a*sqrt(C-n*alfa)/n
       sita=n*(lanbuda-lanbuda0)

  Eing=EF+rou*sin(sita)
       Ning=NF+rou0-rou*cos(sita)
  return
       end
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  该子程序采用距离方向加权平均法进行空间插值
c  x0,y0---------气象站点坐标,n0-----站数（站号）
c      x1,y1---------待求点坐标
c  P-------------气象站实测值
c      P-out---------点插值结果
c  dis(j)--------与第j个站的距离
cc
  subroutine DDWA(x0,y0,p,n0,x1,y1,n1,P_out)
  dimension x0(n0),y0(n0),x1(n1),y1(n1)
  dimension P(n0),P_out(n1)
  dimension dis(n0),num(n0),w(n0),alf(n0)

c  通过最近8个站计算每个目标点的值
  do i=1,n1
      do j=1,n0
        dis(j)=sqrt((x1(i)-x0(j))*(x1(i)-x0(j))
     $     +(y1(i)-y0(j))*(y1(i)-y0(j)))
   if(dis(j).eq.0) then
     P_out(i)=p(j)
     goto 33
   endif
      enddo

     
    !找出最近的8个站及其与待求点的距离
       do j=1,n0
        num(j)=j
       enddo

   do ii=1,8
         s=1000000000
        do j=ii+1,n0 
         if(dis(j).lt.s) then
    s=dis(j)
          k=j
         endif
        enddo
    if(s.lt.dis(ii))then
        dis(k)=dis(ii)
     dis(ii)=s
   tempk=num(ii)
   num(ii)=num(k)
   num(k)=tempk
     endif  
      enddo
    
       !计算权重距离
       dis0=dis(8)             !取8个最远的点的距离（经验衰减距离）
       do j=1,8
   w(j)=(exp(-dis(j)/dis0))**4.0
    enddo

    !计算方向修正系数
   do j=1,8
   s1=0
   s2=0
     do l=1,8
   s1=s1+w(l)
   d=SQRT((x0(num(j))-x0(num(l)))**2
     $  +(y0(num(j))-y0(num(l)))**2)
   cos=(dis(j)*dis(j)+dis(l)*dis(l)-d*d)/(2*dis(j)*dis(l))
   s2=s2+w(l)-w(l)*cos
        enddo
   alf(j)=s2/s1
   enddo

   !计算总修正权重
   do j=1,8
    w(j)=w(j)*(1+alf(j))
   enddo

   !计算目标点值
   s1=0
   s2=0
   do j=1,8
      s1=s1+w(j)*P(num(j))
    s2=s2+w(j)
   enddo
   P_out(i)=s1/s2
33  enddo
  return
  end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       该子程序采用修正距离倒数平方法
c  x0,y0---------气象站点坐标,n0-----站数（站号）
c      x1,y1---------待求点坐标
c  dem-----------待求点高程 
c  P-------------气象站实测值
c      P-out---------点插值结果
c  dis(j)--------与第j个站的距离
cc
  subroutine MIDW(x0,y0,elev_sta,p,n0,x1,y1,dem,n1,P_out)
  dimension x0(n0),y0(n0),elev_sta(n0),x1(n1),y1(n1),dem(n1)
  dimension P(n0),P_out(n1)
  dimension dis(n0),alf(n0)

  b=-2                   ! b是权重指数，一般情况下取负值
  do i=1,n1
 
  ss=0.0
     do j=1,n0
    de=abs(dem(i)-elev_sta(j))
    if(de.eq.0) de=1
c    de=1       !如果de＝1，则变为距离倒数平方法
       dis(j)=SQRT((x1(i)-x0(j))**2+(y1(i)-y0(j))**2)
       if(dis(j).eq.0)then
         P_out(i)=p(j)
   goto 440
     endif
    dis(j)=(dis(j)/de)**b
    ss=ss+dis(j)
   enddo
 
  do j=1,n0
      alf(j)=dis(j)/ss
 write(*,*)alf(j)
   p_out(i)=p_out(i)+alf(j)*p(j)
  enddo
440 enddo
 return
 end