bug: Definition and use of "unnormalized" scale height in cam-fv is inconsistent

[kdraeder]

Describe the bug

There are several, related bugs, which I've tried to separate into managable units.
This is the first I'd like to resolve, even though I filed one (#296) before this one.

The vertical coordinate "scale height" is defined as log(ps/p), ps = surface pressure (or some other reference)
and log = natural log ln. This can be written as log(ps) - log(p).
When the difference of 2 scale heights is calculated, the 2 log(ps) terms cancel,
so we could save some computation by defining an "unnormalized" scale height = log(1/p) = - log(p)
and calculating SH1 - SH2 = log(p2) - log(p1) (note the numbers).
But when one of those unnormalized SH is compared against some threshhold height,
like ramping_end, that threshhold must be calculated the same way - unnormalized - which it was not.

This topic was further confused by the misdefinition of unnormalized scale height as "log(p)" instead of "- log(p)"
and the resulting, unnecessary complexity in defining higher_is_smaller, and unintended consequences
in routines that use higher_is_smaller.

1. List the steps someone needs to take to reproduce the bug.
   Set up a cam-fv assimilation using

• vertical_localization_coord = SCALEHEIGHT
• model_damping_ends_at_level = 5
• an obs_seq.out file which has obs at heights above model level 6

2. What was the expected outcome?
   There should be no increments at and above layer 5.

3. What actually happened?
   There are non-0 increments in those layers.
   The unnormalized scale heights of the state variables have values ranging from
   5.88 (log(360 Pa)) to 11.5 (log(99200 Pa)), while the ramp_start (normalized scale height)
   had a value of ~ 3.0 (log(10^5/3500) - .3 ramp_depth).
   It would appear that all of those state variables are in the damped layer,
   but the choice of using unnormalized SH caused higher_is_smaller to be .true. (which it should not have),
   and that caused v_above to test whether the state variables have values less than ramp_start.
   They do not, so above_ramp_start decided that none of them were in the damped layer.
   and so their increments were not damped.

Which model(s) are you working with?

Cam-fv only (all versions of CAM and WACCM(-X))

Version of DART

Which version of DART are you using? Manhattan, all of the versions

Have you modified the DART code?

Only cam-fv/model_mod.f90
I'll issue a pull request shortly on a branch of Manhattan.

Build information

Please describe:

1. The machine you are running on (e.g. windows laptop, NCAR supercomputer Cheyenne).
   Any machine.
2. The compiler you are using (e.g. gnu, intel).
   Any compiler

[nancycollins]

hi kevin - apparently the version of the wrf model_mod which supported unnormalized scale heights is not on the main branch. that's a problem to solve for another day. can you see this file?

https://github.com/nancycollins/DART_development/blob/all_nsc_changes/models/wrf/model_mod.f90

if so, that has code developed with craig schwartz that does not normalize by surface pressure, and it returns the positive log. in most cases the scale height value is used to compute distances, so the separation between 2 locations. if they're both positive or both negative you get the same distance. but i remember we decided that scale heights should in fact be the positive log of pressure.

wikipedia says scale height is a distance (like in kilometers) across which the atmospheric pressure drops by a factor of e. so asking for the scale height of a single location seems to be the value of that location relative to the surface of the earth? so that's the natural log of the pressure, as a positive value, i'd think.

[kdraeder]

Even if the SH is defined as + log(p), it still has positive and negative values.
WACCM extends up to pressures that are < 1 Pa, so for example log(0.1 Pa) = - 2.3.
That's why the distance between 2 heights is (hopefully) calculated as abs(SH1 - SH2)
and the signs will be taken into account properly to give the distance.

