Adding a wind gustiness parameterization

[rljacob]

The merge of PR #1065 proposed to introduce a gustiness parameterization. This issue discusses how its currently done and how best to do it.

In module seq_flux_mct:
In seq_flux_init_mct: allocate prec_gust and set it to 0.

In seq_flux_atmocnexch_mct: define gust_fac, set it to HUGE(1.0) and then overwrite with seq_infodata_GetData. pass prec_gust and gust_fac to shr_flux_atmocn;

In seq_flux_atmocn_mct: set prec_gust to a2x%rAttr(index_a2x_Faxa_rainc,n); define gust_fac, set it to HUGE(1.0) and then overwrite with seq_infodata_GetData (from drv_in). pass prec_gust and gust_fac to shr_flux_atmocn

In module seq_infodata_mod
Add gust_fac to seq_infodata_type

In seq_infodata_Init, declare gust_fac and add to /seq_infodata_inparm/ namelist block. Set initial value to huge(1.0_SHR_KIND_R8), read value from seq_infodata_inparm namelist; copy value to infodata%gust_fac,

In seq_infodata_GetData_explicit: If present, gust_fac = infodata%gust_fac
In seq_infodata_PutData_explicit: if present, infodata%gust_fac = gust_fac
In seq_infodata_bcast: broadcast infodata%gust_fac
In seq_infodata_print: print out gust_fac value.

In module shr_flux_mod
In shr_flux_atmOcn: Add gust_fac and prec_gust to arguments.
Define new function:
ugust(gprec) = gust_faclog(1._R8+57801.6_R8gprec-3.55332096e7_R8*(gprec**2.0_R8))
Alter vmag calculation to include gustiness:
vmag = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) + ugust(min(prec_gust(n),6.94444e-4_R8)))

[mvertens]

The key issue that I am concerned about is that cesm never use the gustiness parameterization - not that its values happen to be 0 by default.

…

On Sat, Feb 4, 2017 at 2:48 PM, Robert Jacob ***@***.***> wrote: The merge of PR #1065 <#1065> proposed to introduce a gustiness parameterization. This issue discusses how its currently done and how best to do it. In seq_flux_init_mct: allocate prec_gust and set it to 0. In — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1105>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AHlxEw2RytucuK1ZqxAwLLkJv3pO8CQ3ks5rZPIXgaJpZM4L3S97> .

[mvertens]

I feel that since this is a new parameterization that ACME has introduced into share science code that also happens to be a key to CESM coupling, there should be an if-block that ensures that CESM does not use it. The way to do this is to pass cime_model to shr_flux_mct.F90 and have the following type of if-block

if (trim(cime_model) == 'acme') then
        vmag   = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) + ugust(min(prec_gust(n),6.94444e-4_R8)))
else
        !--- vmag+ugust (convective gustiness) Limit to a max precip 6 cm/day = 0.00069444 mm/s.
        !--- reverts to original formula if gust_fac=0 
        vmag   = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2)
end if

[gold2718]

Since this is a new issue, I will reiterate my concern with the proposed solution and my suggestion:

Every other science option is bracketed by logic, I do not see why this one is different. If folks do not like having model names (e.g., cesm & acme), then come up with a namelist parameter which indicates the scientific reason for the change and use that.

For instance, following usage in physics packages, introduce an XML variable, DO_GUSTINESS and an associated namelist parameter, do_gustiness, which will be always off in CESM (set by default in config_component_cesm.xml) and could always be on for ACME (config_component_acme.xml).

This makes the science choices explicit, just like other parameterization choices in the models.

[rljacob]

There already is logic to turn this off; if gust_fac is 0, its off. gust_fac=0 will make the ugust function evaluate to 0. So it functions as both a DO_GUSTNISS logical and a value in the parameterization.

[gold2718]

@rljacob, why are you insisting on changing how parameterizations have always been introduced? To take just one example, in the atmosphere physics package, the CLUBB SGS parameterization is bracketed by if(do_clubb_sgs), not by magic numbers.
It appears to me that you are upending many years of best practice by hiding a parameterization in math rather than by explicit code logic. Why do you feel this is important in this case? What SE principle is involved?

[rljacob]

If you're choosing between two large blocks of code or calling a different subroutine, then yes something like if(do_gustiness) would be appropriate. But this is just one term in one line so gust_fac=0 or not is the more elegant and better solution. We can clarify with some inline documentation and the documentation for gust_fac.

[gold2718]

Elegance is in the eye of the beholder and is not a software-engineering principle.
While the code is a bit shorter in this case, logically and scientifically, it is very similar to the introduction of the flux_diurnal parameterization.
Why do you feel that a departure from the standard protocol is "better" in this case?

[mvertens]

