计算子流域面雨量：黄河mask边界nc文件，和全国的多年月降水nc文件，编译有问题！

mzhang

      program mPRCP
      use netcdf
      integer,parameter::clon1=105,clat1=65,clon2= 280,clat2=160,&
      year=64,bmon=12,dims=3,dims1=2,m=12
      character*255 ext1,ext2,ext3
      integer::ncid1,ncid2,ncid3,varid1,varid2,varid3,varid(4),&
               lon_dimid, lat_dimid,bmon_dimid,year_dimid,mask_dimid
      integer::dimids(dims),dimids1(dims1),start1(3),con1(3)
      real::lon1(clon1),lat1(clat1),lon2(clon2),lat2(clat2)
      character*4 cyr
      real::mask(clon1,clat1)
      real::q(m)=(/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0/)
      real::PRCPm(clon2,clat2,bmon), var1(bmon,year),&
      var2(m,bmon,year)
      integer::i,j,ii,jj,k,k1,k2,cont,status
      ext1='/mnt/lustre/zhang/pub/'
      call check(nf90_open(trim(ext1)//'Yellow_0.25_VIC.nc',&
           NF90_NOWRITE,ncid1))
      call check(nf90_inq_varid(ncid1,'lon',varid1))
      call check(nf90_get_var(ncid1,varid1,lon1))
      call check(nf90_inq_varid(ncid1,'lat',varid2))
      call check(nf90_get_var(ncid1,varid2,lat1))
      call check(nf90_inq_varid(ncid1,'mask',varid3))
      call check(nf90_get_var(ncid1,varid3,mask))
      call check(nf90_close(ncid1))
      ext2='/mnt/lustre/zhang/'
      call check(nf90_create(trim(ext2)//'monthly PRCP19512014.nc',&
           nf90_clobber,ncid2))
      call check(nf90_def_dim(ncid2,'year',year,year_dimid))
      call check(nf90_def_dim(ncid2,'bmon',bmon,bmon_dimid))
      call check(nf90_def_dim(ncid2,'lon',clon1,lon_dimid))
      call check(nf90_def_dim(ncid2,'lat',clat1,lat_dimid))
      call check(nf90_def_dim(ncid2,'mask1',m,mask_dimid))
         dimids=(/mask_dimid,bmon_dimid,year_dimid/)
      call check(nf90_def_var(ncid2,'PRCP',NF90_REAL,dimids,varid(1)))
      call check(nf90_put_att(ncid2,varid(1),'long_name','monthly rain'))
      call check(nf90_put_att(ncid2,varid(1),'units','mm/day'))
      call check(nf90_put_att(ncid2,varid(1),'_FillValue',-999.))
      call check(nf90_put_att(ncid2,varid(1),'missing_value',-999.))
         dimids(1)=lon_dimid
      call check(nf90_def_var(ncid2,'lon',NF_90REAL,dimids(1),varid(2)))
      call check(nf90_put_att(ncid2,varid(2),'units','degrees_east'))
         dimids(1)=lat_dimid
      call check(nf90_def_var(ncid2,'lat',NF_90REAL,dimids(1),varid(3)))
      call check(nf90_put_att(ncid2,varid(3),'units','degrees_north'))
         dimids(1)=mask_dimid
      call check(nf90_def_var(ncid2,'mask1',NF_90REAL,dimids(1),varid(4)))
      call check(nf90_put_att(ncid2,varid(4),'_FillValue',-999.))
      call check(nf90_put_att(ncid2,varid(4),'missing_value',-999.))
      call check(nf90_enddef(ncid2))
      ext3='/mnt/lustre/zhang/monthly/'
      do k1=1951,2014,1
       write(cyr,'(i4)') k1
      call check(nf90_open(trim(ext3)//''//cyr//'m.nc',NF90_NOWRITE,&
           ncid3))
      call check(nf90_inq_varid(ncid3,'lon',varid1))
      call check(nf90_get_var(ncid3,varid1,lon2))
      call check(nf90_inq_varid(ncid3,'lat',varid2))
      call check(nf90_get_var(ncid3,varid2,lat2))
      call check(nf90_inq_varid(ncid3,'PRCP',varid3))
      start1=(/1,1,1/);con1=(/clon2,clat2,bmon/)
      call check(nf90_get_var(ncid3,varid3,PRCPm,start1,con1))

      do k2=1,bmon
      do k=1,m
      var2=-999.0
      cont=0
      var1=0
         do i=1,clon1
          do j=1,clat1
            if(abs(mask(i,j)-q(k))<=0.01)then
               do ii=1,clon2
                do jj=1,clat2
                    if(abs(lon2(ii)-lon1(i))<=0.01.and.&
                       abs(lat2(jj)-lat1(j))<=0.01.and.&
                       PRCPm(ii,jj,k2)>-900) then
                        cont=cont+1
                        var1(k2,k1)=var1(k2,k1)+PRCPm(ii,jj,k2)
                    end if
                       if(cont>0)then
                         var2(q(k),k2,k1)=var1(k2,k1)/cont
                       end if
                end do
               end do
            end if
         end do
        end do
      end do
      end do
       call check(nf90_close(ncid3))
      end do
     print*,'ok'
        call check(nf90_put_var(ncid2,varid(1),var2))
        call check(nf90_put_var(ncid2,varid(2),lon1))
        call check(nf90_put_var(ncid2,varid(3),lat1))
        call check(nf90_put_var(ncid2,varid(4),q))
        call check(nf90_close(ncid2))
!      end program

     contains
       subroutine check(status)
      integer, intent ( in) :: status
      if(status /= nf90_noerr) then
     print *, trim(nf90_strerror(status))
     stop "Stopped"
      end if
      end subroutine check