What's really happening with the "unnormalized" scale height is that it's _re_normalized
with a reference pressure of 1 Pa. The definition can be written +log(1/p) (or -log(p/1).
1 is the new dividing line between positive and negative log values.
log(1/0.9) = +0.1, log(1/1.1) = -.09

When the reference pressure is 10^5 Pa (~Earth's surface), that's the dividing line
between positive and negative logs: log(10^5/99000) = + 0.01, log(10^5/101000) = - 0.01.

So the conventional definition of scale height, +log(ps/p) (== -log(p/ps)) gives scale heights
that are 0 at the surface and positive all the way to the moon.

Arthur(?) brought up the issue of terminology. The definition of scale height "H" is actually
the one you found on wikipedia. Arthur pointed out that we're calculating a "scaled height".
To see what he means we have to start with the hypsometric equation. The difference in geometric height is
z2 - z1 = - H (log(p2) - log(p1)) = + H log(p1/p2)
If we rescale the geometric heights by the "scale height" we get the quantity we've been calling "scale height":
(z2 - z1)/H = log(p1/p2). At the (reference) surface z1 = 0 and p1 = ps;
z2/H = log(ps/p2) is the "scaled height".

[nancycollins]

thanks for that. what i should have said is that i don't see any need to take -log(p) in our code (compared to just log(p)) since as you point out, the actual value returned can be positive or negative. the vertical separation eventually needs to be converted to radians in the location module so it can be combined with the horizontal separation. here's the code that computes it in the threed_sphere locations mod:

vert_dist = abs(loc1%vloc - loc2%vloc) / vert_normalization(loc1%which_vert)

so it is the abs of the distance divided by the normalization factor set by the user in the namelist. (note these can be set on a per-type basis in the threed_sphere locations mod.)

the original implementation did normalize the values by the surface pressure. but since we are computing the vertical separation between two points which may be quite far apart in the horizontal, changes in surface pressure near steep topography at the two points changes the vertical distance returned. in revisiting the code later the consensus was that this wasn't what we wanted.

arthur is right, we're computing a 'scaled distance' but at this point the terminology is too ingrained, i think, to change. we can certainly document this better in the threed_sphere locations module -- what we're expecting a model_mod to do if the vertical localization type is set to scale height.

so the bottom line, in dart, i think is: scale height computations should not be scaled by surface pressure, and there's no reason to take the negative of the log before returning it. it certainly does have to be treated the same everywhere and if that's not the case, that is a bug to be fixed.

i confess that i added code in mpas_atm, wrf (uncommitted) and in cam-fv that allowed the user to set a namelist value to do the old computation that scaled by surface pressure for regression testing. it may be at this point that the additional complexity is not worth it and that code could be removed and the logic simplified.

[kdraeder]

If we exclusively use +log(p), then what we're calling "scale height"
will decrease (change in the negative direction) with height,
which is the opposite of what everyone expects.
When we print vertical locations (cam-fv does this, and others may in the future),
we start with a pressure, and then say "which is equivalent to scale height "
and for reference we give the scale height of the top of the model.
The top of the model will have a smaller "scale height" than the other vertical locations,
which makes it look lower to non-DART people.

And it contributes to the confusion between "larger or smaller" and "increasing or decreasing"
and figuring out whether a point is above or below a threshhold.

I agree that normalizing by the surface pressure is dangerous.
That's why the normalization should be done relative to a fixed reference pressure; 10^5 Pa
(or 1 Pa, if it's "unnormalized").

It can be made to work, but it will take some explaining to dispel the appearance
that we don't know what we're doing.

[nancycollins]

"it can be made to work"? it is already working. maybe a bug fix is needed for code that is now being called that wasn't tested before, but it's not like it's a big development job. i'm not sure i agree with your argument that it's going to confuse people because mostly it is completely invisible. it's used in the distance calculations only and not written out in the obs_seq.final file. my thoughts on changing it now are: 1) we can call it "scaled height" in our documentation; 2) it's this way in the mpas_atm and this will be inconsistent; 3) no one from MMM has complained about this definition; 4) pressure also decreases with increasing height so it's not like we don't already handle this; 5) WACCM users have been localizing in scale height for several years, i believe, and if this really matters it will be a change for them which is in itself unexpected; and 6) the code has been there for several years without any complaints from users. maybe ask nick whether he has had problems with the current way the code works?

[kdraeder]

I've taken another look at the original code and the first fixes I made for Gio's problem
and I see now that the original code handled the log(p) in a way that works, but at the time it didn't occur to me
that scale height would decrease with height. So my fixes for Gio were incorrect, and they led me to look at
what scale height actually means in cam-fv. What I found was not consistent with my understanding,
seemed more complicated than it needed to be, and seemed to be contributing to the problem
I was trying to solve. Vertical locations that were being printed were not the ones being used,
which was ... unhelpful.
So, yes, +log(p) worked, but this seems like an opportunity to bring this flagship clean code
into allignment with commonly accepted definitions. If we're going to define a new vertical coordinate system
that's different from scale height and pressure and geometric height, we should call it something
besides "scale height". I nominate "logp". A coordinate that increases in the opposite direction and has a different reference point than the conventional scaled height is a different coordinate.
I'm probably taking this too seriously, as a result of spending too much time figuring out what was going on
in Gio's case.

2)-3) Those models don't have vertical location (not distance) tests yet, but Chris reported at a standup recently
that there are some problems in mpas-atm at the top of the model, so it may need code with vertical threshholds. I'm in favor of consistency too, within DART and with the wider community, and I favor it over
continuing to do something suboptimal, just because we've done it that way for a while.

