McPhee implementation for sea ice

[cbegeman]

@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 https://github.com/xylar/palm_les_lanl/pull/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.

[katsmith133]

@cbegeman

Please forgive my ignorance of the code structure, since I haven't been using it that long, but for issues 1) and 2), would it not be sufficient to use an if statement for the McPhee parameterization to set usurf and vsurf on these lines: https://github.com/xylar/palm_les_lanl/blob/de6b0b7c39d7e7503391b7307460dacae3fb0cf9/trunk/SOURCE/surface_layer_fluxes_mod.f90#L1081-L1082 to a constant drift velocity that can be specified in the input file? Perhaps that's more of a inelegant hack though...

[cbegeman]

No apologies necessary!

I wouldn't go that route, mostly because you want the output to reflect the boundary conditions that you impose. You might want to actually test out my hypothesis that those lines from init_3d_model cause a problem because the boundary conditions appear to be properly set here: https://github.com/xylar/palm_les_lanl/blob/de6b0b7c39d7e7503391b7307460dacae3fb0cf9/trunk/SOURCE/boundary_conds.f90#L260-L267 that is, as long as u_init reflects u_surface in the input file you may be ok.

[katsmith133]

@cbegeman Alright, after a little bit of fumbling around, I was able to test those lines in init_3d_model.f90 and confirm that they do not affect the boundary conditions and u_init does indeed reflect the value of ug_surface in the input file. So that's good to go.

The fumbling around came from not realizing that these lines: https://github.com/xylar/palm_les_lanl/blob/1e43aa3af2c95d0eb5a71895e7baa73e598457f9/trunk/SOURCE/prognostic_equations.f90#L564-L568 will overwrite u_init(nzt+1) = ug_surface before it gets to boundary_conditions.f90 if rayleigh_damping_geostrophic = .False. (the default is .True.) is not specified in the input file. Is this the intended result of these lines? Just making sure that, if .True., those lines are meant to supersede the ug_surface value set in the input file.

[cbegeman]

@katsmith133 Good catch. We should probably change the default in the modules file to rayleigh_damping_geostrophic = False and add a warning in check_parameters.f90

if (rayleigh_damping_geostrophic .AND. (ug_surface != 0 .or. vg_surface != 0) ) then print('warning, surface velocities will be overwritten with geostrophic velocities') endif

or similar.

I could pretty quickly submit this as another bugfix, or you could tackle it. Up to you, just let me know.

[katsmith133]

@cbegeman As I am out of practice with all of that, I think maybe it would be good for me to take this simple task on and shake off the cobwebs. That being said, bear with me as I fumble through it, haha.

[cbegeman]

@katsmith133 Great, but don't hesitate to ping me if you have questions. w.r.t. the warning message, decide whether you want to terminate the simulation or continue, which is specified in an argument to the message routine.

[cbegeman]

@katsmith133 Perhaps this is obvious but I thought I should point out that the ug_surface namelist option generates geostrophic velocities through the coriolis term and is probably not how you would want to set surface velocities for sea ice drift that is generated from surface stress. You could do this by prescribing momentum flux boundary conditions, but there isn't currently an appropriate namelist parameter that I can find for setting velocity dirichlet boundary conditions (ending up as surf%u_surface) to non-zero, arbitrary (i.e., nongeostrophic) values.

[katsmith133]

@cbegeman Ah yes, I am seeing that now. At first glance, I only followed what was happening with ug(nzt+1) (which is used to set u_init(nzt+1), if I understand that correctly) and not how the rest of the profile was set or used (which I now see gets used in coriolis.f90). It should have been obvious, given that its called ug, haha. Thanks for mentioning that.

I feel as though creating a namelist parameter for setting non-geostrophic, non-zero velocity dirichlet boundary conditions is what we want for sea ice drift, as all of the hardware to compute the resultant spatially varying momentum flux is already there (wrapped up in the calc_uvw_abs subroutine), there just needs to be a way to set just u_init(nzt+1), without setting ug as well. Would you agree?

From what I can gather, u_init is set in check_parameters.f90 here: https://github.com/xylar/palm_les_lanl/blob/1e43aa3af2c95d0eb5a71895e7baa73e598457f9/trunk/SOURCE/check_parameters.f90#L1712-L1766 So perhaps an if statement here concerning a non-geostrophic, non-zero surface velocity that would set u_init(nzt+1) with a specified velocity value instead of ug, u_profile, or the pressure gradient would do it?

[cbegeman]

@katsmith133 Yes, that's right. I agree that dirichlet conditions are what you want and then use the MOST approach to get momentum fluxes.

An if-statement in the suggested spot would do it. You'll alsowant a new namelist parameter. I'd be pretty tempted to just change the namelist parameters ug_surface and accompanying vertical_gradient settings to just u_surface for consistency with initialization of scalars, and then then add a separate namelist option to determine whether to use the initial profile in the coriolis term. You'll probably also want to add u_top for your convenience as setting it up with vertical gradients would be cumbersome.