Since log law of the wall holds in a spatially-averaged sense, Bou-Zeid et al. (2015) filter the resolved velocity near the wall by 2 times the grid size. The current version of the code has some averaging but perhaps not enough:
https://github.com/xylar/palm_les_lanl/blob/c293074463bb0bda5a6796182b56793a1d163d91/trunk/SOURCE/surface_layer_fluxes_mod.f90#L1054-L1055
(these terms are used to define the surface stress). What do you think about averaging u(i,j,k),u(i+1,j,k),u(i,j+1,k),u(i+1,j+1,k) for the surface stress calculation instead (and likewise for v)? It's probably a minor effect but I don't think it would hurt to implement.

when multiple _xy or _xz or _yz variables are defined, the section check (from #11) fails. Reproduced here

!--       Make sure that section_nn is defined
          IF ( (data_output(i)(ilen-2:ilen) == '_xy') .AND. ( (section_xy(i) == -9999) ) ) THEN
             message_string = 'to output _xy variables, the depth of at least one xy_section must be specified via the namelist parameter section_xy'
             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( (data_output(i)(ilen-2:ilen) == '_xz') .AND. ( (section_xz(i) == -9999) ) ) THEN
             message_string = 'to output _xz variables, the depth of at least one xz_section must be specified via the namelist parameter section_xz'
             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( (data_output(i)(ilen-2:ilen) == '_yz') .AND. ( (section_yz(i) == -9999) ) ) THEN
             message_string = 'to output _yz variables, the depth of at least one yz_section must be specified via the namelist parameter section_yz'
             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
          ENDIF

With the section_xy(i) check, I think when you have more than one slice, or a variable not in the first position of the output list this will cause a failure incorrectly. I think what we want is to make sure there is at least one value in section_xy, section_xz, and/or section_yz that is not -9999.

Currently, the code uses a density profile for the atmosphere as rho_0 in the governing equations. I had been doing something analogous for the ocean, but @qingli411 pointed out that technically this term should be a constant in the ocean as a result of the boussinesq approx. I just wanted to run this past you all before I change it to a constant for the ocean case. The depth-averaged density profile will stick around in the numerator of the buoyancy term g*(rho-rho_ref(z))/rho_0.

I would like to add the ability to for PALM to solve for multiple passive scalars, instead of just one. As far as I see it, there are two primary ways to approach this and I wanted to get input on whether anyone thinks one way is better (or preferred).

1. Add an additional indicies to the s field, the size of which can be set in the input file by some integer namelist parameter called nscl or something like that. This requires I make quite a few changes to the code that might also affect things that others are doing, but is the most generic way to do this.

2. Follow what was done in the atmospheric part of the code and make a new variable called something like biochem_species which is wholly separate from the passive scalar s and can be solved for and messed with in a completely isolated part of the code. This is less generic and adds a whole new set of variables associated with it, but changes less of the code others may be using and will be more like an isolated module.

@cbegeman , @vanroekel , @xylar I know I have briefly spoken to you all about this a while ago, but just wondering if you had any stronger feelings, concerns, or suggestions...

On the master branch, the following error arises in combine_plot_fields when:

• there is any section_xy designated (I tried section_xy = 0 and section_xy=nz+1) AND
• the run is restarted

I haven't been able to locate any faults in the master branch code base so I'm wondering whether it might be the combine_plot_fields file.

*** combine_plot_fields ***
uncoupled run

NetCDF output enabled
XY-section: 8 file(s) found

forrtl: severe (67): input statement requires too much data, unit 110, file /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_master_restart/RUN_ifort.grizzly_hdf5_mpirun_test_oceanml/PLOT2D_XY_000000
Image PC Routine Line Source
combine_plot_fiel 000000000041D14E for__io_return Unknown Unknown
combine_plot_fiel 000000000043F571 for_read_seq_xmit Unknown Unknown
combine_plot_fiel 000000000040A918 Unknown Unknown Unknown
combine_plot_fiel 0000000000408FAE Unknown Unknown Unknown
libc-2.17.so 00002B7B2C4E23D5 __libc_start_main Unknown Unknown
combine_plot_fiel 0000000000408EA9 Unknown Unknown Unknown

@vanroekel Your thoughts?

@vanroekel @lconlon @qingli411 @xylar

As far as I can tell, a wall at the bottom boundary is hard-coded here:
https://github.com/xylar/palm_les_lanl/blob/42c59e36d9fffe82b0f4b3159be04f1860408ef7/trunk/SOURCE/init_grid.f90#L1911-L1912
because the nzb+1 index is used instead of nzb. topo is then used to define wall_flags_0.

I propose a new namelist parameter that designates whether the bottom or the top boundary is a wall, rather than hard-coding a wall at the bottom boundary or implementing a wall only if the boundary is associated with a constant flux layer.

Thoughts?

@xylar, When you have a chance, can you look over this version of the 3 equations? I think this is compatible with what we discussed in our last meeting.

After doing the tedious algebra to get the quadratic coefficients and solving it, I'm getting interface salinities that are way too low. I haven't found my mistake yet.

It seems to me that init_coupling, which is called from palm.f90, does not need to be called for uncoupled runs. Do you see any reason why I shouldn't add a line to only call this function if the coupling mode is not uncoupled?

I am motivated to do so because when the totalview debugger is used, it gets hung up on the standard input method of getting the coupling mode. Another solution to this problem (I think) would be to uncomment the method "get_environment_variable".

@qingli411 @vanroekel - perhaps you have thoughts about this since you may be more familiar with coupled runs?

@vanroekel @xylar I'm noticing that the modeled SGS momentum fluxes do not smoothly vary below the surface layer (SGS model is Moeng-Wyngaard). Saiki et al. (2010) discuss an issue with SGS models in the vicinity of boundaries, which is that the model can fail to produce a smooth transition between surface fluxes (prescribed) and vertical fluxes though grid cells near the boundary. As the resolved fluxes tend toward 0 near the boundary, there should be a concomitant increase in SGS fluxes to the prescribed surface flux. The fix that Saiki et al. propose is supposed to be included in PALM (Maronga et al., 2015) but I've noticed that it is not implemented in the version we've been working with. Do you mind taking a look at Saiki et al. and letting me know if you think this approach would be worth implementing? I'd be using the MOST shape functions from McPhee's work in equations 11 and 14. Interestingly, I'm not seeing the same discontinuities in SGS fluxes of pt and sa in the cases I'm running. I'm not sure why this is yet.
Thanks!

Does anyone have experience modifying the include path for one of the modules?

I've made modifications to "surface_layer_fluxes_mod.f90" such that it now uses the module eqn_state_seawater_mod (and the functions contained there).

The modification to the include path that I presumed to be necessary was adding a line to the Makefile
https://github.com/xylar/palm_les_lanl/blob/08bc2609da11ee9616f513c6013ca2eb28537330/trunk/SOURCE/Makefile#L1429-L1436
eqn_state_seawater.o \

However, I get the compiler error:
error #7002: Error in opening the compiled module file. Check INCLUDE paths.

Any thoughts? This isn't urgent, so enjoy the Thanksgiving holiday.

@katsmith133

I think that the McPhee implementation is mostly ready for the sea ice case. As I see it, there are basically 3 components to be concerned about:

1. Setting up the boundary conditions to allow you to prescribe sea ice drift velocities
2. Momentum fluxes
3. Scalar fluxes

(1) is the main component needing some tweaking. The issue that I see off-the-bat is these lines:
https://github.com/xylar/palm_les_lanl/blob/de6b0b7c39d7e7503391b7307460dacae3fb0cf9/trunk/SOURCE/init_3d_model.f90#L1470-L1474
which I think will set u(nzt+1,:,:) and v(nzt+1,:,:) to 0 when constant_flux_layer = 'top'. I haven't thought about this thoroughly but I think place to start is to evaluate whether you need these lines in the code at all, given that boundary_conds should be doing the work. I'm happy to chat more about this given that I've been working with/around wall_flags_0 for a while now ;)

(2) should basically be what you want. You're correct in stating that the relative velocity between the boundary and the first cell below the boundary is computed in calc_uvw_abs. I did make some minor changes to momentum flux computations in this commit that you may want to (re)base any further modifications on.

(3) should also be good to go. The main thing to figure out is what settings you want. You might want to refer to #55 and figure out whether you want gamma_z_dependent = True. I did notice a small bug if you want 2-d averaging of input pt,sa when gamma_z_dependent = False and most_xy_av = True. I'll submit this bugfix shortly.

I just noticed that there might be a bug in the diff terms of the RANS equations. The discretization of this term is missing a factor of 1/2

Lines 3820-3830 and 3930-3940, subroutines diffusion_e and diffusion_e_ij

@cbegeman It seems I may have spoken too soon when saying I was ok with merging #59. After adding commit 266ccaf to the PR to fix the non-AMD closure issue, I only re-tested the non-AMD closure bit-for-bit matching and not whether the AMD closure part of the code still worked. Turns out that commit, in particular I've narrowed it down to this part of the commit:
https://github.com/xylar/palm_les_lanl/blob/266ccaf93e950b432e41a36188a67bf9de12faaf/trunk/SOURCE/init_grid.f90#L1895-L1905
broke the AMD closure part of the code. What I mean by that is now when I try to run that same test case that you provided, the one in:
/lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_amd_mcphee_onesided,
PALM runs with the smallest possible time step and pretty much doesn't get anywhere. Here is beginning of the RUN_CONTROL output:

opposed to before that commit when AMD was working:

If I revert back to before that commit, then PALM runs appropriately with the AMD closure, but we get the same bit-for-bit mismatch between what were the master and amd_closure branches before. Hence the whole reason for this commit.

Not sure exactly what the problem is, but I have a question to hopefully help me better understand what is going on in this part of the code: what is the reason for applying the current code implementation? In the original code before any of the AMD commits, there was just this line for flat topography:
topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
Then before commit 266ccaf, it was:
topo(nzb+1:nzt,:,:) = IBSET( topo(nzb+1:nzt,:,:), 0 )
and then after the commit it was:
https://github.com/xylar/palm_les_lanl/blob/266ccaf93e950b432e41a36188a67bf9de12faaf/trunk/SOURCE/init_grid.f90#L1895-L1905
Is there a reason that the original code with just topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 ) would not work with AMD or was incorrect?