4. People expect pressure to decrease with height. They expect scale height to increase with height.
   As Nancy pointed out, if they never look, they won't notice it. If they look, they'll be baffled for a while.

5)-6) As I wrote above, the code did work in the ways people asked it to. They didn't need to know
the details of how scale height was defined. That would be true for a new model_mod too, so there wouldn't
be any surprise.

I've spent too much time on this, so I'm going to leave it to Nancy and Helen to decide which parts of this PR
to approve and I won't advocate any more. The candidates are, in order of increasing controversy:

1. Changing the definition of no_assim_above_pressure, ramp_end, and model_top to use no_normalization_of_scale_heights instead of .false. The way it is, unnormalized variables are compared against normalized vertical locations, if no_normalization_of_scale_heights = .true.
2. Applying nint() to query_locations which are just getting the which_vert.
3. Removing the vert_only_distance section from above_ramp_start. Vert_only_dist is the vertical distance between a state variable (or ob), and the bottom of the increment damping ramp, while total_dist is between the state variable and the base_ob, which has nothing to do with increment damping. That comparison is done just to decide whether to print a warning and calculate horiz_dist, which subracts the vertical distance to the ramp from the total distance to the base_ob, which I think is not meaningful.
4. Changing higher_is_smaller to higher_is_decreasing
5. Changing unnormalized "scale height" from +log(p) to -log(p), and related changes.

[nancycollins]

would you be able to close your existing pull request and open another one with just the proposed changes to the model_mod.f90 in it, relative to the main branch? the current pull request has 1000+ changed files from trying to merge main onto manhattan, and i don't think it includes your proposed version of the model_mod.f90 because the global variable for "scale height too high" isn't there.

[jlaucar]

We'll pick up this discussion briefly in the stand-up. One of the problems is that Kevin is trying to fix a terminology problem that is out of our hands. The use of "scale-height" for correlation distances in the vertical in the DA community evolved elsewhere and was picked up by our WRF collaborators when they started implementing things. If I recall correctly, that community used the difference in log(p/ps) to define the scale-height "distance". If the ps is the same, it would disappear in the difference. However, when a difference was computed with significantly different ps, this ended up playing a significant role in this distance definition which caused significant problems in areas of varying orography as Nancy pointed out. I believe this led us to remove the ps from our definition so that the DART scale-height distance became log(p1) - log(p2). This would be better characterized as a log(p) vertical coordinate or log(p) distance. However, the term "scale-height" had already been widely used in the peer-reviewed literature for the original ps normalized differences. I believe that the only difference comes back to Kevin's concern about the sign of the scaled height as a function of geometric height. Kevin's argument about the sign now becomes less obvious. The way the term is being used in DA is not being thought of as a "scaled-height", but as a log(p) but called a scale height. I would have to confirm that when used with WRF our sign behavior is such that the "scale-height" gets less as one moves up from the earth's surface. Sounds like Nancy already knows that answer. So, how to proceed. In the short-term, the behavior in CAM needs to stay the same unless there was a fundamental bug. My understanding is that the bug had nothing to do with the computation of distances in "scale-height" but rather in the way in which things were computed in columns by a software expediency. Second, we do not have time to be fixing the more general terminology problem right now because of the widespread adoption in the DA community. In particular, whatever we do in CAM must remain consistent (or become consistent) with the definitions in WRF and MPAS that have been used a lot. Finally, a conversation with people like Ryan and maybe some NCEP DA folks would need to take place to assess whether changing terminology to log(p) or something similar would be useful in the long term. I have not had time to verify all the things I say here, so I expect the two of you will correct me if I am mistaken about any of these issues. Jeff

