[2015年 12月 4日 金曜日 12:38:01 JST]
[~/2015.Ibnu.Indian.Ocean/ECCO.dr080g.Fortran/make_budget]
[am@aofd165]
$ csh make_hbudget.csh
$ csh make_sbudget.csh
$ ls -lh ecco-jpl-h* |head -2 ;echo "....."; ls -lh ecco-jpl-s* |head -2
-rw-rw-r-- 1 am am 124M 12月 4 11:47 ecco-jpl-hbudget.1993.nc
-rw-rw-r-- 1 am am 124M 12月 4 11:48 ecco-jpl-hbudget.1994.nc
.....
-rw-rw-r-- 1 am am 112M 12月 4 12:00 ecco-jpl-sbudget.1993.nc
-rw-rw-r-- 1 am am 112M 12月 4 12:00 ecco-jpl-sbudget.1994.nc
$ more make_hbudget.csh make_sbudget.csh>tmp.txt
::::::::::::::
make_hbudget.csh
::::::::::::::
#!/bin/csh
cat > name.txt << EOF
uadv gtUave temperature tendency due to zonal advecton
vadv gtVave temperature tendency due to meridonal advection
wadv gtWave temperature tendency due to vertical advection
xdif gtXdiffave temperature tendency due to zonal diffusion
ydif gtYdiffave temperature tendency due to meridonal diffusion
zdif gtZdiffave temperature tendency due to vertical diffusion
gmdif gtZdiffGMave temperature tendency due to vertical component of GM diffusion
kppmix gtKPPave temperature tendency due to KPP mixing
kf_cor gtKFave temperature tendency due to Kalman filter correction
sflux gtExtave temperature tendency due to surface heat flux include SWR
sdamp gtRelaxave temperature tendency due to surface restoring
tice gtIceave temperature tendency due to ice formation
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) :: long(12)
character(len=120) :: var1(12)
character(len=120) :: var2(12)
integer, parameter :: nlon = 360
integer, parameter :: nlat = 224
integer, parameter :: nlev = 46
character(len=120) :: dir
character(len=120) :: ifile
character(len=120) :: ofile
character(len=120) :: var
character(len=120) :: subdir, stamp
logical :: yes
integer :: ier
real :: buf (nlon,nlat)
real :: tmp (nlon,nlat)
real :: rlon(nlon)
real :: rlat(nlat)
real :: rlev(nlev)
real, allocatable :: out (:,:)
real, allocatable :: xdim(:)
real, allocatable :: ydim(:)
real, allocatable :: zdim(:)
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
integer :: i, j, k
integer :: ivar
real :: fac1, fac2
open (10,file='name.txt')
do i = 1, 12
read(10,'(a6,1x,a12,1x,a70)') var1(i),var2(i),long(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
fac1 = 30.0
fac2 = 5.0
if (mod(iyr,4) == 0) fac2 = 6.0
write(dir,'(a,i4.4,a)') '/work4/data/ECCO.dr080g/dr080g_',iyr,'/'
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
subdir = 'n10day_01_09/'
var = 'gtUave'
stamp = '_08_08.00001_02160_720.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
call read_nc_dim(ifile,'lon',nlon,rlon)
call read_nc_dim(ifile,'lat',nlat,rlat)
call read_nc_dim(ifile,'depth',nlev,rlev)
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))
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
allocate(out(im,jm))
print *,'',is,rlon(is),js,rlat(js),im,jm
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
write(ofile,'(a,i4.4,a)')
& 'ecco-jpl-hbudget.',iyr,'.nc'
! & 'dat/ecco-jpl-hbudget.',iyr,'.nc'
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, 10
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& 'deg C/s',
& 'xyzt','real4', 0., 1.0, undef)
end do
do ivar = 11, 12
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& 'deg C/s',
& 'xyt','real4', 0., 1.0, undef)
end do
do ivar = 1, 12
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
mon = 0
do ntim = 1, 4
if (ntim == 1) then
subdir = 'n10day_01_09/'
stamp = '_08_08.00001_02160_720.cdf'
end if
if (ntim == 2) then
subdir = 'n10day_10_18/'
stamp = '_08_08.02160_04320_720.cdf'
end if
if (ntim == 3) then
subdir = 'n10day_19_27/'
stamp = '_08_08.04320_06480_720.cdf'
end if
if (ntim == 4) then
subdir = 'n10day_28_37/'
stamp = '_08_08.06480_08880_720.cdf'
end if
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)
stop
end if
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do itim = 1, 3
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)
if (ivar <= 10) then
do k = 1, km
if (ntim == 4 .and. itim == 3) then
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp,k,itim+1,undef)
where (tmp > undef+1)
buf = (buf*fac1 + tmp*fac2)/(fac1+fac2)
end where
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
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,buf,k,itim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp,k,itim+1,undef)
where (tmp > undef+1)
buf = (buf*fac1 + tmp*fac2)/(fac1+fac2)
end where
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
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)
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 = "gfortran"
#set OPT = "-fconvert=big-endian"
#set NCLIB = "$home/local/netcdf/*.o -L/usr/local/netcdf-3.6.2/lib -lnetcdf"
set FC = "ifort"
set OPT = "-assume byterecl -ftrapuv -fpe0 -traceback -CB"
set NCLIB = "../netcdf/*.o -L/usr/local/lib -lnetcdf"
$FC $OPT a.f $NCLIB && ./a.out
#/bin/rm a.f a.out
::::::::::::::
make_sbudget.csh
::::::::::::::
#!/bin/csh
cat > name.txt << EOF
uadv gsUave salinity tendency due to zonal advecton
vadv gsVave salinity tendency due to meridonal advection
wadv gsWave salinity tendency due to vertical advection
xdif gsXdiffave salinity tendency due to zonal diffusion
ydif gsYdiffave salinity tendency due to meridonal diffusion
zdif gsZdiffave salinity tendency due to vertical diffusion
gmdif gsZdiffGMave salinity tendency due to vertical component of GM diffusion
kppmix gsKPPave salinity tendency due to KPP mixing
kf_cor gsKFave salinity tendency due to Kalman filter correction
sflux gsExtave salinity tendency due to evaporation - precipitation
sdamp gsRelaxave salinity tendency due to surface restoring
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) :: long(11)
character(len=120) :: var1(11)
character(len=120) :: var2(11)
integer, parameter :: nlon = 360
integer, parameter :: nlat = 224
integer, parameter :: nlev = 46
character(len=120) :: dir
character(len=120) :: ifile
character(len=120) :: ofile
character(len=120) :: var
character(len=120) :: subdir, stamp
logical :: yes
integer :: ier
real :: buf (nlon,nlat)
real :: tmp (nlon,nlat)
real :: rlon(nlon)
real :: rlat(nlat)
real :: rlev(nlev)
real, allocatable :: out (:,:)
real, allocatable :: xdim(:)
real, allocatable :: ydim(:)
real, allocatable :: zdim(:)
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
integer :: i, j, k
integer :: ivar
real :: fac1, fac2
open (10,file='name.txt')
do i = 1, 11
read(10,'(a6,1x,a12,1x,a70)') var1(i),var2(i),long(i)
! print *,'',trim(var1(i)),' ',trim(var2(i)),' ',trim(long(i))
end do
! stop
!-----------------------------------------------------------------------
! 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
fac1 = 30.0
fac2 = 5.0
if (mod(iyr,4) == 0) fac2 = 6.0
write(dir,'(a,i4.4,a)') '/work4/data/ECCO.dr080g/dr080g_',iyr,'/'
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
subdir = 'n10day_01_09/'
var = 'gsUave'
stamp = '_08_08.00001_02160_720.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
call read_nc_dim(ifile,'lon',nlon,rlon)
call read_nc_dim(ifile,'lat',nlat,rlat)
call read_nc_dim(ifile,'depth',nlev,rlev)
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))
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
allocate(out(im,jm))
print *,'',is,rlon(is),js,rlat(js),im,jm
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
write(ofile,'(a,i4.4,a)')
& 'ecco-jpl-sbudget.',iyr,'.nc'
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, 9
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& 'psu/s',
& 'xyzt','real4', 0., 1.0, undef)
end do
do ivar = 10, 11
call def_nc_var(ofile, trim(var1(ivar)), trim(long(ivar)),
& 'psu/s',
& 'xyt','real4', 0., 1.0, undef)
end do
do ivar = 1, 11
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
mon = 0
do ntim = 1, 4
if (ntim == 1) then
subdir = 'n10day_01_09/'
stamp = '_08_08.00001_02160_720.cdf'
end if
if (ntim == 2) then
subdir = 'n10day_10_18/'
stamp = '_08_08.02160_04320_720.cdf'
end if
if (ntim == 3) then
subdir = 'n10day_19_27/'
stamp = '_08_08.04320_06480_720.cdf'
end if
if (ntim == 4) then
subdir = 'n10day_28_37/'
stamp = '_08_08.06480_08880_720.cdf'
end if
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)
stop
end if
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do itim = 1, 3
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)
if (ivar <= 9) then
do k = 1, km
if (ntim == 4 .and. itim == 3) then
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp,k,itim+1,undef)
where (tmp > undef+1)
buf = (buf*fac1 + tmp*fac2)/(fac1+fac2)
end where
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
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,buf,k,itim,undef)
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,tmp,k,itim+1,undef)
where (tmp > undef+1)
buf = (buf*fac1 + tmp*fac2)/(fac1+fac2)
end where
else
call read_nc_data(ifile,trim(var2(ivar)),
& nlon,nlat,buf,k,itim,undef)
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)
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 = "gfortran"
#set OPT = "-fconvert=big-endian"
#set NCLIB = "$home/local/netcdf/*.o -L/usr/local/netcdf-3.6.2/lib -lnetcdf"
set FC = "ifort"
set OPT = "-assume byterecl -ftrapuv -fpe0 -traceback -CB"
set NCLIB = "../netcdf/*.o -L/usr/local/lib -lnetcdf"
$FC $OPT a.f $NCLIB && ./a.out
#/bin/rm a.f a.out