From b53740a7bd1cc416f634589075b8c8b89f0ef761 Mon Sep 17 00:00:00 2001 From: Xu Lu Date: Tue, 26 Mar 2024 21:46:08 -0400 Subject: [PATCH] Fix HAFS GSI debug build and run issues (#679) **DUE DATE for merger of this PR into `develop` is 2/19/2024 (six weeks after PR creation).** **DUE DATE for this PR is extended to 3/19/2024 because @XuLu-NOAA is on leave.** **Description** Xu Lu (xu.lu@noaa.gov) and Biju Thomas (biju.thomas@noaa.gov) fixed bugs regarding HAFS GSI debug build and run issues. This is in corresponding to issue #661 Fixes #661 1. In read_radar.f90, uninitialized toff is making all the ground-based radar observations be placed at -3h instead of 0h, which creates wrong increments for FGAT and 4DEnVar. 2. In read_radar.f90, uninitialized zsges will crash the debug mode. 3. In read_radar.f90, t4dvo should be used instead of t4dv in the read_radar_l2rw_novadqc subroutine. 4. In radinfo.90, maxscan should be increased to at least 252 to allow more scans, or it will crash the debug mode. 5. In read_fl_hdob.f90, dlnpsob is replaced with 1000. since the SFMR does not sample surface pressure, and the uninitialized dlnpsob creates issues later in setupspd.f90 in the debug mode. 6. In mod_fv3_lola.f90, (i,j+1) should be used instead of (i+1,j) in searching for V edges. 7. In stpcalc.f90, when tried to find the best stepsize from outpen around L838-864, the minimum outstp(i) is stored in stp(ii), but the istp_use is asigned with i instead of ii. Create inconsistency when assigning stp(istp_use) to stpinout at L872. Should use istp_use=ii instead. **Type of change** - [Yes] Bug fix (non-breaking change which fixes an issue) **How Has This Been Tested?** Regression test on Orion: ``` Test project /work/noaa/hwrf/save/xulu/mergeversions/GSI/build CMake Warning (dev) at CTestTestfile.cmake:9 (subdirs): Syntax Warning in cmake code at /work/noaa/hwrf/save/xulu/mergeversions/GSI/build/regression/CTestTestfile.cmake:7:10 1/7 Test #4: [=[netcdf_fv3_regional]=] ........ Passed 365.11 sec 2/7 Test #7: [=[global_enkf]=] ................ Passed 430.29 sec 3/7 Test #3: [=[rrfs_3denvar_glbens]=] ........ Passed 605.35 sec 4/7 Test #2: [=[rtma]=] ....................... Passed 969.78 sec 5/7 Test #6: [=[hafs_3denvar_hybens]=] ........***Failed 1455.47 sec 6/7 Test #1: [=[global_4denvar]=] ............. Passed 1682.40 sec 7/7 Test #5: [=[hafs_4denvar_glbens]=] ........***Failed 1758.90 sec ``` The failed hafs_3denvar and 4denvar are within expectation due to the fix for toff. As demonstrated in the single observation tests in the following figure, the uninitialized toff can result in increment degradations due to wrongly assigned observation times: ![image](https://github.com/NOAA-EMC/GSI/assets/26603014/0de870e1-f8c8-4b6d-8039-57f417b76367) --- src/gsi/mod_fv3_lola.f90 | 8 +++----- src/gsi/radinfo.f90 | 2 +- src/gsi/read_fl_hdob.f90 | 4 ++-- src/gsi/read_radar.f90 | 24 +++++++++++++----------- src/gsi/stpcalc.f90 | 4 ++-- 5 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/gsi/mod_fv3_lola.f90 b/src/gsi/mod_fv3_lola.f90 index e8df85068e..11bb3b6e37 100644 --- a/src/gsi/mod_fv3_lola.f90 +++ b/src/gsi/mod_fv3_lola.f90 @@ -951,7 +951,6 @@ subroutine definecoef_regular_grids(nxen,nyen,grid_lon,grid_lont,grid_lat,grid_l do i=1,nxen ! center lat/lon of the edge rlat=half*(grid_lat(i,j)+grid_lat(i+1,j)) -! rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) diff=(grid_lon(i,j)-grid_lon(i+1,j))**2 if(diff < sq180)then rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) @@ -979,12 +978,11 @@ subroutine definecoef_regular_grids(nxen,nyen,grid_lon,grid_lont,grid_lat,grid_l do j=1,nyen do i=1,nxen+1 rlat=half*(grid_lat(i,j)+grid_lat(i,j+1)) -! rlon=half*(grid_lon(i,j)+grid_lon(i,j+1)) - diff=(grid_lon(i,j)-grid_lon(i+1,j))**2 + diff=(grid_lon(i,j)-grid_lon(i,j+1))**2 if(diff < sq180)then - rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) + rlon=half*(grid_lon(i,j)+grid_lon(i,j+1)) else - rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)-360._r_kind) + rlon=half*(grid_lon(i,j)+grid_lon(i,j+1)-360._r_kind) endif xr=cos(rlat*deg2rad)*cos(rlon*deg2rad) yr=cos(rlat*deg2rad)*sin(rlon*deg2rad) diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90 index ede58b9bca..4ad17626e6 100644 --- a/src/gsi/radinfo.f90 +++ b/src/gsi/radinfo.f90 @@ -897,7 +897,7 @@ subroutine radinfo_read ! Allocate arrays to receive angle dependent bias information. ! Open file to bias file (satang=satbias_angle). Read data. - maxscan=250 + maxscan=252 if (.not.adp_anglebc) maxscan = 90 ! default value for old files if (adp_anglebc) then diff --git a/src/gsi/read_fl_hdob.f90 b/src/gsi/read_fl_hdob.f90 index 1ef3d8617f..4041740d52 100644 --- a/src/gsi/read_fl_hdob.f90 +++ b/src/gsi/read_fl_hdob.f90 @@ -48,7 +48,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si use kinds, only: r_single,r_kind,r_double,i_kind use constants, only: zero,one_tenth,one,two,ten,deg2rad,t0c,half,& three,four,rad2deg,tiny_r_kind,huge_r_kind,r0_01,& - r60inv,r10,r100,r2000,hvap,eps,omeps,rv,grav + r60inv,r10,r100,r2000,hvap,eps,omeps,rv,grav,r_missing use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,& tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& rlats,rlons,twodvar_regional,fv3_regional @@ -1133,7 +1133,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si cdata_all( 1,iout)=woe ! wind error cdata_all( 2,iout)=dlon ! grid relative longitude cdata_all( 3,iout)=dlat ! grid relative latitude - cdata_all( 4,iout)=dlnpsob ! ln(surface pressure in cb) + cdata_all( 4,iout)=r_missing ! ln(surface pressure in cb) !Since dlnpsob is not provided by SFMR, force it to be r_missing. Not used in setupspd.f90 cdata_all( 5,iout)=spdob*sqrt(two)*half ! u obs cdata_all( 6,iout)=spdob*sqrt(two)*half ! v obs cdata_all( 7,iout)=rstation_id ! station id diff --git a/src/gsi/read_radar.f90 b/src/gsi/read_radar.f90 index 5b1cffbf0c..a824bbbe4e 100644 --- a/src/gsi/read_radar.f90 +++ b/src/gsi/read_radar.f90 @@ -3282,7 +3282,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) use convinfo, only: nconvtype,ncmiter,ncgroup,ncnumgrp,icuse,ioctype,pmot_conv use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model use mpimod, only: npe - use obsmod, only: reduce_diag + use obsmod, only: reduce_diag,time_offset implicit none @@ -3323,7 +3323,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) integer(i_kind) iret,kx0 integer(i_kind) nreal,nchanl,ilat,ilon,ikx integer(i_kind) idomsfc - real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff + real(r_kind) usage,ff10,sfcr,skint,t4dvo real(r_kind) eradkm,dlat_earth,dlon_earth real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist @@ -3467,7 +3467,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) staheight=this_stahgt !station elevation tiltangle=corrected_tilt*deg2rad - t4dvo=toff+thistime + t4dvo=thistime+time_offset timemax=max(timemax,t4dvo) timemin=min(timemin,t4dvo) @@ -3586,7 +3586,8 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) end if if(usage >= 100._r_kind)rusage(ndata)=.true. - call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) + call deter_sfc2(dlat_earth,dlon_earth,t4dvo,idomsfc,skint,ff10,sfcr) + call deter_zsfc_model(dlat,dlon,zsges) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude @@ -3594,7 +3595,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) - cdata(7) = t4dv ! obs time (hour) + cdata(7) = t4dvo ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) @@ -3699,6 +3700,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) use obsmod, only: doradaroneob,oneobradid,time_offset,reduce_diag use mpeu_util, only: gettablesize,gettable use convinfo, only: nconvtype,icuse,ioctype + use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model use mpimod, only: npe use read_l2bufr_mod, only: radar_sites,radar_rmesh,radar_zmesh,elev_angle_max,del_time,range_max,radar_pmot use constants, only: eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,grav_equator @@ -3744,7 +3746,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) integer(i_kind) iret,kx0 integer(i_kind) nreal,nchanl,ilat,ilon,ikx integer(i_kind) idomsfc - real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff + real(r_kind) usage,ff10,sfcr,skint,t4dvo real(r_kind) eradkm,dlat_earth,dlon_earth real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist @@ -4071,7 +4073,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) slat=sin(dlat_earth) staheight=this_stahgt !station elevation tiltangle=corrected_tilt*deg2rad - t4dvo=toff+thistime + t4dvo=time_offset+thistime timemax=max(timemax,t4dvo) timemin=min(timemin,t4dvo) ! Exclude data if it does not fall within time window @@ -4166,7 +4168,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) if (l4dvar) then timedif = zero else - timedif=abs(t4dvo-toff) + timedif=abs(t4dvo-time_offset) endif crit1 = timedif/r6+half call map3grids_m(1,save_all,zflag,zl_thin,nlevz, & @@ -4233,8 +4235,8 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) end if ! Get information from surface file necessary for conventional data here -! call deter_zsfc_model(dlat,dlon,zsges) -! call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) + call deter_zsfc_model(dlat,dlon,zsges) + call deter_sfc2(dlat_earth,dlon_earth,t4dvo,idomsfc,skint,ff10,sfcr) nsuper2_kept=nsuper2_kept+1 cdata(1) = error ! wind obs error (m/s) @@ -4243,7 +4245,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) - cdata(7) = t4dvo+time_offset ! obs time (hour) + cdata(7) = t4dvo ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90 index c66bb58291..34030763db 100644 --- a/src/gsi/stpcalc.f90 +++ b/src/gsi/stpcalc.f90 @@ -844,10 +844,10 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & if(outpen(i) < outpensave)then stp(ii)=outstp(i) outpensave=outpen(i) - istp_use=i + istp_use=ii end if end do - if(istp_use /= istp_iter) then + if(istp_use /= nsteptot) then final_ii=ii exit stepsize end if