リストに記載された日付・時刻のデータから集合平均を求める
リストの例：filelist.Jun.intense-0.356.txt
tadv.125.Jun.w-0.356.eps
tadv.500.Jun.intense.w-0.356.eps
wvf.conv.950hPa.Jun.weak-0.356.eps
[2016年 1月 28日 木曜日 08:41:49 JST]
[/work4/am/Kuroshio.Convection/Composite.Analysis/Create.Composite]
[am@aofd165]
::::::::::::::
create_composite.run.sh
::::::::::::::
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.fish.nagasaki-u.ac.jp
# Directory: /work4/am/Kuroshio.Convection/Set.Study.Area/Omega.Climatology
#
# Revision history:
# This file is created by /usr/local/mybin/nbscr.sh at 13:54 on 03-13-2014.
usage(){
cat <<EOF
Usage : $0 <month> <name>
month: month for which the composite mean will be created (e.g., Jun)
name : conditions for computing composite (e.g., intense-0.356)
EOF
}
echo "Shell script, $(basename $0) starts."
echo
# Default values
exe=create_composite
month=Jun
name=intense-0.356
if [ $# -ne 2 ]; then
echo Error in $0 : Wrong number of arguments.
usage
exit 1
fi
month=$1
name=$2
namelist=${exe}.${month}.${name}.namelist
outdir=output
if [ ! -d $outdir ]; then
mkdir -vp ${outdir}
fi
indir_nc=/work3/satsuki/MSM/data/MSM-P
dir="/work4/am/Kuroshio.Convection/Composite.Analysis/Date.Time.of.Each.Case/output"
date_list=${dir}/date.time.${name}.txt
filelist=filelist.${month}.${name}.txt
echo "filelist= $filelist"
lend=${#indir_nc}
lines=$(awk '{ if($1 !="#") print $0}' $date_list |wc -l | awk '{print $1}' )
echo $lines
echo $lines > $filelist
awk -v d="$indir_nc" -v l=$lend '{if ($1 != "#") \
printf "%29s/%04d/%02d%02d.nc\n %04d %02d %02d %02d\n", d,$3,$4,$5,$3,$4,$5,$6}' $date_list \
>> $filelist
cat <<EOF > $namelist
¶
filelist="${filelist}",
month="${month}",
ctl="${outdir}/composite.${month}.${name}.ctl",
ofle="${outdir}/composite.${month}.${name}.bin",
&end
EOF
$exe < $namelist
echo "Done $(basename $0)"
echo
::::::::::::::
run.all.sh
::::::::::::::
#!/bin/sh
date_s=$(date)
#exe=create_composite.run.sh
#month=Jun
#name=intense-0.356
#(time $exe $month $name)
#echo
#echo
#echo
#echo
#echo
exe=create_composite.run.sh
month=Jun
name=weak-0.356
(time $exe $month $name)
date_e=$(date)
echo "$exe started at $date_s"
echo "$exe finished at $date_e"
::::::::::::::
create_composite.Jun.intense-0.356.namelist
::::::::::::::
¶
filelist="filelist.Jun.intense-0.356.txt",
month="Jun",
ctl="output/composite.Jun.intense-0.356.ctl",
ofle="output/composite.Jun.intense-0.356.bin",
&end
::::::::::::::
create_composite.Jun.weak-0.356.namelist
::::::::::::::
¶
filelist="filelist.Jun.weak-0.356.txt",
month="Jun",
ctl="output/composite.Jun.weak-0.356.ctl",
ofle="output/composite.Jun.weak-0.356.bin",
&end
::::::::::::::
create_composite.test.namelist
::::::::::::::
¶
filelist="filelist.test.intense-0.356.txt",
month="Jun",
ctl="output/composite.test.intense-0.356.ctl",
ofle="output/composite.test.intense-0.356.bin",
&end
::::::::::::::
makefile
::::::::::::::
#
#
#j Fortran77とFortran90/95の混成プログラムのコンパイル用makefile
#
#j 作成 Sun Feb 15 15:16:02 2009
#
#j 注：ごく簡単な動作チェックしか行っていない
#j サフィックスルールがこれで正しいかどうかチェックする必要あり
#
#j 参考文献
#j makefile の書き方
#j http://www.eis.osakafu-u.ac.jp/~yabu/soft/makefile.html
#j Makefile小技と新バージョンのDelFEM(FEMとUIの日記@New York)
#j http://d.hatena.ne.jp/etopirika5/20091207/1260207955
#
#-----------------------------------------
#j マクロの定義
#-----------------------------------------
#j コンパイラの指定
FC=ifort
#FC=gfortran
#FC=g77
#FC=f77
#j OBJDIR (オブジェクトファイルをおくディレクトリ)
OBJDIR=./obj
#-----------------------------------------
#j コンパイルオプション(ifort)
#-----------------------------------------
#j デバッグ用オプション(はじめて実行するときには必ずこのオプションをつけてコンパイルする)
FDFLAGS= -check -CB -traceback -fpe0
#j コンパイルオプション
FFLAGS=-module ${OBJDIR}
#j -module : モジュールファイル(.mod)を置く場所を指定する
#j 最適化用(他にも何種類かある. ifort -helpかマニュアルを見て調べる)
#FFLAGS = -O3 -module ${OBJDIR}
#j 倍精度
#FFLAGS = -r8 -module ${OBJDIR}
#j 入力データのバイト・スワップ(大型計算機のバイナリデータを読む時などにつかう)
#FFLAGS = #-convert big_endian -module ${OBJDIR}
#-----------------------------------------
#j リンク用のオプション
#-----------------------------------------
LDFLAGS=-module ${OBJDIR} -L/usr/local/lib -lnetcdf
#j ターゲット名(最終的に作りたい実行ファイル名)
TARGET1=create_composite
TARGET2=
TARGETS=$(TARGET1) $(TARGET2)
#j MOD (Fortran90以降のモジュールファイル)
MOD=
#j SRC7 (Fortran77のソースファイル)
SRC7=
#j SRC9 (Fortran90のソースファイル)
SRC9=create_composite.f90 print_header.f90 print_ioinfo.f90 \
error_message.f90 basename.f90
#j OBJ1 (ターゲットを作るのに必要なオブジェクトファイル名)
OBJM=$(MOD:.f90=.o)
OBJ7=$(SRC7:.f=.o)
OBJ9=$(SRC9:.f90=.o)
OBJTMP=$(OBJM) $(OBJ7) $(OBJ9)
OBJ1=$(OBJTMP:%=${OBJDIR}/%)
MOD2=
#j SRC7 (Fortran77のソースファイル)
SRC72=
#j SRC9 (Fortran90のソースファイル)
SRC92=
#j OBJ (ターゲットを作るのに必要なオブジェクトファイル名)
OBJM2=$(MOD2:.f90=.o)
OBJ72=$(SRC72:.f=.o)
OBJ92=$(SRC92:.f90=.o)
OBJTMP2=$(OBJM2) $(OBJ72) $(OBJ92)
OBJ2=$(OBJTMP2:%=${OBJDIR}/%)
OBJ=$(OBJ1) $(OBJ2)
#-----------------------------------------
#j ここからコンパイルのルールの記述
#-----------------------------------------
all: mkobjd $(TARGETS)
$(TARGET1): $(OBJ1)
$(FC) -o $@ $(OBJ1) ${LDFLAGS}
$(TARGET2): $(OBJ2)
$(FC) -o $@ $(OBJ2) ${LDFLAGS}
#j 暗黙のサフィックスルールを無効にする
.SUFFIXES:
#j サフィックスの追加
.SUFFIXES: .o .f .f90
#j f77のソースファイルのコンパイル(.f用のサフィックスルール)
${OBJDIR}/%.o : %.f
$(FC) -c ${FDFLAGS} ${FFLAGS} -c -o $@ $<
#j f90のソースファイルのコンパイル(.f90用のサフィックスルール)
${OBJDIR}/%.o : %.f90
$(FC) -c ${FDFLAGS} ${FFLAGS} -c -o $@ $<
#j オブジェクトファイルを削除
clean:
rm -rf $(OBJ) ${OBJDIR}/*.mod
#j オブジェクトファイルと実行ファイルを削除
distclean:
rm -rf $(OBJ) ${OBJDIR}/*.mod $(TARGETS) \
rm -rf ${OBJDIR}
#j オブジェクトファイルを置くディレクトリを作成
mkobjd:
if [ ! -d ${OBJDIR} ]; then \
mkdir -p ${OBJDIR} ; \
fi
run:
cd .. ; \
./run.sh 2>&1 > log.txt ; \
cd src
#j tarファイルを作る
tar:
tar cvf $(TARGET1).tar ./*.f90 ./*.f ./*.c ./*.h ./*.txt makefile
::::::::::::::
basename.f90
::::::::::::::
!***********************************************************************
! .f
! NOTE
! Reference
!***********************************************************************
subroutine basename(infle, bifle, idxbi)
!------------------------- DECLARATION ---------------------------------
! implicit none
! implicit double precision(a-h,o-z)
!------------------------- INCLUDE FILE --------------------------------
! include 'comblk.h'
!-------------------------- PARAMETERS ---------------------------------
! parameter()
!------------------------- COMMON VARIABLES ----------------------------
! common //
!----------------------------ARGUMENTS ---------------------------------
! double precision
! real
! integer
! dimension
! Input
character infle*(*)
! Output
character bifle*(*)
integer idxbi
! Input/output
!------------------------- LOCAL VARIABLES -----------------------------
character testchar*1
!-----------------------------------------------------------------------
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
last=index(infle,' ')-1
do i=last, 1, -1
testchar=infle(i:i)
id=index(testchar,'/')
if(id.ne.0)then
idx=i+1
goto 1000
end if
end do !i
! If there is no slash in infle
bifle=infle(1:last)
idxbi=last
! write(*,'(A)')infle
! write(*,'(A)')bifle
! write(*,*)idxbi
return
! If there is a slash in infle
1000 bifle=infle(idx:last)
idxbi=last-idx+1
! write(*,'(A)')infle
! write(*,'(A)')bifle
! write(*,*)idxbi
return
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
end
::::::::::::::
create_composite.f90
::::::::::::::
program create_composite
!
! Author: am
!
! Host: aofd165
! Directory: /work4/am/Kuroshio.Convection/Set.Study.Area/Omega.Climatology
! [2014年 3月 13日 木曜日 13:16:31 JST]
! Original program:
! Directory: /work3/satsuki/ASCAT/AVE/large_omega.f90
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 17:24 on 08-03-2012.
! Revision history:
! 2011-02-22 16:38
! Initial Version
! References
! http://www.rain.hyarc.nagoya-u.ac.jp/laboratory/OB/shimizu/HTML/netcdf.html
! Tips: NetCDF. http://iprc.soest.hawaii.edu/users/sasakiyo/Japanese/misc/tips_netcdf.html
! NetCDF Tips. http://www.sci.hokudai.ac.jp/grp/poc/top/software/other/netcdf_tips/index.htm
! use
! implicit none
include '/usr/local/include/netcdf.inc'
! ncid:ファイルのID番号、 varid: 変 数のID番号
integer ncid,varid
character(len=500):: filelist,infle
character ctl*500
character ofle*500, ofle_base*500
character var*50
real,allocatable :: lon(:),lat(:),p(:),time(:)
integer(2),allocatable :: uin(:,:,:,:), vin(:,:,:,:),tempin(:,:,:,:),rhin(:,:,:,:)
real(8),allocatable :: zin(:,:,:,:),win(:,:,:,:)
! real,allocatable :: sdata(:,:,:)
integer,allocatable :: dimids(:),nshape(:)
character*70 err_message
character*100 dimname !dimnameは各次元の名前。
real,parameter::RMISS=-999.9 !欠損値
real,allocatable::olon(:),olat(:)
integer,allocatable::idx(:),jdx(:)
real,allocatable::u(:,:,:,:),v(:,:,:,:),temp(:,:,:,:),rh(:,:,:,:)
! win, zinは倍精度だが、w,zは単精度にしている！
real,allocatable ::w(:,:,:,:),z(:,:,:,:)
character month*3
integer stat,sta
character strm*1000
! composite mean
real,allocatable :: wave_c(:,:,:),zave_c(:,:,:)
real,allocatable :: uave_c(:,:,:),vave_c(:,:,:),tempave_c(:,:,:),&
& rhave_c(:,:,:)
integer :: ndat_c=0
namelist /para/filelist,month,ctl,ofle
im=241 ! Number of grids in longitudinal direction 192
jm=253 ! Number of grids in meridional direction !94
km=16 ! Number of vertical levels
lm=8 ! NUmber of 3 hourly data in a day (00:00, 03:00, ..., 21:00)
allocate(lon(im),lat(jm),p(km),time(lm))
allocate(w(im,jm,km,lm),win(im,jm,km,lm),z(im,jm,km,lm),zin(im,jm,km,lm))
allocate(temp(im,jm,km,lm),tempin(im,jm,km,lm),rh(im,jm,km,lm),&
& rhin(im,jm,km,lm))
allocate(u(im,jm,km,lm),uin(im,jm,km,lm),v(im,jm,km,lm),vin(im,jm,km,lm))
! composite mean
allocate(wave_c(im,jm,km),zave_c(im,jm,km))
allocate(tempave_c(im,jm,km),rhave_c(im,jm,km))
allocate(uave_c(im,jm,km),vave_c(im,jm,km))
uave_c(:,:,:)=0.0
vave_c(:,:,:)=0.0
tempave_c(:,:,:)=0.0
rhave_c(:,:,:)=0.0
zave_c(:,:,:)=0.0
wave_c(:,:,:)=0.0
! Read namelist
read(*,para)
! Check I/O files
print '(A)','filelist = ',trim(filelist)
print '(A,A)','ofle = ',trim(ofle)
print '(A,A)','ctl = ',trim(ctl)
! Open file list of MSM data
open(10,file=filelist,action="read")
read(10,*)nm !nm=Number of cases
loop_time: do n=1,nm
read(10,'(A)')infle
read(10,*)iyyyy,imm,idd,ihh
print '(a,i5,a,i5,a,i4.4,1x,i2.2,1x,i2.2,1x,i2.2)',&
"File ",n," of ",nm, " Date/Time ",iyyyy,imm,idd,ihh
print '(A,A)','Input: ',trim(infle)
sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得
if(sta /= 0)call error_message("nf_open",sta)
var='w'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。
sta = nf_inq_varid(ncid,var,varid)
if(sta /= 0)call error_message("nf_inq_varid",sta)
! スケールファクターを得る。(*1)
! sta = nf_get_att_real(ncid,varid,'scale_factor',scale_w)
! if(sta /= 0)call error_message("nf_inq_varid",sta)
! オフセットの 値を得る(*2)
! sta = nf_get_att_real(ncid,varid,'add_offset',offset_w)
! if(sta /= 0)call error_message("nf_get_att_real",sta)
! varidからndimsを得る。ndimsとは次元の数。二次元データなら2
sta = nf_inq_varndims(ncid, varid, ndims)
if(sta /= 0)call error_message("nf_inq_varndims",sta)
! print *,'ndims= ', ndims
! print *
! varidからdimidsを得る。dimidsとはそれぞれの次元に対するID
allocate(dimids(ndims))
sta = nf_inq_vardimid(ncid, varid, dimids)
! do i=1,ndims
! write(6,'("dimids(",i2,")=",i9 )') i,dimids(i)
! enddo
! print *
allocate(nshape(ndims))
do i=1,ndims
sta = nf_inq_dimlen(ncid, dimids(i), nshape(i))
if(sta /= 0)call error_message("nf_inq_dimlen",sta)
! dimidsからnshapeを得る。
! nshapeとはそれぞれの次元に対する最大データ数 (格子数)
! write(6,'("nshape(",i2,")=",i9 )') i,nshape(i)
enddo !i
! print *
do i=1,ndims
sta = nf_inq_dimname(ncid, dimids(i), dimname)
if(sta /= 0)call error_message("nf_inq_dimlen",sta)
! write(6,'("dimname(",i2,")=",1x,a23)') i,dimname
enddo
sta = nf_get_var_double(ncid, varid, win)
if(sta /= 0)call error_message("nf_get_var_double",sta)
deallocate(dimids,nshape)
var='u'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
sta = nf_get_att_real(ncid,varid,'scale_factor',scale_u)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_u)
stat = nf_inq_varndims(ncid, varid, ndims)
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
allocate(nshape(ndims))
do i=1,ndims
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
enddo !i
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
enddo
stat = nf_get_var_INT2(ncid, varid, uin)
deallocate(dimids,nshape)
var='v'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
sta = nf_get_att_real(ncid,varid,'scale_factor',scale_v)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_v)
stat = nf_inq_varndims(ncid, varid, ndims)
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
allocate(nshape(ndims))
do i=1,ndims
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
enddo !i
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
enddo
stat = nf_get_var_INT2(ncid, varid, vin)
deallocate(dimids,nshape)
var='temp'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
sta = nf_get_att_real(ncid,varid,'scale_factor',scale_temp)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_temp)
stat = nf_inq_varndims(ncid, varid, ndims)
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
allocate(nshape(ndims))
do i=1,ndims
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
enddo !i
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
enddo
stat = nf_get_var_INT2(ncid, varid, tempin)
deallocate(dimids,nshape)
var='rh'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
sta = nf_get_att_real(ncid,varid,'scale_factor',scale_rh)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_rh)
stat = nf_inq_varndims(ncid, varid, ndims)
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
allocate(nshape(ndims))
do i=1,ndims
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
enddo !i
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
enddo
stat = nf_get_var_INT2(ncid, varid, rhin)
deallocate(dimids,nshape)
var='z'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
stat = nf_inq_varndims(ncid, varid, ndims)
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
allocate(nshape(ndims))
do i=1,ndims
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
enddo !i
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
enddo
stat = nf_get_var_double(ncid, varid, zin)
deallocate(dimids,nshape)
var='time'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
if(sta /= 0)call error_message("nf_inq_varid",sta)
stat = nf_get_var_real(ncid, varid, time)
if(sta /= 0)call error_message("nf_get_var_real",sta)
u(:,:,:,:)=uin(:,:,:,:)*scale_u + offset_u
v(:,:,:,:)=vin(:,:,:,:)*scale_v + offset_v
temp(:,:,:,:)=tempin(:,:,:,:)*scale_temp + offset_temp
rh(:,:,:,:)=rhin(:,:,:,:)*scale_rh + offset_rh
z(:,:,:,:)=sngl(zin(:,:,:,:)) !*scale_z + offset_z
w(:,:,:,:)=sngl(win(:,:,:,:)) !*scale_z + offset_z
ldx=ihh/3+1
uave_c(:,:,:)=uave_c(:,:,:)+u(:,:,:,ldx)
vave_c(:,:,:)=vave_c(:,:,:)+v(:,:,:,ldx)
tempave_c(:,:,:)=tempave_c(:,:,:)+temp(:,:,:,ldx)
rhave_c(:,:,:)=rhave_c(:,:,:)+rh(:,:,:,ldx)
zave_c(:,:,:)=zave_c(:,:,:)+z(:,:,:,ldx)
wave_c(:,:,:)=wave_c(:,:,:)+w(:,:,:,ldx)
STA = NF_CLOSE(NCID)
if(sta /= 0)call error_message("NF_CLOSE",sta)
enddo loop_time
close(20)
close(10)
! Compute compoiste means
uave_c(:,:,:) =uave_c(:,:,:)/float(nm)
vave_c(:,:,:) =vave_c(:,:,:)/float(nm)
tempave_c(:,:,:)=tempave_c(:,:,:)/float(nm)
rhave_c(:,:,:) =rhave_c(:,:,:)/float(nm)
zave_c(:,:,:) =zave_c(:,:,:)/float(nm)
wave_c(:,:,:) =wave_c(:,:,:)/float(nm)
! Read p
open(10,file=filelist,action="read")
read(10,*)nm
read(10,'(A)')infle
print '(A,A)','Input: ',trim(infle)
sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得
if(sta /= 0)call error_message("nf_open",sta)
var='p'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
if(sta /= 0)call error_message("nf_inq_varid",sta)
sta = nf_get_var_real(ncid, varid, p)
if(sta /= 0)call error_message("nf_get_var_real",sta)
! Read lon, lat
var='lon'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
if(sta /= 0)call error_message("nf_inq_varid",sta)
sta = nf_get_var_real(ncid, varid, lon)
if(sta /= 0)call error_message("nf_get_var_real",sta)
var='lat'
print '(A,A,A)','Reading ',var(1:lnblnk(var)),' ...'
sta = nf_inq_varid(ncid,var,varid)
if(sta /= 0)call error_message("nf_inq_varid",sta)
sta = nf_get_var_real(ncid, varid, lat)
if(sta /= 0)call error_message("nf_get_var_real",sta)
close(10)
! Print GrADS data
print '(A,A)','Output: ',trim(ctl)
open(20,file=ctl)
call basename(ofle, ofle_base, idxbi)
write(20,'(A,A)')'DSET ^',ofle_base
write(20,'(A)')'OPTIONS sequential template yrev'
write(20,'(A)')'UNDEF -999.9'
write(20,'(A,i5,A,500f12.3)')'XDEF ',im,' LINEAR ', 120.0, 0.125
write(20,'(A,i5,A,500f12.6)')'YDEF ',jm,' LINEAR ', 22.4, 0.1
write(20,'(A,i5,A,500f10.3)')'ZDEF ',km,' LEVELS 1000 975 950 925 900 850 800 700 600 500 400 300 250 200 150 100'
write(20,'(A,i3,A,A3,i4,A)')'TDEF ',1,' LINEAR 00z01',month,iyyyy ,' 1dy'
write(20,'(A)')'VARS 6'
write(20,'(A,i5,A)')'u ',16,' 99 eastward_wind [m/s]'
write(20,'(A,i5,A)')'v ',16,' 99 northward_wind [m/s]'
write(20,'(A,i5,A)')'temp ',16,' 99 temperature [K]'
write(20,'(A,i5,A)')'rh ',16,' 99 relative humidity [%]'
write(20,'(A,i5,A)')'z ',16,' 99 geopotential height [m]'
write(20,'(A,i5,A)')'w ',16,' 99 vertical velocity in p [Pa/s]'
write(20,'(A)')'ENDVARS'
close(20)
open(20,file=ofle,form="unformatted")
do k=1,km
write(20)((uave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
do k=1,km
write(20)((vave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
do k=1,km
write(20)((tempave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
do k=1,km
write(20)((rhave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
do k=1,km
write(20)((zave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
do k=1,km
write(20)((wave_c(i,j,k),i=1,im),j=1,jm)!,k=1,km)
enddo
close(20)
write(*,'(a)')'Done program create_composite.'
write(*,*)
end program create_composite
::::::::::::::
error_message.f90
::::::::::::::
subroutine error_message(subr_name, stat)
! Description:
!
! Author: am
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work4/am/Kuroshio.Convection/Set.Study.Area/Omega.Climatology
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 13:16 on 03-13-2014.
include '/usr/local/include/netcdf.inc'
character(len=*),intent(in)::subr_name
integer,intent(in)::stat
character*70 err_message
err_message = nf_strerror(sta)
print *
print '(A,A)', '*** Error in ',trim(subr_name)
print '(A,i5)','*** Stat = ',stat
print '(A)', '*** Error mesasge: '
print '(A,A)', '*** ',err_message
print '(A)', '*** Program Terminated.'
print *
stop
end subroutine error_message
::::::::::::::
print_header.f90
::::::::::::::
!------------------------------------------------------------------------------
! This source code can be translated into TeX documentation using a Perl script
! called ProTex.
! ProTex is available via:
! http://gmao.gsfc.nasa.gov/software/protex/ .
! On-line manual of ProTex can be found at:
! http://gmao.gsfc.nasa.gov/software/protex/doc/protex/ .
!
! Usage: protex src > out.tex
! src: source file name(s)
! out: output tex file name
! Type "protex -h" for help.
!
!------------------------------------------------------------------------------
!BOP
!
! !ROUTINE: print_header
!
! !DESCRIPTION: 出力ファイルに以下のヘッダ情報を印字する\\
! 作成日時 \\
! 現在作業しているディレクトリ名 \\
! ユーザー名 \\
! ホスト名 \\
! プログラム名 \\
! \\
! 使用法 \\
! call print_header(iunit) \\
! iunit: 出力ファイルの機番(整数) \\
!
! !INTERFACE:
subroutine print_header(iunit)
!
! !USES:
! include
! use
!
!
! !INPUT PARAMETERS:
integer iunit
! !OUTPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
character(len=500) :: dirinfo, userinfo, hostinfo
INTEGER date_time(8)
CHARACTER REAL_CLOCK(3)*12
character(len=500) :: progname,shellinfo
!
! !BUGS:
! None known at this time.
!
! !SEE ALSO:
!
! !SYSTEM ROUTINES:
! None
!
! !FILES USED:
!
!
! !REVISION HISTORY:
! [Thu Aug 20 14:17:59 JST 2009] - A. Manda - Initial version
!
! !REMARKS:
! このサブルーチンは, g77, ifortで動作することを確認している.
!
! 作成したデータが、何時、何処で、誰が、どのプログラムを使って作成した
! か記録しておくと、後々作業をし直す必要が生じたときに役に立つ.
! また、以前使ったプログラムを一部改変して利用するときなど、そのプログ
! ラムを探すのに役に立つことが多い
! !TO DO:
!
!EOP
!------------------------------------------------------------------------------
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
!BOC
! Get date, time, current directory, username, and hostname
! p.345 of for_lang.pdf (intel fortran user manual)
! DATE_AND_TIME is an intrisic S/R and works on most of Fotran 90
! Compilers. g77 also supports this S/R.
CALL DATE_AND_TIME (REAL_CLOCK (1), REAL_CLOCK (2), &
& REAL_CLOCK (3), date_time)
call getenv("PWD", dirinfo)
call getenv("USER", userinfo)
!j 使用しているシェルのチェック
!j 現在, sh, bash, csh, tcshのみに対応
call getenv("SHELL",shellinfo)
idxcsh=index(shellinfo,"/csh")
idxtcsh=index(shellinfo,"/tcsh")
idxsh=index(shellinfo,"/sh")
idxbash=index(shellinfo,"/bash")
if(idxcsh.ne.0.or.idxtcsh.ne.0)then
call getenv("HOST", hostinfo) !csh, tcsh, etc.
else if(idxsh.ne.0.or.idxbash.ne.0)then
call getenv("HOSTNAME", hostinfo) !sh, bash, etc.
endif
write(iunit,'(a, &
& i2.2,a,i2.2,a,i4.4, a,i2.2,a,i2.2,a,i2.2, a,i3.2,a)') &
& '# Date and time: ', &
& date_time(2),'/',date_time(3),'/',date_time(1), &
& ' at ',date_time(5),':',date_time(6),':',date_time(7), &
& ' ',-date_time(4)/60,':00'
write(iunit,'(A,A)')'# hostname: ',hostinfo(1:lnblnk(hostinfo))
write(iunit,'(A,A)')'# cwd: ',dirinfo(1:lnblnk(dirinfo))
write(iunit,'(A,A)')'# user: ',userinfo(1:lnblnk(userinfo))
! Get program name
call getarg( 0, progname)
write(iunit,'(A,A)')'# program name: ',progname(1:lnblnk(progname))
return
end subroutine print_header
!EOC
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
::::::::::::::
print_ioinfo.f90
::::::::::::::
!------------------------------------------------------------------------------
! This source code can be translated into TeX documentation using a Perl script
! called ProTex.
! ProTex is available via:
! http://gmao.gsfc.nasa.gov/software/protex/ .
! On-line manual of ProTex can be found at:
! http://gmao.gsfc.nasa.gov/software/protex/doc/protex/ .
!
! Usage: protex src > out.tex
! src: source file name(s)
! out: output tex file name
! Type "protex -h" for help.
!
!------------------------------------------------------------------------------
!BOP
!
! !ROUTINE: print_ioinfo
!
! !DESCRIPTION: 出力ファイルにヘッダ情報として、入出力ファイル名を印字する\\
! 使用法 \\
! call print_ioinfo(iunit,infle,iokind) \\
! iunit: 出力ファイルの機番(整数) \\
! fle: 入出力ファイル名(文字型) \\
! iokind: 入出力ファイルの種類(文字型): input or output \\
! \\
! !INTERFACE:
subroutine print_ioinfo(iunit,fle,iokind)
!
! !USES:
!
! !INPUT PARAMETERS:
integer,intent(in):: iunit
character,intent(in):: fle*(*), iokind*(*)
! !OUTPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
!
! !BUGS:
!
! !SEE ALSO:
!
! !SYSTEM ROUTINES:
!
! !FILES USED:
!
! !REVISION HISTORY:
! [Thu Aug 20 14:17:59 JST 2009] - A. Manda - Initial version
!
! !REMARKS:
! このサブルーチンは, g77, ifortで動作することを確認している.
!
! 入出力に使ったファイル名を記録しておくと、後々作業をし直す必要
! が生じたときに役に立つ.
!
! !TO DO:
!
!EOP
!------------------------------------------------------------------------------
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
!BOC
write(iunit,'(A,A,A,A)')'# ',iokind,': ',fle(1:lnblnk(fle))
return
end subroutine print_ioinfo
!EOC
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
!------------------------------------------------------------------------------
::::::::::::::
advt_crossec.gs
::::::::::::::
'reinit'
* Default values
lon=125
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/tadv.'lon'.Jun.w-0.356.eps'
figfile2='../Fig/tadv.'lon'.Jun.diff.w-0.356.eps'
'set display color white'
'c'
'set parea 0.5 6 1 8'
title='Strong Updraft 'lon'E'
*'set map 1 1 6'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm ;*msm201206.ctl'
'q file'
*'set strsiz 0.1 0.15'
*'draw string .5 9 Input: 'inputmsm
*'draw string .5 8.8 Input: 'inputmsm2
'set z 1 16'
'dynamic -var z temp u v'
*'define dtx = cdiff(temp,x)'
*'define dty = cdiff(temp,y)'
*'define dx = cdiff(lon,x)*3.1416/180'
*'define dy = cdiff(lat,y)*3.1416/180'
*tadv'=-1*( ('u'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v'*'dty'/'dy' )/6.37e6*84600'
*'set z 10'
*'set gxout shade2'
*'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
*'rgbset2'
*'d 'tadv
*'cbarn'
*'q dims'
*say result
*'allclose'
*return
'set mpdset hires'
'set lon 'lon
'set lat 22 40'
'set grads off'
'mul 2 1 1 1 -novpage'
*'set vpage 1 6 2 8'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
*'set zlog on'
*'set z 1 16'
* Temperature advection
'set gxout shade2'
'rgbset2'
'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
'd 'tadv'*86400'
*'cbarn 0.8 0 4.3 0.12'
'cbarn 1 1 8 5'
'q dims'
say result
'set gxout contour'
'set cthick 10'
'set ccolor 15'
'set clopts 1 10 0.2'
'set clevs -3 0 3 6 9 12 15 18 21 24 27 30 33'
'd u'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
*'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
*'d u'
* Meridional Circulation
'set gxout vector'
'set cthick 10'
'set ccolor 1'
'set arrscl 0.5 5'
'd skip(va,10,1);skip(-w*10,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'close 1'
title='Weak Updraft 'lon'E'
'open 'inputmsm2
'set z 1 16'
'dynamic -var z temp u v'
*'define dtx = cdiff(temp,x)'
*'define dty = cdiff(temp,y)'
*'define dx = cdiff(lon,x)*3.1416/180'
*'define dy = cdiff(lat,y)*3.1416/180'
*tadv'=-1*( ('u'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v'*'dty'/'dy' )/6.37e6*84600'
'set lon 'lon
'set lat 22 40'
*'set vpage 1 6 2 8'
'mul 2 1 2 1 -novpage'
*'set vpage 1 6 2 8'
'set grads off'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
* Temperature advection
'set gxout shade2'
'rgbset2'
'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
'd 'tadv'*86400'
* Zonnal velocity
*'set gxout shade2'
*'color -levs -kind blue->white->red'
*'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
*'rgbset2'
*'set z 1 16'
*'set zlog on'
*'d u'
'set gxout contour'
'set cthick 10'
'set ccolor 15'
'set clopts 1 10 0.2'
'set clevs -3 0 3 6 9 12 15 18 21 24 27 30 33'
'd u'
* Meridional Circulation
'set gxout vector'
'set cthick 10'
'set ccolor 1'
'set arrscl 0.5 5'
'd skip(va,10,1);skip(-w*10,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'print 'figfile
say 'Fig file = 'figfile
'allclose'
::::::::::::::
advt_crossec.legend.gs
::::::::::::::
'reinit'
* Default values
lon=125
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/tadv.'lon'.Jun.w-0.356.eps'
figfile2='../Fig/tadv.'lon'.Jun.diff.w-0.356.eps'
'set display color white'
'c'
'set parea 0.5 6 1 8'
title='Strong Updraft 'lon'E'
*'set map 1 1 6'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm ;*msm201206.ctl'
'q file'
*'set strsiz 0.1 0.15'
*'draw string .5 9 Input: 'inputmsm
*'draw string .5 8.8 Input: 'inputmsm2
'set z 1 16'
'dynamic -var z temp u v'
*'define dtx = cdiff(temp,x)'
*'define dty = cdiff(temp,y)'
*'define dx = cdiff(lon,x)*3.1416/180'
*'define dy = cdiff(lat,y)*3.1416/180'
*tadv'=-1*( ('u'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v'*'dty'/'dy' )/6.37e6*84600'
*'set z 10'
*'set gxout shade2'
*'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
*'rgbset2'
*'d 'tadv
*'cbarn'
*'q dims'
*say result
*'allclose'
*return
'set mpdset hires'
'set lon 'lon
'set lat 22 40'
'set grads off'
'mul 2 1 1 1 -novpage'
*'set vpage 1 6 2 8'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
*'set zlog on'
*'set z 1 16'
* Temperature advection
'set gxout shade2'
'rgbset2'
'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
'd 'tadv'*86400'
*'cbarn 0.8 0 4.3 0.12'
'cbarn 1 1 7 5'
'q dims'
say result
'set gxout contour'
'set cthick 10'
'set ccolor 15'
'set clopts 1 10 0.2'
'set clevs -3 0 3 6 9 12 15 18 21 24 27 30 33'
'd u'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
*'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
*'d u'
* Meridional Circulation
'set gxout vector'
'set cthick 10'
'set ccolor 1'
'set arrscl 0.5 5'
'd skip(va,10,1);skip(-w*10,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'close 1'
'print 'figfile
say 'Fig file = 'figfile
'allclose'
::::::::::::::
advt_horiz.gs
::::::::::::::
'reinit'
* Default values
press=500
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/tadv.'press'.Jun.intense.w-0.356.eps'
figfile2='../Fig/tadv.'press'.Jun.weak.w-0.356.eps'
figfile3='../Fig/tadv.'press'.Jun.diff.w-0.356.eps'
'set display color white'
'c'
*'set parea 0.5 6 1 8'
title='Strong Updraft 'press'hPa'
*'set map 1 1 6'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm ;*msm201206.ctl'
'q file'
'set strsiz 0.1 0.15'
'draw string .5 10.5 Input: 'inputmsm
*'draw string .5 10.5 Input: 'inputmsm2
'set z 10' ;* 16'
'define dtx = cdiff(temp,x)'
'define dty = cdiff(temp,y)'
'define dx = cdiff(lon,x)*3.1416/180'
'define dy = cdiff(lat,y)*3.1416/180'
tadv'=-1*( ('u'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v'*'dty'/'dy' )/6.37e6*84600'
*'mul 2 1 1 1 -novpage'
'set mpdset hires'
'set gxout shade2'
'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
'rgbset2'
'd 'tadv
'q dims'
say result
'cbarn 1 0 4.3 0.4'
'set gxout contour'
'set cthick 10'
'set ccolor 1'
'set clopts 1 10 0.2'
'set clevs 3 6 9 12 15 18 21 24 27 30 33'
'd mag(u,v)'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
*'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
*'d u'
* Meridional Circulation
*'set cthick 10'
*'set arrscl 0.3 20'
*'d skip(v,10,1);skip(-w*200,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'print 'figfile
say 'Fig file = 'figfile
'allclose'
'c'
'set strsiz 0.1 0.15'
'draw string .5 10.5 Input: 'inputmsm2
title='Weak Updraft 'press'hPa'
'open 'inputmsm2
'set z 10' ;* 16'
'define dtx = cdiff(temp,x)'
'define dty = cdiff(temp,y)'
'define dx = cdiff(lon,x)*3.1416/180'
'define dy = cdiff(lat,y)*3.1416/180'
tadv'=-1*( ('u'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v'*'dty'/'dy' )/6.37e6*84600'
*'set vpage 1 6 2 8'
*'mul 2 1 2 1 -novpage'
*'set vpage 1 6 2 8'
'set grads off'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
* Temperature advection
'set gxout shade2'
'rgbset2'
'set clevs -10 -5 -2 -1 -0.5 0.5 1 2 5 10'
'd 'tadv
* Zonnal velocity
*'set gxout shade2'
*'color -levs -kind blue->white->red'
*'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
*'rgbset2'
*'set z 1 16'
*'set zlog on'
*'d u'
'set gxout contour'
'set cthick 10'
'set ccolor 1'
'set clopts 1 10 0.2'
'set clevs 6 9 12 15 18 21 24 27 30 33'
'd mag(u,v)'
* Meridional Circulation
*'set cthick 10'
*'set arrscl 0.3 20'
*'d skip(v,10,1);skip(-w*200,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'print 'figfile2
say 'Fig file = 'figfile2
'allclose'
title='d.tadv & d.U 'lon'E'
'c'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm
say; say 'Open ctl file: 'inputmsm2; say
'open ' inputmsm2
'set z 10' ;* 16'
'define dtx = cdiff(temp.1,x)'
'define dty = cdiff(temp.1,y)'
'define dx = cdiff(lon,x)*3.1416/180'
'define dy = cdiff(lat,y)*3.1416/180'
tadv1'=-1*( ('u.1'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v.1'*'dty'/'dy' )/6.37e6*84600'
'define dtx = cdiff(temp.2,x)'
'define dty = cdiff(temp.2,y)'
'define dx = cdiff(lon,x)*3.1416/180'
'define dy = cdiff(lat,y)*3.1416/180'
tadv2'=-1*( ('u.2'*'dtx')/(cos('lat'*3.1416/180)*'dx') + 'v.2'*'dty'/'dy' )/6.37e6*84600'
*'set lon 'lon
*'set lat 22 40'
*'set vpage 1 6 2 8'
*'mul 2 1 1 1 -novpage'
'set xlint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
*'set zlog on'
*'set z 1 16'
* Zonnal velocity
'set gxout shade2'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
'set clevs -3 -2 -1 -0.5 0.5 1 2 3'
'rgbset2'
'd tadv1-tadv2'
'cbarn 0.8 0 4.3 0.35'
'set gxout contour'
'set clevs 1 2 3 4 5'
'set cstyle 1'
'set ccolor 1'
'set cthick 7'
'set clopts 1 10 0.2'
'd u.1-u.2'
'set clevs -5 -4 -3 -2 -1'
'set cstyle 2'
'set ccolor 1'
'set cthick 7'
'd u.1-u.2'
'draw title 'title
'print 'figfile3
say 'Fig file = 'figfile3
'allclose'
::::::::::::::
convergence.gs
::::::::::::::
'reinit'
* Default values
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/convergence.Jun.w-0.356.eps'
figfile2='../Fig/convergence.diff.Jun.w-0.356.eps'
'set display color white'
'c'
'set strsiz 0.1 0.15'
'draw string 1 9 Input: 'inputmsm
'draw string 1 8.8 Input: 'inputmsm2
'open ' inputmsm
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
**Convergence field in intense updraft case
'mul 2 1 1 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set z 1'
'q dims'
linetmp=sublin(result,4)
press=subwrd(linetmp,6)
'set gxout shaded'
'color 2 20 2 -kind white->rainbow'
*'rgbset2'
'd -hdivg(u,v)*1.e6'
'set arrscl 0.3 10'
*'d skip(u,25,25);skip(v,25,25)'
*'cbarn 0.55 0 2.2 3.7'
'draw title Strong Updraft 'press'hPa'
'close 1'
**Convergence field in weak updraft case
'open ' inputmsm2
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
'mul 2 1 2 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'q dims'
linetmp=sublin(result,4)
press=subwrd(linetmp,6)
'set gxout shaded'
'color 2 20 2 -kind white->rainbow'
'd -hdivg(u,v)*1.e6'
'set arrscl 0.3 10'
*'d skip(u,25,25);skip(v,25,25)'
'cbarn 0.7 0 5.5 3.5';*横幅、縦幅、x軸、y軸
'draw title Weak Updraft 'press'hPa'
say 'Fig. file = 'figfile
'print 'figfile
'close 1'
'c'
say inputmsm
'open ' inputmsm
'q file'
say inputmsm2
'open ' inputmsm2
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
* 'mul 2 1 2 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'q dims'
linetmp=sublin(result,4)
press=subwrd(linetmp,6)
*'color 2 20 2 -kind white->rainbow'
'set z 1'
dconv'=-hdivg('u.1','v.1')+hdivg('u.2','v.2')'
dconv6'='dconv'*1.e6'
'set gxout shaded'
'set clevs -20 -15 -10 -5 5 10 15 20'
'rgbset2'
'd 'dconv6
'set arrscl 0.3 10'
*'d skip(u,25,25);skip(v,25,25)'
'cbarn 0.7 0 5.5 3.5';*横幅、縦幅、x軸、y軸
'draw title Strong - Weak 'press'hPa'
say 'Fig. file = 'figfile2
'print 'figfile2
'allclose'
::::::::::::::
diff_grad.ept.700.gs
::::::::::::::
'reinit'
* Default values
pres=700
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/diff.merid.grad.EPT.'pres'hPa.Jun.weak-0.356.eps'
'set display color white'
'c'
'open ' inputmsm
'q file'
'set mpdset hires'
*'set lon 121 137.5'
*'set lat 24 37.5'
'set grads off'
**largeomega veroticity
*'mul 2 1 1 1'
'set xlint 5'
'set ylint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'L=2.5*1.e6'
'Rv=461.0'
'p=700'
'set z 8'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp*pow((1000/p),po)'
'ept=temp*pow((1000/p),po) * exp1'
'set clskip 5'
'set cint 2'
'epty=cdiff(ept,y)'
'eptx=cdiff(ept,x)'
'xdiff=cdiff(lon,x)'
'ydiff=cdiff(lat,y)'
'eptdiff=sqrt(pow(epty/ydiff,2)+pow(eptx/xdiff,2))'
*'set clevs 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5'
'rgbset2'
*'d eptdiff'
*'d -epty'
*'cbarn'
*'set gxout contour'
*'d ept'
*'draw title Largeomega'
*'close 1'
**smallomega voroticity
'open ' inputmsm2
'q file'
'set mpdset hires'
*'set lon 121 137.5'
*'set lat 24 37.5'
'set grads off'
*'mul 2 1 2 1'
'set xlint 5'
'set ylint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
'L=2.5*1.e6'
'Rv=461.0'
'p=700'
'set z 8'
*mix
'denom=1.0/(temp.2-55.0)-log(rh.2/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp.2- 273.2)/(273.2*temp.2)'
'es = 6.11*exp(A*L/Rv)'
'e = rh.2 * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp.2*pow((1000/p),po)'
'ept2=temp.2*pow((1000/p),po) * exp1'
'set clskip 5'
'set cint 2'
'epty2=cdiff(ept2,y)'
'eptx2=cdiff(ept2,x)'
'xdiff2=cdiff(lon,x)'
'ydiff2=cdiff(lat,y)'
'eptdiff2=sqrt(pow(epty2/ydiff2,2)+pow(eptx2/xdiff2,2))'
'rgbset2'
'set clevs 0 0.5 1 1.5 2 2.5 3 3.5'
'set clevs 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5'
'd eptdiff-eptdiff2'
'set gxout contour'
'd ept-ept2'
*'cbarn 0.9 0 5.2 3.3';*横幅、縦幅、x軸、y軸
'cbarn'
*'draw title Smallomega'
'print 'figfile
::::::::::::::
dynamic.gs
::::::::::::::
function dyanmic(args)
* http://gradsaddict.blogspot.jp/2013/05/script-dynamic-calculates-temperature.html
'q dims'
xline=sublin(result,2)
yline=sublin(result,3)
zline=sublin(result,4)
tline=sublin(result,5)
lons=subwrd(xline,6)' 'subwrd(xline,8)
lats=subwrd(yline,6)' 'subwrd(yline,8)
if(subwrd(zline,7)='Z');levs=subwrd(zline,6);else;levs=subwrd(zline,6)' 'subwrd(zline,8);endif
time=subwrd(tline,6)
check=1
a=1
v=1
while(check=1)
line=subwrd(args,a)
if(line='');break;endif
if(line='-help')
help()
return
endif
if(line='-var');hgt=subwrd(args,a+1);tmp=subwrd(args,a+2);uwnd=subwrd(args,a+3);vwnd=subwrd(args,a+4);v=0;endif
a=a+1
endwhile
if(v=1)
say ''
say 'You have Chosen the Default variables for Height, Temperature, and Wind'
say 'To specify your variables, type "-var" then your variable names in the order of Height, Temperature, Uwind, Vwind'
say ''
'hgt=hgtprs'
'tmp=tmpprs'
'vwnd=vgrdprs'
'uwnd=ugrdprs'
endif
if(v!=1)
say ''
say 'You have Chosen to specify the variables for Height, Temperature, and Wind'
say 'You have specified the variables as:'
say 'Height: 'hgt
say 'Temp: 'tmp
say 'U-wind: 'uwnd
say 'V-wind: 'vwnd
'hgt='hgt
'tmp='tmp
'vwnd='vwnd
'uwnd='uwnd
endif
say 'Please wait while I calculate your variables'
'pi=3.14159'
'dtr=pi/180'
'a=6.37122e6'
'omega=7.2921e-5'
'g=9.8'
'R=287'
'define f=2*omega*sin(lat*dtr)'
'define p=lev*100'
say '...'
'dy=cdiff(lat,y)*dtr*a'
'dx=cdiff(lon,x)*dtr*a*cos(lat*dtr)'
'dhgtx=cdiff(hgt,x)'
'dhgty=cdiff(hgt,y)'
'define ug=-1*(g/f)*(dhgty/dy)'
'define vg=(g/f)*(dhgtx/dx)'
'define ua=uwnd-ug'
'define va=vwnd-vg'
say '....'
'dugx=cdiff(ug,x)'
'dugy=cdiff(ug,y)'
'dvgx=cdiff(vg,x)'
'dvgy=cdiff(vg,y)'
'dtdx=cdiff(tmp,x)/dx'
'dtdy=cdiff(tmp,y)/dy'
say '.....'
'define Q1=-1*(R/p)*(dugx/dx*dtdx + dvgx/dx*dtdy)'
'define Q2=-1*(R/p)*(dugy/dy*dtdx + dvgy/dy*dtdy)'
'define divq=hdivg(Q1,Q2)'
'define tadv=(-uwnd*dtdx-vwnd*dtdy)'
say '......'
'divg=hdivg(uwnd,vwnd)'
'vort=hcurl(uwnd,vwnd)'
'dvdx=cdiff(vort,x)/dx'
'dvdy=cdiff(vort,y)/dy'
'define vadv=(-uwnd*dvdx-vwnd*dvdy)'
say '.......'
'def1=cdiff(uwnd,x)/dx-cdiff(vwnd,y)/dy-vwnd*tan(dtr*lat)/a'
'def2=cdiff(vwnd,x)/dx+cdiff(uwnd,y)/dy+uwnd*tan(dtr*lat)/a'
'f1=-(dtdx*(divg+def1)+dtdy*(vort+def2))/2'
'f2=(dtdx*(vort-def2)-dtdy*(divg-def1))/2'
'fn=(dtdx*f1+dtdy*f2)/mag(dtdx,dtdy)'
'fs=(dtdx*f2-dtdy*f1)/mag(dtdx,dtdy)'
'fnx=-1*((dtdx*f1)/mag(dtdx,dtdy))*10e9'
'fny=-1*((dtdy*f2)/mag(dtdx,dtdy))*10e9'
say '........'
'define F=(fn+fs)*10e9'
say 'DONE!'
say ''
say 'The following variables have been defined for the dimensions:'
say 'Longitude: 'lons
say 'Latitude: 'lats
say 'Pressure Levels: 'levs
say 'Time: 'time
say '--------------------------------------------------------------------'
say 'Defined Variables: Variable Name Units '
say ' -Geostrophic Wind : ug,vg [m/s] '
say ' -Ageostrophic Wind: ua,va [m/s] '
say ' -Q-Vectors : Q1,Q1 [pa/m2/s] '
say ' -Temp Advection : tadv [K/s] '
say ' -Vort Advection : vadv [-] '
say ' -Frontogenesis : F [K/m/s]x10^9 '
say ' -Fn Vector : fnx,fny [K/m/s]x10^9 '
say ' -Deformation : def1,def2 [m] '
say '--------------------------------------------------------------------'
say 'To plot a variable simply type "d " then pick a Name from the list '
return
function help()
say '---------------------------------------------------'
say ' Dynamic v2.0 '
say '---------------------------------------------------'
say 'Usage:'
say 'Calculates advanced dynamic variables:'
say ' -Qvectors, Frontogenesis, advection, Geostrophic winds, etc.'
say ' -ALL variables must have MKS units, convert before running the script'
say ''
say 'Required: No required arguments'
say '---------------------------------------------------'
say ''
say 'Optional: -help - Pulls up Help Page'
say ' -var - Allows you to specify variable names'
say ' - Defaults: hgt=hgtprs, tmp=tmpprs, uwnd=ugrdprs, vwnd=vgrdprs'
say ''
say '---------------------------------------------------'
say ''
say 'Example: dynamic -var height tempk u v'
say ' calculates dynamic variables using the variables, height, tempk, u, v'
say
say 'Version 2.0: Orginially Developed Sept 2012'
say ' -Help page and arguments added May 2013 as part of v2.0'
**
***::::::::::::::
grad.ept.700.gs
::::::::::::::
'reinit'
* Default values
pres=700
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/merid.grad.EPT.'pres'hPa.Jun.weak-0.356.eps'
'set display color white'
'c'
'open ' inputmsm
'q file'
'set mpdset hires'
*'set lon 121 137.5'
*'set lat 24 37.5'
'set grads off'
**largeomega veroticity
'mul 2 1 1 1'
'set xlint 5'
'set ylint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'L=2.5*1.e6'
'Rv=461.0'
'p=700'
'set z 8'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp*pow((1000/p),po)'
'ept=temp*pow((1000/p),po) * exp1'
'set clskip 5'
'set cint 2'
'epty=cdiff(ept,y)'
'eptx=cdiff(ept,x)'
'xdiff=cdiff(lon,x)'
'ydiff=cdiff(lat,y)'
'eptdiff=sqrt(pow(epty/ydiff,2)+pow(eptx/xdiff,2))'
'set clevs 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5'
*'set clevs 0 1 2 3 4 5 6 7 8 9 10'
'rgbset2'
'd eptdiff'
*'d -epty'
*'cbarn'
'set gxout contour'
'd ept'
'draw title Intense Updraft'
'close 1'
**smallomega voroticity
'open ' inputmsm2
'q file'
'set mpdset hires'
*'set lon 121 137.5'
*'set lat 24 37.5'
'set grads off'
'mul 2 1 2 1'
'set xlint 5'
'set ylint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
'L=2.5*1.e6'
'Rv=461.0'
'p=700'
'set z 8'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp*pow((1000/p),po)'
'ept=temp*pow((1000/p),po) * exp1'
'set clskip 5'
'set cint 2'
'epty=cdiff(ept,y)'
'eptx=cdiff(ept,x)'
'xdiff=cdiff(lon,x)'
'ydiff=cdiff(lat,y)'
'eptdiff=sqrt(pow(epty/ydiff,2)+pow(eptx/xdiff,2))'
'rgbset2'
'set clevs 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5'
'd eptdiff'
'set gxout contour'
'd ept'
'cbarn 0.9 0 5.2 3.3';*横幅、縦幅、x軸、y軸
'draw title Weak Updraft'
'print 'figfile
::::::::::::::
u_crossec.gs
::::::::::::::
'reinit'
* Default values
lon=120
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/u_crossec.'lon'.Jun.w-0.356.eps'
figfile2='../Fig/u_crossec.'lon'.Jun.diff.w-0.356.eps'
'set display color white'
'c'
'set parea 0.5 6 1 8'
title='Strong Updraft 'lon'E'
*'set map 1 1 6'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm ;*msm201206.ctl'
'q file'
'set mpdset hires'
'set lon 'lon
'set lat 25 40'
'set grads off'
'mul 2 1 1 1 -novpage'
*'set vpage 1 6 2 8'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
*'set zlog on'
'set z 1 16'
* Zonnal velocity
'set gxout shade2'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
'rgbset2'
'd u'
* Meridional Circulation
*'set cthick 10'
*'set arrscl 0.3 20'
*'d skip(v,10,1);skip(-w*200,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'close 1'
title='Weak Updraft 'lon'E'
'open 'inputmsm2
'set lon 'lon
'set lat 25 40'
*'set vpage 1 6 2 8'
'mul 2 1 2 1 -novpage'
*'set vpage 1 6 2 8'
'set grads off'
'set xlint 5'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
* Zonnal velocity
'set gxout shade2'
*'color -levs -kind blue->white->red'
'set clevs 0 3 6 9 12 15 18 21 24 27 30 33'
'rgbset2'
'set z 1 16'
*'set zlog on'
'd u'
'cbarn 0.8 0 4.3 0.35'
* Meridional Circulation
*'set cthick 10'
*'set arrscl 0.3 20'
*'d skip(v,10,1);skip(-w*200,10,1)'
*EPT
*'L=2.5*1.e6'
*'Rv=461.0'
*'p=850'
*'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
*'tl=1.0/denom+55.0'
*'A = (temp- 273.2)/(273.2*temp)'
*'es = 6.11*exp(A*L/Rv)'
*'e = rh * es * 0.01'
*'r = (0.622 * e/ (p - e))*1.E3'
*'po=0.2854*(1.0 - 0.28*0.001*r)'
*'arg4 = 3.376/tl - 0.00254'
*'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
*'exp1=exp( arg4 * arg5 )'
*'pt=temp*pow((1000/p),po)'
*'ept=temp*pow((1000/p),po) * exp1'
*'set gxout contour'
*'set clevs 240 250 260 270 280 290 300 310 320 330 340 350 360 370'
*'d ept'
'draw title 'title
'print 'figfile
say 'Fig file = 'figfile
'allclose'
title='d.U 'lon'E'
'c'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm
say; say 'Open ctl file: 'inputmsm2; say
'open ' inputmsm2
'set lon 'lon
'set lat 25 40'
*'set vpage 1 6 2 8'
'mul 2 1 1 1 -novpage'
'set xlint 5'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
*'set zlog on'
'set z 1 16'
* Zonnal velocity
'set gxout shade2'
*'color -levs 0 5 15 20 25 30 -kind blue->white->red'
'set clevs -10 -5 -3 -2 -1 1 2 3 5 10'
'rgbset2'
'd u.1-u.2'
'cbarn 0.8 0 4.3 0.35'
'draw title 'title
'print 'figfile2
say 'Fig file = 'figfile2
'allclose'
::::::::::::::
uv.gs
::::::::::::::
'reinit'
* Default values
press=600
casename='w-0.356'
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/uv.'press'hPa.Jun.'casename'.eps'
'set display color white'
'c'
'open ' inputmsm
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
title='Strong Updraft 'press'hPa'
'mul 2 1 1 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set z 9' ;* 600 hPa
'set gxout shade2'
'set clevs 4 6 8 10 12 14 16 18 20 25'
'rgbset2'
'd mag(u,v)'
'set gxout vector'
'set arrscl 0.3 30'
'd skip(u,10,10);skip(v,10,10)'
'draw title 'title
'close 1'
title='Weak Updraft 'press'hPa'
'open ' inputmsm2
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
'mul 2 1 2 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set z 9' ;* 600 hPa
'set gxout shade2'
'set clevs 4 6 8 10 12 14 16 18 20 25'
'rgbset2'
'd mag(u,v)'
'set arrscl 0.3 200'
'set gxout vector'
'set arrscl 0.3 30'
'd skip(u,10,10);skip(v,10,10)'
'draw title 'title
'cbarn 0.8 0 4.3 0.35'
'print 'figfile
say
say 'Fig file = 'figfile; say
'allclose'
::::::::::::::
w_crossec.gs
::::::::::::::
'reinit'
* Default values
lon=125
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/w_crossec.'lon'E.Jun.w-0.356.eps'
'set display color white'
'c'
'set map 1 1 6'
say; say 'Open ctl file: 'inputmsm; say
'open ' inputmsm ;*msm201206.ctl'
'q file'
'set mpdset hires'
'set lon 'lon
'set lat 27 31'
'set grads off'
'mul 2 1 1 1'
'set xlint 1'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set z 1 12'
'L=2.5*1.e6'
'Rv=461.0'
'p=850'
*ept
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp*pow((1000/p),po)'
'ept=temp*pow((1000/p),po) * exp1'
*'set clevs 270 280 290 300 310 320 330 340 350 360 370'
*'set clevs 340 345 350 355 360 365 370 375 380 385'
'set clevs 0 1 2 3 4 5 6 7 8 9 10 11 12'
*'rgbset2'
'set gxout shaded'
'rgbset2'
'd -w*10'
*'cbarn'
'set cthick 5.5'
'set arrscl 0.3 20'
*'d skip(u,10,10);skip(v,10,10)'
'draw title Strong Updraft 'lon'E'
'close 1'
**smallomega
'open 'inputmsm2
'set lon 'lon
'set lat 27 31'
'mul 2 1 2 1'
'set xlint 1'
*'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'set clevs 270 280 290 300 310 320 330 340 350 360 370'
*'set clevs 340 345 350 355 360 365 370 375 380 385'
'set clevs 0 1 2 3 4 5 6 7 8 9 10 11 12'
'set z 1 12'
'L=2.5*1.e6'
'Rv=461.0'
'p=850'
*ept
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'po=0.2854*(1.0 - 0.28*0.001*r)'
'arg4 = 3.376/tl - 0.00254'
'arg5= r*(1.0 + 0.81 * 1.0E-3*r)'
'exp1=exp( arg4 * arg5 )'
'pt=temp*pow((1000/p),po)'
'ept=temp*pow((1000/p),po) * exp1'
'set gxout shaded'
'rgbset2'
'd -w*10'
'set cthick 5.5'
'set arrscl 0.3 20'
*'d skip(u,10,10);skip(v,10,10)'
'draw title Weak Updraft 'lon'E'
'cbarn'
'print 'figfile
say
say 'Fig file = 'figfile ; say
'allclose'
::::::::::::::
wvf.conv.gs
::::::::::::::
'reinit'
* Default values
pres=950
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/wvf.conv.'pres'hPa.Jun.weak-0.356.eps'
'set display color white'
'c'
'open ' inputmsm
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
**largeomega veroticity
'mul 2 1 1 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'set z 3'
'L=2.5*1.e6'
'Rv=461.0'
'p=950'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
*'color 17 20 -kind white->blue'
'set clevs 0 10 20 30 40 50'
*'set clevs -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30'
'rgbset2'
*'d v*r'
'd -hdivg(u*r,v*r)*1.e5'
'set gxout contour'
*'set clskip 10'
'set clskip 2'
'set cint 0.5'
'set cmin 17'
*'d v*r'
*'cbarn'
*'set gxout contour'
'set arrscl 0.3 200'
'd skip(r*u,10,10);skip(r*v,10,10)'
*'d v'
*'set arrscl 0.3 10'
*'d skip(u,25,25);skip(v,25,25)'
*'cbarn 0.55 0 2.2 3.7'
'draw title Intense Updraft'
'close 1'
**smallomega voroticity
'open ' inputmsm2
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
'mul 2 1 2 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'set z 3'
'L=2.5*1.e6'
'Rv=461.0'
'p=950'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'rgbset2'
*'set clevs'
'set clevs 0 10 20 30 40 50'
*'d v*r'
'd -hdivg(u*r,v*r)*1.e5'
*'color 17 20 -kind white->blue'
'set gxout contour'
*'set clab '
'set clskip 2'
'set cint 0.5'
'set cmin 17'
'set arrscl 0.3 200'
'd skip(r*u,10,10);skip(r*v,10,10)'
*'set gxout contour'
*'d v'
'cbarn 0.9 0 5.2 3.0';*横幅、縦幅、x軸、y軸
'draw title Weak Updraft'
'print 'figfile
::::::::::::::
wvf950.gs
::::::::::::::
'reinit'
* Default values
pres=950
inputmsm='../output/composite.Jun.intense-0.356.ctl'
inputmsm2='../output/composite.Jun.weak-0.356.ctl'
figfile='../Fig/wvf.'pres'hPa.Jun.weak-0.356.eps'
'set display color white'
'c'
'open ' inputmsm
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
**largeomega veroticity
'mul 2 1 1 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'L=2.5*1.e6'
'Rv=461.0'
'p=950'
'set z 3'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
*'color 17 20 -kind white->blue'
'set arrscl 0.3 200'
'set clevs 13 13.6 14.2 14.8 15.4 16 16.6 17.2 17.8 18.4 19'
'rgbset2'
*'color 13 19 -kind white->blue'
'ur=u*r'
'vr=v*r'
'd r'
'set arrscl 0.3 200'
'd skip(ur,10,10);skip(vr,10,10)'
'set gxout contour'
*'set clskip 10'
'set clskip 2'
'set cint 0.5'
'set cmin 17'
*'d v*r'
*'cbarn'
*'set gxout contour'
*'d v'
*'set arrscl 0.3 10'
*'d skip(u,25,25);skip(v,25,25)'
*'cbarn 0.55 0 2.2 3.7'
'draw title Intense Updraft'
'close 1'
**smallomega voroticity
'open ' inputmsm2
'q file'
'set mpdset hires'
'set lon 121 137.5'
'set lat 24 37.5'
'set grads off'
'mul 2 1 2 1'
'set xlint 3'
'set ylint 2'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set gxout shaded'
*'d v'
*'set gxout contour'
'L=2.5*1.e6'
'Rv=461.0'
'p=950'
'set z 3'
*mix
'denom=1.0/(temp-55.0)-log(rh/100.0)/2840.0'
'tl=1.0/denom+55.0'
'A = (temp- 273.2)/(273.2*temp)'
'es = 6.11*exp(A*L/Rv)'
'e = rh * es * 0.01'
'r = (0.622 * e/ (p - e))*1.E3'
'rgbset2'
*'set clevs'
'set arrscl 0.3 200'
*'color 13 19 -kind white->blue'
'set clevs 13 13.6 14.2 14.8 15.4 16 16.6 17.2 17.8 18.4 19'
'rgbset2'
'd r'
'ur=u*r'
'vr=v*r'
'set arrscl 0.3 200'
'd skip(ur,10,10);skip(vr,10,10)'
*'color 10 20 -kind white->blue'
'set clskip 2'
'set cint 0.5'
'set cmin 17'
*'d r*v'
*'set gxout contour'
*'d v'
'cbarn 0.9 0 5.2 3.3';*横幅、縦幅、x軸、y軸
'draw title Weak Updraft'
'print 'figfile