背景风场与台风风场插值耦合提供给swan程序

limingxing

1. 制作tfparameter.dat
收集到台风信息后，用坐标转换工具LonLat2xy.for将坐标转换为xy坐标
台风参数计算.xls 表中输入相应的数据
将获得的结果制成tfparameter.dat格式
2. 提取NCEP风场数据
(1)编辑经纬度.txt，格式为：
经度起始
经度结束
纬度起始
纬度结束
(2)用readnecpwind2.f90读取uv.yearmon.bln中的风场文件
结果为yearmon.dat
需要对此文件进行的处理是：截取首、尾多余时刻的风场数据
（NCEP每天的数据是以世界时6,12,18,0排序的；ECMWF是以0,6,12,18排序的）
注意下载的台风数据是地方时，-8为世界时
3. 坐标转换
需要制作lonlat.dat文件，包含风场数据的经纬度信息，格式为
经度 纬度
(如果各场台风提取的风场网格一样，lonlat.dat文件和xy.dat文件各场台风都一样，可略过此步)
利用LonLat2xy.for计算出xy.dat为xy坐标（软件运行选择0-北京54坐标）
4. 插值和耦合等
需要tfparameter.dat
截取后的yearmon.dat
xy.dat
应用combineprogram.f90（应用中需要自己修改一些参数，若整个风场网格一定，各场台风的这些参数是一致的）
可获得插值耦合结果：spaceinterresult.dat
      耦合风场：windpressurefield.dat      
      风和气压数据：wind.dat
      用于swan计算格式的风数据：swanwind.dat
插值耦合程序如下：
!将typhooninterpolation.f90，Typhoon field calculation inter.f90和creatwindforswan.f90合到一起
program combinedate
implicit none
!===============================需设定的部分========================================================================
real, parameter ::x00=126.8026775,y00=194.5496492                   !x00：风场网格原点x坐标(km)；y00：风场网格原点y坐标(km)
integer, parameter ::num_notes_x=656,num_notes_y=401  !num_notes_x：风场网格x向网格结点数；num_notes_y：风场网格y向网格结点数
real, parameter ::dis_unitx=0.2,dis_unity=0.2     !x00/(num_notes_x-1);y00/(num_notes_y-1)
!===================================================================================================================
integer i,j,n1,n2,num_input,inp
real lon,lat
real,allocatable ::u(:),v(:)
double precision, allocatable ::x(:),y(:)
character file_lluv*10,date*6
inp=num_notes_x*num_notes_y


! write(*,*)"输入读取数据的年月(如200310)"
! read(*,*)date
! file_lluv=date//'.dat'
!open(5,file=file_lluv,status="old",action='read')
open(6,file='xy.dat',status="old",action='read')
open(7,file='xyuv.dat',status="old",action='read')
!open(8,file='uv.dat',status="replace",action='write')
!open(9,file='ll.dat',status="replace",action='write')   
   n1=0
   do while(.not.eof(7))
     read(7,*)
         n1=n1+1
   enddo         
   rewind(7)
   n2=0
   do while(.not.eof(6))
     read(6,*)
         n2=n2+1
   enddo         
   rewind(6)
num_input=n1/n2
write(*,*)"num_input=",num_input
allocate (u(n2),v(n2),x(n2),y(n2))
do i=1,n2
read(6,*)x(i),y(i)
end do
!do i=1,num_input
!   do j=1,n2
!   read(5,*)lon,lat,u(j),v(j)
!   write(7,*)x(j),y(j),u(j),v(j)
!   end do
!end do
close(5)
close(6)
close(7)
call spaceinterplotation(num_input,x00,y00,num_notes_x,num_notes_y,inp,dis_unitx,dis_unity)
call Typhoon_field_calculation(num_input,x00,y00,num_notes_x,num_notes_y,inp,dis_unitx,dis_unity)
call output_tecplotANDswan(num_input,num_notes_x,num_notes_y,inp)

end