Here is another summary of my understanding:

ugust in shr_flux_atmOcn is the following function:

ugust(gprec) = gust_fac*log(1._R8+57801.6_R8*gprec-3.55332096e7_R8*(gprec**2.0_R8))
gust_fac is being sent to shr_flux_atmOcn and this is being used to evaluate the function ugust.

By default gust_fac is hug(1._r8) - int value in shr_infodata_mod.F90). However, it is zero by default in the namelist input for ugust_fac.

The issue I am concerned about is whether CESM can accept a non-default value that a user could set other than zero in a CESM release.

[rljacob]

I agree setting the initial value to HUGE(1.0) can be removed. That seems like an overly cautious way to make sure the namelist setting is taking effect.

[rljacob]

"flux_dirunal" is different because that logical controls if a completely different subroutine is called, one that has dozens of lines different from the non-diurnal version. This is one term in one line of code.

[jgfouca]

@rljacob @mvertens @gold2718 , looking at this conversation, it's not clear to me what the resolution is.

[rljacob]

Go ahead and merge it as-is. @mvertens will then do a PR with an implementation for CESM.

[mvertens]

It would help if just the following two minor changes were made to this PR. This would ensure that gust_fac is always being initialized as zero. I think that for the purposes of this PR all that needs to be done is: 1) Change line 486 of seq_infodata_mod.F90 gust_fac = huge(1.0_SHR_KIND_R8) <==== current gust_fac = 0._r8 <==== new at line 486 of seq_infodata_mod.F90. 2) Change line 910 of seq_flux_mct.F90: real(r8) :: gust_fac = huge(1.0_r8) !wind gust factor <==== current real(r8) :: gust_fac <==== new since you are obtaining gust_fac from seq_infodata in line 910 of seq_flux_mct.F90. call seq_infodata_GetData(infodata, & read_restart=read_restart, & dead_comps=dead_comps, & atm_nx=atm_nx, atm_ny=atm_ny, & ocn_nx=ocn_nx, ocn_ny=ocn_ny, & ocn_prognostic=ocn_prognostic, & flux_diurnal=flux_diurnal, & gust_fac = gust_fac ) Moving forwards, we will be getting input from CESM scientists as to how they want to see this science change implemented to satisfy CESM requirements. However, that can be put into a upcoming PR.

…

On Mon, Feb 6, 2017 at 11:50 AM, Robert Jacob ***@***.***> wrote: Go ahead and merge it as-is. @mvertens <https://github.com/mvertens> will then do a PR with an implementation for CESM. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1105 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AHlxE1vnZ0CG6V_jGBtamohPsMh1g56Vks5rZ2t_gaJpZM4L3S97> .

[rljacob]

Sure I'll make those changes.

[singhbalwinder]

Sorry, just saw this conversation. Implementation is not complete yet as currently I don't know a way to turn this on based on a particular compset. Originally, I wanted to turn this on (gust_fac>0) only for a new compset but at that time, I couldn't quite figure out how to do it. I think I talked with @rljacob at that time but I don't remember if we were able to find a solution.

This was my first time to make changes in the coupler codes and it sounds like I might have broken some existing rules. I will try to answer based on how I understand it works.

1. real(r8) :: gust_fac = huge(1.0_r8): I generally do that so that if this variable is used before it is initialized the simulation breaks (to be extra cautious). In this case, it is initialized right after so we can remove this. I also realize that it will be an automatic save variable (it is not required to be save and was not intentional) so I think it is a okay to just declare it without initializing. I followed flux_diurnal variable to add this variable.

2. In seq_flux_init_mct: allocate prec_gust and set it to 0.: I followed prec variable to add this variable. prec was also allocated and set to zero just after that, so I just followed suit. This can also be changed if necessary.

Rest of the code is just to read gust_fac from the namelist and pass it along so that it can be used in shr_flux_atmOcn. I can re-organize the code if it doesn't fit the current way of adding code to the coupler.

After implementing this code, I ran tests to make sure the code stays BFB with the old code.

[rljacob]

@mvertens what do you about not initializing gust_fac at all in the declaration? As @singhbalwinder said, that gives it an implied "save" attribute.

[mvertens]

@jacob - given all of these questions - lets just proceed with this PR as and clarify the best approach in a following PR.

…

On Mon, Feb 6, 2017 at 1:38 PM, Robert Jacob ***@***.***> wrote: @mvertens <https://github.com/mvertens> what do you about not initializing gust_fac at all in the declaration? As @singhbalwinder <https://github.com/singhbalwinder> said, that gives it an implied "save" attribute. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1105 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AHlxE8wnU7AG0IT1yBDzujoaFldvxR42ks5rZ4TMgaJpZM4L3S97> .