WAM Physics: Technically include Russell's Empirical Poynting flux model

[ZhuxiaoLi66]

After the outline test on the comparison of the integral poynting flux from both WAM and Russell's Empirical Poyntingflux model, it has been verified that the ratio of the integral flux is applicable in WAM as the the Joule Heating factor. The code work to include Russell's Empirical Poyntingflux model into WAM is under construction. After this work done, we need to verify the Joule heating, neutral temperature and composition output from WAM, especially during storm time.

[akubaryk]

The new branch is epf, commit here: 5f739ad3

• Everything epf related has been moved into phys instead of gsm/phys.
• I performed a minor rework of read_model_parameters so that the string length was not a necessary part of the function call, and renamed those files appropriately. I also had to format the file properly for the idea_phys library compile, which really just meant removing & at the end of lines for line continuation. I also removed all trailing whitespace for cleanliness.
• read_model_parameters is called during efield_init and never again (much like the read_coef from Weimer).
• call Pflux_NS_Bx occurs after the tilt calculation in prep_weimer so that only has to be calculated once, and the result is immediately assigned into pf_nh_integral.
• The st_fac calculation was moved outside a do loop because it will never change per-timestep.

If you do a test run, you'll see that jh_nh_integral and pf_nh_integral are printed out at each timestep -- I tried to make only the lead proc print it out, but there are still a lot of prints because of how physics functions. That will need to be removed before we consider this code for pull. We can also send pf_nh_integral out of phys for write into the jh.nc4 file that the previous commit in poynting_flux introduced -- let me know if that's something you want.

[ZhuxiaoLi66]

@akubaryk good works, thanks for the great help. it will be very helpful if we can get pf_nh_integral output in jh.nc4.

[timfullerrowell]

Zhuxiao,

Thinking more about what you said. You are right, a better correlation would be a scatter plot of the current ratio with a ratio calculated from the PF(t)/JH(t-1).

Tim

On Tue, Jan 7, 2020 at 5:36 PM ZhuxiaoLi notifications@github.com wrote:

@akubaryk https://github.com/akubaryk good works, thanks for the great help. it will be very helpful if we can get pf_nh_integral output in jh.nc4.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/NOAA-SWPC/WAM/issues/14?email_source=notifications&email_token=AH5BFOFSRPY7LE74Y3M5TA3Q4UN2RA5CNFSM4KDJJ4CKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIKY3WI#issuecomment-571837913, or unsubscribe https://github.com/notifications/unsubscribe-auth/AH5BFOBRVUSZCF57LPUQK5TQ4UN2RANCNFSM4KDJJ4CA .

[akubaryk]

epf branch has been updated to remove the print and output in jh.nc4 instead.

[ZhuxiaoLi66]

The following is the 1stp time lag corr of flux ratio for both PF(t)/JH(t-1) and PF(t-1)/JH(t-1) cases. we have got more noise in the former case.

[ZhuxiaoLi66]

The following plots are the corresponding ones for 2015 St. Patrick case.

[akubaryk]

Are these values with the JH factor applied? My quick glance at one hour of output showed PF values that were ~double the (unscaled) northern hemisphere JH.

[ZhuxiaoLi66]

it is not the future JH factor. Since we have to use the former time step ratio factor (pf/jh), Tim and I want to see how much the factor changes after 1 time step. it is the ratio of the factors with 1 stp time lag.

[ZhuxiaoLi66]

[ZhuxiaoLi66]

