adds reduced formulation (down gradient) option for horiz momentum flux terms

[katsmith133]

Adds option to use a down gradient formulation for the horizontal momentum flux terms. Set config_use_grad_diff_horiz_mom_flux = .true. to use and adjust coefficient C_mom_flux.

Running a "wind2" test case (tau = 0.02, no cooling), with the following namelist parameter values:
config_adc_enable = .true.
config_adc_timestep = 1.0
config_adc_use_old_length_scale = .true.
config_adc_truncate_tend = .true.
config_adc_decimals_to_keep = 12
config_adc_tau_o = 100
config_adc_length_multiple = 2.0
config_adc_epsilon = 0.0
config_adc_use_single_column = .true.
config_adc_sigmat = 0.72
config_adc_Ko = 4.574296
config_adc_c_b_tracer = 0.33
config_adc_c_b = 0.5
config_adc_alpha_0 = 0.8
config_adc_alpha_1 = 0.984
config_adc_alpha_2 = 0.568
config_adc_alpha_tracer1 = 0.2175
config_adc_alpha_tracer2 = 0.3625
config_adc_c11 = 0.1
config_adc_Cmom = 0.5
config_adc_Ctherm = 0.5
config_adc_Cmom_w3 = 5.0
config_adc_c_slow = 2.5
config_adc_slow_w_factor = 0.121
config_adc_c_slow_tracer = 5.936
config_adc_dissipation_constant = 10.0
config_adc_CwwE = 1.0
config_adc_CwwD = 1.5
config_adc_kappaFL = 0.0
config_adc_kappaVAR = 0.0
config_adc_kappaW3 = 0.0
config_adc_bc_wstar = 0.3
config_adc_frictionVelocityMin = 1.0e-5
config_adc_bc_const = 3.0
config_adc_bc_const_wp2 = 0.0
config_adc_use_splat_parameterization = .true.
config_adc_splat_tend_max = 1.0e-5
config_adc_splat_wp2_val = 3.0
config_adc_up2_vp2_factor = 1.0

Below is a comparison between PALM LES (black), the default horiz. momentum flux formulation (blue), and the new down gradient formulation with c_mom_flux = 0.05 (red) at 24 hrs (top row) and 120 hrs (bottom row).

The new formulation appears to not be unstable like the default. Below is a hovmoller plot of temperature for the default formulation (left) and the new down gradient formulation with c_mom_flux = 0.05 (right).

[vanroekel]

@katsmith133 this is great, very encouraging results as well. I have just one comment to improve computational efficiency.

Also, can you run a few more wind related cases (especially mixed convection/wind?). Thanks for working on this!

[katsmith133]

Here is a mixed convection/wind case ("t2wind" 0.01N/m2 and -100W/m2) using the same parameter values as before, comparing PALM LES (black), the default horiz. momentum flux formulation (blue), and the new down gradient formulation with c_mom_flux = 0.05 (red) at 24 hrs (top row) and 120 hrs (bottom row)

And a hovmoller plot of temperature for the LES (left), default formulation (middle), and the new down gradient formulation with c_mom_flux = 0.05 (right):

Again the new formulation appears to be more stable. I think with some additional tuning (I have done none here) we could get it to match the LES even more. I will run a few more wind cases and post as I finish them...

[BrodiePearson]

@katsmith133 It looks like the current down-gradient momentum flux assumption in your PR looks like:
uw(i2,k,iCell) = -C_mom_flux*sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell)))*length(k,iCell)*Uz

I think that the diagnostic formulation for the momentum fluxes that we discussed earlier, calculated by rearranging the tendency equations and assuming that the time derivative is zero, would result in the equation below. There might also be a need to adjust the first indices, for example I'm not sure if tauVel should have an i index here. Note the - w2(i1,k,iCell)*Uz/tauVel(k,iCell) term below has a similar form to the down-gradient model you coded. I don't know if you will need to add all the terms below, but you might get an improvement in the convective case simply by adding the third-order flux divergence and/or the (non-pressure) buoyancy terms :

uw(i2,k,iCell) = -( 1.0_RKIND/tauVel(k,iCell) ) * ( &
-(uw2(i1,k-1,iCell) - uw2(i1,k,iCell)) / dzmid &
+ 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND &
+ (alpha_1 - alpha_2)*u2(i1,k,iCell) + (alpha_1 + alpha_2)*w2(i1,k,iCell))*Uz &
+ 0.5_RKIND*(alpha_1 - alpha_2)*uv(i1,k,iCell)*Vz &
- C_b*grav*(alphaT(k,iCell)*ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell)) &
- w2(i1,k,iCell)*Uz &
+ grav*(alphaT(k,iCell)*ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell)) &
+ kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND )

The analogous equation for vw would look like:

vw(i2,k,iCell) = -( 1.0_RKIND/tauVel(k,iCell) ) * ( &
-(vw2(i1,k-1,iCell) - vw2(i1,k,iCell)) / dzmid &
+ 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND &
+ (alpha_1 - alpha_2)*v2(i1,k,iCell) + (alpha_1 + alpha_2)*w2(i1,k,iCell))*Vz &
+ 0.5_RKIND*(alpha_1 - alpha_2)*uv(i1,k,iCell)*Uz &
- C_b*grav*(alphaT(k,iCell)*vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell)) &
- w2(i1,k,iCell)*Vz &
+ grav*(alphaT(k,iCell)*vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell)) &
+ kappa_FL*(vw(i1,k-1,iCell) - vw(i1,k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND )

[katsmith133]

Great, thanks @BrodiePearson! This makes sense now that I see it in writing. I will work on putting that in and seeing if it helps in the mixed convection/wind case.

On another note, @amrapallig asked for plots of sigma for the cases. Here is the purely wind case:

And the mixed convection/wind case: