changes for reanalysis runs #591
Conversation
|
Since I couldn't edit Jeff's PR fork, I made my own fork and installed the changes mentioned by reviewers of the first PR. I probably should've had my own fork for this to begin with because I had made all the initial code changes, and tested them as well. So I think this PR should supersede the first one from this point on. In any case I would request final reviews for this PR from @BrettHoover-NOAA , @HaixiaLiu-NOAA , @ilianagenkova , @jswhit , and @RussTreadon-NOAA. Thanks all. |
RussTreadon-NOAA
requested review from
BrettHoover-NOAA,
ilianagenkova,
HaixiaLiu-NOAA,
jswhit2 and
RussTreadon-NOAA
src/gsi/read_satwnd.f90
Outdated
| obsdat(4) > 100000000.0_r_kind) cycle loop_readsb | ||
| if(ppb >r10000) ppb=ppb/r100 | ||
| if (ppb>rmiss .or. hdrdat(3)>rmiss .or. obsdat(4)>rmiss) cycle loop_readsb | ||
| ppb=ppb/r100 |
There was a problem hiding this comment.
A pressure value range check was removed here: if(ppb >r10000)
I think it's a good idea to keep it, in case a data provider changes the pressure units.
There was a problem hiding this comment.
@jack-woollen I agree with Iliana, I think this is the only remaining issue for satwnd and I can approve once this check is reinstated (or if it was included elsewhere and I missed it, let me know)
There was a problem hiding this comment.
@ilianagenkova @BrettHoover-NOAA How about if I add a comment such as
if(ppb >r10000) ppb=ppb/r100 ! ppb<10000 may indicate data reported in hPa
Although it seems like a realistic cutoff might be ppb>r1000
There was a problem hiding this comment.
I think this would work, the original limit is probably better for keeping the check completely unambiguous about realistic Pa/hPa values
There was a problem hiding this comment.
@BrettHoover-NOAA Okay! I'll leave it at that. Though I have to say, haven't seen much 10K hPa data! Let's say the data can be Pa/daPa/hpa.
There was a problem hiding this comment.
@jack-woollen , adding that comment is very helpful, especially for users who don't have experience with retrievals and what goes on data provider's end. Thanks!
src/gsi/setupoz.f90
Outdated
| @@ -538,11 +554,12 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& | |||
| if (ozone_diagsave .and. luse(i)) then | |||
| rdiagbuf(1,k,ii) = ozobs(k) | |||
| rdiagbuf(2,k,ii) = ozone_inv(k) ! obs-ges | |||
| errorinv = sqrt(varinv4diag(k)*rat_err4diag) | |||
| errorinv = sqrt(varinv4diag(k)*rat_err2) | |||
There was a problem hiding this comment.
The variable rat_err2 should be reverted back to rat_err4diag. A few lines before this, rat_err2 was set to zero for any monitored data. However, in the diagnostic data files, we want to save real error inverse values rather than zeros. Therefore, the new variable rat_err4diag was introduced to correctly store the values in the diagnostic data file.
There was a problem hiding this comment.
@HaixiaLiu-NOAA I'm looking again at subroutine setupozlev where the above change is.
Grepping for rat_err4diag (and adding definitions of rat_err2) gets:
real(r_kind),dimension(nlevs):: ratio_errors,error,rat_err4diag
real(r_kind) omg,rat_err2
rat_err2 = ratio_errors(k)**2
rat_err4diag=rat_err2
errorinv = sqrt(varinv4diag(k)*rat_err4diag(k))
So really, rat_err4diag(k)=ratio_errors(k)**2 , is that right? If so I can simplify.
There was a problem hiding this comment.
both ratio_errors and rat_err2 are set to 0s in the following section of the code:
! If not assimilating this observation, reset inverse variance to zero
if (iouse(k)<1) then
varinv3(k)=zero
ratio_errors(k)=zero
rat_err2 = zero
end if
Then a few lines down, we have the following:
! Optionally save data for diagnostics
if (ozone_diagsave .and. luse(i)) then
rdiagbuf(1,k,ii) = ozobs(k)
rdiagbuf(2,k,ii) = ozone_inv(k) ! obs-ges
errorinv = sqrt(varinv4diag(k)*rat_err4diag)
There was a problem hiding this comment.
@HaixiaLiu-NOAA Got it. I see what your saying. Thx!
src/gsi/m_extOzone.F90
Outdated
| @@ -1021,7 +1039,7 @@ subroutine ozlev_ncread_(dfile,dtype,ozout,nmrecs,ndata,nodata, gstime,twind) | |||
| if (ozone(ilev, iprof) < -900.0_r_kind) cycle ! undefined | |||
| if (err(ilev, iprof) < -900.0_r_kind) cycle ! undefined | |||
| if (iuse_oz(ipos(ilev)) < 0) then | |||
| usage = 100._r_kind | |||
| usage = 10000._r_kind | |||
There was a problem hiding this comment.
what is the reason for changing the usage from 100 to 10000? This section is in the subroutine ozlev_ncread_. There is a similar subroutine ozlev_bufrread_ where the usage is still set to 100. I am unsure the impact of changing the usage from 100 to 10000, but it is important to ensure consistency for the same variables in both subroutines.
There was a problem hiding this comment.
@HaixiaLiu-NOAA Thanks again for checking. The first one is my blooper, I made half the required change. The second one is from NASA I guess. There's a couple other places there where usage is set to 10000 but I agree the one you pointed out should probably be 100. Pushing thos changes now.
|
@BrettHoover-NOAA I agree with you to merge the PR with minimal changes to read_satwnd.f90, to get things moving. I've tested a more complete refactoring to address some of the issues we discussed and got identical gsistat output from observer runs. That version is in cactus:/lfs/h2/emc/global/noscrub/Jack.Woollen/gfsda.gsidev/src/gsi. If you want to go through it when you have time and put it another PR later on that's good with me. Cheers. |
|
Thanks @jack-woollen, I'm fine with @RussTreadon-NOAA approving these changes. |
|
Hera ctests NOAA-EMC/GSI The Comparison of the initial Jo terms between the update (
Comparison of the initial wind fit-to-obs statistics from the update and control show differences for subtype 257 of prepbufr report types 245, 246, and 247. Below is the difference in the assimilated observation counts for these obs: update control The update code is allowing more observations to pass quality control and be assimilated. Why? FYI, the above values were taken from
The This PR can not be merged into |
|
Tagging @jack-woollen , @BrettHoover-NOAA , and @ilianagenkova for awareness. For some reason the update code processes more satwnd obs for report types 245, 246, and 247 from subtype 0257. We can't merge PR #591 into |
src/gsi/read_satwnd.f90
Outdated
| @@ -1609,6 +1666,11 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis | |||
| deallocate(lmsg,tab,nrep) | |||
| ! Close unit to bufr file | |||
| call closbf(lunin) | |||
|
|
|||
|
|
|||
| do i=1,1000 | |||
There was a problem hiding this comment.
Is this do block with a print a left over section of debug code? Can we remove it?
|
Yes, I'll get rid of.
…On Mon, Aug 7, 2023 at 10:40 AM RussTreadon-NOAA ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/gsi/read_satwnd.f90
<#591 (comment)>:
> @@ -1609,6 +1666,11 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
deallocate(lmsg,tab,nrep)
! Close unit to bufr file
call closbf(lunin)
+
+
+do i=1,1000
Is this do block with a print a left over section of debug code? Can we
remove it?
—
Reply to this email directly, view it on GitHub
<#591 (review)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AO3XO6NT43D6ZPPWVOEFVLLXUD44ZANCNFSM6AAAAAA2R6APBY>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
src/gsi/read_satwnd.f90
Outdated
| do_qc = subset(1:7)=='NC00503'.and.nint(hdrdat(1))>=270 | ||
| do_qc = do_qc.or.subset=='NC005081'.or.subset=='NC005091' | ||
| do_qc = do_qc.or.qcdat(1,1)<rmiss | ||
| if(.not.do_qc) goto 99 |
There was a problem hiding this comment.
WCOSS implementation standards state that we should not use goto statements. NCO has a long-standing bugzilla (#216) requiring that GFS applications, including, GSI, reduce the number of goto statements in every implementation. Can we handle the .not.do_qc condition without using a goto?
There was a problem hiding this comment.
@RussTreadon-NOAA I could add a do_qc check in every extra block .or. make a gigantic if(do_qc) block.
There was a problem hiding this comment.
Which approach yields more readable & logical code?
There was a problem hiding this comment.
I like adding do_qc as a condition in each extra block.
There was a problem hiding this comment.
Works for me, but sorry for the extra work for you.
There was a problem hiding this comment.
I found the bug. It has to do with not doing QC on NC00501*, where goes13 can be found. It'll take a couple days to get the changes worked out. Thanks for finding the problem!
There was a problem hiding this comment.
Right, the NC00501* parts of the code is the culprit.
|
The difference between "update" and "control" - o-g 01 uv asm 245 0257 count 0 46 131 0 0 120 38 75 64 2 0 476 o-g 01 uv asm 245 0257 count 0 29 91 0 0 59 13 17 24 2 0 235 appears to be coming from an older GOES instrument (SAID=257 is GOES-13). May be a SAID value (or range) check in the code needs to be revisited to include these earlier GOES instruments. |
|
@jack-woollen , appreciate your quickly identifying the problem. If you want me to test changes, just let me know. It's easy to drop in updated code, recompile, and rerun. |
|
@RussTreadon-NOAA @ilianagenkova Here's the bulk of the new bits. I think its ready to test again. Thanks! |
|
@jack-woollen , your changes work for I'll rerun all 9 ctests later today. Thank you for the quick fix! |
|
Rerun Hera ctests Update working copy of The global var failures are true failures. Examination of The counts from the control code are The update is assimilating more satwnd observations. The above are taken from
|
|
@RussTreadon-NOAA @ilianagenkova Here's the bulk of the new bits. |
|
Hera ctest using Recompile The A check of the update and control The |
src/gsi/m_extOzone.F90
Outdated
| nodata = 0 | ||
| read_loop1: do iprof = 1, nprofs | ||
| do ilev = 1, levs | ||
| if (ozone(ilev, iprof) .lt. -900.0) cycle ! undefined |
There was a problem hiding this comment.
Add _r_kind suffix to -900.0 (GSI coding standard).
src/gsi/m_extOzone.F90
Outdated
| ozout(4,ndata)=dlat ! grid relative latitude | ||
| ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) | ||
| ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) | ||
| ozout(7,ndata)=0. ! total ozone error flag |
There was a problem hiding this comment.
Add _r_kind suffix to 0.0 or use zero from constants module.
(GSI coding standard).
src/gsi/m_extOzone.F90
Outdated
| ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) | ||
| ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) | ||
| ozout(7,ndata)=0. ! total ozone error flag | ||
| ozout(8,ndata)=0. ! profile ozone error flag |
There was a problem hiding this comment.
Add _r_kind suffix to 0.0 or use zero from constants module.
(GSI coding standard).
src/gsi/read_ssmi.f90
Outdated
|
|
||
|
|
||
| !--- Extract brightness temperature data. Apply gross check to data. | ||
| ! If obs fails gross check, reset to missing obs value. | ||
|
|
||
| call ufbrep(lnbufr,mirad,1,nchanl*nscan,iret,str2) | ||
|
|
||
| !print*,mirad |
There was a problem hiding this comment.
Since this print and abort are commented out, I guess we don't need them. Is there any reason to keep them in the code as inactive lines?
There was a problem hiding this comment.
Thank you @jack-woollen for making these changes and also catching ones I missed. Recompile and rerun Hera ctests. All tests pass except netcdf_fv3_regional. This failure was due to the scalability test. Nothing anomalous found in gsi.x wall times for netcdf_fv3_regional test.
Approve.
|
@jack-woollen , thanks for handling the satwnd code with grace! I wished I was of more help. |
|
Did anyone try to cycle in the workflow with this PR? In the process of updating the hash for global-workflow, I'm running into an issue that might be related. See NOAA-EMC/global-workflow#1835 for more information. |
) **Description** PR #591 removed jacobian information from the netcdf ozone diagnostic file. This caused `enkf.x` to crash. This PR adds the removed ozone jacobian arrays back to the netcdf ozone diagnostic file. Fixes #618 **Type of change** - [x] Bug fix (non-breaking change which fixes an issue) **How Has This Been Tested?** The revised code was tested in the 20210814 18 gdas cycle of a C192L127 enkf parallel. The updated `gsi.x` created an oznstat file which was successfully processed by `enkf.x`. **Checklist** - [x] My code follows the style guidelines of this project - [x] I have performed a self-review of my own code - [x] New and existing tests pass with my changes
Description
Fixes #565
NOTE: close PR #564 when this PR is merged into
developType of change
Please delete options that are not relevant.
How Has This Been Tested?
Checklist
DUE DATE for this PR is 8/31/2023. If this PR is not merged into
developby this date, the PR will be closed and returned to the developer.