Skip to content

Adds print statements for Murphy Koop and added new tests #159

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Adds print statements for Murphy Koop and added new tests #159

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c595410
Added print statements for Murphy Koop and added new tests
singhbalwinder Jul 7, 2020
f5a1e9b
fixes a comment
singhbalwinder Jul 7, 2020
f74099a
data arrays for MK now has all the vals
singhbalwinder Jul 7, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
131 changes: 100 additions & 31 deletions analysis-scripts/test_qv_sat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,14 @@
https://journals.ametsoc.org/doi/abs/10.1175/1520-0450%281992%29031%3C1507%3APFTSVP%3E2.0.CO%3B2)
in python to make sure it matches our result in C++.

[BSINGH: Extended the original code to add new test cases along with Murphy and Koop specific
print atements. Test cases are obtained from Murphy and Koop (2005) paper (accessed 07/05/2020):
https://pdfs.semanticscholar.org/519c/bb52a54c5abadff7263bf98678fccce5b594.pdf

Added control variables "print_flatau" and "print_murphyKoop" to decide whether to
print Flatau specific or MurphyKoop specific (or both) print statements]


#Notes:
1. Units of output from Flatau are never given. They end up being hPa, which can be computed by
looking at output when T=T0 => es=a1 = 6.12... which is 100x smaller than the 612 Pa we know is the
Expand Down Expand Up @@ -36,6 +44,12 @@
import numpy as np
import pylab as pl

#Setting print_* variables to "True" would print values of that scheme
#currently implemented schemes are "Flatau" and "Murphy and Koop"
print_flatau = False
print_murphyKoop = True


#DATA USED IN CALCULATIONS:
#=============================
T0=273.15 #"conversion between Celsius an Kelvins" from Flateau. Note this was what was wrong in initial P3 impl.
Expand Down Expand Up @@ -360,8 +374,8 @@ def MurphyKoop_esi(t):
#allowing us to read the expected value directly off the table. Multiplying by 100 b/c
#table values are in hPa.

print 'When T=T0, polysvp should exactly match 1st Flatau coef. Diff = ',polysvp_liq(T0) - a_liq[0]*100.
print '\n'
print('When T=T0, polysvp should exactly match 1st Flatau coef. Diff = '+str(polysvp_liq(T0) - a_liq[0]*100.))
print('\n')

#unfortunately, polysvp uses liquid saturation at exactly T=T0, so this test is expected to fail.
#print('polysvp_ice(T0) - a_ice[0] = ',polysvp_ice(T0) - a_ice[0]*100.)
Expand All @@ -388,15 +402,41 @@ def MurphyKoop_esi(t):
T243P500_data=[243.15, 5e4, 37.98530141245404, 50.98455924912173,\
0.00023634717905493638, 0.0003172707211143376]

#Following Test values are directly coming from Murphy and Koop (2005) paper
#Table C1, titled "VALUES RECOMMENDED FOR CHECKING COMPUTER CODES"

#DEFINITIONS:ice vapor pressure(ivp); liquid vapor pressure(lvp)
#"expected values" are from the paper

#Test values @150 @ 1e5 Pa (expected values-> ivp: 6.106e-6; lvp:1.562e-5)
T150_data = [150, 1e5]

#Test values @180 @ 1e5 Pa (expected values-> ivp: 0.0053975; lvp: 0.011239)
T180_data = [180, 1e5]

#Test values @210 @ 1e5 Pa (expected values-> ivp: 0.70202; lvp: 1.2335)
T210_data = [210, 1e5]

#Test values @240 @ 1e5 Pa (expected values-> ivp: 27.272; lvp: 37.667)
T240_data = [240, 1e5]

#Test values @273.16 @ 1e5 Pa (expected values-> ivp: 611.657; lvp: 611.657)
T273_16_data = [273.16, 1e5]

#Test values @300 @ 1e5 Pa (expected values-> lvp:3536.8, no ivp mentioned as temp is 300K)
T300_data = [300, 1e5]

Copy link
Contributor

@PeterCaldwell PeterCaldwell Jul 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One last question - why include the MK expected values in a comment rather than in the data array? Not only does it break the paradigm started by the Flatau data, but it would be nice to print out the difference between the values we're getting from MK calculations and what's expected from their table...

Copy link
Contributor Author

@singhbalwinder singhbalwinder Jul 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will add these to the arrays. I noticed that no other field from the "data" arrays (except for the first two) were used in the script. Therefore, I decided not to put these values in the array (as they are not used) but specify them in comments so that we know what should we expect without referring to the paper.

Also, I noticed that, the dt calculation in polysvp looks like the following:
dt = T-T0
I modified it to:
dt = max(T-T0, -80.0)
so that it is exactly the same as in the code. Should I push this change as well?

Copy link
Contributor

@PeterCaldwell PeterCaldwell Jul 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