…

On Wed, Oct 13, 2021 at 4:23 PM nancy collins ***@***.***> wrote: would you be able to close your existing pull request and open another one with just the proposed changes to the model_mod.f90 in it, relative to the main branch? the current pull request has 1000+ changed files from trying to merge main onto manhattan, and i don't think it includes your proposed version of the model_mod.f90 because the global variable for "scale height too high" isn't there. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#298 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANDHUIUBPWMZUD7YQBFM2RDUGYBGLANCNFSM5FZGMP7A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[kdraeder]

would you be able to close your existing pull request and open another one with just the proposed changes to the model_mod.f90 in it, relative to the main branch? the current pull request has 1000+ changed files from trying to merge main onto manhattan, and i don't think it includes your proposed version of the model_mod.f90 because the global variable for "scale height too high" isn't there.

I followed Helens instructions for changing the base on the issue, and then pushed the commit again.
My PR page shows that I'm trying to merge from main. Am I missing something?
If I need to close it, I'd prefer to do that after I have a decision about what should be in the new version.

[nancycollins]

great - i now see only a single changed file. thank you!! i can review that.

but i still don't see a global variable around line 260 in the proposed model_mod.f90 for no_assim_above_scale_height, so this can't be the final version you want us to review? also obs_too_high() is still being called in model_interpolate()...

[kdraeder]

Lots of potential changes came up while debugging Gio's problem.
I've tried to divide them into magagable, but cohesive Issues,
so this Issue does't have the things that are more directly related to obs_too_high.
I planned submit that issue next. Jeff, in general, supports this strategy,
but if it's causing problems we can develop a new strategy.

Based on Jeff's response to this discussion, 5) in my list of candidate changes is rejected.
It also seems that 4) should be rejected in the name of minimizing changes.
Any thoughts about 1)-3)?

[nancycollins]

kevin - maybe you talked about this in the standup - if so, my apologies.

i'm getting 'issues' and 'pull requests' confused. if you want to deal with other changes in another issue, that's great. but i'm assuming that in the end there will be a single pull request and release that deals with all the cam top of model issues together. committing updates to your branch incrementally is still a great idea because it helps you isolate changes for testing.

is your plan that you will be keeping none of the changes that are currently in model_mod.f90 in the open pull request? some are sign changes which we aren't doing, and some are debugging code that shouldn't stay in the release once you've tested things. (if that's not the case i will review them and give you specific feedback.)

as far as thoughts about issues 1,2 and 3 - it would help to see the code changes for issues 1 and 3. for 1) i don't know what 'no_assim_above_pressure' has to do with ramping and scale height. if i see the changes to fix a bug then it probably makes perfect sense. for 3) you say total distance doesn't matter, only the vertical separation for the top ramp. so 'vert_only_distance' seems useful? but i probably don't remember how the code works so seeing your proposed changes would sort that out for me.

for 2) adding nint() probably falls into the 'nuisance change' for me since that's what the fortran standard does with or without that call when assigning a real to an integer.

[kdraeder]

My (mis)understanding is that we open an issue about a manageable topic
and then have possibly multiple PRs related to it, as solutions are developed and tested.
I can't say much about the end game; merging changes into main. I'm happy for it to happen
however makes sense to the person who's doing it; several smaller merges from a series of issues,
or one merge from a combination that we assemble at the end.

To me 1)-3) are still potentially going in.
1) I believe this is required to be in the final package.
2) There are places in the existing code which use nint(query_location), which is why I felt comfortable
adding more of those, even though I suspected that they are not necessary. Do we want to remove the
existing ones, or go with inertia; it exists and is not causing a problem, so we'll leave it in?
3) When I do
git difftool --tool=diffuse --no-prompt origin/main -- model_mod.f90
on the model_mod which I pushed, I see the vert_only_dist section highlighted.
Are you not able to see it?
Vert_only_dist is used only in a comparison with total_dist, which makes no sense to me,
and in the calculation of horiz_dist, which is used only in a (commented out) diagnostic statement.
The calculation, which is actually used by the subroutine, of the vertical distance
between the variable and the ramp bottom happens in
ramp_dist = v_distance(test_value,ramp_start).
Now that I look even more closely, the 3 lines before that section are also unneeded;
this_loc =
ramp_start_loc =
It looks to me like this whole section falls in the "unneeded diagnostic statements" category.