!===================================================================
!对NCEP数据进行插值
subroutine spaceinterplotation(num_input,x00,y00,num_notes_x,num_notes_y,inp,dis_unitx,dis_unity)
integer np,nx,ny,nt,nrec
integer ine
integer num1,num2,num3
character cha*1,file_grid*15,string*100,cgrid*5
integer i,j,k,numtime
real ri,rj,rk,rdm,xp,yp
real,allocatable :: x(:),y(:),u(:),v(:)
real,allocatable :: px(:),py(:),puinter(:),pvinter(:),pvelinter(:)
real,allocatable :: dm1(:),dm2(:),dm3(:),w1(:),w2(:),w3(:)
integer,allocatable :: lnodes(:,:),ndm(:,:)
character fileout*30,file_ext*3,cha_tmp*1

!--read Wind grid data file
write(*,*)'Wind grid data files imported……'
open(11,file='xy.dat',status="old",action='read')     !QuikSCAT风场网格
   np=0
   do while(.not.eof(11))
     read(11,*)
         np=np+1
   enddo
   write(*,*)np
   allocate(x(np),y(np),u(np),v(np))
   rewind(11)
   do i=1,np
     read(11,*)x(i),y(i)
   enddo
   close(11)

!-------------------------------------------------------------

! open(14,file='huangbodonghai_wind.dat',status='old',action='read')  !由台风模型生成的风场网格，使用Typhoon field calculation.f90

allocate(px(inp),py(inp),puinter(inp),pvinter(inp),pvelinter(inp))
do i=1,inp
   px(i)= (x00 + dis_unitx * (mod(i-1,num_notes_x)))*1000
   py(i)= (y00 + dis_unity * ((i-1)/num_notes_x))*1000
   write(35,*) px(i),py(i)
enddo
! rewind(14)

!--compute the weight
!--compute three nearest nodes
  write(*,*)'compute the weight……'
  allocate(ndm(inp,3))
  allocate(dm1(inp),dm2(inp),dm3(inp))
  allocate(w1(inp),w2(inp),w3(inp))
  !--start compute the weight--
  do i=1,inp
    dm1(i)=10.e12
    dm2(i)=10.e12
    dm3(i)=10.e12
        do j=1,np
          rdm=sqrt((x(j)-px(i))**2+(y(j)-py(i))**2)
          if(rdm<dm1(i))then
            dm3(i)=dm2(i)
                ndm(i,3)=ndm(i,2)
                dm2(i)=dm1(i)
                ndm(i,2)=ndm(i,1)
                dm1(i)=rdm
                ndm(i,1)=j
          elseif(rdm<dm2(i))then
            dm3(i)=dm2(i)
                ndm(i,3)=ndm(i,2)
                dm2(i)=rdm
                ndm(i,2)=j
          elseif(rdm<dm3(i))then
            dm3(i)=rdm
                ndm(i,3)=j
          else
          endif
        enddo
  enddo
  do i=1,inp
    if(dm1(i)<=0.00001)dm1(i)=0.00001
        if(dm2(i)<=0.00001)dm2(i)=0.00001
        if(dm3(i)<=0.00001)dm3(i)=0.00001
    w1(i)=1./dm1(i)
        w2(i)=1./dm2(i)
        w3(i)=1./dm3(i)
  enddo

!--read the swan result file
open(101,file='xyuv.dat',status='old',action='read')
open(119,file='spaceinterresult.dat',status='replace',action='write')
! open(4,file='windpressurefield.dat',status='unknown')
! write(4,*)'TITLE= "fjhd"'
! write(4,*)'VARIABLES = "X", "Y", "U", "V", "Vel"'
! open(201,file='wavedata.dat',status='replace',action='write')
! open(105,file='tecplot_adcirc_wave.dat',status='replace',action='write')

!  write(105,*)'TITLE= "SWAN TO ADCIRC"'
!  write(105,*)'VARIABLES = "X.", "Y.", "HS"'
!do i = 1, 12
!   write(*,*)'step:   ',i,'interplation……'
!   do j=1,np
!     read(101,*)
!   enddo
!enddo  
do i=1,num_input
!    WRITE(4,110)i,num_notes_x,num_notes_y
  !--read information from the result file
!  write(*,*)'step:   ',i,'interplation……'
!    do i=1,inp
!      read(14,*)xp,yp,pu(i),pv(i),pvel(i),ppr(i)
!         if(hs(j)<0.0)hs(j)=0.0
!         if(rtp(j)<0.0)rtp(j)=0.0
!    enddo
    do j=1,np
      read(101,*)xp,yp,u(j),v(j)
