Changes to density and buoyancy treatment in ocean cases

[cbegeman]

• new rho_reference variable on u- and c-grids used consistently in atmosphere and ocean cases
• define buoyancy terms using density for ocean cases
• removed TKE production from surface buoyancy fluxes, inconsistent with SGS scheme
• in ocean only cases, surface pressure is input in decibars
• boundary fluxes are input in "natural" units (W/m^2 for heat flux, Pa for momentum flux, [kg/m3][m/s][PSU] for salt flux)

[cbegeman]

I'm still testing, just making it available for viewing.

[qingli411]

@vanroekel I think the pressure solver changes here look good to me. For the atmosphere with anelastic approximation rho_ref is a function of z so it should be there to compute the divergence term. For the ocean it doesn't matter as we always assume Boussinesq so rho_ref is a constant. But since the pressure is used in unit of pascal elsewhere, the density should be there. Otherwise we need to do a conversion.

[cbegeman]

@qingli411 Thanks for taking a look at the code. I've added some more changes to this PR. I'm having difficulty writing restart data in coupled runs in this branch. @vanroekel suggested that it could be an issue of having too few grid cells per processor, but I increased the number of grid cells and the issue persists. If you have any ideas on how to investigate this issue, please let me know. Thanks!

The error message is ... writing restart data Error in `./palm': double free or corruption (!prev): 0x0000000002b640d0

[qingli411]

@cbegeman I've never seen this error message. What are the number of processors you use for the atmosphere and ocean? I once tried to use more processors for the ocean than the atmosphere and palm failed to write restart files for the ocean. I'm not sure what was wrong. But it seems to correctly generate the restart files as long as I use the same number of processors for both.

[cbegeman]

@qingli411 I'm using the same number of processors for each: 1 node, 16 or 32 processors each for atmosphere and ocean. The only substantive changes I made to the surface_coupler deal with surface fluxes, but the error persists even when surface fluxes are zero.

[vanroekel]

@cbegeman, confirming that I get your error in coupled that you show above. I'll try see if something jumps out.

[cbegeman]

@vanroekel I was able to track down the issue with resolved TKE. It turns out there was an inconsistency in units ascribed to flow_statistics variables u*2, v*2, and w*2. Now that this fix has been made, I think we're ready to merge.

[vanroekel]

@cbegeman could you explain your last comment a bit more? It's not clear from the PR what you fixed and the plots seem similar to what you showed 1.5 weeks ago, but I do see the resolved TKE here now. It looks like your last commit fixed the issue for coupled runs maybe? Was the fix just in your plotting?

[cbegeman]

@vanroekel There are changes in the code, but they pertain to outputting. This is the relevant commit. I realized that there was an inconsistency in units between the "flow statistics" variables representing u*2, v*2, and w*2. w*2 was in [Pa] while u*2 and v*2 were in [m^2/s^2]. The "master" branch doesn't really see these differences because rho_air was used instead of rho_ocean. Now u*2, v*2 are also in [Pa] and converted to [m^2/s^2] for output using momentumflux_output_conversion factor of rho_ref_zw.

Adding factor of rho_ref_zw, and similarly for sums_vs2_ws_l: https://github.com/xylar/palm_les_lanl/blob/70415ae476f383f6a48851686f855abfd802134a/trunk/SOURCE/advec_ws.f90#L2071-L2079

Dividing by factor of rho_ref_zw: https://github.com/xylar/palm_les_lanl/blob/70415ae476f383f6a48851686f855abfd802134a/trunk/SOURCE/flow_statistics.f90#L445-L450

[vanroekel]

Thanks for the explanation @cbegeman, this makes sense now. I think this is ready to go.

[xylar]

@qingli411, do you want to test this out and approve it or ask for changes? If not, please remove yourself as a reviewer.

[qingli411]

@cbegeman I would like to do a few tests with my cases. If I understand if correctly I will need to modify my namelist to use "natural" units for the surface fluxes, right? If so I think it would be helpful to change the example namelists in this PR to make them consistent with the code.

[cbegeman]

@qingli411 Sounds great. You don't need to make any changes to the namelist file. Let me know if you have any questions.

[qingli411]

@cbegeman I did two quick tests with parameters from McWilliams et al., 1997: a shear case (MSM97-ST) and a Langmuir case (MSM97-LT). Both run for 24 hours. The namelists I used are /usr/projects/climate/qingli/PALM_ocean_MSM97-ST_p3d /usr/projects/climate/qingli/PALM_ocean_MSM97-LT_p3d

These namelist are the same as I used for my previous tests. I'm a little confused here: I thought you modified the code to input the surface fluxes using natural units like W/m^2 (your last bullet point).

Anyway, I have no problem running case MSM97-ST, but got segmentation fault for MSM97-LT. I haven't figured out why but I will keep looking. Could you also try the MSM97-LT namelist to see if you can reproduce this error?

For the MSM97-ST case, the results are a little bit different from the master branch, especially for the temperature variance and temperature flux. See here. These are mean profiles averaged over the last inertial period (~16 hours). Solid lines are the rho_ocean branch, dashed lines are the master branch. In particular, the SGS temperature flux profile doesn't look right to me...

Any thoughts?

[cbegeman]

@qingli411 You can input with natural units if you set flux_input_mode to 'dynamic'. I kept this original feature of the code because it didn't seem necessary to change it. Yes, I reproduce the seg fault. The only difference between the two runs is Stokes forcing, correct? It looks like the only place where rho_ref_zu appears is https://github.com/xylar/palm_les_lanl/blob/70415ae476f383f6a48851686f855abfd802134a/trunk/SOURCE/stokes_force_mod.f90#L550-L554 It's unclear to me what the issue is. I'll look into it, and since you're more familiar with the stokes part of the code, perhaps you can look into it as well? I see similar pt*2 profiles in the NCAR free convection test case, so I'm not yet sure why the difference is greater in your case. @vanroekel, your thoughts? With respect to the resolved temperature flux, I think that can be attributed to the correction I made in buoyancy terms (now dependent on rho rather than pt).

[qingli411]

@cbegeman OK. I see.

Yes, I will look in to the seg fault problem for the Langmuir case.

[vanroekel]

@qingli411 thanks for the testing. I'll take a look at the sea fault as well. In the meantime, could you also plot the NCAR LES line on at least some of these figures for a baseline?

@cbegeman and @qingli411, what is interesting is that the new changes seem to better resolve the turbulence at the boundary layer base. This is suggested by the weaker SGS temperature flux and stronger temperature variance. Given that I'm sure these are the same resolution this makes very little sense. I'll dig into it.

[vanroekel]

@qingli411 can you retry your tests with the changes I suggest here?

[cbegeman]

Thanks @vanroekel! I'll make these changes now.

[cbegeman]

@qingli411 I just ran the code with these changes and compared it with the NCAR free-convection test case. I also ran the stokes case, results at /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_MSM97-LT

[qingli411]

@cbegeman After the above fix the results compared with the master branch are here: MSM97-ST http://portal.nersc.gov/project/acme/qingli/plot_palm_mean_profile_rho_ocean_SGS_flux_MSM97-ST.html

MSM97-LT http://portal.nersc.gov/project/acme/qingli/plot_palm_mean_profile_rho_ocean_SGS_flux_MSM97-LT.html

I haven't test the coupled case but will do that next.

I notice another small bug not related to this PR while doing the tests, but it would be nice to fix it here. I think L239 of eqn_state_seawater.f90 where rho_ocean is computed again can be removed. It's already computed on L234. Thanks!

[cbegeman]

@qingli411 @vanroekel Any surface momentum flux on the master branch should be off by a factor of rho_air_zw(nzt+1)/rho_air = rho_air_zw(nzt+1), which can be about 1.2. All links below are to code on the master branch. The difference should be reduced when 'dynamic' fluxes are used. https://github.com/xylar/palm_les_lanl/blob/4a029194fe7e9e91e1718ed13cc63caad0818a79/trunk/SOURCE/init_3d_model.f90#L883-L902 https://github.com/xylar/palm_les_lanl/blob/4a029194fe7e9e91e1718ed13cc63caad0818a79/trunk/SOURCE/init_3d_model.f90#L917-L925 https://github.com/xylar/palm_les_lanl/blob/4a029194fe7e9e91e1718ed13cc63caad0818a79/trunk/SOURCE/surface_mod.f90#L2318-L2320 https://github.com/xylar/palm_les_lanl/blob/4a029194fe7e9e91e1718ed13cc63caad0818a79/trunk/SOURCE/diffusion_u.f90#L346-L347

[vanroekel]

@cbegeman thanks for digging into the density differences! Once @qingli411 finishes the longer term comparison, I'm ready to sign off if all goes well with that test.

[qingli411]

@cbegeman Thanks for checking this! Here are the comparison with master for the longer run. Case MSM97-ST Case MSM97-LT

These profiles are averaged over an inertial period much later-long enough after the spin-up. I also included the profiles for NCARLES in the comparison so the legend is a bit different: Black: this PR Red: master Blue: ncarles

I think this PR is good to go.

@vanroekel I notice the differences between PALM (both this PR and master) and NCARLES are greater than my previous tests. I'm still testing and trying to figure out why but it seems that the results of PALM have greater sensitivity to horizontal resolution. My previous tests were on a higher resolution.

[cbegeman]

@qingli411 Thanks for running these test cases!

[vanroekel]

Thanks for the test cases an extra review @qingli411. I agree that PALM is more resolution dependent than NCAR LES, I'm fairly sure this due to PALM not having spectral derivatives.