[kdraeder]

no_assim_above_pressure is in this issue and PR because its calculation uses a normalization argument,
which has been set to false, but should use no_normalization_of_scale_heights.

[nancycollins]

thanks. where i got confused was how no_assim_above_pressure related to missing that normalization flag which only applies to scale height. it's needed for initializing no_assim_above_scaleh, not pressure. it does look like it was missing. good catch.

[nancycollins]

kevin - i'm going to post my code review in just a bit. i hope it doesn't come across as harsh - it's hard for me to be clear but soften the tone on some of those items. based on the discussion about the misuse of scale height in DA, and the fact that other modules (wrf, mpas) use positive log() in their code, i think the negative logs need to be removed. i also prefer to make no changes that aren't required to fix the bugs. so i'm objecting to the variable rename and adding the nints(). i am happy to be disagreed with by you or anyone else who is reviewing the changes. i'm retired and the ongoing responsibility for the code is yours, so in the end current group members have to make the final decisions. i have stronger opinions on this code than others because i was part of the re-factoring with johnny when that was done several years ago, so i sometimes can reconstruct the intent, even if there are bugs in the implementation.

[kdraeder]

Nancy, thanks for the code review. I'll incorporate the appropriate parts in the next PR,
since Helen closed the first PR.
Helen wrote

Lets step back from reviewing the code, and start with:

1. Documentation of what the code is supposed to do.
2. A link to restart files + obs_seq that lets someone reproduce the issue (on Cheyenne or wherever).
3. If this is going to change people's results, who does it affect. e.g does this affect everyone who localizes in SCALEHEIGHT? All WACCM users? Is this a non-bitwise change for CAM users?

---

This issue was intended to cover just some of the changes I believed were needed to get Gio running.
Since it's been decided that the previous implementation of "scale height" will be kept, there's less reason
(but not none) to continue with that strategy. Or we can close this issue and open a new one that
addresses all of the remaining issues, after we set up the process outlined by Helen.
What's the preference?

1. I'll work on the description this afternoon, focusing on just the parts that need fixing for rttov radiance obs.
2. I'll ask Gio for data files and any modified code and scripts for using them. I doubt that I'll have anything to work with until Monday, and then it will probably take some time to get it all working here.
3. Since we're not changing the definition of scale height, and the non-cam-fv models don't even use no_normalization_of_scale_heights, there should be no impact on them. In cam-fv models it appears that assimilations, that did not use rttov, worked correctly. That should be all of them, since using rttov is what brought all this up.

It's a mystery to me, so far, how they excluded obs in the upper layers. My understanding is that get_close_YYY is called before model_interpolate. Get_close_ converts vertical locations to (in this case) VERTISSCALEHEIGHT. Are those locations discarded after get_close_? If not, those locations are passed to model_interpolate, which passes them to obs_too_high, which does NOT check for scale height. But it also does not error out, so model_interpolate is giving it a which_vert which is not VERTISSCALEHEIGHT.
Maybe we don't need to answer this, since we're moving obs_too_high out of model_interpolate and into get_close_. But I'd feel a little queasy about not understanding it.

[nancycollins]

starting a new issue is probably better but we can do it either way.

here is the calling sequence in filter:

model_interpolate() is called during forward operator computation, for all obs, which happens before any assimilation starts. those obs have the original locations in the obs_seq.out file.

at the start of the assimilation phase there are options to either convert all the obs and state into the requested vertical location, or they can be converted one by one in the get_close routines. there are efficiency and communication bottleneck tradeoffs so it's a namelist option which happens. for a global model with global obs, doing all the conversions up front is probably the better strategy. when they're converted, the converted verticals are saved for the rest of the assimilation phase.

then if posterior forward operators are called, they use model_interpolate on the updated state but they use the original obs locations.

does this help?

[nancycollins]

i have a suggestion to ease kevin's mind about how DA treats scale height.

