From d84a18ffb7702c75e73434929ee189a3cc7132ed Mon Sep 17 00:00:00 2001 From: JingCheng-NOAA <135154465+JingCheng-NOAA@users.noreply.github.com> Date: Wed, 20 Sep 2023 11:04:19 -0400 Subject: [PATCH 1/7] Bring in HAFSv1 Related Changes (#608) **Description** Bring in HAFSv1 related maxobs changes and the capability of assimilating GOES-R high-resolution AMVs. Resolved #599 **Type of change** Use "maxobs" as a condition to check whether the number of observations exceeds the limit, to avoid the out of bound/dimension issue in read_anowbufr.f90 read_dbz_nc.f90 read_gmi.f90 read_goesglm.f90 read_radar.f90 read_radar_wind_ascii.f90 This update also added the capability of assimilating the CIMSS enhanced GOES-R AMVs in a new "satwhr" bufr file. Please delete options that are not relevant. - [x] Bug fix (non-breaking change which fixes an issue) - [x] New feature (non-breaking change which adds functionality) **How Has This Been Tested?** This updates passed the HAFS related regression tests. All tests are performed on Orion. **Checklist** - [x] My code follows the style guidelines of this project - [x] I have performed a self-review of my own code - [x] I have commented my code, particularly in hard-to-understand areas - [x] New and existing tests pass with my changes - [x] Any dependent changes have been merged and published **DUE DATE for this PR is 9/21/2023.** If this PR is not merged into develop by this date, the PR will be closed and returned to the developer. --------- Co-authored-by: jswhit Co-authored-by: Ting.Lei@noaa.gov Co-authored-by: Jili Dong Co-authored-by: Bin Liu Co-authored-by: daprediction Co-authored-by: TingLei-NOAA Co-authored-by: TingLei-NOAA <63461756+TingLei-NOAA@users.noreply.github.com> Co-authored-by: xulu Co-authored-by: Li Bi Co-authored-by: edward.safford Co-authored-by: MichaelLueken-NOAA <63728921+MichaelLueken-NOAA@users.noreply.github.com> Co-authored-by: RussTreadon-NOAA Co-authored-by: Michael Lueken Co-authored-by: AndrewEichmann-NOAA Co-authored-by: Rahul Mahajan Co-authored-by: Emily Liu Co-authored-by: emilyhcliu <36091766+emilyhcliu@users.noreply.github.com> Co-authored-by: Andrew Collard <40322596+ADCollard@users.noreply.github.com> Co-authored-by: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> --- src/gsi/read_anowbufr.f90 | 1 + src/gsi/read_dbz_nc.f90 | 1 + src/gsi/read_gmi.f90 | 1 + src/gsi/read_goesglm.f90 | 1 + src/gsi/read_obs.F90 | 7 ++--- src/gsi/read_radar.f90 | 2 ++ src/gsi/read_radar_wind_ascii.f90 | 1 + src/gsi/read_satwnd.f90 | 45 ++++++++++++++++++++++++------- src/gsi/setupuwnd10m.f90 | 2 +- src/gsi/setupvwnd10m.f90 | 2 +- src/gsi/setupw.f90 | 4 +-- src/gsi/setupwspd10m.f90 | 2 +- 12 files changed, 51 insertions(+), 18 deletions(-) diff --git a/src/gsi/read_anowbufr.f90 b/src/gsi/read_anowbufr.f90 index e2b744eb6a..1873d0b877 100644 --- a/src/gsi/read_anowbufr.f90 +++ b/src/gsi/read_anowbufr.f90 @@ -307,6 +307,7 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& ndata=ndata+1 nodata=nodata+1 + if(ndata>maxobs) exit cdata_all(iconc,ndata) = conc ! pm2_5 obs cdata_all(ierror,ndata) = obserror ! pm2_5 obs error diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index cddbd14de4..f6ac9aa112 100644 --- a/src/gsi/read_dbz_nc.f90 +++ b/src/gsi/read_dbz_nc.f90 @@ -417,6 +417,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index 0f45aa7e28..f59529662a 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -528,6 +528,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& flgch = 0 iobs=iobs+1 + if(iobs>maxobs) exit end do read_loop end do read_subset 690 continue diff --git a/src/gsi/read_goesglm.f90 b/src/gsi/read_goesglm.f90 index e0124abbf2..bf8639c72d 100644 --- a/src/gsi/read_goesglm.f90 +++ b/src/gsi/read_goesglm.f90 @@ -276,6 +276,7 @@ subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis) icntpnt=icntpnt+1 ndata=ndata+1 + if(ndata>maxobs) exit nodata=nodata+1 iout=ndata isort(icntpnt)=iout diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 8e451c75be..cb4a7c4b8f 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -436,10 +436,10 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) end do nread = nread + 1 end do airploop - else if(trim(filename) == 'satwndbufr')then + else if(index(filename,'satwnd') /=0 .or. index(filename,'satwhr') /=0) then lexist = .false. loop: do while(ireadmg(lnbufr,subset,idate2) >= 0) -! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034 and NC005039) +! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034, NC005039, NC005099) ! are added as the GOES-R bufr file provide do not contain other winds. ! May not be necessary with the operational satwnd BUFR if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or.& @@ -450,6 +450,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or.& trim(subset) == 'NC005032' .or. trim(subset) == 'NC005034' .or.& trim(subset) == 'NC005039' .or. & + trim(subset) == 'NC005099' .or. & trim(subset) == 'NC005090' .or. trim(subset) == 'NC005091' .or.& trim(subset) == 'NC005067' .or. trim(subset) == 'NC005068' .or. trim(subset) == 'NC005069' .or.& trim(subset) == 'NC005047' .or. trim(subset) == 'NC005048' .or. trim(subset) == 'NC005049' .or.& @@ -1503,7 +1504,7 @@ subroutine read_obs(ndata,mype) else if(obstype == 'uv' .or. obstype == 'wspd10m' .or. & obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then ! Process satellite winds which seperate from prepbufr - if ( index(infile,'satwnd') /=0 ) then + if ( index(infile,'satwnd') /=0 .or. index(infile,'satwhr') /=0 ) then call read_satwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs_sub1(1,i)) string='READ_SATWND' diff --git a/src/gsi/read_radar.f90 b/src/gsi/read_radar.f90 index 40c77a7ee2..9ce156e736 100644 --- a/src/gsi/read_radar.f90 +++ b/src/gsi/read_radar.f90 @@ -2911,6 +2911,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then if(zflag == 0)then @@ -4031,6 +4032,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) end if !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit ithin=1 !number of obs to keep per grid box if(radar_no_thinning) then ithin=-1 diff --git a/src/gsi/read_radar_wind_ascii.f90 b/src/gsi/read_radar_wind_ascii.f90 index 2e1b06a50c..604d9d0eca 100644 --- a/src/gsi/read_radar_wind_ascii.f90 +++ b/src/gsi/read_radar_wind_ascii.f90 @@ -509,6 +509,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg !#################### Data thinning ################### icntpnt=icntpnt+1 + if(icntpnt>maxobs) exit if(ithin > 0)then if(zflag == 0)then diff --git a/src/gsi/read_satwnd.f90 b/src/gsi/read_satwnd.f90 index 1679708787..943cf4d47b 100644 --- a/src/gsi/read_satwnd.f90 +++ b/src/gsi/read_satwnd.f90 @@ -18,6 +18,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 253: EUMETSAT IR winds, 254: EUMETSAT WV deep layer winds ! 257,258,259: MODIS IR,WV cloud top, WV deep layer winds ! 260: VIIR IR winds +! 241: CIMSS enhanced AMV winds ! respectively ! For satellite subtype: 50-80 from EUMETSAT geostationary satellites(METEOSAT) ! 100-199 from JMA geostationary satellites(MTSAT) @@ -77,6 +78,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 2021-07-25 Genkova - added code for Metop-B/C winds in new BUFR,NC005081 ! ! 2022-01-20 Genkova - added missing station_id for polar winds ! 2022-01-20 Genkova - added code for Meteosat and Himawari AMVs in new BUFR +! 2022-12-10 Bi - added code for CIMSS enhanced AMVs in new BUFR ! ! ! input argument list: @@ -155,7 +157,6 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis real(r_kind),parameter:: r799=799.0_r_kind real(r_kind),parameter:: r1200= 1200.0_r_kind real(r_kind),parameter:: r10000= 10000.0_r_kind - real(r_double),parameter:: rmiss=10d7 ! Declare local variables @@ -212,7 +213,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis real(r_double),dimension(13):: hdrdat real(r_double),dimension(4):: obsdat - real(r_double),dimension(2) :: hdrdat_test + real(r_double),dimension(2) :: hdrdat_test,hdrdat_005099 real(r_double),dimension(3,5) :: heightdat real(r_double),dimension(6,4) :: derdwdat real(r_double),dimension(3,12) :: qcdat @@ -509,7 +510,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis !GOES-R section of the 'if' statement over 'subsets' else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039') then + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then ! Commented out, because we need clarification for SWCM/hdrdat(9) from Yi Song ! NOTE: Once it is confirmed that SWCM values are sensible, apply this logic and replace lines 685-702 ! if(hdrdat(9) == one) then @@ -537,6 +538,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis itype=246 else if(trim(subset) == 'NC005031') then ! WV clear sky/deep layer itype=247 + else if(trim(subset) == 'NC005099') then + itype=241 endif else ! wind is not recognised and itype is not assigned cycle loop_report @@ -735,7 +738,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis do_qc = subset(1:7)=='NC00503'.and.nint(hdrdat(1))>=270 do_qc = do_qc.or.subset(1:7)=='NC00501' do_qc = do_qc.or.subset=='NC005081'.or.subset=='NC005091' - do_qc = do_qc.or.qcret>0 + do_qc = do_qc.or.qcret>0 ! assign types and get quality info: start @@ -1051,7 +1054,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif ! get quality information THIS SECTION NEEDS TO BE TESTED!!! call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}') - irep_array = int(rep_array) + irep_array = max(1,int(rep_array)) allocate( amvivr(2,irep_array)) call ufbrep(lunin,amvivr,2,irep_array,iret, 'TCOV CVWD') pct1 = amvivr(2,1) ! use of pct1 (a new variable in the BUFR) is introduced by Nebuda/Genkova @@ -1175,9 +1178,12 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif ! Extra block for VIIRS NOAA20: End ! Extra block for GOES-R winds: Start - else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & !IR(LW) / CS WV / VIS GOES-R like winds - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' ) then !CT WV / IR(SW) GOES-R like winds + else if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & !IR(LW) / CS WV / VIS GOES-R like winds + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then !CT WV / IR(SW) GOES-R like winds + if ( trim(subset) == 'NC005099' ) then + hdrdat(10)=61.23 ! set zenith angle for CIMSS AMVs to 67 to pass QC, no value in origional data + end if if(hdrdat(1) >=r250 .and. hdrdat(1) <=r299 ) then ! the range of NESDIS satellite IDs ! The sample newBUFR has SAID=259 (GOES-15) ! When GOES-R SAID is assigned, pls check @@ -1209,6 +1215,10 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis c_station_id='WV'//stationid c_sprvstg='WV' !write(6,*)'itype= ',itype + else if(trim(subset) == 'NC005099') then ! WV clear sky/deep layer + itype=241 + c_station_id='IR'//stationid + c_sprvstg='IR' endif ! call ufbint(lunin,rep_array,1,1,iret, '{AMVAHA}') @@ -1223,6 +1233,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! call ufbrep(lunin,amviii,12,irep_array,iret, 'LTDS SCLF SAID SIID CHNM SCCF ORBN SAZA BEARAZ EHAM PRLC TMDBST') ! deallocate( amviii ) + if (itype /= 241) then + call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}') irep_array = int(rep_array) allocate( amvivr(2,irep_array)) @@ -1253,7 +1265,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(wrf_nmm_regional) then ! type 251 has been determine not suitable to be subjected to pct1 range check - if(itype==240 .or. itype==245 .or. itype==246) then + if(itype==240 .or. itype==245 .or. itype==246 .or. itype==241) then if (pct1 < 0.04_r_kind) qm=15 if (pct1 > 0.50_r_kind) qm=15 elseif (itype==251) then @@ -1279,6 +1291,16 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if (isflg == 1 .and. ppb > 850.0_r_kind) qm=15 ! low over land endif + else ! Assign values for the mnemonics/variables missing in original datafile for type 241 + + call ufbint(lunin,hdrdat_005099,2,1,iret, 'GNAPS PCCF'); + qifn=hdrdat_005099(2); + qm=2.0 ! do not reject the wind + pct1=0.4 ! do not reject the wind + ee=1.0 ! do not reject the wind + + endif + ! winds rejected by qc dont get used if (qm == 15) usage=r100 if (qm == 3 .or. qm ==7) woe=woe*r1_2 @@ -1288,9 +1310,12 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(itype==246 ) then; c_prvstg='GOESR' ; c_sprvstg='WVCT' ; endif if(itype==247 ) then; c_prvstg='GOESR' ; c_sprvstg='WVCS' ; endif if(itype==251 ) then; c_prvstg='GOESR' ; c_sprvstg='VIS' ; endif + if(itype==241 ) then; c_prvstg='GOESR' ; c_sprvstg='IR' ; endif !to be revisited I.Genkova + endif ! Extra block for GOES-R winds: End else ! wind is not recognised and itype is not assigned + write(6,*) 'read_satwnd: WIND IS NOT RECOGNIZEd and we are in hell' cycle loop_readsb endif @@ -1338,7 +1363,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! 3 snow ! 4 mixed if( .not. twodvar_regional) then - if(itype ==245 .or. itype ==252 .or. itype ==253 .or. itype ==240) then + if(itype ==245 .or. itype ==252 .or. itype ==253 .or. itype ==240 .or. itype ==241) then if(hdrdat(2) >20.0_r_kind) then call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg) if(isflg /= 0) cycle loop_readsb @@ -1465,7 +1490,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! GOES-R wind are identified/recognised here by subset, but it could be done by itype or SAID ! After completing the evaluation of GOES-R winds, REVISE this section!!! if(trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or. trim(subset) == 'NC005032' .or. & - trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' ) then + trim(subset) == 'NC005034' .or. trim(subset) == 'NC005039' .or. trim(subset) == 'NC005099') then obserr=obserr/two endif diff --git a/src/gsi/setupuwnd10m.f90 b/src/gsi/setupuwnd10m.f90 index 4552a7e81a..24a4e3d4f7 100644 --- a/src/gsi/setupuwnd10m.f90 +++ b/src/gsi/setupuwnd10m.f90 @@ -428,7 +428,7 @@ subroutine setupuwnd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) diff --git a/src/gsi/setupvwnd10m.f90 b/src/gsi/setupvwnd10m.f90 index 0c601e716b..0f5b46900a 100644 --- a/src/gsi/setupvwnd10m.f90 +++ b/src/gsi/setupvwnd10m.f90 @@ -428,7 +428,7 @@ subroutine setupvwnd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index 087c3c34ab..784df1dfbe 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -877,7 +877,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) @@ -1146,7 +1146,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(itype ==244) then ! AVHRR, use same as MODIS qcgross=r0_7*cgross(ikx) endif - if( itype == 245 .or. itype ==246) then + if( itype == 245 .or. itype ==246 .or. itype ==241) then if(presw <400.0_r_kind .and. presw >300.0_r_kind ) qcgross=r0_7*cgross(ikx) endif if(itype == 253 .or. itype ==254) then diff --git a/src/gsi/setupwspd10m.f90 b/src/gsi/setupwspd10m.f90 index 22618fbf9e..ad50c5b0c1 100644 --- a/src/gsi/setupwspd10m.f90 +++ b/src/gsi/setupwspd10m.f90 @@ -635,7 +635,7 @@ subroutine setupwspd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d lowlevelsat=itype==242.or.itype==243.or.itype==245.or.itype==246.or. & itype==247.or.itype==250.or.itype==251.or.itype==252.or. & itype==253.or.itype==254.or.itype==257.or.itype==258.or. & - itype==259 + itype==259.or.itype==241 if (lowlevelsat .and. twodvar_regional) then call windfactor(presw,factw) data(iuob,i)=factw*data(iuob,i) From aa656c6ee23d9a8d253be3dad95615b29ac0d033 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Fri, 22 Sep 2023 08:05:11 -0600 Subject: [PATCH 2/7] fix for zero total ozone pressure in EnKF (issue #625) (#626) The fix involves setting the pressure to 0.001Pa for ozone obs that have zero pressure (to avoid Inf when log(p) calculated), and turning off ob-space vertical localization. Has no effect on current operational setup which uses model-space vertical localization (modelspace_vloc=T). --------- Co-authored-by: jswhit2 --- src/enkf/enkf_obsmod.f90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/enkf/enkf_obsmod.f90 b/src/enkf/enkf_obsmod.f90 index ba4b2946b1..ea8f6446fb 100644 --- a/src/enkf/enkf_obsmod.f90 +++ b/src/enkf/enkf_obsmod.f90 @@ -262,7 +262,6 @@ subroutine readobs() allocate(corrlengthsq(nobstot),lnsigl(nobstot),obtimel(nobstot)) lnsigl=1.e10 do nob=1,nobstot - oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units if (obloclon(nob) < zero) obloclon(nob) = obloclon(nob) + 360._r_single radlon=deg2rad*obloclon(nob) radlat=deg2rad*obloclat(nob) @@ -283,6 +282,13 @@ subroutine readobs() lnsigl(nob)=latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh) end if endif + ! total column ozone has pressure set to zero, set to 0.001Pa + ! and turn vertical localization off (no effect if modelspace_vloc=T) + if (obpress(nob) < 0.001 .and. obtype(nob)(1:3) .eq. ' oz') then + lnsigl(nob) = 1.e30 ! turn ob-space vert localization off + obpress(nob) = 0.001 ! set to a non-zero value + endif + oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units corrlengthsq(nob)=latval(deglat,corrlengthnh,corrlengthtr,corrlengthsh)**2 if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then corrlengthsq(nob)=latval(deglat,corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh)**2 From 2f4e7fe8124603c6ac8c14ed2cf1aa393d6ec07d Mon Sep 17 00:00:00 2001 From: jderber-NOAA <75998838+jderber-NOAA@users.noreply.github.com> Date: Fri, 22 Sep 2023 12:17:22 -0400 Subject: [PATCH 3/7] Optimize the reading of ensembles and setup for global multiscale runs (#594) This update improves the efficiency of the GSI, especially for multiscale runs. Details can be found in issue #585 --- src/gsi/apply_scaledepwgts.f90 | 152 ++++++--- src/gsi/control_vectors.f90 | 8 +- src/gsi/cplr_gfs_ensmod.f90 | 263 ++++++++-------- src/gsi/general_specmod.f90 | 31 -- src/gsi/general_spectral_transforms.f90 | 4 +- src/gsi/general_sub2grid_mod.f90 | 167 ++++++++++ src/gsi/get_gefs_ensperts_dualres.f90 | 402 +++++++++++++----------- src/gsi/hybrid_ensemble_isotropic.F90 | 265 ++++++++-------- src/gsi/read_iasi.f90 | 6 +- src/gsi/read_prepbufr.f90 | 6 +- src/gsi/setupcldtot.F90 | 19 +- src/gsi/setuprad.f90 | 10 +- src/gsi/stpcalc.f90 | 3 +- 13 files changed, 790 insertions(+), 546 deletions(-) diff --git a/src/gsi/apply_scaledepwgts.f90 b/src/gsi/apply_scaledepwgts.f90 index e97b6fb614..e4952b28fa 100644 --- a/src/gsi/apply_scaledepwgts.f90 +++ b/src/gsi/apply_scaledepwgts.f90 @@ -16,18 +16,18 @@ function fwgtofwvlen (rvlft,rvrgt,rcons,rlen,rinput) ! !$$$ end documentation block - use kinds, only: r_kind,i_kind,r_single + use kinds, only: r_kind implicit none real(r_kind),intent(in) :: rvlft,rvrgt,rcons,rlen,rinput real(r_kind) :: fwgtofwvlen real(r_kind) :: rlen1,rtem1,rconshalf - rlen1=rlen/10.0_r_kind ! rlen corresponds to a (-5,5) region - rconshalf=0.5_r_kind*rcons if(rinput > rvlft .and. rinput < rvrgt) then fwgtofwvlen=rcons else + rlen1=rlen/10.0_r_kind ! rlen corresponds to a (-5,5) region + rconshalf=0.5_r_kind*rcons rtem1=min(abs(rinput-rvlft),abs(rinput-rvrgt)) fwgtofwvlen=rconshalf*(1.0_r_kind+tanh(5.0_r_kind-rtem1/rlen1)) endif @@ -41,23 +41,21 @@ subroutine init_mult_spc_wgts(jcap_in) ! !$$$ end documentation block - use kinds, only: r_kind,i_kind,r_single - use constants, only: zero,half,one,two,three,rearth,pi,tiny_r_kind + use kinds, only: r_kind,i_kind + use constants, only: zero,half,one,rearth,pi,tiny_r_kind use mpimod, only: mype - use general_sub2grid_mod, only: general_sub2grid_create_info - use egrid2agrid_mod,only: g_create_egrid2agrid - use general_sub2grid_mod, only: sub2grid_info use hybrid_ensemble_parameters, only: nsclgrp use hybrid_ensemble_parameters, only: spc_multwgt,spcwgt_params,r_ensloccov4scl implicit none integer(i_kind),intent(in ) :: jcap_in - integer(i_kind) i + integer(i_kind) i,l,ks integer(i_kind) ig - real(r_kind) rwv0,rtem1,rtem2 - real (r_kind):: fwgtofwvlen + real(r_kind) :: rwv0,rtem1,rtem2 + real(r_kind) :: fwgtofwvlen real(r_kind) :: totwvlength + real(r_kind),dimension(0:jcap_in,nsclgrp) :: spcwgt logical :: l_sum_spc_weights ! Spectral scale decomposition is differernt between SDL-cross and SDL-nocross @@ -67,9 +65,9 @@ subroutine init_mult_spc_wgts(jcap_in) l_sum_spc_weights = .true. end if - spc_multwgt(0,1)=one + spcwgt(0,1)=one do ig=2,nsclgrp - spc_multwgt(0,ig)=zero + spcwgt(0,ig)=zero end do @@ -79,13 +77,13 @@ subroutine init_mult_spc_wgts(jcap_in) rtem1=zero do ig=1,nsclgrp if(ig /= 2) then - spc_multwgt(i,ig)=fwgtofwvlen(spcwgt_params(1,ig),spcwgt_params(2,ig),& - spcwgt_params(3,ig),spcwgt_params(4,ig),totwvlength) - spc_multwgt(i,ig)=min(max(spc_multwgt(i,ig),zero),one) + spcwgt(i,ig)=fwgtofwvlen(spcwgt_params(1,ig),spcwgt_params(2,ig),& + spcwgt_params(3,ig),spcwgt_params(4,ig),totwvlength) + spcwgt(i,ig)=min(max(spcwgt(i,ig),zero),one) if(l_sum_spc_weights) then - rtem1=rtem1+spc_multwgt(i,ig) + rtem1=rtem1+spcwgt(i,ig) else - rtem1=rtem1+spc_multwgt(i,ig)*spc_multwgt(i,ig) + rtem1=rtem1+spcwgt(i,ig)*spcwgt(i,ig) endif endif enddo @@ -93,20 +91,52 @@ subroutine init_mult_spc_wgts(jcap_in) if(rtem2 >= zero) then if(l_sum_spc_weights) then - spc_multwgt(i,2)=rtem2 + spcwgt(i,2)=rtem2 else - spc_multwgt(i,2)=sqrt(rtem2) + spcwgt(i,2)=sqrt(rtem2) endif else - if(mype == 0)write(6,*) ' rtem2 < zero ',i,rtem2,(spc_multwgt(i,ig),ig=1,nsclgrp) - spc_multwgt(i,2)=zero + if(mype == 0)write(6,*) ' rtem2 < zero ',i,rtem2,(spcwgt(i,ig),ig=1,nsclgrp) + spcwgt(i,2)=zero endif enddo +!! Code borrowed from spvar in splib + + spc_multwgt = zero + do ig=1,nsclgrp + do i=0,jcap_in + ks=2*i + spc_multwgt(ks+1,ig)=spcwgt(i,ig) + end do + do i=0,jcap_in + do l=MAX(1,i-jcap_in),MIN(i,jcap_in) + ks=l*(2*jcap_in+1-l)+2*i + spc_multwgt(ks+1,ig) = spcwgt(i,ig) + spc_multwgt(ks+2,ig) = spcwgt(i,ig) + end do + end do + end do + return end subroutine init_mult_spc_wgts +subroutine destroy_mult_spc_wgts +!$$$ subprogram documentation block +! +! subprogram: destroy_mult_spc_wgts +! +!$$$ end documentation block + + use hybrid_ensemble_parameters, only: spc_multwgt,spcwgt_params + implicit none + + deallocate(spc_multwgt,spcwgt_params) -subroutine apply_scaledepwgts(grd_in,sp_in,wbundle,spwgts,wbundle2) + return +end subroutine destroy_mult_spc_wgts + + +subroutine apply_scaledepwgts(m,grd_in,sp_in) ! ! Program history log: ! 2017-03-30 J. Kay, X. Wang - copied from Kleist's apply_scaledepwgts and @@ -115,49 +145,67 @@ subroutine apply_scaledepwgts(grd_in,sp_in,wbundle,spwgts,wbundle2) ! use constants, only: one use control_vectors, only: control_vector - use kinds, only: r_kind,i_kind - use kinds, only: r_single - use general_specmod, only: general_spec_multwgt + use kinds, only: r_kind,i_kind,r_single use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: general_sub2grid,general_grid2sub use general_specmod, only: spec_vars use general_sub2grid_mod, only: sub2grid_info - use mpimod, only: mpi_comm_world,mype + use hybrid_ensemble_parameters, only: spc_multwgt,en_perts,nsclgrp,n_ens + use mpimod, only: mype implicit none ! Declare passed variables - type(gsi_bundle),intent(in) :: wbundle - type(gsi_bundle),intent(inout) :: wbundle2 + integer,intent(in) :: m type(spec_vars),intent (in):: sp_in type(sub2grid_info),intent(in)::grd_in - real(r_kind),dimension(0:sp_in%jcap),intent(in):: spwgts ! Declare local variables - integer(i_kind) kk + integer(i_kind) kk,ig,n,ig2,i,j - real(r_kind),dimension(grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc) :: hwork - real(r_kind),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work - real(r_kind),dimension(sp_in%nc):: spc1 - -! Beta1 first -! Get from subdomains to - call general_sub2grid(grd_in,wbundle%values,hwork) - work=reshape(hwork,(/grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc/)) - - do kk=1,grd_in%nlevs_alloc -! Transform from physical space to spectral space - call general_g2s0(grd_in,sp_in,spc1,work(:,:,kk)) - -! Apply spectral weights - call general_spec_multwgt(sp_in,spc1,spwgts) -! Transform back to physical space - call general_s2g0(grd_in,sp_in,spc1,work(:,:,kk)) + real(r_single),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc,nsclgrp) :: hwork2 + real(r_kind),dimension(grd_in%nlat,grd_in%nlon) :: work + real(r_kind),dimension(sp_in%nc,grd_in%nlevs_alloc):: spc1 + real(r_kind),dimension(sp_in%nc):: spc2 + + do n=1,n_ens +! Get from subdomains to full grid + call general_sub2grid(grd_in,en_perts(n,1,m)%valuesr4(:),hwork2(:,:,:,1)) + +!$omp parallel do schedule(static,1) private(i,j,kk,work) + do kk=1,grd_in%nlevs_loc + do j=1,grd_in%nlon + do i=1,grd_in%nlat + work(i,j)=hwork2(i,j,kk,1) + end do + end do +! Transform from physical space to spectral space + call general_g2s0(grd_in,sp_in,spc1(1,kk),work) + + end do +!$omp parallel do schedule(static,1) private(kk,ig,ig2,i,j,work,spc2) + do ig2=1,nsclgrp*grd_in%nlevs_loc + ig=(ig2-1)/grd_in%nlevs_loc+1 + kk=ig2-(ig-1)*grd_in%nlevs_loc + + do i=1,sp_in%nc + spc2(i)=spc1(i,kk)*spc_multwgt(i,ig) + end do +! Apply spectral weights +! Transform back to physical space + call general_s2g0(grd_in,sp_in,spc2,work) + + do j=1,grd_in%nlon + do i=1,grd_in%nlat + hwork2(i,j,kk,ig)=work(i,j) + end do + end do + end do + do ig=1,nsclgrp +! Transfer work back to subdomains + call general_grid2sub(grd_in,hwork2(:,:,:,ig),en_perts(n,ig,m)%valuesr4(:)) + end do end do -! Transfer work back to subdomains - hwork=reshape(work,(/grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc/)) - call general_grid2sub(grd_in,hwork,wbundle2%values) - return end subroutine apply_scaledepwgts diff --git a/src/gsi/control_vectors.f90 b/src/gsi/control_vectors.f90 index 97578124d2..73f605b95f 100644 --- a/src/gsi/control_vectors.f90 +++ b/src/gsi/control_vectors.f90 @@ -897,21 +897,23 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) itot=max(m3d,0)+max(m2d,0) if(l_hyb_ens)itot=itot+n_ens*naensgrp allocate(partsum(itot)) + partsum=zero_quad do ii=1,nsubwin !$omp parallel do schedule(dynamic,1) private(i) do i = 1,m3d - partsum(i) = dplevs(xcv%step(ii)%r3(i)%q,ycv%step(ii)%r3(i)%q,ihalo=1) + partsum(i) = partsum(i)+dplevs(xcv%step(ii)%r3(i)%q,ycv%step(ii)%r3(i)%q,ihalo=1) enddo !$omp parallel do schedule(dynamic,1) private(i) do i = 1,m2d - partsum(m3d+i) = dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1) + partsum(m3d+i) = partsum(m3d+i)+dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1) enddo if(l_hyb_ens) then do ig=1,naensgrp nigtmp=n_ens*(ig-1) !$omp parallel do schedule(dynamic,1) private(i) do i = 1,n_ens - partsum(m3d+m2d+nigtmp+i) = dplevs(xcv%aens(ii,ig,i)%r3(1)%q,ycv%aens(ii,ig,i)%r3(1)%q,ihalo=1) + partsum(m3d+m2d+nigtmp+i) = partsum(m3d+m2d+nigtmp+i) + & + dplevs(xcv%aens(ii,ig,i)%r3(1)%q,ycv%aens(ii,ig,i)%r3(1)%q,ihalo=1) end do end do end if diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index 550f85d209..6f7d1184c5 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -29,6 +29,7 @@ module get_gfs_ensmod_mod integer(i_kind) :: kas,kae,kasm,kaem,kaemz,mas,mae,masm,maem,maemz integer(i_kind) :: ibs,ibe,ibsm,ibem,ibemz,jbs,jbe,jbsm,jbem,jbemz integer(i_kind) :: kbs,kbe,kbsm,kbem,kbemz,mbs,mbe,mbsm,mbem,mbemz + integer(i_kind) :: icw,iql,iqi,iqr,iqs,iqg integer(i_kind) :: n2d type(genex_info) :: s_a2b @@ -180,13 +181,13 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & use gsi_4dvar, only: ens_fhrlevs use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only : assignment(=) - use hybrid_ensemble_parameters, only: n_ens,grd_ens + use hybrid_ensemble_parameters, only: n_ens,grd_ens,ntlevs_ens use hybrid_ensemble_parameters, only: ensemble_path - use control_vectors, only: nc2d,nc3d - !use control_vectors, only: cvars2d,cvars3d + use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d use genex_mod, only: genex_create_info,genex,genex_destroy_info use gridmod, only: use_gfs_nemsio use jfunc, only: cnvw_option + use mpeu_util, only: getindex implicit none @@ -206,9 +207,9 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & integer(i_kind) :: io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,i_ens integer(i_kind) :: ip integer(i_kind) :: nlon,nlat,nsig - integer(i_kind),dimension(n_ens) :: io_pe0 + integer(i_kind),dimension(n_ens) :: io_pe0,iretx real(r_single),allocatable,dimension(:,:,:,:) :: en_full,en_loc - real(r_kind),allocatable,dimension(:) :: sloc + real(r_single),allocatable,dimension(:,:,:) :: sloc integer(i_kind),allocatable,dimension(:) :: m_cvars2dw,m_cvars3dw integer(i_kind) :: m_cvars2d(nc2d),m_cvars3d(nc3d) type(sub2grid_info) :: grd3d @@ -234,10 +235,10 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & !!!!!!!!for example, n2d = nc3d*nsig + nc2d n2d=nc3d*grd_ens%nsig+nc2d - ias=1 ; iae=0 ; jas=1 ; jae=0 ; kas=1 ; kae=0 ; mas=1 ; mae=0 + ias=0 ; iae=0 ; jas=0 ; jae=0 ; kas=1 ; kae=0 ; mas=1 ; mae=0 if(mype==io_pe) then - iae=nlat - jae=nlon + iae=nlat+1 + jae=nlon+1 kae=n2d mas=n_io_pe_s ; mae=n_io_pe_em endif @@ -248,8 +249,13 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & ibe=ibs+grd_ens%lat1-1 jbs=grd_ens%jstart(mype+1) jbe=jbs+grd_ens%lon1-1 + jbs=jbs-1 + jbe=jbe+1 + ibs=ibs-1 + ibe=ibe+1 + + ibsm=ibs ; ibem=ibe ; jbsm=jbs ; jbem=jbe - ibsm=ibs-ip ; ibem=ibe+ip ; jbsm=jbs-ip ; jbem=jbe+ip kbs =1 ; kbe =n2d ; mbs =1 ; mbe =n_ens kbsm=kbs ; kbem=kbe ; mbsm=mbs ; mbem=mbe iaemz=max(iasm,iaem) ; jaemz=max(jasm,jaem) @@ -273,8 +279,6 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & m_cvars2dw=-999 m_cvars3dw=-999 - - !! read ensembles if ( mas == mae ) then @@ -305,45 +309,56 @@ subroutine get_user_ens_gfs_fastread_(ntindex,atm_bundle, & allocate(en_full(1,1,1,1)) end if +! scatter to subdomains: + call mpi_allreduce(m_cvars2dw,m_cvars2d,nc2d,mpi_integer4,mpi_max,mpi_comm_world,ierror) call mpi_allreduce(m_cvars3dw,m_cvars3d,nc3d,mpi_integer4,mpi_max,mpi_comm_world,ierror) - deallocate(m_cvars2dw,m_cvars3dw) -! scatter to subdomains: + ! Check hydrometeors in control variables + icw=getindex(cvars3d,'cw') + iql=getindex(cvars3d,'ql') + iqi=getindex(cvars3d,'qi') + iqr=getindex(cvars3d,'qr') + iqs=getindex(cvars3d,'qs') + iqg=getindex(cvars3d,'qg') ! en_loc=zero + allocate(en_loc(ibsm:ibemz,jbsm:jbemz,kbsm:kbemz,mbsm:mbemz)) call genex(s_a2b,en_full,en_loc) deallocate(en_full) + if(ntindex == ntlevs_ens)call genex_destroy_info(s_a2b) -! call genex_destroy_info(s_a2b) ! check on actual routine name - - allocate(sloc(lat2in*lon2in*(nc2d+nc3d*nsig))) call create_grd23d_(grd3d,nc2d+nc3d*grd%nsig) - iret=0 + + allocate(sloc(grd3d%lat2,grd3d%lon2,grd3d%num_fields)) + iretx=0 +!$omp parallel do schedule(dynamic,1) private(n,k,j,i,sloc) do n=1,n_ens - ii=0 - do k=1,nc2d+nc3d*nsig - do j=jbsm,jbem - do i=ibsm,ibem - ii=ii+1 - sloc(ii)=en_loc(i,j,k,n) + do k=1,grd3d%num_fields + do j=1,grd3d%lon2 + do i=1,grd3d%lat2 + sloc(i,j,k)=en_loc(i+ibsm-1,j+jbsm-1,k,n) enddo enddo enddo - call move2bundle_(grd3d,sloc,atm_bundle(n),m_cvars2d,m_cvars3d,iret) + call move2bundle_(grd3d,sloc,atm_bundle(n),m_cvars2d,m_cvars3d,iretx(n)) enddo + iret=iretx(1) + do n=2,n_ens + iret=iret+iretx(n) + end do deallocate(en_loc,sloc) call general_sub2grid_destroy_info(grd3d,grd) end subroutine get_user_ens_gfs_fastread_ -subroutine move2bundle_(grd3d,sloc,atm_bundle,m_cvars2d,m_cvars3d,iret) +subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret) !$$$ subprogram documentation block ! . . . . @@ -373,58 +388,35 @@ subroutine move2bundle_(grd3d,sloc,atm_bundle,m_cvars2d,m_cvars3d,iret) ! !$$$ - use constants, only: zero,one,two,fv use general_sub2grid_mod, only: sub2grid_info use hybrid_ensemble_parameters, only: en_perts use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d - use mpeu_util, only: getindex implicit none ! Declare passed variables type(sub2grid_info), intent(in ) :: grd3d type(gsi_bundle), intent(inout) :: atm_bundle - real(r_kind), intent(inout) :: sloc(grd3d%lat2*grd3d%lon2*(nc2d+nc3d*grd3d%nsig)) + real(r_single), intent(inout) :: en_loc3(grd3d%lat2,grd3d%lon2,grd3d%num_fields) integer(i_kind), intent(in ) :: m_cvars2d(nc2d),m_cvars3d(nc3d) integer(i_kind), intent(inout) :: iret ! Declare internal variables character(len=*),parameter :: myname_='move2bundle_' - character(len=70) :: filename - integer(i_kind) :: ierr + integer(i_kind) :: ierr,i,j integer(i_kind) :: km1,m - integer(i_kind) :: icw,iql,iqi,iqr,iqs,iqg - real(r_kind),pointer,dimension(:,:) :: ps + real(r_single),pointer,dimension(:,:) :: ps !real(r_kind),pointer,dimension(:,:) :: sst - real(r_kind),dimension(grd3d%lat2,grd3d%lon2,nc2d+nc3d*grd3d%nsig)::en_loc3 - real(r_kind),pointer,dimension(:,:,:) :: u,v,tv,q,oz,cwmr - real(r_kind),pointer,dimension(:,:,:) :: qlmr,qimr,qrmr,qsmr,qgmr - real(r_kind),parameter :: r0_001 = 0.001_r_kind - - -!--- now update halo values of all variables using general_sub2grid - call update_halos_(grd3d,sloc,en_loc3) - - ! Check hydrometeors in control variables - icw=getindex(cvars3d,'cw') - iql=getindex(cvars3d,'ql') - iqi=getindex(cvars3d,'qi') - iqr=getindex(cvars3d,'qr') - iqs=getindex(cvars3d,'qs') - iqg=getindex(cvars3d,'qg') + real(r_single),pointer,dimension(:,:,:) :: u,v,tv,q,oz,cwmr + real(r_single),pointer,dimension(:,:,:) :: qlmr,qimr,qrmr,qsmr,qgmr ! atm_bundle to zero done earlier - call gsi_bundlegetpointer(atm_bundle,'ps',ps, ierr); iret = iret+ierr + call gsi_bundlegetpointer(atm_bundle,'ps',ps, ierr); iret = ierr !call gsi_bundlegetpointer(atm_bundle,'sst',sst, ierr); iret = iret+ierr - do m=1,nc2d -! convert ps from Pa to cb - if(trim(cvars2d(m))=='ps') ps=r0_001*en_loc3(:,:,m_cvars2d(m)) -! if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now - enddo call gsi_bundlegetpointer(atm_bundle,'sf',u , ierr); iret = ierr + iret call gsi_bundlegetpointer(atm_bundle,'vp',v , ierr); iret = ierr + iret @@ -443,19 +435,25 @@ subroutine move2bundle_(grd3d,sloc,atm_bundle,m_cvars2d,m_cvars3d,iret) write(6,'(A)') trim(myname_) // ': For now, GFS requires all MetFields: ps,u,v,(sf,vp)tv,q,oz' write(6,'(A)') trim(myname_) // ': but some have not been found. Aborting ... ' write(6,'(A)') trim(myname_) // ': WARNING!' - write(6,'(3A,I5)') trim(myname_) // ': Trouble reading ensemble file : ', trim(filename), ', IRET = ', iret endif return endif + do m=1,nc2d +! convert ps from Pa to cb + if(trim(cvars2d(m))=='ps') ps=en_loc3(:,:,m_cvars2d(m)) +! if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now + enddo + km1 = en_perts(1,1,1)%grid%km - 1 -!$omp parallel do schedule(dynamic,1) private(m) do m=1,nc3d if(trim(cvars3d(m))=='sf')then u = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1) else if(trim(cvars3d(m))=='vp') then v = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1) else if(trim(cvars3d(m))=='t') then +! Note tv here is sensible temperature. Converted to virtual temperature +! later. tv = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1) else if(trim(cvars3d(m))=='q') then q = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1) @@ -476,9 +474,6 @@ subroutine move2bundle_(grd3d,sloc,atm_bundle,m_cvars2d,m_cvars3d,iret) end if enddo -! convert t to virtual temperature - tv=tv*(one+fv*q) - return end subroutine move2bundle_ @@ -505,51 +500,6 @@ subroutine create_grd23d_(grd23d,nvert) end subroutine create_grd23d_ -subroutine update_halos_(grd,sloc,s) - - use general_sub2grid_mod, only: sub2grid_info,general_sub2grid,general_grid2sub - - implicit none - - ! Declare passed variables - type(sub2grid_info), intent(in ) :: grd - real(r_kind), intent( out) :: s(grd%lat2,grd%lon2,grd%num_fields) - real(r_kind), intent(inout) :: sloc(grd%lat2*grd%lon2*grd%num_fields) - - ! Declare local variables - integer(i_kind) inner_vars,lat2,lon2,nlat,nlon,nvert,kbegin_loc,kend_alloc - integer(i_kind) ii,i,j,k - real(r_kind),allocatable,dimension(:,:,:,:) :: work - - lat2=grd%lat2 - lon2=grd%lon2 - nlat=grd%nlat - nlon=grd%nlon - nvert=grd%num_fields - inner_vars=grd%inner_vars - kbegin_loc=grd%kbegin_loc - kend_alloc=grd%kend_alloc - - - - allocate(work(inner_vars,nlat,nlon,kbegin_loc:kend_alloc)) - call general_sub2grid(grd,sloc,work) - - call general_grid2sub(grd,work,sloc) - deallocate(work) - - ii=0 - do k=1,nvert - do j=1,lon2 - do i=1,lat2 - ii=ii+1 - s(i,j,k)=sloc(ii) - enddo - enddo - enddo - -end subroutine update_halos_ - subroutine ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,io_pe0,i_ens) ! do computation on all processors, then assign final local processor @@ -605,7 +555,7 @@ subroutine ens_io_partition_(n_ens,io_pe,n_io_pe_s,n_io_pe_e,n_io_pe_em,io_pe0,i end subroutine ens_io_partition_ -subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig, & +subroutine parallel_read_nemsio_state_(en_full,m_cvars2dw,m_cvars3d,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & filename,init_head,filenamesfc) @@ -627,7 +577,7 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi integer(i_kind), intent(in ) :: nlon,nlat,nsig integer(i_kind), intent(in ) :: ias,jas,mas integer(i_kind), intent(in ) :: iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz - integer(i_kind), intent(inout) :: m_cvars2d(nc2d),m_cvars3d(nc3d) + integer(i_kind), intent(inout) :: m_cvars2dw(nc2d),m_cvars3d(nc3d) real(r_single), intent(inout) :: en_full(iasm:iaemz,jasm:jaemz,kasm:kaemz,masm:maemz) character(len=*), intent(in ) :: filename character(len=*), optional, intent(in) :: filenamesfc @@ -816,15 +766,27 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi m_cvars3d(k3)=kf+1 do k=1,nsig kf=kf+1 - jj=jas-1 + jj=jas do j=1,nlon jj=jj+1 - ii=ias-1 + ii=ias do i=1,nlat ii=ii+1 en_full(ii,jj,kf,mas)=temp3(i,j,k,k3) enddo enddo + ii=ias + do i=1,nlat + ii=ii+1 + en_full(ii,jasm,kf,mas)=en_full(ii,jaem-1,kf,mas) + en_full(ii,jaem,kf,mas)=en_full(ii,jasm+1,kf,mas) + enddo + jj=jas-1 + do j=jasm,jaem + jj=jj+1 + en_full(iasm,jj,kf,mas)=en_full(iasm+1,jj,kf,mas) + en_full(iaem,jj,kf,mas)=en_full(iaem-1,jj,kf,mas) + end do enddo enddo deallocate(temp3) @@ -853,24 +815,36 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi ! move temp2 to en_full do k2=1,nc2d - m_cvars2d(k2)=kf+1 + m_cvars2dw(k2)=kf+1 kf=kf+1 - jj=jas-1 + jj=jas do j=1,nlon jj=jj+1 - ii=ias-1 + ii=ias do i=1,nlat ii=ii+1 en_full(ii,jj,kf,mas)=temp2(i,j,k2) enddo enddo + ii=ias + do i=1,nlat + ii=ii+1 + en_full(ii,jasm,kf,mas)=en_full(ii,jaem-1,kf,mas) + en_full(ii,jaem,kf,mas)=en_full(ii,jasm+1,kf,mas) + enddo + jj=jas-1 + do j=jasm,jaem + jj=jj+1 + en_full(iasm,jj,kf,mas)=en_full(iasm+1,jj,kf,mas) + en_full(iaem,jj,kf,mas)=en_full(iaem-1,jj,kf,mas) + end do enddo deallocate(temp2) end subroutine parallel_read_nemsio_state_ -subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig, & +subroutine parallel_read_gfsnc_state_(en_full,m_cvars2dw,m_cvars3d,nlon,nlat,nsig, & ias,jas,mas, & iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz, & filename) @@ -884,7 +858,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig ! !$$$ - use constants, only: r60,r3600,zero,one,half,deg2rad + use constants, only: r60,r3600,zero,one,half,deg2rad,zero_single use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d use general_sub2grid_mod, only: sub2grid_info use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -898,7 +872,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig integer(i_kind), intent(in ) :: nlon,nlat,nsig integer(i_kind), intent(in ) :: ias,jas,mas integer(i_kind), intent(in ) :: iasm,iaemz,jasm,jaemz,kasm,kaemz,masm,maemz - integer(i_kind), intent(inout) :: m_cvars2d(nc2d),m_cvars3d(nc3d) + integer(i_kind), intent(inout) :: m_cvars2dw(nc2d),m_cvars3d(nc3d) real(r_single), intent(inout) :: en_full(iasm:iaemz,jasm:jaemz,kasm:kaemz,masm:maemz) character(len=*), intent(in ) :: filename @@ -1021,15 +995,27 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig m_cvars3d(k3)=kf+1 do k=1,nsig kf=kf+1 - jj=jas-1 + jj=jas do j=1,nlon jj=jj+1 - ii=ias-1 + ii=ias do i=1,nlat ii=ii+1 en_full(ii,jj,kf,mas)=temp3(i,j,k,k3) enddo enddo + ii=ias + do i=1,nlat + ii=ii+1 + en_full(ii,jasm,kf,mas)=en_full(ii,jaem-1,kf,mas) + en_full(ii,jaem,kf,mas)=en_full(ii,jasm+1,kf,mas) + enddo + jj=jas-1 + do j=jasm,jaem + jj=jj+1 + en_full(iasm,jj,kf,mas)=en_full(iasm+1,jj,kf,mas) + en_full(iaem,jj,kf,mas)=en_full(iaem-1,jj,kf,mas) + end do enddo enddo @@ -1042,24 +1028,34 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig call read_vardata(atmges, 'pressfc', rwork2d) call move1_(rwork2d,temp2,nlon,nlat) call fillpoles_ss_(temp2,nlon,nlat) - else - temp2=zero - endif ! move temp2 to en_full - kf=kf+1 - m_cvars2d(k2)=kf - jj=jas-1 - do j=1,nlon - jj=jj+1 - ii=ias-1 + kf=kf+1 + m_cvars2dw(k2)=kf + jj=jas + do j=1,nlon + jj=jj+1 + ii=ias + do i=1,nlat + ii=ii+1 + en_full(ii,jj,kf,mas)=temp2(i,j) + enddo + enddo + ii=ias do i=1,nlat ii=ii+1 - en_full(ii,jj,kf,mas)=temp2(i,j) + en_full(ii,jasm,kf,mas)=en_full(ii,jaem-1,kf,mas) + en_full(ii,jaem,kf,mas)=en_full(ii,jasm+1,kf,mas) enddo - enddo + jj=jas-1 + do j=jasm,jaem + jj=jj+1 + en_full(iasm,jj,kf,mas)=en_full(iasm+1,jj,kf,mas) + en_full(iaem,jj,kf,mas)=en_full(iaem-1,jj,kf,mas) + end do + end if enddo -! call close_dataset(atmges) + call close_dataset(atmges) deallocate(rwork2d) deallocate(temp2) @@ -1158,7 +1154,7 @@ subroutine fillpoles_sv_(tempu,tempv,nlon,nlat,clons,slons) real(r_single), intent(inout) :: tempu(nlat,nlon),tempv(nlat,nlon) real(r_kind), intent(in ) :: clons(nlon),slons(nlon) - integer(i_kind) i + integer(i_kind) i,nlatm real(r_kind) polnu,polnv,polsu,polsv ! Compute mean along southern and northern latitudes @@ -1166,11 +1162,12 @@ subroutine fillpoles_sv_(tempu,tempv,nlon,nlat,clons,slons) polnv=zero polsu=zero polsv=zero + nlatm=nlat-1 do i=1,nlon - polnu=polnu+tempu(nlat-1,i)*clons(i)-tempv(nlat-1,i)*slons(i) - polnv=polnv+tempu(nlat-1,i)*slons(i)+tempv(nlat-1,i)*clons(i) - polsu=polsu+tempu(2,i )*clons(i)+tempv(2,i )*slons(i) - polsv=polsv+tempu(2,i )*slons(i)-tempv(2,i )*clons(i) + polnu=polnu+tempu(nlatm,i)*clons(i)-tempv(nlatm,i)*slons(i) + polnv=polnv+tempu(nlatm,i)*slons(i)+tempv(nlatm,i)*clons(i) + polsu=polsu+tempu(2,i )*clons(i)+tempv(2,i )*slons(i) + polsv=polsv+tempu(2,i )*slons(i)-tempv(2,i )*clons(i) end do polnu=polnu/float(nlon) polnv=polnv/float(nlon) diff --git a/src/gsi/general_specmod.f90 b/src/gsi/general_specmod.f90 index c90187bf70..439e26e431 100644 --- a/src/gsi/general_specmod.f90 +++ b/src/gsi/general_specmod.f90 @@ -64,7 +64,6 @@ module general_specmod ! set subroutines to public public :: general_init_spec_vars public :: general_destroy_spec_vars - public :: general_spec_multwgt ! set passed variables to public public :: spec_vars public :: spec_cut @@ -307,36 +306,6 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) return end subroutine general_init_spec_vars - subroutine general_spec_multwgt(sp,spcwrk,spcwgt) -! 2017-03-30 J. Kay, X. Wang - add general_spec_multwgt for scale-dependent -! weighting of mixed resolution ensemble, -! POC: xuguang.wang@ou.edu -! - implicit none - type(spec_vars),intent(in) :: sp - real(r_kind),dimension(sp%nc),intent(inout) :: spcwrk - real(r_kind),dimension(0:sp%jcap),intent(in) :: spcwgt - - integer(i_kind) l,jmax,ks,n - -!! Code borrowed from spvar in splib - jmax=sp%jcap - - do n=0,jmax - ks=2*n - spcwrk(ks+1)=spcwrk(ks+1)*spcwgt(n) - end do - do n=0,jmax - do l=MAX(1,n-jmax),MIN(n,jmax) - ks=l*(2*jmax+(-1)*(l-1))+2*n - spcwrk(ks+1) = spcwrk(ks+1)*spcwgt(n) - spcwrk(ks+2) = spcwrk(ks+2)*spcwgt(n) - end do - end do - - return - end subroutine general_spec_multwgt - subroutine general_destroy_spec_vars(sp) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/general_spectral_transforms.f90 b/src/gsi/general_spectral_transforms.f90 index d4f0959489..541923e450 100644 --- a/src/gsi/general_spectral_transforms.f90 +++ b/src/gsi/general_spectral_transforms.f90 @@ -368,8 +368,8 @@ subroutine sfilter(grd,sp,filter,grid) call general_sptez_s(sp,spec_work,work,-1) - gnlon=float(grd%nlon) -! gnlon=real(grd%nlon,r_kind) +! gnlon=float(grd%nlon) + gnlon=real(grd%nlon,r_kind) do i=1,sp%nc spec_work(i)=spec_work(i)*gnlon end do diff --git a/src/gsi/general_sub2grid_mod.f90 b/src/gsi/general_sub2grid_mod.f90 index f0548643bb..bacea4688b 100644 --- a/src/gsi/general_sub2grid_mod.f90 +++ b/src/gsi/general_sub2grid_mod.f90 @@ -87,6 +87,7 @@ module general_sub2grid_mod interface general_sub2grid module procedure general_sub2grid_r_single_rank11 module procedure general_sub2grid_r_single_rank14 + module procedure general_sub2grid_r_single_rank13 module procedure general_sub2grid_r_single_rank4 module procedure general_sub2grid_r_double_rank11 module procedure general_sub2grid_r_double_rank14 @@ -97,6 +98,7 @@ module general_sub2grid_mod module procedure general_grid2sub_r_single_rank11 module procedure general_grid2sub_r_single_rank41 module procedure general_grid2sub_r_single_rank4 + module procedure general_grid2sub_r_single_rank31 module procedure general_grid2sub_r_double_rank11 module procedure general_grid2sub_r_double_rank41 module procedure general_grid2sub_r_double_rank4 @@ -1019,6 +1021,93 @@ subroutine general_sub2grid_r_single_rank14(s,sub_vars,grid_vars) end subroutine general_sub2grid_r_single_rank14 + subroutine general_sub2grid_r_single_rank13(s,sub_vars,grid_vars) +!$$$ subprogram documentation block +! . . . . +! subprogram: general_sub2grid_r_single_rank4 convert from subdomains to full horizontal grid +! prgmmr: parrish org: np22 date: 2010-02-11 +! +! abstract: generalized version of sub2grid--uses only gsi module kinds. +! All information needed is contained in the structure variable +! "s", instead of various modules. This allows +! for easy adaptation for any collection/ordering of variables +! defined on subdomains, which need to be made available on +! full horizontal grid for horizontal operations. +! The structure variable is specified by subroutine general_sub2grid_setup. +! This version works with single precision (4-byte) real variables. +! Input sub_vars, the desired arrays on horizontal subdomains, has one +! halo row, for now, which is filled with zero, since for ensemble use, +! there is no need for a halo, but is easiest for now to keep it. +! A later version will have variable number of halo rows, filled with proper values. +! +! program history log: +! 2010-02-11 parrish, initial documentation +! +! input argument list: +! s - structure variable, contains all necessary information for +! moving this set of subdomain variables sub_vars to +! the corresponding set of full horizontal grid variables. +! sub_vars - input grid values in vertical subdomain mode (contains one halo row) +! +! output argument list: +! grid_vars - output grid values in horizontal slab mode. +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use mpimod, only: mpi_comm_world,mpi_real4 + implicit none + + type(sub2grid_info),intent(in ) :: s + real(r_single), intent(in ) :: sub_vars(s%lat2*s%lon2*s%num_fields) + real(r_single), intent( out) :: grid_vars(s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc) + + real(r_single) :: sub_vars_r4(s%lat2,s%lon2,s%num_fields) + real(r_single) :: sub_vars0(s%lat1,s%lon1,s%num_fields) + real(r_single) :: work(s%itotsub*(s%kend_alloc-s%kbegin_loc+1)) + integer(i_kind) iloc,iskip,i,i0,j,j0,k,n,k_in,ilat,jlon,ierror,ioffset + + sub_vars_r4 = reshape(sub_vars,(/s%lat2,s%lon2,s%num_fields/)) +! remove halo row +!$omp parallel do schedule(dynamic,1) private(k,j,j0,i0,i) + do k=1,s%num_fields + do j=2,s%lon2-1 + j0=j-1 + do i=2,s%lat2-1 + i0=i-1 + sub_vars0(i0,j0,k)=sub_vars_r4(i,j,k) + end do + end do + end do + + call mpi_alltoallv(sub_vars0,s%recvcounts,s%rdispls,mpi_real4, & + work,s%sendcounts,s%sdispls,mpi_real4,mpi_comm_world,ierror) + + + k_in=s%kend_loc-s%kbegin_loc+1 + +! Load grid_vars array in desired order +!$omp parallel do schedule(dynamic,1) private(k,iskip,iloc,n,i,ilat,jlon,ioffset) + do k=s%kbegin_loc,s%kend_loc + iskip=0 + iloc=0 + do n=1,s%npe + if (n/=1) then + iskip=iskip+s%ijn(n-1)*k_in + end if + ioffset=iskip+(k-s%kbegin_loc)*s%ijn(n) + do i=1,s%ijn(n) + iloc=iloc+1 + ilat=s%ltosi(iloc) + jlon=s%ltosj(iloc) + grid_vars(ilat,jlon,k)=work(i + ioffset) + end do + end do + end do + + end subroutine general_sub2grid_r_single_rank13 subroutine general_sub2grid_r_single_rank4(s,sub_vars,grid_vars) !$$$ subprogram documentation block ! . . . . @@ -1199,6 +1288,84 @@ subroutine general_grid2sub_r_single_rank41(s,grid_vars,sub_vars) end subroutine general_grid2sub_r_single_rank41 + subroutine general_grid2sub_r_single_rank31(s,grid_vars,sub_vars) +!$$$ subprogram documentation block +! . . . . +! subprogram: general_sub2grid convert from subdomains to full horizontal grid +! prgmmr: parrish org: np22 date: 2010-02-11 +! +! abstract: generalized version of grid2sub--uses only gsi module kinds. +! All information needed is contained in the structure variable +! "s", instead of various modules. This allows +! for easy adaptation for any collection/ordering of variables +! defined on subdomains, which need to be made available on +! full horizontal grid for horizontal operations. +! The structure variable is specified by subroutine general_sub2grid_setup. +! This version works with single precision (4-byte) real variables. +! Output sub_vars, the desired arrays on horizontal subdomains, has one +! halo row, for now, which is filled with zero, since for ensemble use, +! there is no need for a halo, but is easiest for now to keep it. +! A later version will have variable number of halo rows, filled with proper values. +! +! program history log: +! 2010-02-11 parrish, initial documentation +! 2010-03-02 parrish - remove setting halo to zero in output +! 2014-12-03 derber - make similar optimization changes already in code for +! double precision. +! +! input argument list: +! s - structure variable, contains all necessary information for +! moving this set of subdomain variables sub_vars to +! the corresponding set of full horizontal grid variables. +! grid_vars - input grid values in horizontal slab mode. +! +! output argument list: +! sub_vars - output grid values in vertical subdomain mode +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use constants, only: zero + use mpimod, only: mpi_comm_world,mpi_real4 + implicit none + + type(sub2grid_info),intent(in ) :: s + real(r_single), intent(in ) :: grid_vars(s%nlat,s%nlon,s%kbegin_loc:s%kend_alloc) + real(r_single), intent( out) :: sub_vars(s%lat2*s%lon2*s%num_fields) + + real(r_single) :: sub_vars_r4(s%lat2,s%lon2,s%num_fields) + real(r_single) :: temp(s%itotsub*(s%kend_loc-s%kbegin_loc+1)) + integer(i_kind) iloc,i,ii,k,n,ilat,jlon,ierror,icount + integer(i_kind),dimension(s%npe) ::iskip + +! reorganize for eventual distribution to local domains + iskip(1)=0 + do n=2,s%npe + iskip(n)=iskip(n-1)+s%ijn_s(n-1)*(s%kend_loc-s%kbegin_loc+1) + end do +!$omp parallel do schedule(dynamic,1) private(n,k,i,jlon,ii,ilat,iloc,icount) + do k=s%kbegin_loc,s%kend_loc + icount=0 + do n=1,s%npe + iloc=iskip(n)+(k-s%kbegin_loc)*s%ijn_s(n) + do i=1,s%ijn_s(n) + iloc=iloc+1 + icount=icount+1 + ilat=s%ltosi_s(icount) + jlon=s%ltosj_s(icount) + temp(iloc)=grid_vars(ilat,jlon,k) + end do + end do + end do + + + call mpi_alltoallv(temp,s%sendcounts_s,s%sdispls_s,mpi_real4, & + sub_vars_r4,s%recvcounts_s,s%rdispls_s,mpi_real4,mpi_comm_world,ierror) + + sub_vars = reshape(sub_vars_r4,(/s%lat2*s%lon2*s%num_fields/)) + end subroutine general_grid2sub_r_single_rank31 subroutine general_grid2sub_r_single_rank4(s,grid_vars,sub_vars) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90 index e244fa9f53..bb5ee374af 100644 --- a/src/gsi/get_gefs_ensperts_dualres.f90 +++ b/src/gsi/get_gefs_ensperts_dualres.f90 @@ -26,7 +26,7 @@ subroutine get_gefs_ensperts_dualres ! ! get_gefs_ensperts_dualres.f90(182): error #6460: This is not a field name that ! is defined in the encompassing structure. [LAT2] -! call genqsat(qs,tsen,prsl,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative) +! call genqsat2(qs,tsen,prsl,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice) ! 2014-11-30 todling - partially generalized to handle any control vector ! (GFS hook needs further attention) ! - also, take SST from members of ensemble @@ -69,7 +69,7 @@ subroutine get_gefs_ensperts_dualres use gsi_enscouplermod, only: gsi_enscoupler_create_sub2grid_info use gsi_enscouplermod, only: gsi_enscoupler_destroy_sub2grid_info use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info - use hybrid_ensemble_parameters, only: nsclgrp,sp_ens,spc_multwgt,global_spectral_filter_sd + use hybrid_ensemble_parameters, only: nsclgrp,sp_ens,global_spectral_filter_sd implicit none real(r_kind),pointer,dimension(:,:) :: ps @@ -78,24 +78,23 @@ subroutine get_gefs_ensperts_dualres ! real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon):: sst_full,dum real(r_kind),pointer,dimension(:,:,:):: p3 real(r_kind),pointer,dimension(:,:):: x2 - type(gsi_bundle),allocatable,dimension(:) :: en_read + type(gsi_bundle),allocatable,dimension(:) :: en_real8 type(gsi_bundle):: en_bar - type(gsi_bundle) :: en_pertstmp1 - type(gsi_bundle) :: en_pertstmp2 ! type(gsi_grid) :: grid_ens - real(r_kind) bar_norm,sig_norm,kapr,kap1 + real(r_kind) bar_norm,sig_norm ! real(r_kind),allocatable,dimension(:,:):: z,sst2 - real(r_kind),allocatable,dimension(:,:,:) :: tsen,prsl,pri,qs + real(r_kind),allocatable,dimension(:,:,:) :: tsen,prsl ! integer(i_kind),dimension(grd_ens%nlat,grd_ens%nlon):: idum - integer(i_kind) istatus,iret,i,ic3,j,k,n,iderivative,im,jm,km,m,ipic + integer(i_kind) istatus,iret,i,ic3,j,k,n,im,jm,km,m,ipic ! integer(i_kind) mm1 integer(i_kind) ipc3d(nc3d),ipc2d(nc2d) integer(i_kind) ier ! integer(i_kind) il,jl logical ice,hydrometeor type(sub2grid_info) :: grd_tmp - integer(i_kind) :: ig + real(r_kind),parameter :: r0_001 = 0.001_r_kind + ! Create perturbations grid and get variable names from perturbations if(en_perts(1,1,1)%grid%im/=grd_ens%lat2.or. & @@ -135,25 +134,15 @@ subroutine get_gefs_ensperts_dualres if ( istatus /= 0 ) & call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus) -! Allocate bundle used for temporary usage - if( nsclgrp > 1 .and. global_spectral_filter_sd )then - call gsi_bundlecreate(en_pertstmp1,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) - call gsi_bundlecreate(en_pertstmp2,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) - if(istatus/=0) then - write(6,*)' get_gefs_ensperts_dualres: trouble creating en_read like tempbundle' - call stop2(999) - endif - end if - - ! Allocate bundle used for reading members - allocate(en_read(n_ens)) + ! Allocate bundle used for real*8 version of members + allocate(en_real8(n_ens)) do n=1,n_ens - call gsi_bundlecreate(en_read(n),en_perts(1,1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d) + call gsi_bundlecreate(en_real8(n),en_perts(1,1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d) if ( istatus /= 0 ) & - call die('get_gefs_ensperts_dualres',': trouble creating en_read bundle, istatus =',istatus) - en_read(n) = zero + call die('get_gefs_ensperts_dualres',': trouble creating en_real8 bundle, istatus =',istatus) end do + ! allocate(z(im,jm)) ! allocate(sst2(im,jm)) @@ -162,7 +151,7 @@ subroutine get_gefs_ensperts_dualres ntlevs_ens_loop: do m=1,ntlevs_ens - call gsi_enscoupler_get_user_Nens(grd_tmp,n_ens,m,en_read,iret) + call gsi_enscoupler_get_user_Nens(grd_tmp,n_ens,m,en_perts(:,1,m),iret) ! Check read return code. Revert to static B if read error detected if ( iret /= 0 ) then @@ -175,67 +164,52 @@ subroutine get_gefs_ensperts_dualres endif en_bar%values=zero + allocate(tsen(im,jm,km)) if (.not.q_hyb_ens) then !use RH - allocate(pri(im,jm,km+1)) - allocate(prsl(im,jm,km),tsen(im,jm,km)) - allocate(qs(im,jm,km)) + allocate(prsl(im,jm,km)) end if do n=1,n_ens + do i=1,nelen + en_real8(n)%values(i)=real(en_perts(n,1,m)%valuesr4(i),r_kind) + end do - call gsi_bundlegetpointer(en_read(n),'q' ,q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(en_real8(n),'q' ,q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(en_real8(n),'t' ,tv,ier);istatus=istatus+ier + call gsi_bundlegetpointer(en_real8(n),'ps',ps,ier);istatus=ier +! Convert ps to correct units + do j=1,jm + do i=1,im + ps(i,j)=r0_001*ps(i,j) + end do + end do +! Convert to real from single and convert tv to virtual temperature do k=1,km do j=1,jm do i=1,im +! Use following 3 lines for results identical to previous version +! tv(i,j,k)= tv(i,j,k)*(one+fv*q(i,j,k)) +! q(i,j,k)=max(q(i,j,k),zero) +! tsen(i,j,k)=tv(i,j,k)/(one+fv*q(i,j,k)) +! Remove following 3 lines for results identical to previous version q(i,j,k)=max(q(i,j,k),zero) + tsen(i,j,k)=tv(i,j,k) + tv(i,j,k)= tsen(i,j,k)*(one+fv*q(i,j,k)) end do end do end do if (.not.q_hyb_ens) then !use RH - call gsi_bundlegetpointer(en_read(n),'ps',ps,ier);istatus=ier - call gsi_bundlegetpointer(en_read(n),'t' ,tv,ier);istatus=istatus+ier + ! Compute RH ! Get 3d pressure field now on interfaces - call general_getprs_glb(ps,tv,pri) -! Get sensible temperature and 3d layer pressure - if (idsl5 /= 2) then - kap1=rd_over_cp+one - kapr=one/rd_over_cp -!$omp parallel do schedule(dynamic,1) private(k,j,i) - do k=1,km - do j=1,jm - do i=1,im - prsl(i,j,k)=((pri(i,j,k)**kap1-pri(i,j,k+1)**kap1)/& - (kap1*(pri(i,j,k)-pri(i,j,k+1))))**kapr - tsen(i,j,k)= tv(i,j,k)/(one+fv*q(i,j,k)) - end do - end do - end do - else -!$omp parallel do schedule(dynamic,1) private(k,j,i) - do k=1,km - do j=1,jm - do i=1,im - prsl(i,j,k)=(pri(i,j,k)+pri(i,j,k+1))*half - tsen(i,j,k)= tv(i,j,k)/(one+fv*q(i,j,k)) - end do - end do - end do - end if + call general_getprs_glb(ps,tv,prsl) ice=.true. - iderivative=0 - call genqsat(qs,tsen,prsl,im,jm,km,ice,iderivative) - do k=1,km - do j=1,jm - do i=1,im - q(i,j,k)=q(i,j,k)/qs(i,j,k) - end do - end do - end do + call genqsat2(q,tsen,prsl,ice) + end if -!$omp parallel do schedule(dynamic,1) private(i,k,j,ic3,hydrometeor,istatus,p3) +! !$omp parallel do schedule(dynamic,1) private(i,k,j,ic3,hydrometeor,istatus,p3) do ic3=1,nc3d hydrometeor = trim(cvars3d(ic3))=='cw' .or. trim(cvars3d(ic3))=='ql' .or. & @@ -246,7 +220,7 @@ subroutine get_gefs_ensperts_dualres if ( hydrometeor ) then - call gsi_bundlegetpointer(en_read(n),trim(cvars3d(ic3)),p3,istatus) + call gsi_bundlegetpointer(en_real8(n),trim(cvars3d(ic3)),p3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' from read in member ',n,m call stop2(999) @@ -260,7 +234,7 @@ subroutine get_gefs_ensperts_dualres end do else if ( trim(cvars3d(ic3)) == 'oz' .and. oz_univ_static ) then - call gsi_bundlegetpointer(en_read(n),trim(cvars3d(ic3)),p3,istatus) + call gsi_bundlegetpointer(en_real8(n),trim(cvars3d(ic3)),p3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' from read in member ',n,m call stop2(999) @@ -270,7 +244,7 @@ subroutine get_gefs_ensperts_dualres end do !c3d do i=1,nelen - en_bar%values(i)=en_bar%values(i)+en_read(n)%values(i)*bar_norm + en_bar%values(i)=en_bar%values(i)+en_real8(n)%values(i)*bar_norm end do @@ -282,10 +256,9 @@ subroutine get_gefs_ensperts_dualres end do ! end do over ensembles if (.not.q_hyb_ens) then !use RH - deallocate(pri) - deallocate(tsen,prsl) - deallocate(qs) + deallocate(prsl) end if + deallocate(tsen) ! Before converting to perturbations, get ensemble spread !!! it is not clear of the next statement is thread/$omp safe. @@ -309,7 +282,7 @@ subroutine get_gefs_ensperts_dualres !$omp parallel do schedule(dynamic,1) private(n,i,ic3,ipic,k,j) do n=1,n_ens do i=1,nelen - en_perts(n,1,m)%valuesr4(i)=en_read(n)%values(i)-en_bar%values(i) + en_perts(n,1,m)%valuesr4(i)=en_real8(n)%values(i)-en_bar%values(i) end do if(.not. q_hyb_ens) then do ic3=1,nc3d @@ -318,8 +291,8 @@ subroutine get_gefs_ensperts_dualres do k=1,km do j=1,jm do i=1,im - en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),limqens) - en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),-limqens) + en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = & + max(min(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),limqens),-limqens) end do end do end do @@ -330,33 +303,22 @@ subroutine get_gefs_ensperts_dualres en_perts(n,1,m)%valuesr4(i)=en_perts(n,1,m)%valuesr4(i)*sig_norm end do end do + if(nsclgrp > 1 .and. global_spectral_filter_sd) then + call apply_scaledepwgts(m,grd_ens,sp_ens) + end if end do ntlevs_ens_loop !end do over bins - call gsi_bundledestroy(en_bar,istatus) - if(nsclgrp > 1 .and. global_spectral_filter_sd) then - do m=1,ntlevs_ens - do n=1,n_ens - en_pertstmp1%values=en_perts(n,1,m)%valuesr4 - do ig=1,nsclgrp - call apply_scaledepwgts(grd_ens,sp_ens,en_pertstmp1,spc_multwgt(:,ig),en_pertstmp2) - en_perts(n,ig,m)%valuesr4=en_pertstmp2%values - enddo - enddo - enddo - endif - do n=n_ens,1,-1 - call gsi_bundledestroy(en_read(n),istatus) + call gsi_bundledestroy(en_real8(n),istatus) if ( istatus /= 0 ) & - call die('get_gefs_ensperts_dualres',': trouble destroying en_read bundle, istatus = ', istatus) + call die('get_gefs_ensperts_dualres',': trouble destroying en_real8 bundle, istatus = ', istatus) end do - deallocate(en_read) - if(nsclgrp > 1 .and. global_spectral_filter_sd) then - call gsi_bundledestroy(en_pertstmp1,istatus) - call gsi_bundledestroy(en_pertstmp2,istatus) - if ( istatus /= 0 ) & - call die('get_gefs_ensperts_dualres',': trouble destroying en_pertstmp2 bundle, istatus = ', istatus) - end if + deallocate(en_real8) + + call gsi_bundledestroy(en_bar,istatus) + + if(nsclgrp > 1 .and. global_spectral_filter_sd) call destroy_mult_spc_wgts + call gsi_enscoupler_destroy_sub2grid_info(grd_tmp) ! mm1=mype+1 @@ -675,7 +637,7 @@ subroutine write_spread_dualres(ibin,bundle) return end subroutine write_spread_dualres -subroutine general_getprs_glb(ps,tv,prs) +subroutine general_getprs_glb(ps,tv,prsl) ! subprogram: getprs get 3d pressure or 3d pressure deriv ! prgmmr: kleist org: np20 date: 2005-09-29 ! @@ -707,107 +669,197 @@ subroutine general_getprs_glb(ps,tv,prs) use kinds,only: r_kind,i_kind use constants,only: zero,half,one_tenth,rd_over_cp,one - use gridmod,only: nsig,ak5,bk5,ck5,tref5,idvc5 - use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta1_ll,eta2_ll,pdtop_ll,pt_ll,& - regional,wrf_mass_regional,twodvar_regional,fv3_regional + use gridmod,only: nsig,ak5,bk5,ck5,tref5,idvc5,idsl5 use hybrid_ensemble_parameters, only: grd_ens implicit none ! Declare passed variables - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2) ,intent(in ) :: ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig) ,intent(in ) :: tv - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig+1),intent( out) :: prs + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2) ,intent(in ) :: ps + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig),intent(in ) :: tv + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig),intent( out) :: prsl ! Declare local variables - real(r_kind) kapr,trk + real(r_kind) kapr,trk,kap1 + real(r_kind),dimension(grd_ens%lat2,nsig+1) :: prs integer(i_kind) i,j,k,k2 ! ,it -! Declare local parameter - real(r_kind),parameter:: ten = 10.0_r_kind - if (regional) then - if(wrf_nmm_regional.or.nems_nmmb_regional) then - do k=1,nsig+1 - do j=1,grd_ens%lon2 + k2=nsig+1 + kap1=rd_over_cp+one + kapr=one/rd_over_cp +!$omp parallel do schedule(dynamic,1) private(k,j,i,trk,prs) + do j=1,grd_ens%lon2 + do i=1,grd_ens%lat2 + prs(i,1)=ps(i,j) + prs(i,k2)=zero + end do + if (idvc5 /= 3) then + do k=2,nsig do i=1,grd_ens%lat2 - prs(i,j,k)=one_tenth* & - (eta1_ll(k)*pdtop_ll + & - eta2_ll(k)*(ten*ps(i,j)-pdtop_ll-pt_ll) + & - pt_ll) + prs(i,k)=ak5(k)+bk5(k)*ps(i,j) end do end do - end do - elseif (fv3_regional) then - do k=1,nsig+1 - do j=1,grd_ens%lon2 + else + do k=1,nsig do i=1,grd_ens%lat2 - prs(i,j,k)=eta1_ll(k)+ eta2_ll(k)*ps(i,j) + trk=(half*(tv(i,j,k-1)+tv(i,j,k))/tref5(k))**kapr + prs(i,k)=ak5(k)+(bk5(k)*ps(i,j))+(ck5(k)*trk) end do end do - end do - - elseif (twodvar_regional) then - do k=1,nsig+1 - do j=1,grd_ens%lon2 + end if +! Get sensible temperature and 3d layer pressure + if (idsl5 /= 2) then + do k=1,nsig do i=1,grd_ens%lat2 - prs(i,j,k)=one_tenth*(eta1_ll(k)*(ten*ps(i,j)-pt_ll) + pt_ll) + prsl(i,j,k)=((prs(i,k)**kap1-prs(i,k+1)**kap1)/& + (kap1*(prs(i,k)-prs(i,k+1))))**kapr end do end do - end do - elseif (wrf_mass_regional) then - do k=1,nsig+1 - do j=1,grd_ens%lon2 + else + do k=1,nsig do i=1,grd_ens%lat2 - prs(i,j,k)=one_tenth*(eta1_ll(k)*(ten*ps(i,j)-pt_ll) + & - eta2_ll(k) + pt_ll) + prsl(i,j,k)=(prs(i,k)+prs(i,k+1))*half end do end do - end do - endif - else - if (idvc5 /= 3) then -!$omp parallel do schedule(dynamic,1) private(k,j,i) - do k=1,nsig - if(k == 1)then - k2=nsig+1 - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - prs(i,j,k)=ps(i,j) - prs(i,j,k2)=zero - end do - end do - else - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - prs(i,j,k)=ak5(k)+bk5(k)*ps(i,j) - end do - end do + end if + end do + + return +end subroutine general_getprs_glb +subroutine genqsat2(q,tsen,prsl,ice) +!$$$ subprogram documentation block +! . . . . +! subprogram: genqsat +! prgmmr: derber org: np23 date: 1998-01-14 +! +! abstract: obtain saturation specific humidity for given temperature. +! +! program history log: +! 1998-01-14 derber +! 1998-04-05 weiyu yang +! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version +! 1903-10-07 Wei Gu, bug fixes,if qs<0,then set qs=0; merge w/ GSI by R Todling +! 2003-12-23 kleist, use guess pressure, adapt module framework +! 2004-05-13 kleist, documentation +! 2004-06-03 treadon, replace ggrid_g3 array with ges_* arrays +! 2005-02-23 wu, output dlnesdtv +! 2005-11-21 kleist, derber add dmax array to decouple moisture from temp and +! pressure for questionable qsat +! 2006-02-02 treadon - rename prsl as ges_prsl +! 2006-09-18 derber - modify to limit saturated values near top +! 2006-11-22 derber - correct bug: es 2._r_kind) .and. & + tsen(i,j,k) < mint(i))then + lmint(i)=k + mint(i)=tsen(i,j,k) end if end do - else - kapr=one/rd_over_cp -!$omp parallel do schedule(dynamic,1) private(k,j,i,trk) - do k=1,nsig - if(k == 1)then - k2=nsig+1 - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - prs(i,j,k)=ps(i,j) - prs(i,j,k2)=zero - end do - end do + end do + do i=1,grd_ens%lat2 + tdry = mint(i) + tr = ttp/tdry + if (tdry >= ttp .or. .not. ice) then + estmax(i) = psat * (tr**xa) * exp(xb*(one-tr)) + elseif (tdry < tmix) then + estmax(i) = psat * (tr**xai) * exp(xbi*(one-tr)) + else + w = (tdry - tmix) / (ttp - tmix) + estmax(i) = w * psat * (tr**xa) * exp(xb*(one-tr)) & + + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr)) + endif + end do + + do k = 1,nsig + do i = 1,grd_ens%lat2 + tdry = tsen(i,j,k) + tr = ttp/tdry + if (tdry >= ttp .or. .not. ice) then + es = psat * (tr**xa) * exp(xb*(one-tr)) + elseif (tdry < tmix) then + es = psat * (tr**xai) * exp(xbi*(one-tr)) else - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - trk=(half*(tv(i,j,k-1)+tv(i,j,k))/tref5(k))**kapr - prs(i,j,k)=ak5(k)+(bk5(k)*ps(i,j))+(ck5(k)*trk) - end do - end do + esw = psat * (tr**xa) * exp(xb*(one-tr)) + esi = psat * (tr**xai) * exp(xbi*(one-tr)) + w = (tdry - tmix) / (ttp - tmix) + es = w * esw + (one-w) * esi +! es = w * psat * (tr**xa) * exp(xb*(one-tr)) & +! + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr)) + + endif + + pw = onep3*prsl(i,j,k) + if(lmint(i) < k)then + esmax=0.1_r_kind*pw + esmax=min(esmax,estmax(i)) + es=min(es,esmax) end if - end do - end if - end if + qs = max(qmin, eps * es / (pw - omeps * es)) + q(i,j,k) = q(i,j,k)/qs + end do + end do + end do return -end subroutine general_getprs_glb +end subroutine genqsat2 + diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index fe2e058dff..fc87026c98 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -591,7 +591,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs,ig) ny=grd_loc%nlat ; nx=grd_loc%nlon ; nz=nlevs if(vvlocal)then -!$omp parallel do schedule(dynamic,1) private(k,j,i,l) +!$omp parallel do schedule(static,1) private(k,j,i,l) do k=1,nz if(iadvance == 1) then @@ -634,7 +634,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs,ig) enddo else -!$omp parallel do schedule(dynamic,1) private(k,j,i,l) +!$omp parallel do schedule(static,1) private(k,j,i,l) do k=1,nz if(iadvance == 1) then @@ -1607,7 +1607,7 @@ subroutine fix_belt(z) real(r_kind) zloc1(ny,nx) integer(i_kind) i,ii,j,jj,k -!$omp parallel do schedule(dynamic,1) private(j,k,i,jj,ii,zloc1) +!$omp parallel do schedule(static,1) private(j,k,i,jj,ii,zloc1) do j=1,nscl do k=1,nnnn1o i=0 @@ -1686,7 +1686,7 @@ subroutine rescale_ensemble_rh_perturbations end if do m=1,ntlevs_ens do ig=1,ntotensgrp -!$omp parallel do schedule(dynamic,1) private(n,i,j,k,w3,istatus) +!$omp parallel do schedule(static,1) private(n,i,j,k,w3,istatus) do n=1,n_ens call gsi_bundlegetpointer(en_perts(n,ig,m),'q',w3,istatus) if(istatus/=0) then @@ -1835,7 +1835,7 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) ipx=1 -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig,iaens) +!$omp parallel do schedule(static,1) private(j,n,ic3,k,i,ipic,ig,iaens) do k=1,km do ic3=1,nc3d ipic=ipc3d(ic3) @@ -1860,7 +1860,6 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) enddo enddo -!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig,iaens) do ic2=1,nc2d ipic=ipc2d(ic2) do j=1,jm @@ -2012,7 +2011,7 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) im=work_ens%grid%im jm=work_ens%grid%jm km=work_ens%grid%km -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig,iaens) +!$omp parallel do schedule(static,1) private(j,n,ic3,k,i,ipic,ig,iaens) do k=1,km do ic3=1,nc3d ipic=ipc3d(ic3) @@ -2036,7 +2035,6 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) enddo enddo enddo -!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig,iaens) do ic2=1,nc2d ipic=ipc2d(ic2) do j=1,jm @@ -2189,7 +2187,7 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin) endif ipx=1 -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens) +!$omp parallel do schedule(static,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens) do n=1,n_ens do ig=1,ntotensgrp do ic3=1,nc3d @@ -2206,6 +2204,7 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin) enddo endif ! iaens>0 enddo + do ic2=1,nc2d iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin) if(iaens>0) then @@ -2354,7 +2353,7 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin) im=a_en(1,1)%grid%im jm=a_en(1,1)%grid%jm km=a_en(1,1)%grid%km -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens) +!$omp parallel do schedule(static,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens) do n=1,n_ens do ig=1,ntotensgrp do ic3=1,nc3d @@ -2686,7 +2685,7 @@ subroutine sqrt_beta_s_mult_cvec(grady) endif ! multiply by sqrt_beta_s -!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i,ii) +!$omp parallel do schedule(static,1) private(ic3,ic2,k,j,i,ii) do j=1,lon2 do ii=1,nsubwin do ic3=1,nc3d @@ -2784,7 +2783,7 @@ subroutine sqrt_beta_s_mult_bundle(grady) endif ! multiply by sqrt_beta_s -!$omp parallel do schedule(dynamic,1) private(ic3,ic2,k,j,i) +!$omp parallel do schedule(static,1) private(ic3,ic2,k,j,i) do j=1,lon2 do ic3=1,nc3d ! check for ozone and skip if oz_univ_static = true @@ -2864,12 +2863,12 @@ subroutine sqrt_beta_e_mult_cvec(grady) call timer_ini('sqrt_beta_e_mult') ! multiply by sqrt_beta_e -!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii,ig) - do j=1,grd_ens%lon2 +!$omp parallel do schedule(static,1) private(nn,k,j,i,ii,ig) + do nn=1,n_ens do ii=1,nsubwin do ig=1,naensgrp - do nn=1,n_ens - do k=1,nsig + do k=1,nsig + do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 grady%aens(ii,ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*grady%aens(ii,ig,nn)%r3(1)%q(i,j,k) enddo @@ -2931,11 +2930,11 @@ subroutine sqrt_beta_e_mult_bundle(aens) call timer_ini('sqrt_beta_e_mult') ! multiply by sqrt_beta_e -!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ig) - do j=1,grd_ens%lon2 +!$omp parallel do schedule(static,1) private(nn,k,j,i,ig) + do nn=1,n_ens do ig=1,naensgrp - do nn=1,n_ens - do k=1,nsig + do k=1,nsig + do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 aens(ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*aens(ig,nn)%r3(1)%q(i,j,k) enddo @@ -2993,19 +2992,18 @@ subroutine init_sf_xy(jcap_in) integer(i_kind),intent(in ) :: jcap_in - integer(i_kind) i,ii,j,k,l,n,jcap,kk,nsigend,ig - real(r_kind),allocatable::g(:),gsave(:) + integer(i_kind) i,ii,j,igg,k,l,n,jcap,kk,nsigend,ig + real(r_kind),allocatable::g(:),gtemp(:) real(r_kind) factor real(r_kind),allocatable::rkm(:),f(:,:),f0(:,:) real(r_kind) ftest(grd_loc%nlat,grd_loc%nlon,grd_loc%kbegin_loc:grd_loc%kend_alloc) real(r_single) out1(grd_ens%nlon,grd_ens%nlat) - real(r_single),allocatable::pn0_npole(:) + real(r_single) pn0_npole real(r_kind) s_ens_h_min real(r_kind) rlats_ens_local(grd_ens%nlat) real(r_kind) rlons_ens_local(grd_ens%nlon) character(5) mapname logical make_test_maps - logical,allocatable,dimension(:)::ksame integer(i_kind) nord_sploc2ens integer(i_kind) nlon_sploc0,nlon_sploc,nlat_sploc,num_fields logical print_verbose @@ -3157,101 +3155,107 @@ subroutine init_sf_xy(jcap_in) if(.not.allocated(spectral_filter)) allocate(spectral_filter(naensloc,sp_loc%nc,grd_sploc%nsig)) if(.not.allocated(sqrt_spectral_filter)) allocate(sqrt_spectral_filter(naensloc,sp_loc%nc,grd_sploc%nsig)) - allocate(g(sp_loc%nc),gsave(sp_loc%nc)) - allocate(pn0_npole(0:sp_loc%jcap)) - allocate(ksame(grd_sploc%nsig)) + allocate(g(sp_loc%nc),gtemp(sp_loc%nc)) do ig=1,naensloc - ksame=.false. - do k=2,grd_sploc%nsig - if(s_ens_hv(k,ig) == s_ens_hv(k-1,ig))ksame(k)=.true. - enddo spectral_filter(ig,:,:)=zero - do k=1,grd_sploc%nsig - if(ksame(k))then - spectral_filter(ig,:,k)=spectral_filter(ig,:,k-1) - else + level_loop: do k=1,grd_sploc%nsig + do kk=1,k-1 + if(s_ens_hv(k,ig) == s_ens_hv(kk,ig))then + spectral_filter(ig,:,k)=spectral_filter(ig,:,k-1) + cycle level_loop + end if + end do + if(ig > 1)then + do igg=1,ig-1 + do kk=1,grd_sploc%nsig + if(s_ens_hv(k,ig) == s_ens_hv(kk,igg))then + spectral_filter(ig,:,k)=spectral_filter(igg,:,kk) + cycle level_loop + end if + end do + end do + end if + + do i=1,grd_sploc%nlat + f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k,ig))**2) + enddo + + + do j=2,grd_sploc%nlon do i=1,grd_sploc%nlat - f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k,ig))**2) + f0(i,j)=f0(i,1) enddo + end do - do j=2,grd_sploc%nlon - do i=1,grd_sploc%nlat - f0(i,j)=f0(i,1) - enddo - enddo - call general_g2s0(grd_sploc,sp_loc,g,f0) + call general_g2s0(grd_sploc,sp_loc,g,f0) - call general_s2g0(grd_sploc,sp_loc,g,f) + call general_s2g0(grd_sploc,sp_loc,g,f) -! adjust so value at np = 1 - f=f/f(grd_sploc%nlat,1) - f0=f - call general_g2s0(grd_sploc,sp_loc,g,f) - call general_s2g0(grd_sploc,sp_loc,g,f) - if(mype == 0)then - nsigend=k - do kk=k+1,grd_sploc%nsig - if(s_ens_hv(kk,ig) /= s_ens_hv(k,ig))exit - nsigend=nsigend+1 - enddo - write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k,ig),maxval(abs(f0-f)) -900 format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', & +! adjust so value at np = 1 + f=f/f(grd_sploc%nlat,1) + f0=f + call general_g2s0(grd_sploc,sp_loc,g,f) + call general_s2g0(grd_sploc,sp_loc,g,f) + if(mype == 0)then + nsigend=k + do kk=k+1,grd_sploc%nsig + if(s_ens_hv(kk,ig) /= s_ens_hv(k,ig))exit + nsigend=nsigend+1 + enddo + write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k,ig),maxval(abs(f0-f)) +900 format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', & i10,f10.2,e20.10) - end if + end if -! correct spectrum by dividing by pn0_npole - gsave=g +! correct spectrum by dividing by pn0_npole -! obtain pn0_npole - do n=0,sp_loc%jcap - g=zero - g(2*n+1)=one - call general_s2g0(grd_sploc,sp_loc,g,f) - pn0_npole(n)=f(grd_sploc%nlat,1) - enddo +! obtain pn0_npole +!$omp parallel do schedule(static,1) private(n,gtemp,f) + do n=0,sp_loc%jcap + gtemp=zero + gtemp(2*n+1)=one + call general_s2g0(grd_sploc,sp_loc,gtemp,f) + pn0_npole=f(grd_sploc%nlat,1) + g(2*n+1)=g(2*n+1)/pn0_npole + enddo - g=zero - do n=0,sp_loc%jcap - g(2*n+1)=gsave(2*n+1)/pn0_npole(n) - enddo -! obtain spectral_filter +! obtain spectral_filter - ii=0 - do l=0,sp_loc%jcap - if(ig>naensgrp) then - factor=one/g(1) + ii=0 + do l=0,sp_loc%jcap + if(ig>naensgrp) then + factor=one/g(1) + else + factor=one + if(l>0) factor=half + end if + do n=l,sp_loc%jcap + ii=ii+1 + if(sp_loc%factsml(ii)) then + spectral_filter(ig,ii,k)=zero else - factor=one - if(l>0) factor=half + spectral_filter(ig,ii,k)=factor*g(2*n+1) + end if + ii=ii+1 + if(l == 0 .or. sp_loc%factsml(ii)) then + spectral_filter(ig,ii,k)=zero + else + spectral_filter(ig,ii,k)=factor*g(2*n+1) end if - do n=l,sp_loc%jcap - ii=ii+1 - if(sp_loc%factsml(ii)) then - spectral_filter(ig,ii,k)=zero - else - spectral_filter(ig,ii,k)=factor*g(2*n+1) - end if - ii=ii+1 - if(l == 0 .or. sp_loc%factsml(ii)) then - spectral_filter(ig,ii,k)=zero - else - spectral_filter(ig,ii,k)=factor*g(2*n+1) - end if - enddo enddo - end if - enddo + enddo + enddo level_loop enddo !ig loop - deallocate(g,gsave,pn0_npole,ksame) + deallocate(g,gtemp) ! Compute sqrt(spectral_filter). Ensure spectral_filter >=0 zero -!$omp parallel do schedule(dynamic,1) private(k,i) +!$omp parallel do schedule(static,1) private(k,i) do ig=1,naensloc do k=1,grd_sploc%nsig do i=1,sp_loc%nc - if (spectral_filter(ig,i,k) < zero) spectral_filter(ig,i,k)=zero + spectral_filter(ig,i,k) = max(spectral_filter(ig,i,k),zero) sqrt_spectral_filter(ig,i,k) = sqrt(spectral_filter(ig,i,k)) end do end do @@ -3337,13 +3341,14 @@ subroutine sf_xy(ig,f,k_start,k_end) if(.not.use_localization_grid) then if(ig>naensgrp) then +!$omp parallel do schedule(static,1) private(k,g) do k=k_start,k_end call general_g2s0(grd_ens,sp_loc,g,f(:,:,k)) g(:)=g(:)*spectral_filter(ig,:,k_index(k)) call general_s2g0(grd_ens,sp_loc,g,f(:,:,k)) enddo else -!$omp parallel do schedule(dynamic,1) private(k) +!$omp parallel do schedule(static,1) private(k) do k=k_start,k_end call sfilter(grd_ens,sp_loc,spectral_filter(ig,:,k_index(k)),f(1,1,k)) enddo @@ -3353,6 +3358,7 @@ subroutine sf_xy(ig,f,k_start,k_end) vector=.false. if(ig>naensgrp) then +!$omp parallel do schedule(static,1) private(k,g,work) do k=k_start,k_end call g_agrid2egrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) call general_g2s0(grd_ens,sp_loc,g,f(:,:,k)) @@ -3361,7 +3367,7 @@ subroutine sf_xy(ig,f,k_start,k_end) call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) enddo else -!$omp parallel do schedule(dynamic,1) private(k,work) +!$omp parallel do schedule(static,1) private(k,work) do k=k_start,k_end call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) call sfilter(grd_ens,sp_loc,spectral_filter(ig,:,k_index(k)),f(1,1,k)) @@ -3421,6 +3427,7 @@ subroutine sqrt_sf_xy(ig,z,f,k_start,k_end) if(.not.use_localization_grid) then +!$omp parallel do schedule(static,1) private(k,g) do k=k_start,k_end g(:)=z(:,k)*sqrt_spectral_filter(ig,:,k_index(k)) call general_s2g0(grd_ens,sp_loc,g,f(:,:,k)) @@ -3429,6 +3436,7 @@ subroutine sqrt_sf_xy(ig,z,f,k_start,k_end) else vector=.false. +!$omp parallel do schedule(static,1) private(k,g,work) do k=k_start,k_end g(:)=z(:,k)*sqrt_spectral_filter(ig,:,k_index(k)) call general_s2g0(grd_sploc,sp_loc,g,work) @@ -3488,6 +3496,7 @@ subroutine sqrt_sf_xy_ad(ig,z,f,k_start,k_end) if(.not.use_localization_grid) then +!$omp parallel do schedule(static,1) private(k,g) do k=k_start,k_end call general_s2g0_ad(grd_ens,sp_loc,g,f(:,:,k)) z(:,k)=g(:)*sqrt_spectral_filter(ig,:,k_index(k)) @@ -3496,6 +3505,7 @@ subroutine sqrt_sf_xy_ad(ig,z,f,k_start,k_end) else vector=.false. +!$omp parallel do schedule(static,1) private(k,g,work) do k=k_start,k_end call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) call general_s2g0_ad(grd_sploc,sp_loc,g,work) @@ -3608,7 +3618,7 @@ subroutine bkerror_a_en(grady) ! Declare local variables integer(i_kind) ii,ip,istatus,k,ig,ig2 real(r_kind),allocatable,dimension(:,:) :: z - real(r_kind),allocatable,dimension(:) :: ztmp + real(r_kind),allocatable,dimension(:) :: z2 ! Initialize timer call timer_ini('bkerror_a_en') @@ -3624,34 +3634,29 @@ subroutine bkerror_a_en(grady) call sqrt_beta_e_mult(grady) ! Apply variances, as well as vertical & horizontal parts of background error -! !$omp parallel do schedule(dynamic,1) private(ii) - do ii=1,nsubwin - if (naensgrp==1) then + if (naensgrp==1) then + do ii=1,nsubwin call bkgcov_a_en_new_factorization(1,grady%aens(ii,1,1:n_ens)) - else - allocate(z(naensgrp,nval_lenz_en)) + end do + else + allocate(z(nval_lenz_en,naensgrp)) + allocate(z2(nval_lenz_en)) + do ii=1,nsubwin do ig=1,naensgrp - call ckgcov_a_en_new_factorization_ad(ig,z(ig,:),grady%aens(ii,ig,1:n_ens)) + call ckgcov_a_en_new_factorization_ad(ig,z(1,ig),grady%aens(ii,ig,1:n_ens)) enddo - allocate(ztmp(naensgrp)) - do k=1,nval_lenz_en - ztmp=zero - do ig=1,naensgrp - do ig2=1,naensgrp - ztmp(ig) = ztmp(ig) + z(ig2,k) * alphacvarsclgrpmat(ig,ig2) + do ig=1,naensgrp + z2=zero + do ig2=1,naensgrp + do k=1,nval_lenz_en + z2(k) = z2(k) + z(k,ig2) * alphacvarsclgrpmat(ig,ig2) enddo enddo - do ig=1,naensgrp - z(ig,k) = ztmp(ig) - enddo - enddo - deallocate(ztmp) - do ig=1,naensgrp - call ckgcov_a_en_new_factorization(ig,z(ig,:),grady%aens(ii,ig,1:n_ens)) + call ckgcov_a_en_new_factorization(ig,z2,grady%aens(ii,ig,1:n_ens)) enddo - deallocate(z) - endif - enddo + enddo + deallocate(z,z2) + endif ! multiply by sqrt_beta_e_mult call sqrt_beta_e_mult(grady) @@ -3710,11 +3715,7 @@ subroutine bkgcov_a_en_new_factorization(ig,a_en) real(r_kind) hwork(grd_loc%inner_vars,grd_loc%nlat,grd_loc%nlon,grd_loc%kbegin_loc:grd_loc%kend_alloc) real(r_kind),allocatable,dimension(:):: a_en_work - call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus) - if(istatus/=0) then - write(6,*)'bkgcov_a_en_new_factorization: trouble getting pointer to ensemble CV' - call stop2(999) - endif + ipnt=1 ! Apply vertical smoother on each ensemble member ! To avoid my having to touch the general sub2grid and grid2sub, @@ -3725,7 +3726,7 @@ subroutine bkgcov_a_en_new_factorization(ig,a_en) call stop2(999) endif iadvance=1 ; iback=2 -!$omp parallel do schedule(dynamic,1) private(k,ii,is,ie) +!$omp parallel do schedule(static,1) private(k,ii,is,ie) do k=1,n_ens call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig) ii=(k-1)*a_en(1)%ndim @@ -3755,7 +3756,7 @@ subroutine bkgcov_a_en_new_factorization(ig,a_en) ! Retrieve ensemble components from long vector ! Apply vertical smoother on each ensemble member iadvance=2 ; iback=1 -!$omp parallel do schedule(dynamic,1) private(k,ii,is,ie) +!$omp parallel do schedule(static,1) private(k,ii,is,ie) do k=1,n_ens ii=(k-1)*a_en(1)%ndim is=ii+1 @@ -3865,9 +3866,10 @@ subroutine ckgcov_a_en_new_factorization(ig,z,a_en) deallocate(a_en_work) ! Apply vertical smoother on each ensemble member + iadvance=2 ; iback=1 +!$omp parallel do schedule(static,1) private(k) do k=1,n_ens - iadvance=2 ; iback=1 call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig) enddo @@ -3937,9 +3939,10 @@ subroutine ckgcov_a_en_new_factorization_ad(ig,z,a_en) endif ! Apply vertical smoother on each ensemble member + iadvance=1 ; iback=2 +!$omp parallel do schedule(static,1) private(k) do k=1,n_ens - iadvance=1 ; iback=2 call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig) enddo @@ -4125,7 +4128,7 @@ subroutine hybens_grid_setup end if if(global_spectral_filter_sd .and. nsclgrp > 1)then - allocate(spc_multwgt(0:jcap_ens,nsclgrp)) + allocate(spc_multwgt(sp_ens%nc,nsclgrp)) allocate(spcwgt_params(4,nsclgrp)) spc_multwgt=1.0 diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index 038188f92a..4e688cd36e 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -829,7 +829,11 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! Prevent out of bounds reference from temperature if ( bufr_index(l) == 0 ) cycle i = bufr_index(l) - data_all(l+nreal,itx) = temperature(i) ! brightness temerature + if(i /= 0)then + data_all(l+nreal,itx) = temperature(i) ! brightness temerature + else + data_all(l+nreal,itx) = tbmin + end if end do nrec(itx)=irec diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index b72e584155..9efd06418c 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -2050,9 +2050,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if(oelev>7000.0_r_kind) cycle loop_k_levs if(abs(diffvv)>5.0_r_kind.and.oelev<5000.0_r_kind) cycle loop_k_levs ! write(6,*)'sliu diffuu,vv::',diffuu, diffvv - uob=0.0 - vob=0.0 - oelev=0.0 + uob=zero + vob=zero + oelev=zero tkk=0 do ikkk=k,klev diffhgt=obsdat(4,ikkk)-obsdat(4,k) diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index a30ef92a90..694c8f1df3 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -45,13 +45,22 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di use mpeu_util, only: die,perr use kinds, only: r_kind,r_single,r_double,i_kind + use constants, only: zero,one,r1000,r10,r100 + use constants, only: huge_single,wgtlim,three + use constants, only: tiny_r_kind,five,half,two,r0_01 + use constants, only: zero,one, h1000 + use obsmod, only: rmiss_single,time_offset + use obsmod, only: netcdf_diag, binary_diag, dirname, ianldate + use obsmod, only: luse_obsdiag + use m_obsLList, only: obsLList + use m_obsdiagNode, only: obs_diags use m_obsNode, only: obsNode use m_qNode, only: qNode use m_qNode, only: qNode_appendto + use m_dtime, only: dtime_setup, dtime_check, dtime_show use gsi_4dvar, only: nobs_bins,hr_obsbin - use obsmod, only: netcdf_diag, binary_diag, dirname, ianldate use nc_diag_write_mod,only: nc_diag_init, nc_diag_header,nc_diag_metadata, & nc_diag_write, nc_diag_data2d use nc_diag_read_mod, only: nc_diag_read_init,nc_diag_read_get_dim, & @@ -60,26 +69,18 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di use guess_grids, only: geop_hgtl,hrdifsig,nfldsig,ges_tsen,ges_prsl use gridmod, only: nsig,get_ijk - use constants, only: zero,one,r1000,r10,r100 - use constants, only: huge_single,wgtlim,three - use constants, only: tiny_r_kind,five,half,two,r0_01 use qcmod, only: npres_print use jfunc, only: jiter use convinfo, only: nconvtype use convinfo, only: icsubtype - use m_dtime, only: dtime_setup, dtime_check, dtime_show use rapidrefresh_cldsurf_mod, only: i_cloud_q_innovation, & cld_bld_hgt,i_ens_mean use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use mpimod, only: mpi_comm_world - use constants, only: zero,one, h1000 use gsdcloudlib_pseudoq_mod, only: cloudLWC_pseudo,cloudCover_Surface_col - use m_obsLList, only: obsLList - use m_obsdiagNode, only: obs_diags - use obsmod, only: luse_obsdiag implicit none diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index ebdd8de52a..f686e19332 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -776,6 +776,12 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& call stop2(275) end if + if (abi2km .and. regional) then + abi2km_bc = zero + abi2km_bc(2) = 233.5_r_kind + abi2km_bc(3) = 241.7_r_kind + abi2km_bc(4) = 250.5_r_kind + end if ! PROCESSING OF SATELLITE DATA ! Loop over data in this block @@ -1085,10 +1091,6 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& endif predbias=zero - abi2km_bc = zero - abi2km_bc(2) = 233.5_r_kind - abi2km_bc(3) = 241.7_r_kind - abi2km_bc(4) = 250.5_r_kind !$omp parallel do schedule(dynamic,1) private(i,mm,j,k,tlap,node,bias) do i=1,nchanl mm=ich(i) diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90 index dd60703ce2..849d2ff5c9 100644 --- a/src/gsi/stpcalc.f90 +++ b/src/gsi/stpcalc.f90 @@ -263,7 +263,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & real(r_quad),parameter:: one_tenth_quad = 0.1_r_quad ! Declare local variables - integer(i_kind) i,j,mm1,ii,iis,final_ii,ibin,ipenloc,it + integer(i_kind) i,j,mm1,ii,final_ii,ibin,ipenloc,it integer(i_kind) istp_use,nstep,nsteptot,kprt real(r_quad),dimension(4,ipen):: pbc real(r_quad),dimension(4,nobs_type):: pbcjo @@ -430,7 +430,6 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc=zero_quad pjcalc=.false. if(iter == 0 .and. kprt >= 2 .and. ii == 1)pjcalc=.true. - iis=ii ! Delta stepsize sges(1)= stp(ii-1) From ca19008b521cccf07a84037cfece3847232025cc Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Mon, 25 Sep 2023 16:12:33 -0400 Subject: [PATCH 4/7] Intel 2022 updates (#629) **Description** This PR fixes two types of bugs discovered when cycling `gsi.x` and `enkf.x` with intel/2022 in the global workflow 1. modify variables written to netcdf diagnostic files by `gsi.x` to be consistent with codes which read netcdf diagnostic files 2. modify `lrun_subdirs=.true.` option of `gsi.x` to properly handle the case in which sub-directories already exist in the run directory Fixes #623 **Type of change** - [x] Bug fix (non-breaking change which fixes an issue) **How Has This Been Tested?** Ctests have been on Hera, Orion, and WCOSS2 (Cactus) with acceptable behavior. A global parallel covering the period 2021073106 through 2021080118 has been run on Hera, Orion, and WCOSS2 (Cactus). All global workflow jobs ran as expected. **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 --- src/gsi/obsmod.F90 | 17 +++++++++++------ src/gsi/setupaod.f90 | 2 +- src/gsi/setupoz.f90 | 5 +++-- src/gsi/setuprad.f90 | 2 +- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 26f8ff1bbf..633bde91ab 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -1022,12 +1022,12 @@ subroutine init_directories(in_pe,num_pe) integer(i_kind),intent(in ) :: in_pe integer(i_kind),intent(in ) :: num_pe - logical :: l_mkdir_stat + logical :: l_mkdir_stat, l_dir_exist character(len=144):: command character(len=8):: pe_name, loc_pe_name character(len=128):: loc_dirname - integer(i_kind) :: i + integer(i_kind) :: i, ierror if (lrun_subdirs) then write(pe_name,'(i4.4)') in_pe @@ -1038,10 +1038,15 @@ subroutine init_directories(in_pe,num_pe) write(loc_pe_name,'(i4.4)') i loc_dirname = 'dir.'//trim(loc_pe_name) #ifdef __INTEL_COMPILER - l_mkdir_stat = MAKEDIRQQ(trim(loc_dirname)) - if(.not. l_mkdir_stat) then - write(6, *) "Failed to create directory ", trim(loc_dirname), " for PE ", loc_pe_name - call stop2(678) + INQUIRE(directory=trim(loc_dirname), exist=l_dir_exist) + if (.not.l_dir_exist) then + l_mkdir_stat = MAKEDIRQQ(trim(loc_dirname)) + if(.not.l_mkdir_stat) then + ierror=GETLASTERRORQQ() + write(6, *) "INIT_DIRECTORIES: ***ERROR** Failed to create directory ", & + trim(loc_dirname)," for PE ", loc_pe_name, ' ierror= ', ierror + call stop2(678) + endif endif #else command = 'mkdir -p -m 755 ' // trim(loc_dirname) diff --git a/src/gsi/setupaod.f90 b/src/gsi/setupaod.f90 index 5fe4233ada..58707acd6a 100644 --- a/src/gsi/setupaod.f90 +++ b/src/gsi/setupaod.f90 @@ -844,7 +844,7 @@ subroutine contents_netcdf_diag_ call nc_diag_metadata("Observation_Class", obsclass) call nc_diag_metadata_to_single("Latitude",(cenlat)) ! observation latitude (degrees) call nc_diag_metadata_to_single("Longitude",(cenlon)) ! observation longitude (degrees) - call nc_diag_metadata_to_single("Obs_Time",(dtime))!-time_offset)) ! observation time (hours relative to analysis time) + call nc_diag_metadata_to_single("Time",(dtime))!-time_offset)) ! observation time (hours relative to analysis time) call nc_diag_metadata_to_single("Sol_Zenith_Angle",(pangs)) ! solar zenith angle (degrees) call nc_diag_metadata_to_single("Sol_Azimuth_Angle",(data_s(isazi_ang,n))) ! solar azimuth angle (degrees) call nc_diag_metadata("Surface_type", nint(data_s(istyp,n))) diff --git a/src/gsi/setupoz.f90 b/src/gsi/setupoz.f90 index d7a85de0b2..7112e967ba 100644 --- a/src/gsi/setupoz.f90 +++ b/src/gsi/setupoz.f90 @@ -1720,6 +1720,7 @@ subroutine contents_netcdf_diag_(odiag) type(obs_diag),pointer,intent(in):: odiag ! Observation class character(7),parameter :: obsclass = ' ozlev' + integer(i_kind),parameter :: ione = 1 real(r_kind),dimension(miter) :: obsdiag_iuse call nc_diag_metadata_to_single("Latitude", data(ilate,i) ) call nc_diag_metadata_to_single("Longitude", data(ilone,i) ) @@ -1731,9 +1732,9 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata_to_single("Obs_Minus_Forecast_unadjusted",ozone_inv ) call nc_diag_metadata_to_single("Reference_Pressure", preso3l*r100 ) ! Pa if(luse(i)) then - call nc_diag_metadata_to_single("Analysis_Use_Flag", one ) + call nc_diag_metadata("Analysis_Use_Flag", ione ) else - call nc_diag_metadata_to_single("Analysis_Use_Flag", -one ) + call nc_diag_metadata("Analysis_Use_Flag", -ione ) endif call nc_diag_metadata_to_single("Input_Observation_Error",obserror ) diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index f686e19332..20ab63456e 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -2577,7 +2577,7 @@ subroutine contents_netcdf_diag_(odiags,idv,iob) call nc_diag_metadata_to_single("Elevation",zsges ) ! model (guess) elevation at observation location - call nc_diag_metadata_to_single("Time",dtime,time_offset,'-') + call nc_diag_metadata_to_single("Obs_Time",dtime,time_offset,'-') call nc_diag_metadata_to_single("Scan_Position",data_s(iscan_pos,n) ) ! sensor scan position call nc_diag_metadata_to_single("Sat_Zenith_Angle", zasat,rad2deg,'*') ! satellite zenith angle (degrees) From 728d006c36ae6b17f08b380ff23b3422b97624fe Mon Sep 17 00:00:00 2001 From: Gang <53267411+GangZhao-NOAA@users.noreply.github.com> Date: Fri, 29 Sep 2023 09:11:39 -0400 Subject: [PATCH 5/7] Adding code to analyze the siginificant wave heigh in GSI 3D Analysis (see issue #601) (#614) Adding code to analyze the siginificant wave heigh in GSI 3D Analysis, esp. for FV3-LAM model based DA, eg. RRFS-DA, RRFS-3DRTMA. (Also see the issue in EMC GSI github repository: #601 Adding I/O for Analysis of Significant Wave Height for 3DRTMA) **Description** Significant Wave Height (hereafter as SWH) is one of the standard products provided by the operational (2D)RTMA. To continuously provide the same products in 3DRTMA, the next-generation RTMA, some efforts in GSI code need to be made in order to analyze the SWH in 3D analysis of GSI. The kernel subroutines to assimilate SWH in GSI (such as stphowv.f90, setuphowv.f90, inthowv.f90, gsi_howvOper.f90 and m_howvNode.f90) already had been added for (2D)RTMA years ago by Manuel Pondeca, so for this issue, the code work mainly focus on adding the I/O of SWH in background and analysis fields for 3DRTMA (esp. RRFS-based 3DRTMA), and some necessary modifications in background error, options, variables related to analysis of SWH, etc. Modified code in GSI: 1. rapidrefresh_cldsurf_mod.f90: adding a few variables related to the analysis of howv in 3D analysis 2. gsimod.F90: adding namelist options used for analysis of howv in 3D analysis 3. m_berror_stats_reg.f90: added some code for the special treatment to the static background error (BE) of howv 4. read_prepbufr.f90: adding code to decode the observation of howv in prepbufr file when howv is available in firstguess 5. setuphowv.f90: adding code to use obs of howv when howv is available in firstguess 6. gsi_rfv3io_mod.f90: adding I/O code to read in howv from firstguess and write out howv into analysis. No dependencies are required for this change. This PR is addressing the issue [#601](https://github.com/NOAA-EMC/GSI/issues/601): Adding code to analyze the siginificant wave heigh in GSI 3D Analysis". Fixes #601 **Type of change** Please delete options that are not relevant. - [*] New feature (non-breaking change which adds functionality) **How Has This Been Tested?** - Brief results from ctest (regression test) with the modified code (on WCOSS2 - Cactus): [gang.zhao@clogin07:build] (feature/3drtma_howv)$ ctest -N Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Test #1: global_3dvar Test #2: global_4dvar Test #3: global_4denvar Test #4: hwrf_nmm_d2 Test #5: hwrf_nmm_d3 Test #6: rtma Test #7: rrfs_3denvar_glbens Test #8: netcdf_fv3_regional Test #9: global_enkf Total Tests: 9 Test #1: global_3dvar [gang.zhao@clogin04:build] (feature/3drtma_howv)$ ctest -R global_3dvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 1: global_3dvar 1/1 Test #1: global_3dvar ..................... Passed 1631.12 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1631.14 sec Test #2: global_4dvar [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R global_4dvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 2: global_4dvar 1/1 Test #2: global_4dvar ..................... Passed 2462.19 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 2462.23 sec Test #3: global_4denvar [gang.zhao@clogin04:build] (feature/3drtma_howv)$ ctest -R global_4denvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 3: global_4denvar 1/1 Test #3: global_4denvar ................... Passed 1922.43 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1922.46 sec Test #4: hwrf_nmm_d2 [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R hwrf_nmm_d2 Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 4: hwrf_nmm_d2 1/1 Test #4: hwrf_nmm_d2 ...................... Passed 1214.10 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1214.20 sec Test #5: hwrf_nmm_d3 [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R hwrf_nmm_d3 Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 5: hwrf_nmm_d3 1/1 Test #5: hwrf_nmm_d3 ...................... Passed 736.38 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 736.50 sec Test #6: rtma [gang.zhao@clogin05:build] (feature/3drtma_howv)$ ctest -R rtma Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 6: rtma 1/1 Test #6: rtma ............................. Passed 1027.01 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1027.01 sec Test #7: rrfs_3denvar_glbens [gang.zhao@clogin06:build] (feature/3drtma_howv)$ ctest -R rrfs_3denvar_glbens Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 7: rrfs_3denvar_glbens 1/1 Test #7: rrfs_3denvar_glbens .............. Passed 484.69 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 484.70 sec Test #8: netcdf_fv3_regional [gang.zhao@clogin03:build] (feature/3drtma_howv)$ ctest -R netcdf_fv3_regional Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 8: netcdf_fv3_regional 1/1 Test #8: netcdf_fv3_regional .............. Passed 483.08 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 483.11 sec Test #9: global_enkf [gang.zhao@clogin03:build] (feature/3drtma_howv)$ ctest -R global_enkf Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 9: global_enkf 1/1 Test #9: global_enkf ...................... Passed 488.50 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 488.57 sec - The modified GSI code passed the regression tests (all 9 tasks) on Hera and WCOSS2 (Cactus). - adding the analysis of howv only has very trial influences on the analyses of other variables. Here is the statistics of the differences of other variables (u/v/t/ps/q/t2m/q2m) from the runs of GSI without howv vs. with howv (from a testing case 2023-07-12_14:00:00 UTC in 3km North-American domain): comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/fv_core.res.tile1.nc fcst_fv3lam_nodata_noinfo/INPUT/fv_core.res.tile1.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev u / 602135550 3926.84 25760.8 -0.1026 0.485788 0.588388 6.52152e-06 0.00115817 v / 620166777 -4891.34 32582.5 -0.835774 0.268402 1.10418 -7.88714e-06 0.00197793 T / 155987083 178.048 6497.51 -0.0246582 0.0384064 0.0630646 1.14143e-06 0.000218737 delp / 19559676 -281.532 3008.29 -0.00292969 0.00219727 0.00512695 -1.43935e-05 0.000183727 comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/fv_tracer.res.tile1.nc fcst_fv3lam_nodata_noinfo/INPUT/fv_tracer.res.tile1.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev sphum / 430707614 0.594287 2.77816 -2.6139e-05 3.1759e-05 5.7898e-05 1.37979e-09 8.03072e-08 comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/sfc_data.nc fcst_fv3lam_nodata_noinfo/INPUT/sfc_data.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev t2m / 10665000 43.3899 135.095 -0.00152825 0.00686629 0.00839454 4.06844e-06 5.02866e-05 q2m / 10665000 0.0192553 0.124707 -3.1476e-06 1.77554e-05 2.0903e-05 1.80547e-09 5.89657e-08 It could be seen that the differences are trivial and ignorable. The regression tests were done by following the instructions of "[GSI Ctests (regression tests)](https://github.com/NOAA-EMC/GSI/wiki/GSI-Ctests-(regression-tests))" in [GSI Wiki](https://github.com/NOAA-EMC/GSI/wiki) The modified code had also been tested with a testing case 2023-07-12_14:00:00 UTC for 3km North-American domain Here is a brief summary of the test results: 1. Here is the analysis increment of Significant Wave Height (aka howv hereafter): pure 3dvar, static background error of howv is 0.42 meters, and the de-correlation length scale is 170km. ![HOWV_var_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/4fdeeb82-7258-4344-be69-cce747474312) 2. The following figure shows the distribution of howv in the analysis (used obs is in green, rejected in red). Obviously the location of used obs of howv match the area of non-zero analysis increments of howv. ![var_obs_2023071214_howv_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/d4ed6013-cfc8-486e-8f47-db07ec0e4e53) 3. The following figure is the analysis increment of howv with hybrid envar analysis (using gdas ensemble 80 members and the ensemble weight is 84%), and the static BE of howv is tuned/inflated. The analysis increments are very similar to the results from pure 3dvar run (see the first figure) ![HOWV_hyb_betas016_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/e6e696e8-932b-42ab-9001-3472e970b21c) 4. The last figure shows the analysis increments of howv with hybrid envar analysis (using gdas ensemble 80 members and the ensemble weight is 84%), but the static BE of howv is NOT tuned. It can be observed that the analysis increments is less than the results from the hybrid run with tuning the static BE of howv. That is because the weight of static BE (16%) reduced the background error of howv (ensemble of howv is not available yet), so the impact of obs is decreased. ![HOWV_hyb_betas016_noTune_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/ca25d068-fc86-4d47-a9d2-46e02ac22dac) **Checklist** - [*] My code follows the style guidelines of this project - [*] I have performed a self-review of my own code - [*] I have commented my code, particularly in hard-to-understand areas - [*] New and existing tests pass with my changes - [*] Any dependent changes have been merged and published **DUE DATE for this PR is 10/5/2023.** If this PR is not merged into `develop` by this date, the PR will be closed and returned to the developer. --- src/gsi/gsi_rfv3io_mod.f90 | 78 +++++++++++++++++++++++++--- src/gsi/gsimod.F90 | 13 ++++- src/gsi/m_berror_stats_reg.f90 | 41 +++++++++++---- src/gsi/rapidrefresh_cldsurf_mod.f90 | 41 ++++++++++++++- src/gsi/read_prepbufr.f90 | 7 +++ 5 files changed, 160 insertions(+), 20 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 62b23ee713..6d16be7c13 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -22,6 +22,8 @@ module gsi_rfv3io_mod ! used as background for surface observation operator ! 2022-04-15 Wang - add IO for regional FV3-CMAQ (RRFS-CMAQ) model ! 2022-08-10 Wang - add IO for regional FV3-SMOKE (RRFS-SMOKE) model +! 2023-07-30 Zhao - add IO for the analysis of the significant wave height +! (SWH, aka howv in GSI) in fv3-lam based DA (eg., RRFS-3DRTMA) ! ! subroutines included: ! sub gsi_rfv3io_get_grid_specs @@ -56,6 +58,7 @@ module gsi_rfv3io_mod use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke + use rapidrefresh_cldsurf_mod, only: i_howv_3dda implicit none public type_fv3regfilenameg @@ -133,7 +136,7 @@ module gsi_rfv3io_mod public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc - public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m + public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv @@ -144,7 +147,7 @@ module gsi_rfv3io_mod integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc - integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m + integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv parameter( & k_f10m =1, & !fact10 k_stype=2, & !soil_type @@ -159,7 +162,8 @@ module gsi_rfv3io_mod k_t2m =11, & ! 2 m T k_q2m =12, & ! 2 m Q k_orog =13, & !terrain - n2d=13 ) + k_howv =14, & ! significant wave height (aka howv in GSI) + n2d=14 ) logical :: grid_reverse_flag character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields @@ -767,6 +771,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ! 2022-04-01 Y. Wang and X. Wang - add capability to read reflectivity ! for direct radar EnVar DA using reflectivity as state ! variable, poc: xuguang.wang@ou.edu +! 2023-07-30 Zhao - added code to read significant wave height (howv) field +! from the 2D fv3-lam firstguess file (fv3_sfcdata). ! attributes: ! language: f90 ! machine: ibm RS/6000 SP @@ -816,6 +822,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,:),pointer::ges_t2m=>NULL() real(r_kind),dimension(:,:),pointer::ges_q2m=>NULL() + real(r_kind),dimension(:,:),pointer::ges_howv=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() @@ -1093,6 +1100,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if(mype == 0) write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' else if(trim(vartem)=='t2m') then else if(trim(vartem)=='q2m') then + else if(trim(vartem)=='howv') then else write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop' call stop2(333) @@ -1110,7 +1118,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) do i=1,size(name_metvars2d) vartem=trim(name_metvars2d(i)) if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z") & - .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m"))) then !z is treated separately + .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m") & + .or.(trim(vartem)=="howv"))) then ! z is treated separately if (ifindstrloc(vardynvars,trim(vartem)) > 0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem) @@ -1365,6 +1374,13 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus if (ier/=0) call die(trim(myname),'cannot get pointers for t2m,ier=',ier) endif + +!--- significant wave height (howv) + if ( i_howv_3dda == 1 ) then + call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus ); ier=ier+istatus + if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier) + endif + if(mype == 0 ) then call check(nf90_open(fv3filenamegin(it)%dynvars,nf90_nowrite,loc_id)) call check(nf90_inquire(loc_id,formatNum=ncfmt)) @@ -1546,7 +1562,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif - call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m) + call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv) if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then ! Convert 2m guess mixing ratio to specific humidity @@ -1782,7 +1798,7 @@ end subroutine gsi_bundlegetpointer_fv3lam_tracerchem_nouv end subroutine read_fv3_netcdf_guess -subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) +subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf2d_read @@ -1792,6 +1808,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ! Scatter the field to each PE ! program history log: ! 2023-02-14 Hu - Bug fix for read in subdomain surface restart files +! 2023-07-30 Zhao - added IO to read significant wave height (howv) from 2D FV3-LAM +! firstguess file (fv3_sfcdata) ! ! input argument list: ! it - time index for 2d fields @@ -1805,7 +1823,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ! !$$$ end documentation block use kinds, only: r_kind,i_kind - use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype + use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype,mpi_itype + use mpeu_util, only: die use guess_grids, only: fact10,soil_type,veg_frac,veg_type,sfc_rough, & sfct,sno,soil_temp,soil_moi,isli use gridmod, only: lat2,lon2,itotsub,ijn_s @@ -1813,8 +1832,11 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_noerr use mod_fv3_lola, only: fv3_h_to_ll,nxa,nya use constants, only: grav + use constants, only: zero implicit none @@ -1822,6 +1844,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) real(r_kind),intent(in),dimension(:,:),pointer::ges_z real(r_kind),intent(in),dimension(:,:),pointer::ges_t2m real(r_kind),intent(in),dimension(:,:),pointer::ges_q2m + real(r_kind),intent(in),dimension(:,:),pointer::ges_howv type (type_fv3regfilenameg),intent(in) :: fv3filenamegin character(len=max_varname_length) :: name integer(i_kind),allocatable,dimension(:):: dim @@ -1835,6 +1858,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) integer(i_kind) kk,n,ns,j,ii,jj,mm1 character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: dynvars !='fv3_dynvars' +! for checking the existence of howv in firstguess file + integer(i_kind) id_howv + integer(i_kind) iret_bcast ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: sfc_fulldomain @@ -1850,6 +1876,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) allocate(work(itotsub*n2d)) allocate( sfcn2d(lat2,lon2,n2d)) +!-- initialisation of the array for howv + sfcn2d(:,:,k_howv) = zero + if(mype==mype_2d ) then allocate(sfc_fulldomain(nx,ny)) @@ -1877,6 +1906,24 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo + +!--- check the existence of significant wave height (howv) in 2D FV3-LAM firstguess file +! if howv is set in anavinfo (as i_howv_3dda=1), then check its existence in firstguess, +! but if it is not found in firstguess, then stop GSI run and set i_howv_3dda = 0. + if ( i_howv_3dda == 1 ) then + iret = nf90_inq_varid(gfile_loc,'howv',id_howv) + if ( iret /= nf90_noerr ) then + iret = nf90_inq_varid(gfile_loc,'HOWV',id_howv) ! double check with name in uppercase + end if + if ( iret /= nf90_noerr ) then + i_howv_3dda = 0 ! howv does not exist in firstguess, then stop GSI run. + call die('gsi_fv3ncdf2d_read','Warning: CANNOT find howv in firstguess, aborting..., iret = ', iret) + else + write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'gsi_fv3ncdf2d_read:: Found howv in firstguess ', & + trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').' + end if + end if + !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!! do i=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,i,name,len) @@ -1904,6 +1951,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) k=k_t2m else if( trim(name)=='Q2M'.or.trim(name)=='q2m' ) then k=k_q2m + else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then + k=k_howv else cycle endif @@ -2036,6 +2085,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain) endif ! mype +!-- broadcast the updated i_howv_3dda to all tasks (!!!!) + call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,& @@ -2058,6 +2109,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ges_t2m(:,:)=sfcn2d(:,:,k_t2m) ges_q2m(:,:)=sfcn2d(:,:,k_q2m) endif + if ( i_howv_3dda == 1 ) then + ges_howv(:,:)=sfcn2d(:,:,k_howv) + endif deallocate (sfcn2d,a) return end subroutine gsi_fv3ncdf2d_read @@ -3192,6 +3246,8 @@ subroutine wrfv3_netcdf(fv3filenamegin) ! 2019-11-22 CAPS(C. Tong) - modify "add_saved" to properly output analyses ! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da ! 2022-04-01 Y. Wang and X. Wang - add code for updating reflectivity +! 2023-07-30 Zhao - added code for the output of the analysis of +! significant wave height (howv) ! ! input argument list: ! @@ -3234,6 +3290,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_t2m =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_q2m =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_howv =>NULL() integer(i_kind) i,k @@ -3350,6 +3407,9 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus endif + if ( i_howv_3dda == 1 ) then + call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus + endif if (ier/=0) call die('wrfv3_netcdf','cannot get pointers for fv3 met-fields, ier =',ier) if (laeroana_fv3cmaq) then @@ -3559,6 +3619,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_fv3ncdf_write_sfc(fv3filenamegin,'t2m',ges_t2m,add_saved) call gsi_fv3ncdf_write_sfc(fv3filenamegin,'q2m',ges_q2m,add_saved) endif +!-- output analysis of howv + if ( i_howv_3dda == 1 ) then + call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved) + endif if(allocated(g_prsi)) deallocate(g_prsi) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 2656a2dce4..70618120d0 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -176,7 +176,8 @@ module gsimod i_coastline,i_gsdqc,qv_max_inc,ioption,l_precip_clear_only,l_fog_off,& cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax,& - i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & + corp_howv, hwllp_howv use gsi_metguess_mod, only: gsi_metguess_init,gsi_metguess_final use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax @@ -503,6 +504,9 @@ module gsimod ! 3. fv3_cmaq_regional = .true. ! 4. berror_fv3_cmaq_regional = .true. ! 09-15-2022 yokota - add scale/variable/time-dependent localization +! 2023-07-30 Zhao - added namelist options for analysis of significant wave height +! (aka howv in GSI code): corp_howv, hwllp_howv +! (in namelist session rapidrefresh_cldsurf) ! !EOP !------------------------------------------------------------------------- @@ -1560,6 +1564,10 @@ module gsimod ! = 2(clean Qg as in 1, and adjustment to the retrieved Qr/Qs/Qnr throughout the whole profile) ! = 3(similar to 2, but adjustment to Qr/Qs/Qnr only below maximum reflectivity level ! and where the dbz_obs is missing); +! corp_howv - real, static background error of howv (stddev error) +! = 0.42 meters (default) +! hwllp_howv - real, background error de-correlation length scale of howv +! = 170,000.0 meters (default 170 km) ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & @@ -1580,7 +1588,8 @@ module gsimod i_coastline,i_gsdqc,qv_max_inc,ioption,l_precip_clear_only,l_fog_off,& cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax, & - i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & + corp_howv, hwllp_howv ! chem(options for gsi chem analysis) : ! berror_chem - .true. when background for chemical species that require diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 2ff8a6aa94..601339e1ac 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -870,16 +870,17 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt hwllp(i,n)=hwllp(i,nrf2_ps) end do else if (n==nrf2_howv) then - call read_howv_stats(mlat,1,2,cov_dum) + call read_howv_stats(mlat,1,2,cov_dum,mype) do i=1,mlat corp(i,n)=cov_dum(i,1,1) !#ww3 hwllp(i,n) = cov_dum(i,1,2) end do hwllp(0,n) = hwllp(1,n) hwllp(mlat+1,n) = hwllp(mlat,n) - - if (mype==0) print*, 'corp(i,n) = ', corp(:,n) - if (mype==0) print*, ' hwllp(i,n) = ', hwllp(:,n) + if (mype==0) then + print*, myname_, ' static BE corp( :,n) (for ', trim(adjustl(cvars2d(n))), ')= ', corp(:,n) + print*, myname_, ' static BE hwllp(:,n) (for ', trim(adjustl(cvars2d(n))), ')= ', hwllp(:,n) + end if ! corp(:,n)=cov_dum(:,1) !do i=1,mlat ! corp(i,n)=0.4_r_kind !#ww3 @@ -1055,7 +1056,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt end subroutine berror_read_wgt_reg !++++ -subroutine read_howv_stats(nlat,nlon,npar,arrout) +subroutine read_howv_stats(nlat,nlon,npar,arrout,mype) !$$$ subprogram documentation block ! . . . . ! subprogram: read_howv_stats @@ -1090,6 +1091,9 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! program history log: ! 2016-08-03 stelios ! 2016-08-26 stelios : Compatible with GSI. +! 2023-07-30 Zhao - added code to set the background error +! standard deviation (corp_howv) and de-correlation +! length scale (hwllp_howv) for non-2DRTMA run ! input argument list: ! filename - The name of the file ! output argument list: @@ -1102,10 +1106,14 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) !$$$ end documentation block ! use kinds,only : r_kind, i_kind + use gridmod, only : twodvar_regional + use rapidrefresh_cldsurf_mod, only : corp_howv, hwllp_howv + use gsi_io, only : verbose ! implicit none ! Declare passed variables integer(i_kind), intent(in )::nlat,nlon,npar + integer(i_kind), intent(in ) :: mype ! "my" processor ID real(r_kind), dimension(nlat ,nlon, npar), intent( out)::arrout ! Declare local variables integer(i_kind) :: reclength,i,j,i_npar @@ -1117,12 +1125,18 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! filename(1) = 'howv_var_berr.bin' filename(2) = 'howv_lng_berr.bin' -! - arrout(:,:,1)=0.42_r_kind - arrout(:,:,2)=50000.0_r_kind +!-- first, assign the pre-defined values to corp and hwllp + if ( twodvar_regional ) then + arrout(:,:,1)=0.42_r_kind ! values were specified by Manuel and Stelio for 2DRTMA + arrout(:,:,2)=50000.0_r_kind ! values were specified by Manuel and Stelio for 2DRTMA + else + arrout(:,:,1) = corp_howv ! 0.42_r_kind used in 3dvar (default) if not set in namelist + arrout(:,:,2) = hwllp_howv ! 17000.0_r_kind used in 3dvar (default) if not set in namelist + end if reclength=nlat*r_kind -! +!-- secondly, if files for corp and hwllp are available, then read them in for +! corp and hwllp. If the files are not found, then use the pre-defined values. do i_npar = 1,npar inquire(file=trim(filename(i_npar)), exist=file_exists) if (file_exists)then @@ -1132,9 +1146,16 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) read(unit=lun34 ,rec=j) (arrout(i,j,i_npar), i=1,nlat) enddo close(unit=lun34) + if (verbose .and. mype .eq. 0) then + write(6,'(1x,A,1x,A2,1x,A)') trim(adjustl(myname)), '::', & + trim(filename(i_npar))//' is used for background error of howv.' + end if else - print*,myname, trim(filename(i_npar)) // ' does not exist' + if (verbose .and. mype .eq. 0) then + write(6,'(1x,A,1x,A2,1x,A)') trim(adjustl(myname)), '::', & + trim(filename(i_npar))//' does not exist for static BE of howv, using pre-defined values.' + end if end if end do end subroutine read_howv_stats diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 1ee35fffba..122d2872d0 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -28,7 +28,11 @@ module rapidrefresh_cldsurf_mod ! option for checking and adjusting the profile of Qr/Qs/Qg/Qnr ! retrieved through cloud analysis to reduce the background ! reflectivity ghost in analysis. (default is 0) -! +! 2023-07-30 Zhao added options for analysis of significant wave height +! (SWH, aka howv in GSI code): +! corp_howv: to set the static background error of howv +! hwllp_howv: to set the de-correlation length scale +! i_howv_3dda: control the analysis of howv in 3D analysis (if howv is in anavinfo) ! ! Subroutines Included: ! sub init_rapidrefresh_cldsurf - initialize RR related variables to default values @@ -181,6 +185,18 @@ module rapidrefresh_cldsurf_mod ! = 2(clean Qg as in 1, and adjustment to the retrieved Qr/Qs/Qnr throughout the whole profile) ! = 3(similar to 2, but adjustment to Qr/Qs/Qnr only below maximum reflectivity level ! and where the dbz_obs is missing); +! corp_howv - namelist real, static BE of howv (standard error deviation) +! hwllp_howv - namelist real, static BE de-correlation length scale of howv +! i_howv_3dda - integer, control the analysis of howv in 3D analysis (either var or hybrid) +! = 0 (howv-off: default) : no analysis of howv in 3D analysis. +! = 1 (howv-on) : if variable name "howv" is found in anavinfo, +! set it to be 1 to turn on analysis of howv; +! note: in hybrid envar run, the static BE is redueced by beta_s (<1.0), +! since there is no ensemble of howv currently yet, then no ensemble +! contribution to the total BE of howv, so the total BE of howv is actually +! just the reduced static BE of howv. If to make the analysis of howv +! in hyrbid run is as similar as the analysis of howv in pure 3dvar run, +! the static BE of howv used in hybrid run needs to be tuned (inflated actually). ! ! attributes: ! language: f90 @@ -252,6 +268,8 @@ module rapidrefresh_cldsurf_mod public :: l_saturate_bkCloud public :: l_rtma3d public :: i_precip_vertical_check + public :: corp_howv, hwllp_howv + public :: i_howv_3dda logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -310,6 +328,8 @@ module rapidrefresh_cldsurf_mod logical l_saturate_bkCloud logical l_rtma3d integer(i_kind) i_precip_vertical_check + real(r_kind) :: corp_howv, hwllp_howv + integer(i_kind) :: i_howv_3dda contains @@ -325,6 +345,8 @@ subroutine init_rapidrefresh_cldsurf ! 2008-06-03 Hu initial build for cloud analysis ! 2010-03-29 Hu change names to init_rapidrefresh_cldsurf ! 2011--5-04 Todling inquire MetGuess for presence of hyrometeors & set default +! 2023-07-30 Zhao added code for initialization of some variables used +! in analysis of significant wave height ! ! input argument list: ! @@ -337,8 +359,12 @@ subroutine init_rapidrefresh_cldsurf !$$$ use kinds, only: i_kind use gsi_metguess_mod, only: gsi_metguess_get + use mpimod, only: mype + use state_vectors, only: ns2d,svars2d + implicit none integer(i_kind) ivar,i,ier + integer(i_kind) i2 logical have_hmeteor(5) character(len=2),parameter :: hydrometeors(5) = (/ 'qi', & 'ql', & @@ -418,6 +444,19 @@ subroutine init_rapidrefresh_cldsurf l_saturate_bkCloud= .true. l_rtma3d = .false. ! turn configuration for rtma3d off i_precip_vertical_check = 0 ! No check and adjustment to retrieved Qr/Qs/Qg (default) + corp_howv = 0.42_r_kind ! 0.42 meters (default) + hwllp_howv = 170000.0_r_kind ! 170,000.0 meters (170km as default for 3DRTMA, 50km is used in 2DRTMA) + i_howv_3dda = 0 ! no analysis of significant wave height (howv) in 3D analysis (default) + +!-- searching for specific variable in state variable list (reading from anavinfo) + do i2=1,ns2d + if ( trim(svars2d(i2))=='howv' .or. trim(svars2d(i2))=='HOWV' ) then + i_howv_3dda = 1 + if ( mype == 0 ) then + write(6,'(1x,A,1x,A8,1x,A,1x,I4)')"init_rapidrefresh_cldsurf: anavinfo svars2d (state variable): ",trim(adjustl(svars2d(i2))), " is found in anavinfo, set i_howv_3dda = ", i_howv_3dda + end if + end if + end do ! i2 : looping over 2-D anasv return end subroutine init_rapidrefresh_cldsurf diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 9efd06418c..304fa62590 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -149,6 +149,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! 2020-05-04 wu - no rotate_wind for fv3_regional ! 2020-09-05 CAPS(C. Tong) - add flag for new vadwind obs to assimilate around the analysis time only ! 2023-03-23 draper - add code for processing T2m and q2m for global system +! 2023-07-30 Zhao - added code to extract obs of significant wave height (howvob) from bufr record +! in prepbufr file for 3D analysis ! input argument list: ! infile - unit from which to read BUFR data @@ -1132,6 +1134,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (howvob) call ufbint(lunin,owave,1,255,levs,owavestr) if (cldchob) call ufbint(lunin,cldceilh,1,255,levs,cldceilhstr) endif +! Extract obs of howv in 3D Analysis +! (if-block is to avoid potential issue if decoding the bufr record twice in 2DRTMA run) + if ( .not. twodvar_regional ) then + if (howvob) call ufbint(lunin,owave,1,255,levs,owavestr) + endif if(kx==224 .and. newvad) then call ufbint(lunin,fcstdat,3,255,levs,'UFC VFC TFC ') end if From ba5a2ca673cc54212dc1c3d1cb2a81f96b3ce075 Mon Sep 17 00:00:00 2001 From: hongli-wang <53354098+hongli-wang@users.noreply.github.com> Date: Fri, 29 Sep 2023 11:51:33 -0600 Subject: [PATCH 6/7] Refine PM2.5 DA for the RRFS_SD model (#609) **Description** Refine the PM2.5 DA for the RRFS_SD model by use of veg_type, which is used to decide whether the obs is in urban area or not. Different thresholds for innovations outside/inside urban areas will be used. Add new namelist parameters, such as threshold for innovations, anowbufr type Read in station terrain height, PM10 et al if extended BUFR format for anow air quality data is used Fixes #606 **Type of change** Please delete options that are not relevant. - [ ] Bug fix (non-breaking change which fixes an issue) - [X ] New feature (non-breaking change which adds functionality) - [ ] Breaking change (fix or feature that would cause existing functionality to not work as expected) - [ ] This change requires a documentation update **How Has This Been Tested?** **Checklist** - [X ] My code follows the style guidelines of this project - [X] I have performed a self-review of my own code - [X] I have commented my code, particularly in hard-to-understand areas - [ ] New and existing tests pass with my changes - [ ] Any dependent changes have been merged and published **DUE DATE for this PR is 9/22/2023.** If this PR is not merged into develop by this date, the PR will be closed and returned to the developer. --- src/gsi/chemmod.f90 | 40 ++++++++++++++++++------- src/gsi/gsimod.F90 | 15 ++++++++-- src/gsi/m_berror_stats_reg.f90 | 6 ++-- src/gsi/read_anowbufr.f90 | 50 +++++++++++++++++++++++++++----- src/gsi/satthin.F90 | 7 ++++- src/gsi/setuppm2_5.f90 | 53 ++++++++++++++++++++++------------ 6 files changed, 128 insertions(+), 43 deletions(-) diff --git a/src/gsi/chemmod.f90 b/src/gsi/chemmod.f90 index 14a90c818c..06bfe6dce6 100644 --- a/src/gsi/chemmod.f90 +++ b/src/gsi/chemmod.f90 @@ -40,21 +40,23 @@ module chemmod public :: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3 ! fv3smoke - public :: naero_smoke_fv3,aeronames_smoke_fv3,pm2_5_innov_threshold + public :: naero_smoke_fv3,aeronames_smoke_fv3 + public :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold + public :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold public :: naero_gocart_wrf,aeronames_gocart_wrf public :: pm2_5_guess,init_pm2_5_guess,& aerotot_guess,init_aerotot_guess public :: init_chem - public :: berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs + public :: berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs public :: oblat_chem,oblon_chem,obpres_chem,diag_incr,oneobschem public :: site_scale,nsites public :: tunable_error public :: in_fname,out_fname,incr_fname,maxstr public :: code_pm25_ncbufr,code_pm25_anowbufr public :: code_pm10_ncbufr,code_pm10_anowbufr - + public :: anowbufr_ext public :: l_aoderr_table public :: laeroana_gocart,laeroana_fv3cmaq,laeroana_fv3smoke,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & @@ -79,7 +81,8 @@ module chemmod integer(i_kind) :: icvt_cmaq_fv3 real(r_kind) :: raod_radius_mean_scale,raod_radius_std_scale real(r_kind) :: ppmv_conv = 96.06_r_kind/28.964_r_kind*1.0e+3_r_kind - real(r_kind) :: pm2_5_innov_threshold + real(r_kind) :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold + real(r_kind) :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold logical :: wrf_pm2_5 @@ -90,7 +93,9 @@ module chemmod logical :: aero_ratios - logical :: oneobtest_chem,diag_incr,berror_chem,berror_fv3_cmaq_regional + logical :: oneobtest_chem,diag_incr,berror_chem + logical :: berror_fv3_cmaq_regional,berror_fv3_sd_regional + logical :: anowbufr_ext character(len=max_varname_length) :: oneob_type_chem integer(i_kind), parameter :: maxstr=256 real(r_kind) :: maginnov_chem,magoberr_chem,conconeobs,& @@ -103,7 +108,7 @@ module chemmod real(r_kind),parameter :: pm2_5_teom_max=900.0_r_kind !ug/m3 !some parameters need to be put here since convinfo file won't !accomodate, stands for maximum realistic value of surface pm2.5 - real(r_kind),parameter :: pm10_teom_max=150.0_r_kind !ug/m3 + real(r_kind),parameter :: pm10_teom_max=3000.0_r_kind !ug/m3 real(r_kind),parameter :: elev_missing=-9999.0_r_kind @@ -157,10 +162,10 @@ module chemmod 'AOLGAJ', 'AISO1J', 'AISO2J', 'AISO3J', 'ATRP1J', 'ATRP2J',& 'ASQTJ', 'AOLGBJ', 'AORGCJ'] ! fv3smoke - integer(i_kind), parameter :: naero_smoke_fv3=2 + integer(i_kind), parameter :: naero_smoke_fv3=3 character(len=max_varname_length), dimension(naero_smoke_fv3), parameter :: & - aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust' ] + aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust','coarsepm'] ! FV3CMAQ integer(i_kind), parameter :: naero_cmaq_fv3=70 ! !number of cmaq aerosol species aero6 @@ -286,8 +291,15 @@ subroutine init_chem !initialiazes default values to &CHEM namelist parameters berror_chem=.false. - berror_fv3_cmaq_regional=.false. ! Set .true. to use berror for fv3_cmaq_regional, whose cv has 10 characters + berror_fv3_cmaq_regional=.false. ! .False. : Dont perform aerosal DA for the online RRFS_CMAQ model so dont need to read in B for RRFS_CMAQ. + ! .true. : Use berror for fv3_cmaq_regional, whose cv has 10 characters + berror_fv3_sd_regional=.false. ! .False. : Dont perform aerosal DA for the RRFS_SD model so dont need to read in B for RRFS_SD. + ! .true. to use berror for rrfs_sd model, whose cv has 10 characters oneobtest_chem=.false. + anowbufr_ext=.false. ! .False. : use default anowbufr data + ! .True. : use the extented bufr data + ! that includes PM10, station elevation + ! etal in addition to pm2.5. maginnov_chem=30_r_kind magoberr_chem=2_r_kind oneob_type_chem='pm2_5' @@ -307,9 +319,15 @@ subroutine init_chem laeroana_gocart = .false. laeroana_fv3cmaq = .false. ! .true. for performing aerosol analysis for regional FV3-CMAQ model(Please other parameters requred in gsimod.F90) laeroana_fv3smoke = .false. - pm2_5_innov_threshold = 20.0_r_kind + pm2_5_innov_threshold = 15.0_r_kind + pm2_5_urban_innov_threshold = 30.0_r_kind + pm2_5_bg_threshold = 2.0_r_kind + pm10_innov_threshold = 15.0_r_kind + pm10_urban_innov_threshold = 30.0_r_kind + pm10_bg_threshold = 2.0_r_kind + pm10_obs_threshold = 140.0_r_kind ! Barry's manuscript l_aoderr_table = .false. - icvt_cmaq_fv3 = 1 ! 1. Control variable is individual aerosol specie; 2: CV is total mass per I,J,K mode + icvt_cmaq_fv3 = 1 ! 1: Control variable is individual aerosol specie; 2: CV is total mass per I,J,K mode raod_radius_mean_scale = 1.0_r_kind ! Tune radius of particles when calculating AOD using CRTM raod_radius_std_scale = 1.0_r_kind ! Tune standard deviation of particles when calculating AOD using CRTM with CMAQ LUTs. aod_qa_limit = 3 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 70618120d0..d0ca1c0fbf 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -182,12 +182,15 @@ module gsimod use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax use chemmod, only : init_chem,berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,& + berror_fv3_sd_regional,& maginnov_chem,magoberr_chem,& oneob_type_chem,oblat_chem,& + anowbufr_ext,& oblon_chem,obpres_chem,diag_incr,elev_tolerance,tunable_error,& in_fname,out_fname,incr_fname, & laeroana_gocart, l_aoderr_table, aod_qa_limit, luse_deepblue, lread_ext_aerosol, & - laeroana_fv3cmaq,laeroana_fv3smoke,pm2_5_innov_threshold,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & + laeroana_fv3cmaq,laeroana_fv3smoke,pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold,& + crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & icvt_cmaq_fv3, raod_radius_mean_scale,raod_radius_std_scale use chemmod, only : wrf_pm2_5,aero_ratios @@ -1594,6 +1597,12 @@ module gsimod ! chem(options for gsi chem analysis) : ! berror_chem - .true. when background for chemical species that require ! conversion to lower case and/or species names longer than 5 chars +! berror_fv3_cmaq_regional - .true. use background error stat for online +! RRFS_CMAQ model. Control variable +! names extended up to 10 chars +! berror_fv3_sd_regional - .true. use background error stat for online +! RRFS_SD model. Control variable +! names extended up to 10 chars ! oneobtest_chem - one-ob trigger for chem constituent analysis ! maginnov_chem - O-B make-believe residual for one-ob chem test ! magoberr_chem - make-believe obs error for one-ob chem test @@ -1615,13 +1624,15 @@ module gsimod ! luse_deepblue - whether to use MODIS AOD from the deepblue algorithm ! lread_ext_aerosol - if true, reads aerfNN file for aerosol arrays rather than sigfNN (NGAC NEMS IO) - namelist/chem/berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,maginnov_chem,magoberr_chem,& + namelist/chem/berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,& + oneobtest_chem,anowbufr_ext,maginnov_chem,magoberr_chem,& oneob_type_chem,oblat_chem,oblon_chem,obpres_chem,& diag_incr,elev_tolerance,tunable_error,& in_fname,out_fname,incr_fname,& laeroana_gocart, laeroana_fv3cmaq,laeroana_fv3smoke,l_aoderr_table, aod_qa_limit, & crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & icvt_cmaq_fv3,pm2_5_innov_threshold, & + pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold,& raod_radius_mean_scale,raod_radius_std_scale, luse_deepblue,& aero_ratios,wrf_pm2_5, lread_ext_aerosol diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 601339e1ac..bf9fb20674 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -12,7 +12,7 @@ module m_berror_stats_reg use kinds,only : i_kind,r_kind use constants, only: zero,one,max_varname_length,half use gridmod, only: nsig - use chemmod, only : berror_chem,berror_fv3_cmaq_regional,upper2lower,lower2upper + use chemmod, only : berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,upper2lower,lower2upper use m_berror_stats, only: usenewgfsberror,berror_stats implicit none @@ -312,7 +312,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt use constants, only: zero,one,ten,three use mpeu_util,only: getindex use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd - use chemmod, only: berror_fv3_cmaq_regional + use chemmod, only: berror_fv3_cmaq_regional,berror_fv3_sd_regional implicit none @@ -466,7 +466,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt var=upper2lower(varshort) if (trim(var) == 'pm25') var = 'pm2_5' else - if ( berror_fv3_cmaq_regional) then + if ( berror_fv3_cmaq_regional .or. berror_fv3_sd_regional) then read(inerr,iostat=istat) varlong, isig var=varlong else diff --git a/src/gsi/read_anowbufr.f90 b/src/gsi/read_anowbufr.f90 index 1873d0b877..449ce5cdf8 100644 --- a/src/gsi/read_anowbufr.f90 +++ b/src/gsi/read_anowbufr.f90 @@ -50,7 +50,9 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& iconc,ierror,ilat,ilon,itime,iid,ielev,isite,iikx,ilate,ilone,& elev_missing,site_scale,tunable_error,& code_pm25_ncbufr,code_pm25_anowbufr,& - code_pm10_ncbufr,code_pm10_anowbufr + code_pm10_ncbufr,code_pm10_anowbufr,& + anowbufr_ext,pm2_5_teom_max,pm10_teom_max + use mpimod, only: npe implicit none @@ -71,10 +73,12 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& nyob=3,ndhr=4,ntyp=5,ncopopm=6 !see headr input format below integer(i_kind), parameter :: nfields=6 + integer(i_kind), parameter :: nfields_b=12 !output format parameters integer(i_kind), parameter:: nchanl=0,nreal=ilone real(r_kind),parameter :: r360 = 360.0_r_kind + real(r_kind),parameter :: r90 = 90.0_r_kind real(r_kind),parameter :: percent=1.e-2_r_kind real(r_kind), parameter :: anow_missing=1.0e11_r_kind,& @@ -96,8 +100,10 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& real(r_kind), dimension(5) :: rinc character(len=8) :: subset character(len=80) :: headr + character(len=80) :: obstr real(r_double), dimension(nfields) :: indata + real(r_double), dimension(nfields_b) :: indata_a,indata_b real(r_kind) :: tdiff,obstime,t4dv real(r_kind) :: dlat,dlon,error_1,error_2,obserror,dlat_earth,dlon_earth @@ -141,7 +147,6 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& ! reading each report from bufr do while (ireadmg(lunin,subset,idate) == 0) - if (trim(obstype)=='pm2_5') then if ( (subset == 'NC008031') .or. (subset == 'NC008032' ) ) then @@ -149,9 +154,16 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& ncbufr=.true. write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset else if (subset == 'ANOWPM') then - headr='SID XOB YOB DHR TYP COPOPM' - anowbufr=.true. - write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset + if (anowbufr_ext) then + headr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG' + obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO' + anowbufr=.true. + write(6,*)'READ_PM2_5_BUFR_EXT: AIRNOW data type, subset=',subset + else ! default ANOWBUFR Table + headr='SID XOB YOB DHR TYP COPOPM' + anowbufr=.true. + write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset + end if else cycle endif @@ -162,6 +174,17 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& headr='PTID CLONH CLATH TPHR TYPO COPOPM' ncbufr=.true. write(6,*)'READ_PM10: AIRNOW data type, subset=',subset + else if (subset == 'ANOWPM') then + if (anowbufr_ext) then + headr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG' + obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO' + anowbufr=.true. + write(6,*)'READ_PM10_BUFR_EXT: AIRNOW data type, subset=',subset + else + headr='SID XOB YOB DHR TYP COPOPM' + anowbufr=.true. + write(6,*)'READ_PM10: AIRNOW data type, subset=',subset + end if else cycle endif @@ -176,8 +199,17 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& imin=0 do while (ireadsb(lunin) == 0) - call ufbint(lunin,indata,nfields,1,iret,headr) - + if (anowbufr_ext) then + call ufbint(lunin,indata_a,nfields_b,1,iret,headr) + indata(1:5) = indata_a(1:5) + call ufbint(lunin,indata_b,nfields_b,1,iret,obstr) + if (trim(obstype)=='pm2_5') indata(ncopopm)=indata_b(3) + if (trim(obstype)=='pm10') indata(ncopopm)=indata_b(5) + site_elev = indata_b(4) + else + call ufbint(lunin,indata,nfields,1,iret,headr) + end if + if (anowbufr) then kx=indata(ntyp) read(sid,'(Z8)')site_id @@ -198,13 +230,15 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& nread = nread + 1 conc=indata(ncopopm) - if ( iret > 0 .and. (conc < conc_missing ) .and. & (conc >= zero)) then if(indata(nxob) >= r360) indata(nxob) = indata(nxob) - r360 if(indata(nxob) < zero) indata(nxob) = indata(nxob) + r360 + if(indata(nxob) > r360)cycle + if(indata(nyob) > r90)cycle + dlon_earth_deg=indata(nxob) dlat_earth_deg=indata(nyob) diff --git a/src/gsi/satthin.F90 b/src/gsi/satthin.F90 index 02d19b198c..2018d80be7 100644 --- a/src/gsi/satthin.F90 +++ b/src/gsi/satthin.F90 @@ -134,6 +134,8 @@ module satthin use obsmod, only: time_window_max use constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,& rad2deg,r1000 + use chemmod, only: laeroana_fv3smoke + implicit none ! set default to private @@ -961,7 +963,10 @@ subroutine getsfc(mype,mype_io,use_sfc,use_sfc_any) end if if (.not.lobserver) then if(allocated(veg_frac)) deallocate(veg_frac) - if(allocated(veg_type)) deallocate(veg_type) +! veg_type will be used in setuppm2_5.f90 for rrfs_sd PM2.5 DA + if(.not. laeroana_fv3smoke )then + if(allocated(veg_type)) deallocate(veg_type) + endif if(allocated(soil_type)) deallocate(soil_type) if(allocated(soil_moi)) deallocate(soil_moi) if(allocated(sfc_rough)) deallocate(sfc_rough) diff --git a/src/gsi/setuppm2_5.f90 b/src/gsi/setuppm2_5.f90 index 79e6129cd8..ad940cce78 100644 --- a/src/gsi/setuppm2_5.f90 +++ b/src/gsi/setuppm2_5.f90 @@ -94,8 +94,8 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) use gsi_4dvar, only: nobs_bins,hr_obsbin use gridmod, only : get_ij,get_ijk - use guess_grids, only : nfldsig,hrdifsig + use guess_grids, only : veg_type use gsi_bundlemod, only : gsi_bundlegetpointer,GSI_BundlePrint use gsi_chemguess_mod, only : gsi_chemguess_get,gsi_chemguess_bundle use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle @@ -115,7 +115,8 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) use chemmod, only: naero_gocart_wrf,aeronames_gocart_wrf,& upper2lower,lower2upper,laeroana_gocart,wrf_pm2_5 use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq - use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke,pm2_5_innov_threshold + use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke + use chemmod, only: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold use gridmod, only : cmaq_regional,wrf_mass_regional,fv3_cmaq_regional implicit none @@ -146,7 +147,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) real(r_kind) :: pm2_5ges real(r_kind) :: ratio_errors,error real(r_kind) :: innov,innov_error2,rwgt,valqc,tfact,innov_error,elevges,& - elevdiff,conc,elevobs,ps_ges,site_id,tv_ges + elevdiff,conc,elevobs,ps_ges,site_id,tv_ges,veg_type_ges real(r_kind) errinv_input,errinv_adjst,errinv_final real(r_kind) err_input,err_adjst,err_final @@ -278,7 +279,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) call stop2(453) endif - do i=2,naero_smoke_fv3 + do i=2,naero_smoke_fv3-1 ! remove contribution from coarsepm aeroname=trim(aeronames_smoke_fv3(i)) call gsi_bundlegetpointer(gsi_chemguess_bundle(1),trim(aeroname),& rank3,ier) @@ -705,13 +706,18 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) if (wrf_mass_regional .or. fv3_cmaq_regional .or. laeroana_fv3smoke) then call tintrp2a11(ges_ps,ps_ges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) - call tintrp2a11(ges_tv(:,:,1,nfldsig),tv_ges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) conc=conc/(ps_ges*r1000/(rd*tv_ges)) endif - - +! + if (laeroana_fv3smoke) then + if (.not. allocated(veg_type)) then + print*,"VEG_TYPE NOT ALLOCATED, WILL NOT BE USED IN PM2.5 DA FOR RRFS_SD",mype + else + call intrp2a11(veg_type(:,:,1),veg_type_ges,dlat,dlon,mype) + endif + endif !if elevobs is known than calculate difference otherwise !assume that difference is acceptable @@ -740,17 +746,20 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) mype,nfldsig) innov = conc - pm2_5ges if (laeroana_fv3smoke) then - if ( -1.0*innov >= pm2_5_innov_threshold .or. & - (innov > pm2_5_innov_threshold .and. pm2_5ges >=1.0_r_kind).or. & - (conc >= 40.0_r_kind .and. pm2_5ges >=1.0_r_kind).or. & - conc >= 100.0_r_kind ) then - innov = innov + if ( veg_type_ges == 13.0_r_kind ) then + if (abs(innov) < pm2_5_urban_innov_threshold) then + muse(i)=.false. + end if else - innov = 0.0_r_kind + if (abs(innov) < pm2_5_innov_threshold) then + muse(i)=.false. + end if + end if + + if (pm2_5ges < pm2_5_bg_threshold) then muse(i)=.false. end if if (tv_ges-273.15_r_kind < 5.0_r_kind) then - innov = 0.0_r_kind muse(i)=.false. end if @@ -770,21 +779,28 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) mype,nfldsig) call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,& mype,nfldsig) - - if (pm25wc_ges(1) >= 1.0_r_kind) then + if (pm25wc_ges(1) >= pm2_5_bg_threshold) then pm25wc_ges(1)=1.0_r_kind else - pm25wc_ges(2)=0.0_r_kind + pm25wc_ges(1)=0.0_r_kind end if - if (pm25wc_ges(2) >= 1.0_r_kind) then + if (pm25wc_ges(2) >= pm2_5_bg_threshold) then pm25wc_ges(2)=1.0_r_kind else pm25wc_ges(2)=0.0_r_kind end if + if ( (pm25wc_ges(1)+pm25wc_ges(2)) < 1.0_r_kind ) then + muse(i) = .false. + end if else pm25wc_ges = 0.0_r_kind end if + if (oneobtest_chem) then + pm25wc_ges=1.0_r_kind + muse(i) = .true. + end if + error=one/data(ierror,i) ratio_errors=one/sqrt(real(dup(i))) innov_error = error*innov @@ -1115,6 +1131,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Latitude", data(ilate,i) ) call nc_diag_metadata("Longitude", data(ilone,i) ) call nc_diag_metadata("Station_Elevation", data(ielev,i) ) + call nc_diag_metadata("Station_Veg_Type", veg_type_ges ) call nc_diag_metadata("Pressure", ps_ges ) call nc_diag_metadata("Height", data(ielev,i) ) call nc_diag_metadata("Time", dtime-time_offset ) From c56d7bc616057054f592653a7b4fd0438deace4a Mon Sep 17 00:00:00 2001 From: shoyokota <103961291+shoyokota@users.noreply.github.com> Date: Fri, 29 Sep 2023 21:27:50 -0400 Subject: [PATCH 7/7] GitHub Issue NOAA-EMC/GSI#604 Undefined values found in radar reflectivity direct DA (#605) **Description** To prevent some undefined values found in the radar reflectivity direct DA (if_model_dbz=T and l_use_dbz_directDA=F), corresponding parts are fixed. It doesn't change the result except for the case of the execution with the debug option. Fixes #604 **Type of change** - [x] Bug fix (non-breaking change which fixes an issue) **How Has This Been Tested?** The radar reflectivity DA test with the RRFS setting was done on Orion. After this modification, EnVar was completed even with the debug option. This modification didn't change the result in the test without the debug option. **Checklist** - [x] My code follows the style guidelines of this project - [x] I have performed a self-review of my own code - [x] I have commented my code, particularly in hard-to-understand areas - [x] New and existing tests pass with my changes - [x] Any dependent changes have been merged and published **Due date for this PR is 9/15/2023.** If this PR is not merged into `develop` by this date, the PR will be closed and returned to the developer. Co-authored-by: Sho Yokota --- src/gsi/read_dbz_nc.f90 | 16 ++++++++++------ src/gsi/setupdbz.f90 | 15 +++++++++------ 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index f6ac9aa112..7f8604b9d2 100644 --- a/src/gsi/read_dbz_nc.f90 +++ b/src/gsi/read_dbz_nc.f90 @@ -69,11 +69,12 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no use kinds, only: r_kind,r_double,i_kind,r_single use constants, only: zero,half,one,two,deg2rad,rad2deg, & one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, & - eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening + eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,r_missing use gridmod, only: tll2xy,nsig,nlat,nlon use obsmod, only: iadate,doradaroneob,oneoblat,oneoblon,oneobheight, & mintiltdbz,maxtiltdbz,minobrangedbz,maxobrangedbz,& static_gsi_nopcp_dbz,rmesh_dbz,zmesh_dbz + use gsi_4dvar, only: iwinbgn use hybrid_ensemble_parameters,only : l_hyb_ens use obsmod,only: radar_no_thinning,missing_to_nopcp use convinfo, only: nconvtype,ctwind,icuse,ioctype @@ -147,7 +148,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no integer(i_kind) :: maxobs,nchanl,ilat,ilon,scount real(r_kind) :: thistiltr,thisrange,this_stahgt,thishgt - real(r_kind) :: thisazimuthr,t4dv, & + real(r_kind) :: thisazimuthr, & dlat,dlon,thiserr,thislon,thislat, & timeb real(r_kind) :: radartwindow @@ -337,6 +338,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978 rmins_an=mins_an !convert to real number + timeb=real(mins_an-iwinbgn,r_kind) !assume all observations are at the analysis time ivar = 1 @@ -453,7 +455,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no ntmp=ndata ! counting moved to map3gridS - timedif=abs(t4dv) !don't know about this + timedif=zero ! assume all observations are at the analysis time crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,& @@ -481,7 +483,10 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no !!end modified for thinning - thisazimuthr=0.0_r_kind + thisazimuthr=r_missing + thistiltr=r_missing + this_stahgt=r_missing + thisrange=r_missing this_staid=radid !Via equivalence in declaration, value is propagated ! to rstation_id used below. cdata_all(1,iout) = thiserr ! reflectivity obs error (dB) - inflated/adjusted @@ -491,7 +496,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no cdata_all(5,iout) = dbzQC(i,j,k) ! radar reflectivity factor cdata_all(6,iout) = thisazimuthr ! 90deg-azimuth angle (radians) - cdata_all(7,iout) = timeb*r60inv ! obs time (analyis relative hour) + cdata_all(7,iout) = timeb*r60inv ! obs time (relative hour from beginning of the DA window) cdata_all(8,iout) = ikx ! type cdata_all(9,iout) = thistiltr ! tilt angle (radians) cdata_all(10,iout)= this_stahgt ! station elevation (m) @@ -521,7 +526,6 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no !---all looping done now print diagnostic output write(6,*)'READ_dBZ: Reached eof on radar reflectivity file' - write(6,*)'READ_dBZ: # volumes in input file =',nvol write(6,*)'READ_dBZ: # read in obs. number =',nread write(6,*)'READ_dBZ: # elevations outside time window =',numbadtime write(6,*)'READ_dBZ: # of noise obs to no precip obs =',num_nopcp diff --git a/src/gsi/setupdbz.f90 b/src/gsi/setupdbz.f90 index 068842cd6b..453c4a5f8d 100644 --- a/src/gsi/setupdbz.f90 +++ b/src/gsi/setupdbz.f90 @@ -590,14 +590,17 @@ subroutine setupdbz(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,radardbz_d ! Compute observation pressure (only used for diagnostics) dz = zges(k2)-zges(k1) dlnp = prsltmp(k2)-prsltmp(k1) - pobl = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1)) - - presw = ten*exp(pobl) - if ( l_use_dbz_directDA ) then - presq = presw + pobl = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1)) + presw = ten*exp(pobl) + presq = presw else - if( (k1 == k2) .and. (k1 == 1) ) presw=ten*exp(prsltmp(k1)) + if( (k1 == k2) .and. (k1 == 1) ) then + presw = ten*exp(prsltmp(k1)) + else + pobl = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1)) + presw = ten*exp(pobl) + end if end if ! solution to Nan in some members only for EnKF which causes problem?