Hi @xylar,

I'm working out the details of the 3-equation parameterization with McPhee exchange velocities. I have a few implementation decisions to make and I wanted to give you (and others) the opportunity to weigh in. See individual comments below.

Cheers

I noticed that somewhere along the line we broke the restart capability. The error is generated at https://github.com/xylar/palm_les_lanl/blob/74b332fd5bd95b45efbca99b17b35ee1b8230805/trunk/SOURCE/netcdf_interface_mod.f90#L2462

I went back to @vanroekel 's old "palm_les_updates" version and verified that restart worked there. It did run successfully, but with a different error message:
errors in local file ENVPAR
some variables for steering may not be properly set

Do any of you know of a version of the code where you had a successful restart? Or do you have ideas about what the source of the issue is?
Thanks!

Do any of you know the mechanism by which palm makes modifications to the namelist file during simulations? Specifically, I've noticed that initializing_actions is truncated after the first "action", e.g.,
initializing_actions = 'set_constant_parameters initialize_vortex',
becomes
initializing_actions = 'set_constant_parameters'
Thanks!

It would be nice to no longer see this warning when compiling PALM
poismg_mod.f90(1210): warning #8093: A do-variable within a DO body shall not appear in a variable definition context. [J]
CALL MPI_ABORT(MPI_COMM_WORLD,i,j)

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

Let's use the existing file structure for storing test cases:
palm_dir/trunk/TESTS/JOBS/INPUT/namelist_file
palm_dir/trunk/TESTS/JOBS/REFERENCE/output_file

There is only one test case in the directory right now. I don't know if this changed with the latest release.

@xylar @vanroekel @katsmith133 I'm running into high km values in the second grid level below melting ice interfaces as a result of high buoyancy gradients. This problem was noted here, at which time the problem seemed limited to just the first grid level below the interface.

I'd like to get your thoughts on possible remedies, whether that's in code changes to AMD, a more restrictive km_max or different simulation configuration. Thanks!


Note, km_max = 0.1.
Namelist file here: /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_amd_mcphee_lowres_disturb5h_cd1e-2_rotation