!         if(hs(j)<0.0)hs(j)=0.0
!         if(rtp(j)<0.0)rtp(j)=0.0
    enddo
  !--interplation   将QuikSCAT数据插值到风场模型生成的网格上，并取二者的较大值作为输出结果
   do k=1,inp
      !-u
            puinter(k)=w1(k)*u(ndm(k,1))+w2(k)*u(ndm(k,2))+w3(k)*u(ndm(k,3))
                puinter(k)=puinter(k)/(w1(k)+w2(k)+w3(k))
          !-v
                pvinter(k)=w1(k)*v(ndm(k,1))+w2(k)*v(ndm(k,2))+w3(k)*v(ndm(k,3))
                pvinter(k)=pvinter(k)/(w1(k)+w2(k)+w3(k))
            pvelinter(k)   =sqrt(puinter(k)*puinter(k)+pvinter(k)*pvinter(k))
              !if(pvelinter(k).gt.pvel(k))then
                   !write(45,*)px(k),py(k),pvelinter(k)-pvel(k)
                   !pu(k)=puinter(k)
                   !pv(k)=pvinter(k)
                   !pvel(k)=pvelinter(k)
                !endif
                write(119,*)px(k),py(k),puinter(k),pvinter(k),pvelinter(k)
!                write(4,*)px(k),py(k),puinter(k),pvinter(k),pvelinter(k)
   enddo
enddo

! close(4)
close(14)
close(101)
close(119)
110        FORMAT('ZONE T="P_',I4,'", I=',I8,',J=',I8,', F=POINT')

end subroutine

!================================================================================================
!台风数据与NCEP背景风场耦合
subroutine Typhoon_field_calculation(num_input,x00,y00,num_notes_x,num_notes_y,inp,dis_unitx,dis_unity)
real, parameter ::pi=3.1415926,g=9.8
real, parameter ::pa=100800.0,rt=40.000,c1=1.0,den_a=1.225,angle=20.0/180.0*3.1415926
real pr,u,v,r,x,y,f,delt_p,xmax,ymax,c2,c3,c4
real,allocatable,dimension(:)::lat,pr_center,dir_wind,speed_typhoon
real,allocatable,dimension(:)::x_center,y_center
integer i,j,k
real px,py,puinter,pvinter,pvelinter,ppr,uinter,vinter
allocate(lat(num_input),pr_center(num_input),dir_wind(num_input),speed_typhoon(num_input))
allocate(x_center(num_input),y_center(num_input))

!读入数据
open(1,file='tfparameter.dat',status='old')
do i=1,num_input
        read(1,*)j,x_center(i),y_center(i),lat(i),pr_center(i),dir_wind(i),speed_typhoon(i)
enddo
close(1)

!将台风中心气压单位由毫巴(mbar)转换为帕(pa)
do i=1,num_input
        pr_center(i)=pr_center(i)*100
enddo
!将台风移向的角度由角度值转换为弧度值
do i=1,num_input
        dir_wind(i)=dir_wind(i)*pi/180
enddo
!将台风移动速度由km/hr转化为m/s
do i=1,num_input    !
        speed_typhoon(i)=speed_typhoon(i)*1000/3600
enddo

open(2,file='wind.dat',status='unknown')       !同时修改下面子程序中相应位置处的文件名
!计算风速x方向分量，风速y方向分量和气压值,并将结果写入输入文件。

  open(120,file='spaceinterresult.dat',status='old')


