- Merge original netCDF files
- Compute monthly mean
make_data.csh
#!/bin/csh
cat > name.txt << EOF
h Have sea surface height m
u Uave zonal velocity m/s
v Vave meridonal velocity m/s
w Wave vertical velocity m/s
t Tave temperature deg C
s Save salinity psu
EOF
cat > a.f << EOF
program make_data
implicit none
!=======================================================================
!
!=======================================================================
real , parameter :: slon = 25.0
real , parameter :: elon = 125.0
real , parameter :: slat = -30.0
real , parameter :: elat = 30.0
character(len=120) :: var1(6)
character(len=120) :: var2(6)
character(len=120) :: long(6)
character(len=120) :: unit(6)
integer, parameter :: nlon = 360
integer, parameter :: nlat = 224
integer, parameter :: nlev = 46
character(len=120) :: dir
character(len=300) :: ifile
character(len=300) :: ofile
character(len=120) :: var
character(len=120) :: subdir, stamp
logical :: yes
integer :: ier
real :: buf (nlon,nlat)
real :: tmp1(nlon,nlat)
real :: tmp2(nlon,nlat)
real :: tmp3(nlon,nlat)
real :: tmp4(nlon,nlat)
real :: rlon(nlon)
real :: rlat(nlat)
real :: rlev(nlev)
real :: rlon_u(nlon)
real :: rlat_v(nlat)
real, allocatable :: out (:,:)
real, allocatable :: xdim(:)
real, allocatable :: ydim(:)
real, allocatable :: zdim(:)
real, allocatable :: xt(:,:), xu(:,:), xv(:,:)
real, allocatable :: yt(:,:), yu(:,:), yv(:,:)
integer :: im, jm, km, is, ie, js, je
real :: undef = -999.9
character(len=72) :: time_unit = 'days'
real(8) :: time
integer :: id(4), kt
integer :: iyr, mon, itim, ntim, ktim
integer :: i, j, k
integer :: ivar
real :: fac1, fac2, fac3, fac4
open (10,file='name.txt')
do i = 1, 6
read(10,'(a1,1x,a4,1x,a18,1x,a5)')
& var1(i),var2(i),long(i),unit(i)
print *,'',trim(var1(i)),' ',trim(var2(i)),' ',trim(long(i))
end do
!-----------------------------------------------------------------------
! EPrelaxAve_08_08.06480_08880_240.cdf
! Have_08_08.06480_08880_240.cdf
! Save_08_08.06480_08880_240.cdf
! Tave_08_08.06480_08880_240.cdf
! Uave_08_08.06480_08880_240.cdf
! Vave_08_08.06480_08880_240.cdf
! Wave_08_08.06480_08880_240.cdf
!-----------------------------------------------------------------------
do iyr = 1993, 2014
! do iyr = 2015, 2015
fac1 = 10.0
fac2 = 10.0
fac3 = 10.0
fac4 = 5.0
if (mod(iyr,4) == 0) fac4 = 6.0
write(dir,'(a,i4.4,a)') '../dr080g_',iyr,'/'
dir='/work4/data/ECCO.dr080g/dr080g_1993/'
! dir = '/mnt/hgfs/Downloads/'
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
subdir = 'n10day_01_09/'
! subdir = 'Ibnu/'
var = 'Wave'
stamp = '_08_08.00001_02160_240.cdf'
! stamp = '_08_08.00001_01920_240.cdf'
write(ifile,'(a,a,a,a)')
& trim(dir),trim(subdir),trim(var),trim(stamp)
inquire(file=trim(ifile),exist=yes)
if (.not. yes) then
print *,' error : no file : ',trim(ifile)
stop
end if
print *,'ifile=',trim(ifile)
call read_nc_dim(ifile,'lon',nlon,rlon)
call read_nc_dim(ifile,'lat',nlat,rlat)
call read_nc_dim(ifile,'depth',nlev,rlev)
call read_nc_dim(ifile,'lon_u',nlon,rlon_u)
call read_nc_dim(ifile,'lat_v',nlat,rlat_v)
do i = 1, nlon
if (rlon(i) <= slon) is = i
if (rlon(i) <= elon) ie = i
end do
im = ie-is+1
do j = 1, nlat
if (rlat(j) <= slat) js = j
if (rlat(j) <= elat) je = j
end do
jm = je-js+1
do k = 1, nlev
if (rlev(k) <= 300.0) km = k
end do
allocate(xdim(im), ydim(jm), zdim(km))
allocate(xt(im,jm), xu(im,jm), xv(im,jm))
allocate(yt(im,jm), yu(im,jm), yv(im,jm))
allocate(out(im,jm))
do i = 1, im
xdim(i) = rlon(i+is-1)
end do
do j = 1, jm
ydim(j) = rlat(j+js-1)
end do
do k = 1, km
zdim(k) = rlev(k)
end do
print *,'',is,rlon(is),js,rlat(js),im,jm
do j = 1, jm
do i = 1, im
xt(i,j) = rlon(i+is-1)
xu(i,j) = rlon_u(i+is-1)
xv(i,j) = rlon(i+is-1)
yt(i,j) = rlat(j+js-1)
yu(i,j) = rlat(j+js-1)
yv(i,j) = rlat_v(j+js-1)
end do
end do
! ofile = '/mnt/hgfs/Downloads/Ibnu/ECCO/ecco-jpl.grid.nc'
ofile = './ecco-jpl.grid.nc'
print *,'ofile=',trim(ofile)
id(1) = 1; id(2) = 1; id(3) = 1; id(4) = 0
time = 0
call def_nc_file(ofile,'ECCO JPL Smoother MITGCM dr080g',
& im,jm,km,xdim,ydim,zdim,time_unit,id)
call def_nc_var(ofile, 'xt', 'xt', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call def_nc_var(ofile, 'xu', 'xu', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call def_nc_var(ofile, 'xv', 'xv', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call def_nc_var(ofile, 'yt', 'yt', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call def_nc_var(ofile, 'yu', 'yu', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call def_nc_var(ofile, 'yv', 'yv', 'degree',
& 'xyt','real4', 0., 1.0, undef)
call write_nc_data(ofile, 'xt',
& 'xyt', im, jm, xt, 1, 1, time, undef)
call write_nc_data(ofile, 'xu',
& 'xyt', im, jm, xu, 1, 1, time, undef)
call write_nc_data(ofile, 'xv',
& 'xyt', im, jm, xv, 1, 1, time, undef)
call write_nc_data(ofile, 'yt',
& 'xyt', im, jm, yt, 1, 1, time, undef)
call write_nc_data(ofile, 'yu',
& 'xyt', im, jm, yu, 1, 1, time, undef)
call write_nc_data(ofile, 'yv',
& 'xyt', im, jm, yv, 1, 1, time, undef)
! stop
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
write(ofile,'(a,i4.4,a)')
! & 'dat2/ecco-jpl.',iyr,'.nc'
! & '/mnt/hgfs/Downloads/Ibnu/ECCO/ecco-jpl.',iyr,'.nc'
& './ecco-jpl.',iyr,'.nc'
print *,'ofile=',trim(ofile)
id(1) = 1; id(2) = 1; id(3) = 1; id(4) = 0
call def_nc_file(ofile,'ECCO JPL Smoother MITGCM dr080g',
& im,jm,km,xdim,ydim,zdim,time_unit,id)
do ivar = 1, 1
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& trim(unit(ivar)),
& 'xyt','real4', 0., 1.0, undef)
end do
do ivar = 2, 6
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& trim(unit(ivar)),
& 'xyzt','real4', 0., 1.0, undef)
end do
do ivar = 1, 6
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
mon = 0
do ntim = 1, 4
! do ntim = 1, 1
if (ntim == 1) then
subdir = 'n10day_01_09/'
stamp = '_08_08.00001_02160_240.cdf'
! stamp = '_08_08.00001_01920_240.cdf'
end if
if (ntim == 2) then
subdir = 'n10day_10_18/'
stamp = '_08_08.02160_04320_240.cdf'
end if
if (ntim == 3) then
subdir = 'n10day_19_27/'
stamp = '_08_08.04320_06480_240.cdf'
end if
if (ntim == 4) then
subdir = 'n10day_28_37/'
stamp = '_08_08.06480_08880_240.cdf'
end if
! subdir = ''
write(ifile,'(a,a,a,a)')
& trim(dir),trim(subdir),trim(var2(ivar)),trim(stamp)
inquire(file=trim(ifile),exist=yes)
if (.not. yes) then
print *,' error : no file : ',trim(ifile)
print *,'ntim =',ntim
print *,'subdir=',trim(subdir)
stop
end if
print *,'ifile=',trim(ifile)
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do itim = 1, 3
! do itim = 1, 1
mon = mon + 1
id(1) = iyr; id(2) = mon; id(3) = 1; id(4) = 0
call gettime(id(1), id(2), id(3), kt)
if (time_unit == 'hours') time = (kt-1)*24 + id(4)
if (time_unit == 'days' ) time = (kt-1)
ktim = (itim-1)*3
if (ivar >= 2) then
do k = 1, km
if (ntim == 4 .and. itim == 3) then
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp1,k,1+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp2,k,2+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp3,k,3+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp4,k,4+ktim,undef)
do j = 1, nlat
do i = 1, nlon
if (tmp1(i,j) > undef+1 .and.
& tmp2(i,j) > undef+1 .and.
& tmp3(i,j) > undef+1 .and.
& tmp4(i,j) > undef+1 ) then
buf(i,j) = tmp1(i,j)*fac1
& + tmp2(i,j)*fac2
& + tmp3(i,j)*fac3
& + tmp4(i,j)*fac4
buf(i,j) = buf(i,j)/(fac1+fac2+fac3+fac4)
else
buf(i,j) = undef
end if
end do
end do
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp1,k,1+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp2,k,2+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp3,k,3+ktim,undef)
do j = 1, nlat
do i = 1, nlon
if (tmp1(i,j) > undef+1 .and.
& tmp2(i,j) > undef+1 .and.
& tmp3(i,j) > undef+1 ) then
buf(i,j) = tmp1(i,j)*fac1
& + tmp2(i,j)*fac2
& + tmp3(i,j)*fac3
buf(i,j) = buf(i,j)/(fac1+fac2+fac3)
else
buf(i,j) = undef
end if
end do
end do
end if
call move_array(nlon,nlat,buf,im,jm,out,is,js,undef)
call write_nc_data(ofile, trim(var1(ivar)),
& 'xyzt', im, jm, out, k, mon, time, undef)
end do
else
do k = 1, 1
if (ntim == 4 .and. itim == 3) then
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp1,k,1+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp2,k,2+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp3,k,3+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp4,k,4+ktim,undef)
do j = 1, nlat
do i = 1, nlon
if (tmp1(i,j) > undef+1 .and.
& tmp2(i,j) > undef+1 .and.
& tmp3(i,j) > undef+1 .and.
& tmp4(i,j) > undef+1 ) then
buf(i,j) = tmp1(i,j)*fac1
& + tmp2(i,j)*fac2
& + tmp3(i,j)*fac3
& + tmp4(i,j)*fac4
buf(i,j) = buf(i,j)/(fac1+fac2+fac3+fac4)
else
buf(i,j) = undef
end if
end do
end do
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp1,k,1+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp2,k,2+ktim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp3,k,3+ktim,undef)
do j = 1, nlat
do i = 1, nlon
if (tmp1(i,j) > undef+1 .and.
& tmp2(i,j) > undef+1 .and.
& tmp3(i,j) > undef+1 ) then
buf(i,j) = tmp1(i,j)*fac1
& + tmp2(i,j)*fac2
& + tmp3(i,j)*fac3
buf(i,j) = buf(i,j)/(fac1+fac2+fac3)
else
buf(i,j) = undef
end if
end do
end do
end if
call move_array(nlon,nlat,buf,im,jm,out,is,js,undef)
call write_nc_data(ofile, trim(var1(ivar)),
& 'xyt', im, jm, out, k, mon, time, undef)
end do
end if
end do ! loop for ivar
end do ! loop for itim
end do ! loop for ntim
deallocate(xdim, ydim, zdim, out)
deallocate(xt, xu, xv)
deallocate(yt, yu, yv)
end do ! loop for year
stop
end
subroutine move_array(nlon,nlat,buf,im,jm,out,ioff,joff,undef)
implicit none
!=======================================================================
!
!=======================================================================
integer :: nlon, nlat
integer :: im, jm, ioff, joff
real :: buf(nlon,nlat)
real :: out(im,jm)
real :: undef
integer :: i, j
do j = 1, jm
do i = 1, im
out(i,j) = buf(i+ioff-1,j+joff-1)
if (abs(out(i,j)) > 9.99e8) out(i,j) = undef
end do
end do
return
end
EOF
set FC = "ifort"
set OPT = "-assume byterecl -ftrapuv -fpe0 -traceback -CB"
set NCLIB = "../netcdf/*.o -L/usr/local/lib -lnetcdf"
#set FC = "gfortran"
#set OPT = "-fconvert=big-endian"
#set NCLIB = "$home/local/netcdf/*.o -L/usr/local/netcdf-3.6.2/lib -lnetcdf"
$FC $OPT a.f $NCLIB && ./a.out
# /bin/rm a.f a.out