Grab Weiyu's former codes into the tip WAM development for the NetCDF output of WAM-IPE interface variables

[ZhuxiaoLi66]

Doing efforts to manually grab Weiyu's former codes (several subroutines in NEMS/src/atmos/gsm/dyn) into the current WAM develop version code for the NetCDF output of several very helpful WAM-IPE interface variables (total number density, mean mass...). The validation idea for Weiyu's codes is that we would like to get the bit by bit regular WAM sig output between the tip WAM version and the one after adding his codes.
The reference former issue is NOAA-SWPC/WAM-IPE#134.

[ZhuxiaoLi66]

Added the subroutines already, and correspondingly added several output namelist variables in namelist and run scripts. Have got the bit by bit sig output, while haven't gotten the expected NetCDF output files yet. Will double check the codes path and the settings again.

[akubaryk]

I can help you with this as soon as there's a branch available from the repository.

[ZhuxiaoLi66]

Thanks, Adam. will do soon.

[ZhuxiaoLi66]

@akubaryk already pushed a branch to the repository, named WAM_NetCDF_output.

[ZhuxiaoLi66]

Adam, Could you please have a look at the 'WAM_NetCDF_output' branch on Github I pushed several weeks ago? Do hope your expertise in dycore and WAM output can help this. I reviewed and checked my former works, still don't get any clue for not getting the NetCDF output, even I feel it is very close to done. Works Summary:

1. Weiyu's codes are mainly sit in NEMS/src/atmos/gsm/dyn/grid_collect_ipe*.f if you diff makefile makefile_ori, you will get < grid_collect_ipe_grads.o \ < grid_collect_ipe_restart.o \ < grid_collect_ipe_netcdf.o \ < wrtout_netcdf.o \ < read_netcdf_test.o \

   grid_collect_ipe.o \

[ZhuxiaoLi66]

the main subroutine for NetCDF diagnostic variable output is: subroutine atmgg_move_ipe
in code of grid_collect_ipe_grads.f I checked its 'print' output in fcst file, it seems that this subroutine doesn't been called.

[ZhuxiaoLi66]

In addition, please check the file 'namelist_dynamics_def.f' in the ./gms/dyn directory. and please diff namelist_dynamics_def.f namelist_dynamics_def.f_ori

2. In ./oldtests I changed the following files correspondingly to add the namelist variables namelist variables: export FHOUT_GRADS=1 export FHOUT_NC=1 export GRADS_OUTPUT=.true. export NC_OUTPUT=.true. ./oldtest files,
   gfs_fcst_run.IN_Linux rt_gfs.sh rt_gfs.sh_restart

about the namelist variable, do you think the following two also matters? I marked them in my run script.

export WAM_IPE_CPL_RST_INPUT=.false.

export WAM_IPE_CPL_RST_OUTPUT=.true.

I will forward the related contact emails between Weiyu and I to see if they were helpful. You can find my run directories in the first e-mail. Thanks in ahead and maybe we can talk more about this issue when you come to office?

[ZhuxiaoLi66]

run script: /scratch3/NCEPDEV/swpc/save/Zhuxiao.Li/git/WAM-IPE_jdel8_ncout/scripts/compsets//scratch3/NCEPDEV/swpc/save/Zhuxiao.Li/git/WAM-IPE_jdel8_ncout/scripts/compsets output: /scratch3/NCEPDEV/swpc/scrub/Zhuxiao.Li/WAM_output/condT3_swin_20150316_jdel8_compiler_60s_1stp_ncout /scratch3/NCEPDEV/swpc/scrub/Zhuxiao.Li/PTMP/condT3_swin_20150316_jdel8_compiler_60s_1stp_ncout

[akubaryk]

At some point I wanted to pull out w from WAM, and I noticed that this was commented out:

!        IF(NC_output .AND.
!     &    MOD(NINT(deltim) * kdt, FHOUT_NC * 3600) == 0) THEN
!         CALL write_NC(kdt, lonf, latg, ngrids_gg_ipe,
!     &                  buff_final)
!        END IF

and that the "write_NC" subroutine didn't exist, so I just wrote one... it was a quick and dirty thing, so I'd like to clean it up if we actually need this code. With respect to the brief discussion in NOAA-SWPC/WAM-IPE#134, restart files can be in whatever format (though I still don't think binary is optimal for these, either), but there's no reason to write binary files for GrADS since nobody on the team uses that software.

[ZhuxiaoLi66]

Agree the GrADS thing.

[ZhuxiaoLi66]

there is a wrtout_netcdf.f. need to double check.

[ZhuxiaoLi66]

while it is nice to know that you write another one and going to refine it.

[ZhuxiaoLi66]

@akubaryk could you point out your write_NC subroutine? thanks.

[akubaryk]

It would be great if, as we discussed yesterday, we could formulate a clear plan for what kind of output we need directly from the model relative to what can (or should) be post-processed before I push a branch with a feature that we may or may not want to support. Any work that goes into the GSM dycore (which this work is) will have to be rewritten mostly from scratch in FV3. We may want to loop in Valery on this discussion since he has a vision for the diagnostics that we need directly from FV3.