do k=1,num_input
        f=2*0.00007272*sin(lat(k)*pi/180)                       !0.00007272rad/s为地球自转角速度
        delt_p=pa-pr_center(k)
        do i=1,num_notes_y
                y=dis_unity*(i-1)+y00
                do j=1,num_notes_x
                        x=dis_unitx*(j-1)+x00                                                                       
                        r=sqrt((x-x_center(k)/1000)**2+(y-y_center(k)/1000)**2)      
                        c2=0.8+(0.8*10**(0.0231*delt_p/100-1.96))*(2*r/rt)**1.5*exp(0.6*(1-(2*r/rt)**2.5))
                    read(120,*)px,py,puinter,pvinter,pvelinter
                        if(r.ge.0.and.r.le.2*rt) then       
                                pr=delt_p*(1-1/sqrt(1+2*(r/rt)**2))+pr_center(k)
                                u=c1*sin(dir_wind(k))*speed_typhoon(k)*exp(-pi*abs(r-rt)/4/rt)-c2*(-f/2+sqrt(f**2/4+2*delt_p*(1+2*(r/rt)**2)**(-3/2)/den_a/1000000/rt**2))*((x*1000-x_center(k))*sin(angle)+(y*1000-y_center(k))*cos(angle))
                                !这里的rt的单位是km，式中除以1000000(关于rt*rt项)，是将表达式中的单位统一成国际制

                               
                                v=c1*cos(dir_wind(k))*speed_typhoon(k)*exp(-pi*abs(r-rt)/4/rt)-c2*(-f/2+sqrt(f**2/4+2*delt_p*(1+2*(r/rt)**2)**(-3/2)/den_a/1000000/rt**2))*(-(x*1000-x_center(k))*cos(angle)+(y*1000-y_center(k))*sin(angle))
                    write(35,*)(-f/2+sqrt(f**2/4+2*delt_p*(1+2*(r/rt)**2)**(-3/2)/den_a/1000000/rt**2))*((x*1000-x_center(k))*sin(angle)+(y*1000-y_center(k))*cos(angle)),(-f/2+sqrt(f**2/4+2*delt_p*(1+2*(r/rt)**2)**(-3/2)/den_a/1000000/rt**2))*(-(x*1000-x_center(k))*cos(angle)+(y*1000-y_center(k))*sin(angle))

                        else
                                pr=delt_p*(1-1/sqrt(1+r/rt))+pr_center(k)       
                                u=c1*sin(dir_wind(k))*speed_typhoon(k)*exp(-pi*abs(r-rt)/4/rt)-c2*(-f/2+sqrt(f**2/4+delt_p/den_a/r/rt/1000000/(1+r/rt)**2))*((x*1000-x_center(k))*sin(angle)+(y*1000-y_center(k))*cos(angle))
                                v=c1*cos(dir_wind(k))*speed_typhoon(k)*exp(-pi*abs(r-rt)/4/rt)-c2*(-f/2+sqrt(f**2/4+delt_p/den_a/r/rt/1000000/(1+r/rt)**2))*(-(x*1000-x_center(k))*cos(angle)+(y*1000-y_center(k))*sin(angle))                               
                        end if
                     c3=r/(9*rt)
                         c4=c3**4/(1+c3**4)
                         uinter=u*(1-c4)+c4*puinter
                         vinter=v*(1-c4)+c4*pvinter
                          write(2,*)x*1000,y*1000,uinter,vinter,pr
                                                 

                enddo
        enddo
enddo

close(2)
close(120)


end subroutine
!*************************************
!Tecplot画图程序和swan风场格式数据输出
subroutine  output_tecplotANDswan(num_input,num_notes_x,num_notes_y,inp)
integer i,j,n,k
real data_output(inp,5)

open(3,file='wind.dat',status='unknown')
open(4,file='windpressurefield.dat',status='unknown')
open(8,file='swanwind.dat',status='unknown')
!写TECPLOT的总HEADER
write(4,*)'TITLE= "fjhd"'
write(4,*)'VARIABLES = "X", "Y", "U", "V", "PR"'
do i=1,num_input
                WRITE(4,11)i,num_notes_x,num_notes_y
                do j=1,inp
                        read(3,*)(data_output(j,n),n=1,5)
                        write(4,*)(data_output(j,n),n=1,5)
                enddo
write(8,*)(data_output(k,3),k=1,inp)
write(8,*)(data_output(k,4),k=1,inp)
enddo
close(3)
close(4)
close(8)
11        FORMAT('ZONE T="P_',I4,'", I=',I8,',J=',I8,', F=POINT')
end subroutine

栗子

感谢分享

mofangbao

多谢分享啦

lasa123

感谢感谢分享

limingxing

下面回复的是什么情况  能看懂么

limingxing

limingxing 发表于 2015-12-15 17:38
下面回复的是什么情况  能看懂么

哈哈 做过类似工作应该可以看懂的

barry

谢谢分享谢谢分享

txy

厉害，谢谢楼主分享

zdcgl

谢谢楼主分享

plizi

特别感谢 一直在找这方面的内容 不知道有没有matlab 做过类似工作的 谢谢