PROGRAM Prep_Hourly_Precip_Input ! ! PURPOSE: Perform long-term budget adjustment on merged Stage 2/4 hourly ! analysis; fill in missing area with NAM hourly forecast. ! ! SUMMARY: ! 1. Use a "water budget history" (cumulative differences between ! hourly precip input and verifying daily gauge analysis) file ! to adjust the analysis. ! 2. Use NAM hourly precipitation forecasts (from the 00Z cycle's ! 12-36h forecast) to fill in missing data ! 3. Output the adjusted precipitation analyses in GRIB. Also write it out ! in simple binary (2005: difficult to get WRF to read GRIB). ! ! INPUT FILES: ! Unit 11 Merged Stage 2/4 analysis on B-grid ! Unit 12 NAM hourly forecast ! Unit 13 Long-term budget history array ! ! OUTPUT FILES: ! Unit 51 Adjusted, filled (O'ConUS and other no-obs area) hourly precip, ! ready to be used by NDAS ! Unit 52 Grib'd version of Unit 51 file ! ! RECORD OF REVISIONS ! Date Programmer Description of Change ! ========== ========== ===================================================== ! 2001-06-11 Ying Lin Original implementation using Stage 2 analysis. ! Mapping to model grid is done using ipolate within ! this program. ! ! 2002-11-15 Y. Lin Modifed to use both Stage II and IV ! 2002-12-23 Y. Lin Use RFC domain mask to exclude NWRFC data in Stage 4 ! domain, because the regional hourly analysis from ! NWRFC has been problematic. ! ! 2003-04-15 Y. Lin Use RFC domain mask to exclude CNRFC data in Stage 4 ! domain. CNRFC seldom sends hourly data, so the part ! of the CNRFC grid that overlaps with the NWRFC grid ! usually carries the (often) faulty NWRFC value. ! ! 2003-06-09 Y. Lin Use RFC domain mask to exclude Stage IV data points ! that ! 1) are not in any RFC's domain proper, and ! 2) have a zero value ! These are generally points in Canada or over water ! (though some Canadian points are included in one or ! more RFC domains). Reason: some RFC's (e.g. SERFC) ! do not mark these points as 'missing', even when they ! are outside of their range, so they end up with ! values of zero, when they really should be designated ! as "missing". We'll use Stage 2 instead of Stage 4 ! for these points. ! ! 2003-06-09 Y. Lin Added adjustment using cumulative h2o budget history ! ! 2005-03-16 Y. Lin Changes made for WRF-NMM: read in model grid number ! from unit 5, so a single program can be used for ! multiple grids. ! ! 2007-01-30 Y. Lin Fill in missing data (usually O'Conus) with hourly ! precip forecast from a 00Z NAM (12-36h forecast, ! broken down into individual files each containing ! one hour) ! ! 2009-04-16 Y. Lin Take out the merging of the Stage 2/4 (spin off to ! "merge2n4.f") and the ipolate mapping to model grid ! (now handled by copygb) ! INPUT FILES: ! Unit 11 Merged Stage 2/4 analysis on model B-grid ! Unit 12 H2O budget history (read into array 'balance') ! Unit 13 NAM precip forecast (from a 00Z cycle) for this hour ! ! OUTPUT FILES: ! Unit 51 Adjusted/filled precipitation analysis, simple binary ! Unit 52 Adjusted precipitation analysis, GRIB ! ! LANGUAGE: Fortran 90/95 ! IMPLICIT NONE ! Maximum 1-d dimension of model grid: INTEGER, PARAMETER :: Lmax=2800000 ! test point: INTEGER, PARAMETER :: i0=194023 ! Merged Stage 2/4 analysis on B-grid (Panl); NAM 1-hour pcp forecast (Pnam); ! final hourly product ready for use by NDAS (P); long-term precipitation ! budget array (balance) REAL, DIMENSION(Lmax) :: Panl, Pnam, P, balance ! Beginning and ending date of the precip budget balance array (for info only, ! Not used in calculation) INTEGER :: iday1, iday2 ! For read status on budget history file: INTEGER :: istat ! Increment (=-balance/12.) added to Panl REAL :: pinc ! Bitmap for Panl (BIT; also used later for Pout); for Pnam (BITnam). BITnam ! is a dummy array, since NAM native grid precip output does not include a ! bitmap. LOGICAL(KIND=1), DIMENSION(Lmax) :: BIT, BITnam ! GDS and PDS for 1) getgb ('jgds, jpds') ! 2) merged ST2/4 analysis ('kgdsa, kpdsa') ! 3) NAM 1-hour precip fcst amount ('kgdsn, kpdsn') ! 5) Final product (kpds, kgds) INTEGER, DIMENSION(200) :: jgds, kgdsa, kgdsn, kgds INTEGER, DIMENSION(200) :: jpds, kpdsa, kpdsn, kpds ! For do loop index: INTEGER :: i ! Misc for baopenr and getgb. kfa and kfn are the lengths of the merged ST2/4 ! array and the NAM 1-hour precip array. They should be the same length (the ! size of the B-grid). INTEGER :: kfa, kfn, k, iret, iretanl, iretnam ! For 'getgb' searches: jpds = -1 ! Read in merged Stage 2/4 analysis on B-grid CALL BAOPENR(11,'fort.11',iret) write(*,*) 'baopenr on unit 11, iret=', iret CALL GETGB(11,0,Lmax,0,jpds,jgds,kfa,k,kpdsa,kgdsa,BIT,Panl,iretanl) write(*,10) 'st2n4', kpdsa(21)-1,kpdsa(8),kpdsa(9), kpdsa(10),kpdsa(11), & & iretanl, kfa 10 format('GETGB for ',a5,x,5i2.2, ' iret=', i2, ' kf=',i8) ! Read in NAM precip forecast for this hour: CALL BAOPENR(12,'fort.12',iret) write(*,*) 'baopenr on unit 12, iret=', iret CALL GETGB(12,0,Lmax,0,jpds,jgds,kfn,k,kpdsn,kgdsn,BITnam,Pnam,iretnam) write(*,10) 'NAM1h', kpdsn(21)-1,kpdsn(8),kpdsn(9),kpdsn(10),kpdsn(11), & & iretnam, kfn ! Read in cumulative precip budget balance history file. If the file does not ! exist, or if the end of file is reached, assume perfect balance ! ("INQUIRE(UNIT=13,EXIST=bexist,ERR=50)" does not work). READ(13,IOSTAT=istat) iday1,iday2,(balance(i),i=1,kfa) IF( istat == 0 ) then WRITE(*,*) 'Precip budget balance is from', iday1,' to', iday2 else WRITE(*,*) 'Error reading precip budget history, make no adjustment.' balance=0. endif ! IF (iretanl == 0 ) THEN ! If the merged ST2/4 and NAM 1-hour forecast exists, adjust the values ! using the long-term budget balance array (if available), then fill in ! with NAM hourly forecast (if available). kgds=kgdsa kpds=kpdsa DO i = 1, kfa IF (BIT(i)) THEN pinc = -balance(i)/12. pinc = min(pinc, 0.2*Panl(i)) pinc = max(pinc,-0.2*Panl(i)) P(i) = Panl(i) + pinc ELSE IF (iretnam == 0 ) THEN BIT(i) = .TRUE. P(i) = Pnam(i) ELSE P(i) = 999. BIT(i) = .FALSE. END IF END DO ELSE IF (iretanl /= 0 .and. iretnam == 0) then ! The merged Stage 2/4 file does not exist, but there is a valid NAM 1hr ! forecast. Use NAM 1hr forecast. write(6,*) 'st2n4 file does not exist, use model hourly only.' P = Pnam kgds=kgdsn kpds=kpdsn kfa = kfn ELSE ! Neither merged Stage 2/4 nor the NAM 1-hour precip file exists. Make up ! a dummy file. No use specifying bitmap or PDS/GDS: we won't know what ! the GDS would be (since neither B-grid file exists) and we want to keep ! this code non-grid-specific write(6,*) 'st2n4 and model 1hr fcst both missing, set input to 999.' ! Set kfa to Lmax to for binary output: kfa = Lmax ENDIF ! Output precip in binary format: write(51) (P(i), i=1, kfa) ! NAM native-grid precip does not come with a bitmap. For some reason, when ! it is putgb'd with a bitmap (all .true., and with kpds(4)=192), things go ! awry. This is the case for both the B-grid and E-grid files. For now, ! Only output bitmap if NAM precip is not involved. ! IF (iretnam == 0) THEN kpds(4) = 128 END IF ! Output a GRIB version, if either the merged ST2/4 or the NAM 1-hourly existed IF (iretanl == 0 .OR. iretnam == 0) THEN CALL BAOPEN(52,'fort.52',iret) write(*,*) 'BAOPEN on unit 52, iret=', iret CALL PUTGB(52,kfa,kpds,kgds,BIT,P,iret) WRITE(*,*) 'FINISHED PUTGB for processed precip, iret=', iret ENDIF STOP END PROGRAM Prep_Hourly_Precip_Input