校内生产实习南方气旋天气过程数据处理.f90文件

liewenjide

program shixi1
implicit none
integer i,j,it,ifile,idv,idy,idh,idz
integer,parameter::nv=5,nx=111,ny=61,nz=11,nt=4,a=6371*1000
integer::height(11)=(/00,92,85,70,50,40,30,25,20,15,10/)
real,parameter::pi=3.1425926,d=1.0/180.0*pi,g=9.8
real lat(ny)
real,dimension(nv,nx,ny,nz,nt)::vars
real,dimension(nx,ny,nz,nt)::vu,vv,vh,vq,vt,vtt,tp,wo,wp,xsqtl,ysqtl,sqrtsan
character u,v,h,q,t
character::var(5)=(/'u','v','h','q','t'/)
character (len=12) filename
vars(:,:,:,:,:)=0
write(filename(1:4),'(A4)')'0004'
write(filename(9:9),'(A1)')'.'

ifile=0
do idv=1,5
write(filename(12:12),'(A1)')var(idv)
it=0
do idy=13,14
write(filename(5:6),'(I2)')idy
  do idh=00,12,12
  write(filename(7:8),'(I2)')idh
  if(idh==0) then
  write(filename(7:8),'(A2)')'00'
  endif
  it=it+1
   do idz=1,nz
   write(filename(10:11),'(I2)')height(idz)
   if(idz==1) then
   write(filename(10:11),'(A2)')'00'
   endif
   if(idv==4.and.(idz==8.or.idz==9.or.idz==10.or.idz==11)) cycle
   ifile=ifile+1
   open(ifile,file='e:\2014shixi\data\'//trim(adjustl(filename))//'',status='old',form='formatted')   
   read(ifile,*)      !(读出第一行资料格式说明)
   read(ifile,*)((vars(idv,i,j,idz,it),i=1,nx),j=1,ny)
   100    format(10f8.2)
    close(ifile)
   enddo
  enddo
enddo
enddo
do it=1,nt
  do idz=1,nz
   do j=1,ny
    do i=1,nx
    vu(i,j,idz,it)=vars(1,i,j,idz,it)
    vv(i,j,idz,it)=vars(2,i,j,idz,it)
    vh(i,j,idz,it)=vars(3,i,j,idz,it)
    vq(i,j,idz,it)=vars(4,i,j,idz,it)
    vt(i,j,idz,it)=vars(5,i,j,idz,it)
    vtt(i,j,idz,it)= vt(i,j,idz,it)+273.15
    enddo
   enddo
  enddo
enddo

!计算温度平流wp、涡度wo、涡度平流wp、水汽通量散度sqrtsan
do it=1,nt
do idz=1,nz
do j=2,ny-1
  lat(j)=10.0+float(j-2)*1.0
  lat(j)=lat(j)/180.0*pi
do i=2,nx-1
tp(i,j,idz,it)=vu(i,j,idz,it)*(vtt(i+1,j,idz,it)-vtt(i-1,j,idz,it))/(2.0*a*cos(lat(j))*d)+vv(i,j,idz,it)*(vtt(i,j+1,idz,it)-vtt(i,j-1,idz,it))/(2.0*a*d)
wo(i,j,idz,it)=1.0/(2.0*a)*((vv(i+1,j,idz,it)-vv(i-1,j,idz,it))/(cos(lat(j))*d)-(vu(i,j+1,idz,it)-vu(i,j-1,idz,it))/d+2.0*vu(i,j,idz,it)*tan(lat(j)))
wp(i,j,idz,it)=vu(i,j,idz,it)*(wo(i+1,j,idz,it)-wo(i-1,j,idz,it))/(2.0*a*cos(lat(j))*d)+vv(i,j,idz,it)*(wo(i,j+1,idz,it)-wo(i,j-1,idz,it))/(2.0*a*d)
sqrtsan(i,j,idz,it)=1.0/g*((vu(i+1,j,idz,it)*vq(i+1,j,idz,it)-vu(i-1,j,idz,it)*vq(i-1,j,idz,it))/(2.0*a*cos(lat(j))*d)+(vv(i,j+1,idz,it)*vq(i,j+1,idz,it)-vv(i,j-1,idz,it)*vq(i,j-1,idz,it))/(2.0*a*d)+vq(i,j,idz,it)*((vu(i+1,j,idz,it)-vu(i-1,j,idz,it))/(2.0*a*cos(lat(j))*d)+(vv(i,j+1,idz,it)-vv(i,j-1,idz,it))/(2.0*a*d)))
enddo
enddo
enddo
enddo
do it=1,nt
do idz=1,nz
tp(1,1,idz,it)=tp(2,2,idz,it)
tp(1,ny,idz,it)=tp(2,ny-1,idz,it)
tp(nx,1,idz,it)=tp(nx-1,2,idz,it)
tp(nx,ny,idz,it)=tp(nx-1,ny-1,idz,it)
wp(1,1,idz,it)=wp(2,2,idz,it)
wp(1,ny,idz,it)=wp(2,ny-1,idz,it)
wp(nx,1,idz,it)=wp(nx-1,2,idz,it)
wp(nx,ny,idz,it)=wp(nx-1,ny-1,idz,it)
sqrtsan(1,1,idz,it)=sqrtsan(2,2,idz,it)
sqrtsan(1,ny,idz,it)=sqrtsan(2,ny-1,idz,it)
sqrtsan(nx,1,idz,it)=sqrtsan(nx-1,2,idz,it)
sqrtsan(nx,ny,idz,it)=sqrtsan(nx-1,ny-1,idz,it)
enddo
enddo
do it=1,nt
do idz=1,nz
do j=2,ny-1
tp(1,j,idz,it)=tp(2,j,idz,it)
tp(nx,j,idz,it)=tp(nx-1,j,idz,it)
wp(1,j,idz,it)=wp(2,j,idz,it)
wp(nx,j,idz,it)=wp(nx-1,j,idz,it)
sqrtsan(1,j,idz,it)=sqrtsan(2,j,idz,it)
sqrtsan(nx,j,idz,it)=sqrtsan(nx-1,j,idz,it)
enddo
enddo
enddo
do it=1,nt
do idz=1,nz
do i=2,nx-1
tp(i,1,idz,it)=tp(i,2,idz,it)
tp(i,ny,idz,it)=tp(i,ny-1,idz,it)
wp(i,1,idz,it)=wp(i,2,idz,it)
wp(i,ny,idz,it)=wp(i,ny-1,idz,it)
sqrtsan(i,1,idz,it)=sqrtsan(i,2,idz,it)
sqrtsan(i,ny,idz,it)=sqrtsan(i,ny-1,idz,it)
enddo
enddo
enddo
do it=1,nt
do idz=1,nz
do j=1,ny
do i=1,nx
tp(i,j,idz,it)= tp(i,j,idz,it)*10**5
wp(i,j,idz,it)= wp(i,j,idz,it)*10**10
sqrtsan(i,j,idz,it)= sqrtsan(i,j,idz,it)*10**7
enddo
enddo
enddo
enddo
!print*,tp(:,:,1,1)
!print*,wp(:,:,1,1)
print*,sqrtsan(:,:,1,1)

!计算水汽通量sqtl（xsqtl和ysqtl）
do it=1,nt
do idz=1,nz
do j=1,ny
do i=1,nx
xsqtl(i,j,idz,it)=1.0/g*vq(i,j,idz,it)*vu(i,j,idz,it)
ysqtl(i,j,idz,it)=1.0/g*vq(i,j,idz,it)*vv(i,j,idz,it)
enddo
enddo
enddo
enddo

!储存数据
open(205,file='e:\2014shixi\data\datas.grd',form='binary')
open(206,file='e:\2014shixi\data\datas.txt')

do it=1,nt
write(205)(((vu(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((vv(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((vh(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((vq(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((vt(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((tp(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((wp(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((xsqtl(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((ysqtl(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(205)(((sqrtsan(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)

write(206,100)(((vu(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((vv(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((vh(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((vq(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((vt(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((tp(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((wp(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((xsqtl(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((ysqtl(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
write(206,100)(((sqrtsan(i,j,idz,it),i=1,nx),j=1,ny),idz=1,nz)
enddo

end program shixi1