#/bin/sh # # 1st release: Jan. 05, 2000 # 2nd release: Mar. 24, 2004 # 3rd release: Feb. 27, 2007 # M. HAYASAKI # prog_path=$PWD my_name=$PWD/$0 exec_prog=SATWATER ofile51=$HOME/fort/moist/sat_water_pres.dat ofile52=$HOME/fort/moist/sat_water_pres_Bolton1980MWR.dat cd $workdir mkdir $$ ; cd $$ cat > prog.f90 << EOF PROGRAM main IMPLICIT NONE INTEGER(4), PARAMETER :: np=18 REAL(4), PARAMETER :: base_kelvin=273.15 INTEGER(4) :: idegC,ip REAL(4) :: spec_hum, mix_ratio, sat_vapor_pres, vapor_pres REAL(4) :: RH, celsius, absdeg, pres, temp, degC, power REAL(4) :: temp_LCL REAL(4) :: q, r, es, e REAL(4) :: theta, theta_e, sat_theta_e CHARACTER(LEN=200) :: my_name, header01, header02, header03 CHARACTER(LEN=200) :: ofile51, ofile52 REAL(4) :: pr(np) DATA pr/1013.,1000.,925.,850.,700.,600.,500.,400.,300., & 250.,200.,150.,100.,70.,50.,30.,20.,10./ my_name='${my_name}' ! +++ open output file ofile51='${ofile51}' OPEN(51, FILE=TRIM(ofile51)) ofile52='${ofile52}' OPEN(52, FILE=TRIM(ofile52)) header01=TRIM(my_name) header02=' RH, Pres, T(K), T(degC), mix_ratio, shum, Es ' header03=' RH, Pres, T(K), T(degC), mix_ratio, shum, Es , Temp_LCL' WRITE(51,'("# This data had made by: ",a)') TRIM(header01) WRITE(52,'("# This data had made by: ",a)') TRIM(header01) WRITE(51,'("# ",a)') TRIM(header02) WRITE(52,'("# ",a)') TRIM(header03) ! +++ set constant value => RH=100% (fully saturated condition) ! RH = 50.0 RH = 100.0 WRITE(6,'(" Constant RH = ", f10.1)') RH DO ip = 1, 10 pres = pr(ip) DO idegC = -10, 10 ! from -50 degC to 50 degC celsius = 5.0 * FLOAT(idegC) absdeg = celsius + base_kelvin CALL moisture(absdeg,RH,pres, spec_hum,mix_ratio,sat_vapor_pres) q = spec_hum * 1.0E+3 ! specific humidity (g/kg) r = mix_ratio * 1.0E+3 es = sat_vapor_pres IF( idegC == 0 ) THEN write( 6,'(2f10.1, 2f10.2,3f10.4)') RH,pres,absdeg,celsius,r,q,es ENDIF CALL robitzsch_etheta(absdeg,RH,pres, theta,theta_e,sat_theta_e,mix_ratio) r = mix_ratio * 1.0E+3 write(51,'(2f10.1, 2f10.2,3f10.4, 10x, 3f10.2)') & RH,pres,absdeg,celsius,r,q,es, & theta, theta_e, sat_theta_e ! +++ Using equations in Bolton (1980, MWR) CALL bolton_etheta(absdeg,RH,pres, theta,theta_e,sat_theta_e, & spec_hum, mix_ratio, sat_vapor_pres, temp_LCL) q = spec_hum * 1.0E+3 ! specific humidity (g/kg) r = mix_ratio * 1.0E+3 es = sat_vapor_pres write(52,'(2f10.1, 2f10.2,3f10.4,f10.2, 3f10.2)') & RH,pres,absdeg,celsius,r,q,es, temp_LCL, & theta, theta_e, sat_theta_e ENDDO ! idegC ENDDO ! ip STOP END PROGRAM main ! ---------------------------------------------------------------------- subroutine moisture(temp,RH,pres, spec_hum,mix_ratio,sat_vapor_pres) IMPLICIT NONE REAL(4), INTENT(IN) :: temp, RH, pres REAL(4), INTENT(OUT) :: spec_hum, mix_ratio, sat_vapor_pres REAL(4) :: q, r, es, e REAL(4) :: degC, power, vapor_pres ! ### copied from "robitzsh_e-theta" ! spec_hum : specific humidity (q) [g/kg] ! q = 1000*(0.622*E/pr) ! mix_ratio : mixing ratio (r) [g/kg] ! r = 1000*(0.622*E/(pr - E)) ! 0.622=ratio of molecular mass; water/dry-air ! vapor_pres : Vapour pressure (E) [hPa] ! E = RH*Es/100 ! ! *** Tetens's experimental equation (1930) *** ! sat_vapor_pres => Es ! Es = 6.1078*10**((T * A)/(T + B)) ! A = 7.5 } for use in ! B = 237.3 } vapor pressure with respect to water ! * A = 9.5 } for use in vapor pressure ! B = 265.5 } with respect to ice REAL(4), PARAMETER :: p00 = 1000.0 degC = temp - 273.15 ! +++ case selection by Temperature ! !! DO NOT proper for applying to real atmosphere ! Because of existence of "super-cooled water" !! if(c.gt.0) then !! power=(c*7.5)/(c+237.3) !! else !! power=(c*9.5)/(c+265.5) !! endif ! For simplicity, ! "power" set to the constant value in the case of ! "water surface" condition. ! (Of cource, this case is NOT proper for applying.) power = (degC*7.5)/(degC+237.3) ! 0.622 = Mw / Md (ratio of molecular weight between water and dry air) sat_vapor_pres = 6.1078 * (10**power) vapor_pres = sat_vapor_pres * ( RH * (1.0E-2) ) mix_ratio = 0.622 * ( vapor_pres / (pres - vapor_pres) ) spec_hum = 0.622 * ( vapor_pres / pres ) RETURN END SUBROUTINE moisture ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- EOF cat $HOME/fort/sub/robitzsch_e-theta.sub >> prog.f90 cat $HOME/fort/moist/bolton_e-theta.sub >> prog.f90 $FC -o ${exec_prog} prog.f90 size ${exec_prog} #cd .. ; \rm -r $$ #exit time nice -10 ./${exec_prog} cd .. ; rm -r $$ exit ### GMT plotting ipath=$HOME/fort/moist opath=$HOME/fort/moist ### CASE: Robitzsch (1928) #ifile=$ipath/sat_water_pres.dat #psfile=$opath/sat_water_pres.ps ### CASE: Bolton (1980, MWR) ifile=$ipath/sat_water_pres_Bolton1980MWR.dat psfile=$opath/sat_water_pres_Bolton1980MWR.ps rr="-R-51/41/-5/50" jj="-JX5.5/4.0" gmtset UNIX_TIME_POS 0.5/0.5 gmtset MEASURE_UNIT inch # 4.25 10.5 15 0 5 2 T & pressure dependency of specific humidity [g/kg] cat > text.in << EOF 4.25 10.5 15 0 5 2 T & pressure dependency of mixing ratio [g/kg] 4.25 10.1 20 0 4 2 Tetens (1930) experimental approximation form 4.25 9.7 15 0 4 2 !! saturated(RH=100%) & water surface case !! 1.5 3.5 15 0 1 1 pressure: 1000, 850, 700, 500hPa 1.5 3.2 12 0 1 1 circle, triangle, square, open circle 1.5 2.5 10 0 0 1 This approximation differs from real values in warmer condition, 1.5 2.3 10 0 0 1 but its differences is very small (0.2 g/kg in 1000hPa,35degC)! 7.7 0.5 10 0 5 3 minus value represents "super cooled water" 7.7 0.2 10 0 5 3 M. HAYASAKI EOF # *** PLOT *** pstext text.in -R0/8.26/0/11.69 -Jx1 -U/0.5/0.2/"${my_name}" -N -Y0.0 -X0.0 -P -K > $psfile psbasemap $rr $jj -Ba10f5g5:"@+\312@+C":/a5f1g5:"[g/kg]":WeSn -O -K -X1.5 -Y5.2 >> $psfile gawk '$2 == 1000.0 {print $4,$5}' $ifile > t.dat gawk '{print $1,$2}' t.dat | psxy -R -JX -O -K >> $psfile gawk '{print $1,$2}' t.dat | psxy -R -JX -Sc0.05 -G0 -L -W2 -O -K >> $psfile # gawk '$2 == 850.0 {print $4,$5}' $ifile > t.dat gawk '{print $1,$2}' t.dat | psxy -R -JX -O -K >> $psfile gawk '{print $1,$2}' t.dat | psxy -R -JX -St0.05 -G0 -L -W2 -O -K >> $psfile # gawk '$2 == 700.0 {print $4,$5}' $ifile > t.dat gawk '{print $1,$2}' t.dat | psxy -R -JX -O -K >> $psfile gawk '{print $1,$2}' t.dat | psxy -R -JX -Ss0.05 -G0 -L -W2 -O -K >> $psfile # gawk '$2 == 500.0 {print $4,$5}' $ifile > t.dat gawk '{print $1,$2}' t.dat | psxy -R -JX -O -K >> $psfile gawk '{print $1,$2}' t.dat | psxy -R -JX -Sc0.05 -G255 -L -W2 -O -K >> $psfile psbasemap -R -JX -B -O >> $psfile \rm text.in t.dat .gmtdefaults* .gmtcommands* exit ### CHANGES 24Mar2004: convert from C-shell into Bourne-shell form 05Jan2000: 1st release of this program