we have a model_nml namelist item to control how scale height is computed called "no_normalization_of_scale_heights". i propose we remove that variable from the namelist but leave it as an internal module global, and add a namelist item called "use_logp_for_scale_height = .true." in the init code we assign that namelist item to no_normalization_of_scale_heights and nothing else changes. then the name of the item is a clue that this is special, and the documentation can explain the differences between the options.

[kdraeder]

I think that Nancy's suggestion of use_logp_for_scale_height would be helpful.
Thanks also for the clarification about when model_interpolate is called.

Going back to the list of things to do from Helen:

2. Gio has provided code and files to help with a test case, and I've pointed to local files, in #296

1. Here's my description of the parts of model_mod at issue, and what they should do:

1) Localize vertically in "scale height", which currently can mean 1 of 2 things:
   no_normalization_of_scale_heights =
      .false.: scale height = log(p_reference/pressure).   
               Note that some people have set p_reference = p_surface, 
               which varies with location and weather, and is not a robust choice.  
               CAM has a reference pressure = 10^5 Pa, which is used in the definition 
               of its vertical coordinates.  This is a good choice for p_reference. 
      .true.:  scale height = +log(pressure).

2) Exclude observations above a user specified level from being assimilated; 
   no_obs_assim_above_level.

3) Prevent observation increments above a user specified model level; 
   model_damping_ends_at_level.  
   The current code uses the idea that there should be a gradual reduction in the increments, 
   instead of a sharp transition to 0, to be gentler to the model state.

4) Code outside of model_mod may want to get information about the whole column, 
   so the implementation of no_obs_assim_above_level needs to allow that.

Implications
1) The issue of the reference pressure in scale height would be clarified 
   by changing "p_surface" to "p_reference".
1) When debugging be careful to keep track of whether a pressure (or log(p)) 
   changes in the positive direction with height and don't confuse positive 
   with a larger size or with higher.
2),3) no_obs_assim_above_level and model_damping_ends_at_level need to be translated 
   into the vertical localization coordinate.  The correctly defined version needs to be used 
   in obs_too_high (see earlier notes about no_assim_above_level).
4) Model_mod code that uses no_obs_assim_above_level should appear in get_close_{obs,state}, 
   not in model_interpolate.  That requires that it handle VERTISSCALEHEIGHT too, 
   since state and obs can have that in the get_close routines.  
   The efficient place for no_obs_assim_above_level to be checked is near the beginning, 
   so that no further calculations are done for obs that are too high.

[nancycollins]

kevin, this looks good with just a few comment about the implication sections.

1. if the code uses p_reference which is constant instead of p_surface (which varies by location), then the resulting distance computation is exactly the same as the unnormalized choice. the differences are divided by the same values so they cancel out. changing the variable name doesn't seem worth the code changes since it doesn't make sense to substitute a constant reference pressure for p_surface, and there is no code to do that in the model_mod now.

2, 3) correct - and i believe that code is already there except for changing the local in the init code to a global for the scale height because the value is already being computed. then obs_too_high() needs to test the global value if scale height is the selected vert localization coordinate.

4. when you say 'near the beginning' you mean in the get_close routines? that's right. after the base obs is converted to the vertical localization type, then it can be tested to decide whether to compute any other distances or not. if the base obs is too high, it can return immediately with 0 close items. if it's not, then it does the standard distance computations without any other height tests. if the user wants to ramp the impact at the top, that's a secondary computation controlled by the ramp settings. the test should be by calling obs_too_high() which handles all the vert types (once you add the scale height test there). it shouldn't be looking at level in the get_close() routines because the vert types have been converted to the vert localization choice already.

does this make sense?

[kdraeder]

Some of this is veering into solutions, instead of Helen's request for "what should the code do",
but that's useful too.

1) I was thinking narrowly of the scale_height function. Sorry that I didn't say that. In my defense, "p_surface" only appears in that function. The function can have either the global ref_surface_pressure passed to it, or the surface pressure, if that was actually useful for some reason. But inside the function they both serve as the reference pressure in the definition of scale height = log(p_ref/p), so I'm suggesting changes only in that function. I should know by now that if it's functioning correctly, it's not likely to be changed.