[ZhuxiaoLi66]

Great, Adam. I am going to tell you that I will talk with Tim again about this (based on our former communication, his idea is that this is an important step forward, just the efforts timing was a tricky issue). I think the timing issue will be based on our future discussion, the point is that how much efforts we need to input for now. I will organize a discussion later about this.

[ZhuxiaoLi66]

First of all, it doesn't make sense that many important variables are already calculated inside WAM, we dump out many binary things and repeat the calculation outside. especially, the frame work is already there, we need not that much effort on this point. we also should to check WAM standard output, don't think all the output are very useful. This feature of WAM is really out of date.

[ZhuxiaoLi66]

what is your idea on this? Sure, I will loop in Valery. you are right, his idea and package are something we need for the next step.

[akubaryk]

The sigio output format is codified and used by GSI component of WDAS. Because of that, and the WAM history files taking up relatively little disk space, I'm not inclined to move to a different format unless there is a high pressing need for it. If there's additional output that we need from the model that's not in the sigio format (and it's still not clear to me what that is), that should be the focus of the discussion with Tim and Valery.

[ZhuxiaoLi66]

As for the better_ouput variables, please go to NEMS/src/atmos/gsm/dyn/get_w_z.f, and check the variables in module 'get_variables_for_WAM_IPE_coupling'. they are w(wwg), z(zzg), rqg(q,oz,clw,O1,O2,number density), n2g(NO2 number density), den(total number density), gmol (mean mass), most of them are very important for modeling performance analysis, especially for WAM-IPE coupling. also most of the variables are already in grid_collect_ipe.f and wrtout_netcdf.f (even not be used directly for now), kind of be prepared for output.
on the other hand, we do need to double check the current WAM output, to discuss what are the necessary ones, what are not.

[ZhuxiaoLi66]

As for the format options, let us discuss with Valery and Tim. It is an real issue that the sigio output format is used by GSI component of WDAS. by the way, what does 'codified' mean here? as for the 'relatively little disk space', would like to ask how much?

[akubaryk]

I'm not understanding how sigio is a problem. It's not worth dumping the effort into rectifying if we have other options (model output, yes, but also post-processing) given the GSM's limited lifespan. As I said, it's fine if we need additional output in a different format, but there's no reason to abandon sigio for the limited benefit it would bring.

A WAM sigma file is roughly 19MB:

$ du -sh /scratch1/NCEPDEV/swpc/WAM-IPE_DATA/ICS/WAM-IPE/2013031600/siganl.2013031600
19M     /scratch1/NCEPDEV/swpc/WAM-IPE_DATA/ICS/WAM-IPE/2013031600/siganl.2013031600

[ZhuxiaoLi66]

Your idea about the sigio sounds reasonable, given the GSM's limited lifespan... Maybe for now, we should focus on the additional output in the NetCDF format. Also we need Valery's input and package to think more about the next step and what we could do for now...

[akubaryk]

@ZhuxiaoLi66 please take a look at the following NetCDF file, created from the latest updates to the netcdf_output WAM branch: /scratch1/NCEPDEV/stmp4/Adam.Kubaryk/ncout_test/wam_fields.nc

It's worth noting that it adds about ~8.5 seconds to the model run per write (3 hour run with hourly writes went from 1m22s to 1m47s), and give or take 170MB per write (compared to 19MB for a history file). Since these are not intended as model history files and are only being used for diagnostic purposes, we could use NetCDF's built-in compression options to cut down the file size, but having never used this option before, I don't know what kind of potential disk savings or additional processing time there might be. We also might consider simply writing fewer history files when this option is turned on, though that's not going to make much of a difference.

[ZhuxiaoLi66]

@akubaryk got the file and your information. going to take a look at the NetCDF file. Thanks and talk later.

[ZhuxiaoLi66]

@akubaryk thanks for adding the Netcdf output option to WAM. I checked the variables in the wam_fields.nc, they looks pretty reasonable (I proposed the units of the variables, would like to confirm with you here). For T, U, V,

For number density of O, O2, N2, I supposed that the units are (1/m3) here, then we can explain the magnitude difference between WAM and MSIS (not the same date comparison, not necessary have same pattern, while we can see the pattern are quite similar.)

[ZhuxiaoLi66]

As for the extra time and disk space needed by the NetCDF writing, my thought is that since the NetCDF output is only an option (not the routine output for operational) for WAM run, maybe we don't need to worry too much about the disk resource. The sig output need to be transformed into NetCDF format for all the diagnostic runs anyway. As for the compression method on NetCDF, we can try it, I don't know how the coupling (WAM-IPE) run handle this for its NetCDF output?

[ZhuxiaoLi66]

Adam, excuse me to ask the following question (do hope it is not a heart-attack one). Do you think it is facile to add the possibility to do the every min or 3 mins NetCDF output for WAM? just like you did for Joule Heating integral output? we don't always do the very timestep output, while sometimes for the storm runs especially for the satellite validation, we need it. what is your thought about this? Thanks very much for the already very great and hard work!!!

