|
PROGRAM FILTERS
|
|
implicit none
|
|
!
|
|
real ignore_value
|
|
parameter(ignore_value=-99999.0)
|
|
!
|
|
real, allocatable :: long(:),lat(:),zcoord(:),tcoord(:)
|
|
real, allocatable :: ice_lon(:,:),ice_lat(:,:),ice_lon_box(:,:),ice_lat_box(:,:),ice_lon_temp(:,:),ice_lat_temp(:,:)
|
|
!
|
|
character fname*180,fname1*180,oname*80
|
|
!
|
|
integer loc_dim,lnsig,ix,iy,iz,it,ifail,i,imon,ip,j
|
|
real alon_min,alon_max,alat_min,alat_max,dx,dy
|
|
logical miss_data1,lexist,lexist1,lexist2,lexist3
|
|
!
|
|
fname1='DP09_lonlat.nc'
|
|
inquire(file=trim(fname1),exist=lexist1)
|
|
!
|
|
oname='setgrid/DP09_lonlat_v3.dat'
|
|
!
|
|
if (lexist1) then
|
|
print*,' Opening file = '//trim(fname1)
|
|
else
|
|
print*,' Relevant files to create ice mask do not exist'
|
|
stop 1
|
|
end if
|
|
!
|
|
! Extract dimensions of field. This uses a subroutine get_dims
|
|
! which is in a library which is part of libancil suite
|
|
!
|
|
call get_dims(trim(fname1),'latitude',ix,iy,iz,it,ifail)
|
|
print '(a,6i5)',' Opened '//trim(fname1)//' ',ix,iy,iz,it,ifail
|
|
if (ifail.ne.0) then
|
|
stop 1
|
|
end if
|
|
!
|
|
! Allocate memory for arrays
|
|
!
|
|
iz=max(iz,1)
|
|
it=max(it,1)
|
|
!
|
|
allocate(long(ix))
|
|
allocate(lat(iy))
|
|
allocate(zcoord(iz))
|
|
allocate(tcoord(it))
|
|
!
|
|
allocate(ice_lon(ix,iy))
|
|
allocate(ice_lat(ix,iy))
|
|
allocate(ice_lon_temp(ix+1,iy))
|
|
allocate(ice_lat_temp(ix+1,iy))
|
|
allocate(ice_lon_box(ix+1,iy+1))
|
|
allocate(ice_lat_box(ix+1,iy+1))
|
|
!
|
|
! Extract longitudes and latitudes
|
|
!
|
|
ifail=-3
|
|
call get_data1(trim(fname1),'longitude',ix,iy,iz,it,long,lat,zcoord,tcoord,ice_lon,ifail)
|
|
if (ifail.ne.0) then
|
|
print '(a,6i5)','Failed to read data from file '//trim(fname1),ix,iy,iz,it,ifail
|
|
stop 1
|
|
else
|
|
print*,' read longitude ',ix,iy,iz,it,ifail
|
|
end if
|
|
!
|
|
ifail=-3
|
|
call get_data1(trim(fname1),'latitude',ix,iy,iz,it,long,lat,zcoord,tcoord,ice_lat,ifail)
|
|
if (ifail.ne.0) then
|
|
print '(a,6i5)','Failed to read data from file '//trim(fname1),ix,iy,iz,it,ifail
|
|
stop 1
|
|
else
|
|
print*,' read latitude ',ix,iy,iz,it,ifail
|
|
end if
|
|
!
|
|
! Do some checks
|
|
!
|
|
alon_max=-1000.0
|
|
alon_min=1000.0
|
|
alat_max=-1000.0
|
|
alat_min=1000.0
|
|
do j=1,iy
|
|
do i=1,ix
|
|
! if (fname1.eq.'Anice_LonLat_greenland.nc') then
|
|
! if (ice_lon(i,j).lt.180.0) ice_lon(i,j)=ice_lon(i,j)+360.0
|
|
! else if (fname1.eq.'Anice_LonLat_eurasia.nc') then
|
|
! if (ice_lon(i,j).gt.210.0) ice_lon(i,j)=ice_lon(i,j)-360.0
|
|
! end if
|
|
if (i.ge.2) then
|
|
if (ice_lon(i,j).lt.ice_lon(i-1,j)) then
|
|
print*,i,j,ice_lon(i-1,j),ice_lon(i,j)
|
|
end if
|
|
end if
|
|
if (ice_lon(i,j).lt.alon_min) alon_min=ice_lon(i,j)
|
|
if (ice_lon(i,j).gt.alon_max) alon_max=ice_lon(i,j)
|
|
if (ice_lat(i,j).lt.alat_min) alat_min=ice_lat(i,j)
|
|
if (ice_lat(i,j).gt.alat_max) alat_max=ice_lat(i,j)
|
|
end do
|
|
end do
|
|
!
|
|
print*,' Longitudes ',alon_min,alon_max
|
|
print*,' Latitudes ',alat_min,alat_max
|
|
!
|
|
! cdo operator needs the corners of the grid boxes, not the centre so create
|
|
! the corners by assuming corners are half way between adjacent centres.
|
|
! At edge assume that they go an equal distance beyond.
|
|
!
|
|
do j=1,iy
|
|
dx=ice_lon(2,j)-ice_lon(1,j)
|
|
ice_lon_temp(1,j)=ice_lon(1,j)-0.5*dx
|
|
do i=2,ix
|
|
ice_lon_temp(i,j)=0.5*(ice_lon(i-1,j)+ice_lon(i,j))
|
|
end do
|
|
dx=ice_lon(ix,j)-ice_lon(ix-1,j)
|
|
ice_lon_temp(ix+1,j)=ice_lon(ix,j)+0.5*dx
|
|
end do
|
|
!
|
|
do i=1,ix+1
|
|
dy=ice_lon_temp(i,2)-ice_lon_temp(i,1)
|
|
ice_lon_box(i,1)=ice_lon_temp(i,1)-0.5*dy
|
|
do j=2,iy
|
|
ice_lon_box(i,j)=0.5*(ice_lon_temp(i,j-1)+ice_lon_temp(i,j))
|
|
end do
|
|
dy=ice_lon_temp(i,iy)-ice_lon_temp(i,iy-1)
|
|
ice_lon_box(i,iy+1)=ice_lon_temp(i,iy)+0.5*dy
|
|
end do
|
|
!
|
|
do j=1,iy
|
|
dx=ice_lat(2,j)-ice_lat(1,j)
|
|
ice_lat_temp(1,j)=ice_lat(1,j)-0.5*dx
|
|
do i=2,ix
|
|
ice_lat_temp(i,j)=0.5*(ice_lat(i-1,j)+ice_lat(i,j))
|
|
end do
|
|
dx=ice_lat(ix,j)-ice_lat(ix-1,j)
|
|
ice_lat_temp(ix+1,j)=ice_lat(ix,j)+0.5*dx
|
|
end do
|
|
!
|
|
do i=1,ix+1
|
|
dy=ice_lat_temp(i,2)-ice_lat_temp(i,1)
|
|
ice_lat_box(i,1)=ice_lat_temp(i,1)-0.5*dy
|
|
do j=2,iy
|
|
ice_lat_box(i,j)=0.5*(ice_lat_temp(i,j-1)+ice_lat_temp(i,j))
|
|
end do
|
|
dy=ice_lat_temp(i,iy)-ice_lat_temp(i,iy-1)
|
|
ice_lat_box(i,iy+1)=ice_lat_temp(i,iy)+0.5*dy
|
|
end do
|
|
!
|
|
print*,ice_lon(1,1)
|
|
print*,ice_lon_box(1,1),ice_lon_box(2,1),ice_lon_box(1,2),ice_lon_box(2,2)
|
|
print*,ice_lon(ix,1)
|
|
print*,ice_lon_box(ix,1),ice_lon_box(ix+1,1),ice_lon_box(ix,2),ice_lon_box(ix+1,2)
|
|
!
|
|
print*,ice_lon(1,iy)
|
|
print*,ice_lon_box(1,iy),ice_lon_box(2,iy),ice_lon_box(1,iy+1),ice_lon_box(2,iy+1)
|
|
print*,ice_lon(ix,iy)
|
|
print*,ice_lon_box(ix,iy),ice_lon_box(ix+1,iy),ice_lon_box(ix,iy+1),ice_lon_box(ix+1,iy+1)
|
|
!
|
|
! Now write out data in appropriate format
|
|
!
|
|
open(unit=10,file=trim(oname),form='formatted',status='unknown')
|
|
write(10,'(a)')'gridtype = curvilinear'
|
|
write(10,'(a,i5)')'gridsize = ',ix*iy
|
|
write(10,'(a,i3)')'xsize = ',ix
|
|
write(10,'(a,i3)')'ysize = ',iy
|
|
!
|
|
write(10,'(a)')'xvals = '
|
|
write(10,'(1x,8f7.2)')((ice_lon(i,j),i=1,ix),j=1,iy)
|
|
write(10,'(a)')'xbounds = '
|
|
do j=1,iy
|
|
do i=1,ix
|
|
write(10,'(1x,4f7.2)')ice_lon_box(i,j),ice_lon_box(i,j+1),ice_lon_box(i+1,j+1),ice_lon_box(i+1,j)
|
|
end do
|
|
end do
|
|
!
|
|
write(10,'(a)')'yvals = '
|
|
write(10,'(1x,8f7.2)')((ice_lat(i,j),i=1,ix),j=1,iy)
|
|
write(10,'(a)')'ybounds = '
|
|
do j=1,iy
|
|
do i=1,ix
|
|
write(10,'(1x,4f7.2)')ice_lat_box(i,j),ice_lat_box(i,j+1),ice_lat_box(i+1,j+1),ice_lat_box(i+1,j)
|
|
end do
|
|
end do
|
|
close(unit=10)
|
|
!
|
|
stop
|
|
end PROGRAM FILTERS
|