Yes, the cancellation of the log(p_ref) parts of a vertical distance calculation is exactly what prompted the definition and use of log(p) instead of log(p_ref/p) (at least in cam-fv); it gives the same distance for less computation. But the cancellation only works if the same p_ref is used for the calculation of the scale height location of both points. If I understood the discussion about other models decision to use only SH = log(p), they decided that because they ran into problems when they used log(p_surf/p) to calculate distances between points at different horizontal locations that had different p_surf. If they had used log(10^5/p) they wouldn't have had those problems. I'm hoping that everywhere in the code that passes p_surf(i,j) to scale_height for 2 locations always does it for the same i,j, so that a differing reference pressure is not causing unnoticed problems.

And as I pointed out before, log(p) implies a reference pressure of 1 Pa. log(p) below that has positive values, log(p) above it has negative. So there's a reference pressure (point where the coordinate = 0) even in the log(p) coordinate. It's just implicit (hidden). Again, I know it doesn't make any difference in distance calculations.

2), 3) And no_normalization_of_scale_height needs to be passed to the subroutines instead of .false.

4) I believe that the check for obs_too_high can be done even before the conversion to the vertical_localization_coord because obs_too_high can handle the original vertical coordinate of the ob.

[nancycollins]

1. sorry helen! my point was, i think, that if you use a constant reference pressure that you will get identical distance results to simple log(p) because everything else in the distance computation cancels out. so the only reason to use the other namelist option is if you want to explore using the real surface pressures and include surface pressure differences in the distance computations.

2,3,4) you're right about passing the normalization flag into the subroutines, and that obs_too_high can deal with any vertical type. good points.

[kdraeder]

Woops! I wrote this 10 days ago and got distracted from sending it.

---

Some of this is veering into solutions, instead of Helen's request for "what should the code do",
but that's useful too.

1) I was thinking narrowly of the scale_height function. Sorry that I didn't say that. In my defense, "p_surface" only appears in that function. The function can have either the global ref_surface_pressure passed to it, or the surface pressure, if that was actually useful for some reason. But inside the function they both serve as the reference pressure in the definition of scale height = log(p_ref/p), so I'm suggesting changes only in that function. I should know by now that if it's functioning correctly, it's not likely to be changed.

Yes, the cancellation of the log(p_ref) parts of a vertical distance calculation is exactly what prompted the definition and use of log(p) instead of log(p_ref/p) (at least in cam-fv); it gives the same distance for less computation. But the cancellation only works if the same p_ref is used for the calculation of the scale height location of both points. If I understood the discussion about other models decision to use only SH = log(p), they decided that because they ran into problems when they used log(p_surf/p) to calculate distances between points at different horizontal locations that had different p_surf. If they had used log(10^5/p) they wouldn't have had those problems. I'm hoping that everywhere in the code that passes p_surf(i,j) to scale_height for 2 locations always does it for the same i,j, so that a differing reference pressure is not causing unnoticed problems.

And as I pointed out before, log(p) implies a reference pressure of 1 Pa. log(p) below that has positive values, log(p) above it has negative. So there's a reference pressure (point where the coordinate = 0) even in the log(p) coordinate. It's just implicit (hidden). Again, I know it doesn't make any difference in distance calculations.

2), 3) And no_normalization_of_scale_height needs to be passed to the subroutines instead of .false.

4) I believe that the check for obs_too_high can be done even before the conversion to the vertical_localization_coord because obs_too_high can handle the original vertical coordinate of the ob.

[nancycollins]

hi kevin, you did post that already. i agree with points 2-4. i'm going to try one last time with point 1. if you still don't agree, please schedule a google meet or zoom call with me?

there are currently 2 ways to compute scale_heights in the code. one way is wrong. wrong, wrong, wrong. we all agree. however, it used to be the way scale_height was computed. if we change it then you can't do regression tests or compare results to old runs. so i'm resisting fixing this code because it needs to compute values the old (wrong) way. the default for the namelist item that controls this is to use the new way and this item should be documented as only set to false if you're trying to reproduce an old run or do regression testing.

the new way to compute scale height is correct and avoids a couple of divides which don't change the answer and so is more computationally efficient and simpler code. so i don't see anything that needs to be fixed here.

your math is completely correct and if we did want to fix the wrong code i'd agree 100%, but it's wrong on purpose. does this make sense? it's not just inertia or not wanting to rename variables.