[akubaryk]

Here are the new units:

    float w(time, levs, lat, lon) ;
        w:units = "m/s" ;
    float z(time, levs, lat, lon) ;
        z:units = "m" ;
    float u(time, levs, lat, lon) ;
        u:units = "m/s" ;
    float v(time, levs, lat, lon) ;
        v:units = "m/s" ;
    float t(time, levs, lat, lon) ;
        t:units = "K" ;
    float q1(time, levs, lat, lon) ;
        q1:units = "kg/kg" ;
    float q2(time, levs, lat, lon) ;
        q2:units = "kg/kg" ;
    float q3(time, levs, lat, lon) ;
        q3:units = "kg/kg" ;
    float q4(time, levs, lat, lon) ;
        q4:units = "m^-3" ;
    float q5(time, levs, lat, lon) ;
        q5:units = "m^-3" ;
    float n2(time, levs, lat, lon) ;
        n2:units = "m^-3" ;
    float dn(time, levs, lat, lon) ;
        dn:units = "kg*m^-3" ;
    float gm(time, levs, lat, lon) ;
        gm:units = "kg/mol" ;

For compression, it turns out that approach is extremely effective. Not only were the three hourly outputs compressed to 2.2MB (~200x compression), but the wallclock dropped back to the 1m25s range; seems most of the time spent was on I/O, which isn't so surprising in retrospect. You can check the outputs here: /scratch1/NCEPDEV/stmp4/Adam.Kubaryk/compress_ncout_test/wam_fields.nc -- but I think you'll find it's very similar to the last set (with the units added). It makes me wonder if we can't utilize this in the IPE repo for smaller diagnostic files with faster writes as well -- although it should be noted that you can't use compression with parallel writes (another argument for asynchronous I/O in IPE).

It is possible to output at every timestep -- should be as simple as setting fhout_nc=0. With the above information, it shouldn't really be a big deal to run for a relatively long time with every-timestep output; up to an extra minute per model hour on wallclock, 45MB/hour of data storage or roughly 1GB/day is extremely tolerable.

I'll wait for your comment on the above units and then commit up these changes.

[ZhuxiaoLi66]

Great to know!! have to do something else now, will double check the units after an hour, seems fine to me though!

[ZhuxiaoLi66]

All the units are correct!

[ZhuxiaoLi66]

While when I apply ncview or NCL to display the compressed wam_fields.nc directly, it seems that only the surface level has valid values for ncview, doesn't have any for NCL. I bet we need to do the decompression before use it?

[akubaryk]

Well that makes a lot more sense for the file size and speed, oops.

[akubaryk]

So with the full field being written, ~75MB/write and with the added time on the compression, it's more like ~10s/cycle (in one single test). So it's 110GB/day and ~4 hours/day wallclock with every-timestep output, not quite as appealing anymore.

It's true that this time comes out of the post-processing, but this diagnostic work is also being done with 39 idle cores rather than in parallel. Kind of indifferent on the best path forward.

[akubaryk]

One option would be an additional namelist parameter by which the user can select which fields they want, e.g.

#         w, z, u, v, t,q1,q2,q3,q4,q5,n2,den,gmol
nc_fields=F, F, T, T, T, F, F, F, T, F, T, F, F

for U, V, T, O, N2 only? I'm sure there's a better way to do this, but this would work...

[ZhuxiaoLi66]

I think it is a good idea to use the 'nc_fields'. can we also output every 3 mins or 6 mins instead of every 1 min (time step length)?

[ZhuxiaoLi66]

so there is no makeup way for the compression method?

[akubaryk]

We can leave the compression method in, just a matter of trading compute time for disk space. Can even make that a namelist parameter, too. NetCDF4 (or HDF5 beneath it) supports compression levels 1-9 -- 0 can simply not enable it?

For the output, I can change FHOUT_NC to something like DELOUT_NC, where that would be the same unit as DELTIM (seconds) -- e.g. 60 for every timestep, 180 for every three, 3600 for every hour, etc.

[ZhuxiaoLi66]

As long as we can get the valid NetCDF files, something like https://www.unidata.ucar.edu/blogs/developer/entry/netcdf_compression ? not sure it is the same method? seems the file size don't reduced dramatically (~200 x) after the compression.

The 'DELOUT_NC', looks totally fine for me.

[ZhuxiaoLi66]

Adam, could you please change the variable names q1, q2, q3, q4, q5 into qr, Oz, clw, O, O2 if not very complicated? thanks.

[akubaryk]

Sure thing.

[akubaryk]

WAM branch netcdf_output has been updated and GSMWAM-IPE branch wam_netcdf_output also updated.

[ZhuxiaoLi66]

Got it, thanks for the great job!

[akubaryk]

Pulls are open for this. wam_fields.nc is now linked into the ROTDIR with the cycle appended for restart-enabled output.