@akubaryk the above plot is based the output in jh.nc4 which has all valid reasonable values for both pf_nh and jh_nh variables. While when I printed out the two variables (pf integral and jh integral) in idea_ion.f, it seems for the first time step, jh_nh_integral is 0. (make sense for the first step, jh_nh_integral hasn't been calculated, while still an issue need to deal with later which will cause the Jh_fac become infinity), jh_nh_integral looks good at the second time step, while become NaN since 3rd time step. any idea on this problem? thanks. the code:/scratch1/NCEPDEV/swpc/Zhuxiao.Li/save/git/GSMWAM-IPE_epf/WAM/src/phys/idea_ion.f printing output: /scratch1/NCEPDEV/swpc/Zhuxiao.Li/scrub/PTMP/GSM_ept_20150311_storm_ori/fcst.2523296

[akubaryk]

Probably has to do with dtdt going to infinity at the first timestep because jh_fac is infinity? The jh_nh_integral output at kdt=2 is actually from kdt=1 (before the physics are applied), but after that you'll have NaN's because the jh_fac is garbage for kdt=1. Try putting a check in for jh_fac (.gt. 10.0?), else use the old method.

This approach is going to break restarts, by the way.

[ZhuxiaoLi66]

Got it, thanks.

[ZhuxiaoLi66]

is it easy to always set the constant value of jh_nh_integral for the first time step? by the way, could we avoid the break of the restart?

[akubaryk]

The values would have to be written out and ingested with the rest of the restart data.

[ZhuxiaoLi66]

I fixed this issue with the following code in idea_ion.f if(kdt_interval.eq.1) then jh_fac = 2.0 else jh_fac = pf_nh_integral/jh_nh_integral endif

don't know how restart codes can handle this?

[akubaryk]

We will have to write out and read in the jh_nh_integral so it's available at the first timestep.

[ZhuxiaoLi66]

got it, thanks.

[ZhuxiaoLi66]

The following plots shows the performance of the epf branch with the new JH_factor (poynting_flux_integral/Jh_integral), the neutral temperature evolution during the 2015 St. Patrick storm behaves every similar with our original JH parameterization.

[akubaryk]

In order to do seasonal runs, we need to have the restart working properly, or find a way to have the Joule heating integral available at the same timestep. The former is going to be far easier for now, so we'll plan to implement a simple write out and read in with the rest of the lsidea=.true. restart functionality. This will be added to the epf branch -- if you can update the remote branch with your latest code changes, @ZhuxiaoLi66, I'll get on this right away.

[ZhuxiaoLi66]

@akubaryk already pushed the code. Thanks for the help. By the way, It will be perfect that you can find a way to have the JH integral available at the same timestep, the JH_factor will change with time more consistently, as showed in the time lag correlation plot about the ratio (JH_factor).

[akubaryk]

I think I have a fairly clean way to do same-timestep. It'll take a bit more time but I'll update when it's ready.

[akubaryk]

epf branch has been updated again, this commit explains the new physics workflow: 9849897. The code was tested on Mars yesterday and seems to work as intended. I did a crude comparison because analysis tools are limited on Phase 3, but it seemed like the 2013 storm peak was ~50K cooler with no net impact before or after the storm. Will produce a 2015 storm run for Zhuxiao for plotting.

[ZhuxiaoLi66]

The following is the output from Adam's new code, kind of colder than the original parameterization way. (compare them with the above plots named 'VBz_25000'.

[ZhuxiaoLi66]

the following is the ratio (pf_integral(t)/jh_integral(t-1)) for 20131120 storm by the epf_jhfac (before the latest revise).

[ZhuxiaoLi66]

@akubaryk how is your debugging going? I checked the WAM codes again this morning, i think one doubtful place is in gsm/phys/do_physics_one_step.f when it called gloopb before doing the extra Joule heating replenish. I mean in gloopb, the calculation of 'dtdt' in idea_phys effected the other phys calculation in gbphys (grep 'dtdt' in it) which called by gloopb later. So it is not that fine and simple that we can do the replenish on 'dtdt' later after the call of gloopb in do_physics_one_step. Just an brief idea, please double check and do some write out if necessary.

[ZhuxiaoLi66]

we did some validation on the jhfac test (apply former jh_integral and current poyntingflux_integral to calculate the JH_factor in idea_ion.f). the results against satellite data and the baseline run with the original JH_factor parameterization seems pretty off, especially during storm time.

[ZhuxiaoLi66]

The following plots are the number density validation against GRACE data along it orbit.

[ZhuxiaoLi66]

this plot is the validation for original JH_factor parameterization.

[ZhuxiaoLi66]

The difference of neutral temp. at 250km between the current test and the original JH_factor looks consistent with their number density difference in Satellite Validation. the temp. is much cooler at the intense storm time in poyntingflux_JHfactor run.

[ZhuxiaoLi66]

The following is the number density validation against CHAMP data for 20031120 storm.

[timfullerrowell]

Thanks Zhuxiao, Is this with the F10 variation included or excluded. If the former can we try the latter as well. Tim

On Fri, Jan 31, 2020 at 6:29 AM ZhuxiaoLi notifications@github.com wrote:

The following is the number density validation against CHAMP data for 20031120 storm. [image: num_den_jhfac_WAM_CHAMP_20031120] https://user-images.githubusercontent.com/22546571/73491929-9369b800-43a7-11ea-8b13-a11ad469eac7.png

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/NOAA-SWPC/WAM/issues/14?email_source=notifications&email_token=AH5BFODW2ONEJVWQCLO2JHLRANBDLA5CNFSM4KDJJ4CKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEKMT7KI#issuecomment-580468649, or unsubscribe https://github.com/notifications/unsubscribe-auth/AH5BFOGTKUZ64J7BURE4EY3RANBDLANCNFSM4KDJJ4CA .

[ZhuxiaoLi66]

will do this after merge the epf and develop branch in GSMWAM-IPE.