About including expected MK values in a comment rather than including them in the data array like Flatau: I made my comment after seeing those later entries used in the print statements... but closer inspection now shows that those later-entry print statements are actually commented out. So I don't care if you fix this or not. Just including the 2 entries that are actually used is actually kind of nice.

Copy link
Contributor

@PeterCaldwell PeterCaldwell Jul 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

About dt = max(T-T0, -80.0): I don't think it matters since we never call the python with T<-80C(?). I guess we might as well make the change for consistency. I didn't realize we were thresholding T to be warmer than -80C in the model. I guess that should prevent the negative T errors we've been seeing in polysvp, huh?

Copy link
Contributor Author

@singhbalwinder singhbalwinder Jul 7, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, for some reason, I didn't read your first comment earlier and went ahead to fix the code. The code now has complete data arrays with None inserted for the missing data. I have also added max for -80C.

MK tables starts from 150K (-123.15 C), so I had to add max(T-T0,-80) to get values for Flatau (similar to the values produced by the C++ code).


#The last line in this loop prints stuff to copy/paste into C++ tests.
#Tell users what that data consists of:
print('Below, "data to copy" is: [sat ice pres, sat liq pres, qsat_ice, qsat_liq]' )


ind=0
for data in [tmelt_data,T243_data,T303_data, T243P500_data]:
for data in [tmelt_data,T243_data,T303_data, T243P500_data, T150_data, T180_data,
T210_data, T240_data, T273_16_data, T300_data]:
ind+=1
print "Case %i: T=%5.2f p=%5.1f mb"%(ind,data[0],data[1]/100.)
print("Case %i: T=%5.2f p=%5.1f mb"%(ind,data[0],data[1]/100.))

#SATURATION VAPOR PRESSURE CHECKS:
#---------------------------------
Expand All @@ -407,15 +447,6 @@ def MurphyKoop_esi(t):
e_ice_GG=GoffGratch_esi(data[0])
e_liq_MK=MurphyKoop_esl(data[0])
e_ice_MK=MurphyKoop_esi(data[0])
#print(' Naive minus old expected e_ice='+str(e_ice_f - data[2]))
#print(' New minus old expected e_ice='+str(e_ice - data[2]))
print(' F90 ice p vs naive ice p='+str(e_ice - e_ice_f))
print(' Flatau vs GoffGratch ice p='+str(e_ice_GG - e_ice))
print(' Flatau vs MurphyKoop ice p='+str(e_liq_MK - e_ice))

print(' Flatau vs GoffGratch liq p='+str(e_liq_GG - e_liq))
print(' Flatau vs MurphyKoop liq p='+str(e_liq_MK - e_liq))
#print(' actual minus old expected e_liq='+str(e_liq - data[3]))

#SATURATION MIXING RATIO CHECKS:
#---------------------------------
Expand All @@ -424,33 +455,71 @@ def MurphyKoop_esi(t):
qsat_l_MK=vap2mix(e_liq_MK,data[1])
qsat_exact_l=qs_exact(data[0],data[1],is_ice=False)
qsat_num_l = qs_numeric(data[0],data[1],is_ice=False)
#print(' actual minus old expected qsat_liq='+str(qsat_l - data[5]))
print(' exact minus numerical qsat_liq='+str(qsat_exact_l - qsat_num_l))
print(' exact minus Flatau qsat_liq='+str(qsat_exact_l - qsat_l))
print(' Percent error in Flatau liq = %e'%((qsat_exact_l - qsat_l)/qsat_exact_l*100.) )
print(' Percent error in GoffGratch liq = %e'%((qsat_exact_l - qsat_l_GG)/qsat_exact_l*100.) )
print(' Percent error in MurphyKoop liq = %e'%((qsat_exact_l - qsat_l_MK)/qsat_exact_l*100.) )

#if above freezing, don't bother with expensive repeat calculations
if data[0]<T0:
if data[0]<T0:
qsat_i=vap2mix(e_ice,data[1])
qsat_i_GG=vap2mix(e_ice_GG,data[1])
qsat_i_MK=vap2mix(e_ice_MK,data[1])
qsat_exact_i=qs_exact(data[0],data[1],is_ice=True)
qsat_num_i = qs_numeric(data[0],data[1],is_ice=True)
#print(' actual minus old expected qsat_ice='+str(qsat_i - data[4]))
print(' exact minus numerical qsat_ice='+str(qsat_exact_i - qsat_num_i))
print(' exact minus Flatau qsat_ice='+str(qsat_exact_i - qsat_i))
print(' Percent error in Flatau ice = %e'%((qsat_exact_i - qsat_i)/qsat_exact_i*100.) )
print(' Percent error in GoffGratch ice = %e'%((qsat_exact_i - qsat_i_GG)/qsat_exact_i*100.) )
print(' Percent error in MurphyKoop ice = %e'%((qsat_exact_i - qsat_i_MK)/qsat_exact_i*100.) )
else:
qsat_i=qsat_l

