Surface fluxes in PALM

[qingli411]

I just took a closer look at the use of surface fluxes surf_def_h(*)%shf, surf_def_h(*)usws and surf_def_h(*)vsws in PALM. I think we probably don't need to explicitly replace rho_air by rho_ocean in the code for the ocean mode. Even for the atmosphere I don't think we are going to run the code in the 'anelastic' approximation mode. With Boussinesq-approximation the reference density is essentially a constant for both atmosphere and ocean. So rho_air and therefore heatflux_input_conversion and momentumflux_input_conversion are simply constants. In addition, the reference density only appears in the pressure solver. The buoyancy term has already been taken care of by subroutine buoyancy() by using rho_ocean for the ocean mode and pt for the atmosphere mode. I think almost everywhere in the code the variable surf_def_h(*)%shf represents $\overline{w'T'}_0 \rho_{air}$ and the variable surf_def_h(*)%usws represents $\overline{u'w'}_0 \rho_{air}$. There variables are divided by rho_air to get the kinematic fluxes in the temperature and momentum equations. So we only need to be careful to make sure the kinematic fluxes $\overline{w'T'}_0$ and $\overline{u'w'}_0$ etc are correct for the atmosphere or the ocean, which are the input parameters anyway. And wherever we add new lines in which these variables are used, they should be divided by rho_air to get the kinematic fluxes.

I think there might be a bug in the coupler. In surface_coupler where the code converts the surf_def_h(*)%shf and surf_def_h(*)%usws etc from the atmosphere values to the ocean values, there seems to be a missing factor rho_air. Without rho_air these fluxes surf_def_h(*)%shf and surf_def_h(*)%usws etc will become kinematic, which is inconsistent with other places. https://github.com/xylar/palm_les_lanl/blob/92d208fd830d01e97301e0f47cdb0f31fdec6152/trunk/SOURCE/surface_coupler.f90#L499-L503

What do you think? @cbegeman @vanroekel

[xylar]

@cbegeman will respond as well, I'm sure. But the issue for ocean-only simulations where the surface is not the atmosphere-ocean interface is that the fact that dynamics fluxes are scaled by rho_air is very confusing. There is no interface with air, so heat and momentum fluxes would be more appropriately scaled by either rho_ocean or rho_ice in our case. If our kinematic surface momentum and heat fluxes have to be scaled to be a factor of 1000 larger than they are in reality because they are conceptually the kinematic fluxes on the atmosphere side of the atmosphere/ocean interface, that is very confusing for our problem.

[cbegeman]

If you give up rho_air as a function of z and replace it with a constant, then my preference would be to proceed as follows: In the namelist file, add variables rho_air_ref for atmosphere runs and rho_ocean_ref for ocean runs and both for coupled runs. In the code, assign rho_air_ref to rho_ref for atmosphere cores and rho_ocean_ref to rho_ref for ocean cores. Define the factor rho_air_ref/rho_ocean_ref for conversions in the surface coupler between ocean and atmosphere. Define fluxes as rho_ref*\overline{w'X'}

This way we don't have the need to keep rho_air around in uncoupled ocean runs.

[vanroekel]

@qingli411 regarding your bug, it certainly appears to be a bug. It seems odd to divide by rho_ocean. For the broader question, I think my preference would be for there to be a general "rho" variable that is specified on a per core bases (either atmosphere or ocean). If you write down the equations for the ocean, the turbulent fluxes really should be rho_ref*w'T', where rho_ref is ocean specific. To me the cleanest solution would be to specify fluxes in "natural units" (e.g. Wm^-2 for heat flux, Pa for momentum flux) and let each core do the appropriate conversion.

[qingli411]

@cbegeman I agree. We could even keep the z-dependence of rho_ref (i.e., simply rename rho_air to rho_ref in most places) to allow 'anelastic' approximation for the atmosphere, but set rho_ref to rho_air_ref or rho_ocean_ref at all levels for 'boussinesq' approximation.

@vanroekel I think I would prefer to Carolyn's approach, to specify the fluxes (shf, usws, etc.) as rho_ref*\overline{w'X'}, instead of in their "natural units". So that we can minimize the code modifications (if we want to merge our code back to the PALM repository later) and avoid using the specific heat of air and water in many places. Though I agree using the "natural units" would be much cleaner...

All I was saying in my previous comments here and in #15 is that we should avoid redundant code for atmosphere and ocean (e.g., if (ocean) then ... else ... end if blocks) in the solver given that the governing equations under 'boussinesq' approximation are essentially the same for the atmosphere and ocean except a few constant conversion factors.

[cbegeman]

@qingli411 It does seem that the cleanest way to retain the z-dependence for the atmosphere cases is to carry around a rho_ref vector for all cases, otherwise we're stuck with lots of if-then clauses.

@qingli411 @vanroekel How about this compromise: we specify fluxes in their "natural units" in the namelist file but store fluxes as rho_ref*\overline{w'X'} in the body of the code.

Let me know when we've settled on an approach you all are happy with. I'd like to start making these changes early this week.

[vanroekel]

@cbegeman I think having a conversion from "natural units" to rho_ref \overline{w'X"} when the nameless file is read in sounds good to me.

[qingli411]

@cbegeman Yes, this sounds good to me as well. Thanks for doing this! This is very helpful.

[cbegeman]

Addressed by PR#30