2018-09-25_11-03
/work3/am/2018.KUROSHIO.CONVECTION/SNAPSHOT/SATURATED.EPT
am@localhost
$ sat.pseudoadiabats.run.sh
-rwxrwxr-x 1 am am 826K 9月 25 11:11 sat.pseudoadiabats.exe
-rw-rw-r-- 1 am am 676 9月 25 11:11 sat.pseudoadiabats.table.txt
TABLE 78 of Smithsonian Meteorological Tables (p. 319)
https://repository.si.edu/handle/10088/23746
この表は飽和相当温位(列)と気温(行)に対し、気圧が何hPaになるかを示している。
30.0 28.0 26.0 24.0 20.0
386.2 373.8 362.6 352.8 335.5
30 1000 1063.2 -999.0 -999. -999.
28 940 1000 1061.5 -999. -999.
26 883 941 1000 1059.2 -999.
24 830 885 942 1000 -999.
22 779 833 888 944 1055.6
20 733 785 837 892 1000
16 650 698 747 798 899
10 544 588 632 678 770
4 459 499 540 581 665
2 435 473 513 553 635
0 412.4 449 488 527 606
-2 391.4 427 464 502 579
-6 353.8 387.8 423 459 531
-10 321.4 353.4 386.6 420 489
-20 257.7 285.5 314.3 342.9 403
この表は、上の表で与えられた気温と気圧の値に対して、飽和相当温位が何Kになるかを示している。各列の飽和相当温位がほぼ同じ値になっていれば正解（計算誤差があるので全く同じ値にはならない）
REPRODUCTION OF TABLE 78 OF Smithsonian Meteorological Tables (p. 319)
30.0 28.0 26.0 24.0 20.0
386.2 373.8 362.6 352.8 335.5
30.0 386.3 373.9 0.0 0.0 0.0
28.0 386.0 373.8 362.7 0.0 0.0
26.0 386.0 373.7 362.7 352.8 0.0
24.0 385.9 373.8 362.7 352.7 0.0
22.0 386.1 373.8 362.7 352.7 335.7
20.0 386.0 373.7 362.8 352.6 335.6
16.0 385.7 373.5 362.6 352.5 335.5
10.0 385.8 373.4 362.6 352.5 335.5
4.0 385.8 373.4 362.3 352.4 335.4
2.0 385.6 373.5 362.3 352.4 335.3
0.0 385.6 373.5 362.2 352.3 335.4
-2.0 385.5 373.4 362.4 352.4 335.4
-6.0 385.4 373.2 362.1 352.1 335.3
-10.0 385.3 373.1 362.0 352.2 335.1
-20.0 385.2 372.9 361.8 352.1 335.0
2018-09-25_12-21
/work3/am/2018.KUROSHIO.CONVECTION/SNAPSHOT/SATURATED.EPT
am@localhost
$ srcdump.sh sat.pseudoadiabats.run.sh sat.pseudoadiabats.f90 ept_sept.f90 >l.txt
#------------------------------
# List of the following files:
#------------------------------
sat.pseudoadiabats.run.sh
sat.pseudoadiabats.f90
ept_sept.f90
#------------------------------
# Machine info
#------------------------------
localhost
/work3/am/2018.KUROSHIO.CONVECTION/SNAPSHOT/SATURATED.EPT
Tue Sep 25 12:22:23 JST 2018
#======================
# sat.pseudoadiabats.run.sh
#======================
#!/bin/bash
fc=ifort
opt=" -CB -traceback -fpe0"
exe=$(basename $0 .run.sh).exe
src=$(basename $0 .run.sh).f90
nml=$(basename $0 .run.sh).nml
listsrc="$src ept_sept.f90"
for tmp in $listsrc; do
if [ ! -f $tmp ]; then
echo
echo "ERROR in $0 : NO SUCH FILE, $tmp"
echo
exit 1
fi
done
$fc -o $exe $listsrc $opt
if [ $? -ne 0 ]; then
echo
echo COMPLIE ERROR!
echo
exit 1
else
echo
ls -lh $exe
echo
fi
# TABLE 78 OF Smithsonian Meteorological Tables (p. 319)
septcheck_file=$(basename $0 .run.sh).table.txt
cat<<END>$septcheck_file
30.0 28.0 26.0 24.0 20.0
386.2 373.8 362.6 352.8 335.5
30 1000 1063.2 -999.0 -999. -999.
28 940 1000 1061.5 -999. -999.
26 883 941 1000 1059.2 -999.
24 830 885 942 1000 -999.
22 779 833 888 944 1055.6
20 733 785 837 892 1000
16 650 698 747 798 899
10 544 588 632 678 770
4 459 499 540 581 665
2 435 473 513 553 635
0 412.4 449 488 527 606
-2 391.4 427 464 502 579
-6 353.8 387.8 423 459 531
-10 321.4 353.4 386.6 420 489
-20 257.7 285.5 314.3 342.9 403
END
echo
ls -lh $septcheck_file
cat $septcheck_file
echo
nl=$(wc -l ${septcheck_file}|awk '{if (NR==1) print $1}')
nl=$(expr $nl - 2)
nc=$(awk '{if(NR==3)print NF}' ${septcheck_file})
nc=$(expr $nc - 1)
echo "nl = $nl nc = $nc"
cat<<EOF>$nml
¶
septcheck_file="${septcheck_file}"
nl=$nl
nc=$nc
&end
EOF
echo
ls -lh $nml
cat $nml
echo
$exe <$nml
if [ $? -ne 0 ];then
echo
echo "ERROR in $exe !"
echo
exit 1
fi
exit 0
#----------------------
# End of sat.pseudoadiabats.run.sh
#----------------------
#======================
# sat.pseudoadiabats.f90
#======================
program table_ept_sept
real ept,sept,tk,p,RH
real,allocatable::tc(:),pha(:,:),septout(:,:)
real,allocatable::tcsfc(:),septsfc(:)
character(len=100)::septcheck_file
namelist /para/septcheck_file,nl,nc
read(*,nml=para)
print '(A)',trim(septcheck_file)
print '(A,i5)','nl = ',nl
print '(A,i5)','nc = ',nc
debug=0
RH=100.0 !DUMMY FOR CALC EPT (EPT IS NOT USED.)
allocate(tc(nl),pha(nl,nc),septout(nl,nc))
allocate(tcsfc(nc),septsfc(nc))
open(11,file=septcheck_file,action="read")
read(11,*)(tcsfc(i),i=1,nc)
read(11,*)(septsfc(i),i=1,nc)
do l=1,nl
read(11,*)tc(l),(pha(l,m),m=1,nc)
enddo
do l=1,nl
do m=1,nc
if(pha(l,m) < 0.0)then
septout(l,m)=0.0
else
tk=tc(l)+273.15
p=pha(l,m)
print *,tk,p
call ept_sept(ept, sept, es, e, r, tk, p, RH, debug)
septout(l,m)=sept
endif
end do !m
end do !n
print *
print '(A)','REPRODUCTION OF TABLE 78 OF Smithsonian Meteorological Tables (p. 319)'
print *
print '(7x,10f7.1)',(tcsfc(i),i=1,nc)
print '(7x,10f7.1)',(septsfc(i),i=1,nc)
do l=1,nl
print '(10f7.1)', tc(l),(septout(l,m),m=1,nc)
end do
print*
end program table_ept_sept
#----------------------
# End of sat.pseudoadiabats.f90
#----------------------
#======================
# ept_sept.f90
#======================
subroutine ept_sept(ept, sept, es, e, r, tk, p, RH, debug)
! Author: am
!
! EPT
! Bolton, D., 1980: The computation of equivalent potential
! Temperature, Monthly Weather Review, 108, 1046-1053.
! url: http://www.rsmas.miami.edu/users/pzuidema/Bolton.pdf
!
! "The most accurate formula to data is from Bolton (1980).
! It computes EPT to accuracy, relative to numerical
! solutions of the governing equation, better than 0.02 K."
! Davies-Jones (MWR, 2009)
! SATURATED EPT
! https://unidata.github.io/MetPy/latest/api/generated/
! metpy.calc.saturation_equivalent_potential_temperature.html
implicit none
real,intent(out)::ept,sept,es,e,r
! ept: equivalent potential temperature (kelvin)
!sept: saturated equivalent potential temperature (kelvin)
! es: saturated water vapor pressure (hPa)
! e: water vapor pressure (hPa)
! r: water-vapor mixing ratio (g/kg)
real,intent(in)::tk,p,RH
! tk: absolute temperature (kelvin)
! p : pressure (hPa)
! RH : relative humidity (%)
integer,intent(in)::debug
real u,tc
! u : relative humidity (%)
! tc: temperature (degree Celsius)
real,parameter::p0=1000.0
! p0: pressue (hPa)
! Shimizus
real,parameter::L=2.5*1.e6
real,parameter::Rv=461.0
! L : latent heat
! Rv: gas constant
real tl
! tl: Absolute temperature at lifting condensation level (LCL)
! (kelvin)
real pt
! pt : potential temparature
real pow,arg11,arg12,exp1, denom
real A
real numer, C
real rs
!rs : saturated mixing ratio (g/kg)
real,parameter::kappa=0.2854
real arg21,arg22
u=RH
denom=1.0/(tk - 55.0) - log(u/100.0)/2840.0
tl = 1.0/denom + 55.0
! CHECK VALUES FOR tl
! 1/(1/(280-55))+55=280 (RH=100%)
! 1/(1/(280-55)-(log(50/100)/2840))+55=274.758 (RH=50%)
!A = (tk- 273.2)/(273.2*tk)
!es = 6.11*exp(A*L/Rv)
tc=tk-273.15
! Eq.(10) of B80 (p.1047)
es=6.112*exp((17.67*tc)/(tc+243.5)) !hPa
e=RH/100.0 * es !hPa
r = (0.622 * e/ (p - e))*1.E3 !g/kg
! Eq.(43) of B80. (p.1052)
pow=0.2854*(1.0 - 0.28*0.001*r)
pt=tk*(p0/p)**pow
arg11 = 3.376/tl - 0.00254
arg12 = r*(1.0 + 0.81 * 1.0E-3*r)
exp1=exp( arg11 * arg12 )
ept=pt * exp1
! SATURATED EPT
! https://unidata.github.io/MetPy/latest/api/generated/
! metpy.calc.saturation_equivalent_potential_temperature.html
!pt=tk*(p0/p)**kappa
!Eq.(39) of B80 (p.1052)
rs=(0.622 * es/ (p - es))*1.E3 !g/kg
pow=0.2854*(1.0 - 0.28*0.001*rs)
pt=tk*(p0/p)**pow
arg21=3.376/tk - 0.00254
arg22 = rs*(1.0 + 0.81 * 1.0E-3*rs)
!arg21=3036./tk-1.78
!arg22=r*1.E-3*(1.0+0.448*r*1.E-3)
sept=pt*exp(arg21 * arg22)
if(debug/=0)then
print *
print '(A)','CHECK:'
print '(A,f10.2)','tc = ',tc
print '(A,f10.2)','es = ',es
print '(A,f10.2)','e = ',e
print '(A,f10.2)','tl = ',tl
print '(A,f10.2)','rs = ',rs
print '(A,f10.2)','r = ',r
print '(A,f10.2)','pt = ',pt
end if
!
! Eq. (6.3) of Davies-Jones (MWR, 2009)
! numer = (2.771*1.E6 - 1109.0*(tl - 273.15))*r*1.E-3
! denom = 1005.7*tl
! ept=pt*exp(numer/denom)
! Eq. (2.5) of Davies-Jones (MWR, 2009)
! ept = pt*exp((2.690*1.E6 * 1.0E-3 * r)/(1005.7*tl) )
end subroutine ept_sept
#----------------------
# End of ept_sept.f90
#----------------------
NCL
2019-04-18_16-42
/work05/manda/WRF.POST/RIP47.POST/SAMPLE.ALONG.TRAJ/SEPT
$ runncl.sh ncl_ept_sept.ncl 30 1000
tc p es rs sept
30.0 C 1000.0 hPa 42.5 hPa 27.58 g/kg 386.3 K
$ runncl.sh ncl_ept_sept.ncl 0 412.4
tc p es rs sept
0.0 C 412.4 hPa 6.1 hPa 9.36 g/kg 385.6 K
$ runncl.sh ncl_ept_sept.ncl 28.0 1000.0
tc p es rs sept
28.0 C 1000.0 hPa 37.8 hPa 24.44 g/kg 373.8 K
$ runncl.sh ncl_ept_sept.ncl 0.0 449
tc p es rs sept
0.0 C 449.0 hPa 6.1 hPa 8.58 g/kg 373.5 K
$ runncl.sh ncl_ept_sept.ncl 20.0 1000
tc p es rs sept
20.0 C 1000.0 hPa 23.4 hPa 14.88 g/kg 335.6 K
$ runncl.sh ncl_ept_sept.ncl 0.0 606
tc p es rs sept
0.0 C 606.0 hPa 6.1 hPa 6.34 g/kg 335.4 K
$ srcdump.sh ncl_ept_sept.ncl /home/manda/mybin/runncl.sh
#------------------------------
# List of the following files:
#------------------------------
ncl_ept_sept.ncl
/home/manda/mybin/runncl.sh
#------------------------------
# Machine info
#------------------------------
/work05/manda/WRF.POST/RIP47.POST/SAMPLE.ALONG.TRAJ/SEPT
Thu Apr 18 16:45:27 JST 2019
#======================
# ncl_ept_sept.ncl
#======================
scriptname_in = getenv("NCL_ARG_1")
tc = tofloat(getenv("NCL_ARG_2"))
p = tofloat(getenv("NCL_ARG_3"))
if (tc .lt. -100.0 .or. tc .gt. 100.0)then
print("ERROR tc = "+tc)
exit
end if
if (p .lt. 1.0 .or. p .gt. 2000.0)then
print("ERROR p = "+p)
exit
end if
;tc=30.0 ;C
;p=1000.0 ;hPa
p0=1000.0 ;hPa
tk=tc+273.15 ;K
; Eq.(10) of B80 (p.1047)
es=6.112*exp((17.67*tc)/(tc+243.5)) ;hPa
;e=RH/100.0 * es ;hPa
;r = (0.622 * e/ (ph - e))*1000.0 ;g/kg
;qv(0,0,0) = r/1000.0 ;kg/kg
; SATURATED EPT
; https://unidata.github.io/MetPy/latest/api/generated/
; metpy.calc.saturation_equivalent_potential_temperature.html
rs=(0.622 * es/ (p - es))*1.E3 ;g/kg
pow=0.2854*(1.0 - 0.28*0.001*rs)
pt=tk*(p0/p)^pow
arg21=3.376/tk - 0.00254
arg22 = rs*(1.0 + 0.81 * 1.0E-3*rs)
sept=pt*exp(arg21 * arg22)
print(" tc p es rs sept")
print(\
sprintf("%7.1f C",tc)+\
sprintf("%7.1f hPa",p)+\
sprintf("%7.1f hPa",es)+\
sprintf("%7.2f g/kg",rs)+\
sprintf("%7.1f K",sept)\
)
#----------------------
# End of ncl_ept_sept.ncl
#----------------------
#======================
# /home/manda/mybin/runncl.sh
#======================
#!/bin/bash
#
# Universal wrapper script for ncl.
# Pass arguments from the command line to environment variables
#
# version 0.1, Thierry Corti, C2SM ETH Zurich
#
E_BADARGS=65
if [ ! -n "$1" ]
then
echo "Usage: `basename $0` script.ncl argument1 argument2 etc."
exit $E_BADARGS
fi
# save number of arguments to environment variable NCL_N_ARG
export NCL_N_ARGS=$#
# save command line arguments to environment variable NCL_ARG_#
for ((index=1; index<=$#; index++))
do
eval export NCL_ARG_$index=\$$index
done
# run ncl
ncl -Q -n $1
#----------------------
# End of /home/manda/mybin/runncl.sh
#----------------------