qsat_i = qsat_l
qsat_i_MK = qsat_l_MK

if(print_flatau):
print('----------------------------------------------------------')
print('-- Flatau --')
print('----------------------------------------------------------')
#print(' Naive minus old expected e_ice='+str(e_ice_f - data[2]))
#print(' New minus old expected e_ice='+str(e_ice - data[2]))
print(' F90 ice p vs naive ice p='+str(e_ice - e_ice_f))
print(' Flatau vs GoffGratch ice p='+str(e_ice_GG - e_ice))
print(' Flatau vs MurphyKoop ice p='+str(e_liq_MK - e_ice))

print(' Flatau vs GoffGratch liq p='+str(e_liq_GG - e_liq))
print(' Flatau vs MurphyKoop liq p='+str(e_liq_MK - e_liq))
#print(' actual minus old expected qsat_liq='+str(qsat_l - data[5]))
print(' exact minus numerical qsat_liq='+str(qsat_exact_l - qsat_num_l))
print(' exact minus Flatau qsat_liq='+str(qsat_exact_l - qsat_l))
print(' Percent error in Flatau liq = %e'%((qsat_exact_l - qsat_l)/qsat_exact_l*100.) )
print(' Percent error in GoffGratch liq = %e'%((qsat_exact_l - qsat_l_GG)/qsat_exact_l*100.) )
print(' Percent error in MurphyKoop liq = %e'%((qsat_exact_l - qsat_l_MK)/qsat_exact_l*100.) )
#print(' actual minus old expected e_liq='+str(e_liq - data[3]))
if data[0]<T0:
#print(' actual minus old expected qsat_ice='+str(qsat_i - data[4]))
print(' exact minus numerical qsat_ice='+str(qsat_exact_i - qsat_num_i))
print(' exact minus Flatau qsat_ice='+str(qsat_exact_i - qsat_i))
print(' Percent error in Flatau ice = %e'%((qsat_exact_i - qsat_i)/qsat_exact_i*100.) )
print(' Percent error in GoffGratch ice = %e'%((qsat_exact_i - qsat_i_GG)/qsat_exact_i*100.) )
print(' Percent error in MurphyKoop ice = %e'%((qsat_exact_i - qsat_i_MK)/qsat_exact_i*100.) )
#PRINT DATA TO COPY
#---------------------------------
print ' Data to copy: ',e_ice,e_liq,qsat_i,qsat_l

print(' Data to copy: ',e_ice,e_liq,qsat_i,qsat_l)
if(print_murphyKoop):
print('----------------------------------------------------------')
print('-- Murphy and Koop --')
print('----------------------------------------------------------')
#print(' F90 ice p vs naive ice p='+str(e_ice - e_ice_f))
print(' MurphyKoop vs GoffGratch ice p='+str(e_ice_GG - e_ice_MK))
print(' MurphyKoop vs Flatau ice p='+str(e_ice - e_ice_MK))

print(' MurphyKoop vs GoffGratch liq p='+str(e_liq_GG - e_liq_MK))
print(' MurphyKoop vs Flatau liq p='+str(e_liq - e_liq_MK))
print(' exact minus numerical qsat_liq='+str(qsat_exact_l - qsat_num_l))
print(' exact minus MurphyKoop qsat_liq='+str(qsat_exact_l - qsat_l_MK))
print(' Percent error in MurphyKoop liq = %e'%((qsat_exact_l - qsat_l_MK)/qsat_exact_l*100.) )
print(' Percent error in GoffGratch liq = %e'%((qsat_exact_l - qsat_l_GG)/qsat_exact_l*100.) )
print(' Percent error in Flatau liq = %e'%((qsat_exact_l - qsat_l)/qsat_exact_l*100.) )
if data[0]<T0:
#print(' actual minus old expected qsat_ice='+str(qsat_i - data[4]))
print(' exact minus numerical qsat_ice='+str(qsat_exact_i - qsat_num_i))
print(' exact minus MurphyKoop qsat_ice='+str(qsat_exact_i - qsat_i_MK))
print(' Percent error in MurphyKoop ice = %e'%((qsat_exact_i - qsat_i_MK)/qsat_exact_i*100.) )
print(' Percent error in GoffGratch ice = %e'%((qsat_exact_i - qsat_i_GG)/qsat_exact_i*100.) )
print(' Percent error in Flatau ice = %e'%((qsat_exact_i - qsat_i)/qsat_exact_i*100.) )
#PRINT DATA TO COPY
#---------------------------------
print(' Data to copy: ',e_ice_MK,e_liq_MK,qsat_i_MK,qsat_l_MK)

#CHECK THAT Pressure CLAUSIUS-CLAPEYRON FORMULATION IS RIGHT
#=============================
Expand Down