!> @file !! @brief Read and qc afwa, nesdis/ims and autosnow snow data. !! @author gayno org: w/np2 @date 2005-dec-16 !> Read and qc afwa, nesdis/ims and autosnow snow data. !! !! program history log: !! - 2005-dec-16 gayno - initial version !! - 2007-aug-10 gayno - Allow program to run with no nesdis/ims data. !! Add 16th mesh afwa grib data. !! - 2008-feb-04 gayno - Add autosnow cover data for sh. !! - 2009-jun-03 gayno - Add qc check for nesdis/ims and afwa data. !! - 2014-feb-07 gayno - Read nesdis/ims data in grib1 or grib 2 !! format. !! - 2014-sep-30 gayno - Convert weekly nh snow climatology - used to !! qc input data - to grib 2. !! variable definitions: !! - bad_afwa_Xh - is afwa data corrupt? !! - bitmap_afwa_Xh - bitmap of afwa grid (false-non land, true-land) !! - kgds_afwa_Xh - afwa grid description section (grib section 2) !! - nesdis_res - resolution of nesdis/ims data in km !! - snow_dep_afwa_Xh - afwa snow depth data (inches*10 on input, !! meters on output) !! - use_xh_afwa - true if afwa data to be used !! module snowdat use program_setup, only : autosnow_file, & nesdis_snow_file, & nesdis_lsmask_file, & afwa_snow_global_file, & afwa_snow_nh_file, & afwa_snow_sh_file, & afwa_lsmask_nh_file, & afwa_lsmask_sh_file use model_grid, only : imdl, & jmdl integer :: iafwa !< i-dimension of afwa grid integer :: iautosnow !< i-dimension of autosnow grid integer :: inesdis !< i-dimension of nesdis grid integer :: jafwa !< j-dimension of afwa grid integer :: jautosnow !< j-dimension of autosnow grid integer :: jnesdis !< j-dimension of nesdis grid integer :: kgds_afwa_global(200) !< grib1 grid description section for !! global afwa data. integer :: kgds_afwa_nh(200) !< grib1 grid description section for northern !! hemisphere 16th mesh afwa data. integer :: kgds_afwa_nh_8th(200) !< grib1 grid description section for !! northern hemisphere 8th mesh afwa data. integer :: kgds_afwa_sh(200) !< grib1 grid description section for southern !! hemisphere 16th mesh afwa data. integer :: kgds_afwa_sh_8th(200) !< grib1 grid description section for !! southern hemisphere 8th mesh afwa data. integer :: kgds_autosnow(200) !< autosnow grid description section (grib section 2) integer :: kgds_nesdis(200) !< nesdis/ims grid description section (grib section 2) integer :: mesh_nesdis !< nesdis/ims data is 96th mesh (or bediant) integer*1, allocatable :: sea_ice_nesdis(:,:) !< nesdis/ims sea ice flag (0-open water, 1-ice) logical :: bad_afwa_nh !< When true, the northern hemisphere afwa data failed !! its quality control check. logical :: bad_afwa_sh !< When true, the southern hemisphere afwa data failed !! its quality control check. logical :: bad_nesdis !< When true, the nesdis data failed its quality !! control check. logical :: bad_afwa_global !< When true, the global afwa data failed its quality !! control check. logical*1, allocatable :: bitmap_afwa_global(:,:) !< The global afwa data grib bitmap. !!(false-non land, true-land). logical*1, allocatable :: bitmap_afwa_nh(:,:) !< The northern hemisphere afwa data grib bitmap. !! (false-non land, true-land). logical*1, allocatable :: bitmap_afwa_sh(:,:) !< The southern hemisphere afwa data grib bitmap. !! (false-non land, true-land). logical*1, allocatable :: bitmap_nesdis(:,:) !< nesdis data grib bitmap (false-non land, true-land). logical*1, allocatable :: bitmap_autosnow(:,:) !< autosnow data grib bitmap (false-non land, !! true-land). logical :: use_nh_afwa !< True if northern hemisphere afwa data to be used. logical :: use_sh_afwa !< True if southern hemisphere afwa data to be used. logical :: use_global_afwa !< True if global hemisphere afwa data to be used. logical :: use_autosnow !< True if autosnow data to be used logical :: use_nesdis !< True if nesdis/ims data to be used real :: autosnow_res !< Resolution of autosnow in km real :: afwa_res !< Resolution of afwa data in km real :: nesdis_res !< Resolution of the nesdis data in km. real, allocatable :: snow_cvr_nesdis(:,:) !< nesdis/ims snow cover flag (0-no, 100-yes) real, allocatable :: snow_cvr_autosnow(:,:) !< autosnow snow cover flag (0-no, 100-yes) real, allocatable :: snow_dep_afwa_global(:,:) !< The global afwa snow depth. real, allocatable :: snow_dep_afwa_nh(:,:) !< Northern hemisphere afwa snow depth. real, allocatable :: snow_dep_afwa_sh(:,:) !< Southern hemisphere afwa snow depth. ! the afwa 8th mesh binary data has no grib header, so set it from these ! data statements. needed for ipolates routines. data kgds_afwa_nh_8th/5,2*512,-20826,145000,8,-80000,2*47625,0, & 9*0,255,180*0/ data kgds_afwa_sh_8th/5,2*512,20826,-125000,8,-80000,2*47625,128, & 9*0,255,180*0/ contains !> Read autosnow snow cover. !! !! program history log: !! 2008-feb-04 gayno - initial version !! !! files: !! input: !! - autosnow data, grib 2, unit=lugb !! !! condition codes: all fatal !! 74 - bad open of autosnow file !! 75 - bad read of autosnow file !! !! @note Autosnow data is available only for southern hemis. !! Autosnow data is in grib 2. !! !! @author George Gayno org: w/np2 @date 2008-Feb-04 subroutine readautosnow use grib_mod ! grib 2 libraries implicit none type(gribfield) :: gfld integer :: iret, j, k, lugb, lugi integer :: jdisc, jgdtn, jpdtn integer :: jids(200), jgdt(200), jpdt(200) logical :: unpack use_autosnow = .true. if ( len_trim(autosnow_file) == 0 ) then print*,"- WILL NOT USE AUTOSNOW DATA." use_autosnow = .false. return end if print*,"- OPEN AND READ AUTOSNOW FILE ", trim(autosnow_file) lugb=12 call baopenr(lugb,autosnow_file,iret) if (iret /= 0) then print*,'- FATAL ERROR: BAD OPEN OF FILE, IRET IS ', iret call w3tage('SNOW2MDL') call errexit(74) endif call grib2_null(gfld) j = 0 ! search at beginning of file lugi = 0 ! no grib index file jdisc = 0 ! search for discipline; 0 - meteorological products jpdtn = 30 ! search for product definition template number; 30 - satellite product jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid jids = -9999 ! array of values in identification section, set to wildcard jgdt = -9999 ! array of values in grid definiation template 3.m jpdt = -9999 ! array of values in product definition template 4.n jpdt(1) = 1 ! search for parameter category - moisture jpdt(2) = 42 ! search for parameter number - snow cover in percent. unpack = .true. ! unpack data call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, iret) if (iret /=0) then print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret call w3tage('SNOW2MDL') call errexit(75) endif print*,"- DATA VALID AT (YYYYMMDDHH): ", gfld%idsect(6),gfld%idsect(7), & gfld%idsect(8),gfld%idsect(9) call baclose (lugb, iret) !----------------------------------------------------------------------- ! set the grib1 kgds array from the g2 grid definition template array. ! the kgds array is used by ipolates. !----------------------------------------------------------------------- call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds_autosnow, & iautosnow, jautosnow, autosnow_res) allocate (bitmap_autosnow(iautosnow,jautosnow)) bitmap_autosnow = reshape (gfld%bmap , (/iautosnow,jautosnow/) ) allocate (snow_cvr_autosnow(iautosnow,jautosnow)) snow_cvr_autosnow = reshape (gfld%fld , (/iautosnow,jautosnow/) ) call grib2_free(gfld) end subroutine readautosnow !> Read nesdis/ims snow cover/ice data. !! !! program history log: !! 2005-dec-16 gayno - initial version !! 2014-feb-07 gayno - Read 4km ims data in either !! grib1 or grib 2 format. !! files: !! input: !! - ims snow cover and ice file, grib 1 or grib 2 !! - 16th-mesh ims land mask, binary !! !! condition codes: all fatal !! 41 - ims file not grib 1 or grib 2 !! 53 - ims data failed quality check !! 70 - bad read of ims snow cover data !! 71 - bad read of ims ice data !! 72 - bad read of ims grib 1 header !! 73 - bad open of ims file !! 87 - bad open ims land mask file !! 88 - bad read ims land mask file !! !! @note Nesdis/ims data available only for n hemis. Ims data used !! to be created by nesdis, hence the references to "nesdis" !! in this routine. Ims data is now created by the national !! ice center. !! !! @author George Gayno org: w/np2 @date 2005-Dec-16 subroutine readnesdis use grib_mod implicit none integer, parameter :: iunit = 13 ! input grib file unit number integer, parameter :: iunit2 = 43 ! input landmask file unit number integer*4, allocatable :: dummy4(:,:) integer :: i, j integer :: iret integer :: jgds(200) integer :: jpds(200) integer :: lskip integer, parameter :: lugi = 0 ! grib index file unit number - not used integer :: jdisc, jgdtn, jpdtn, k integer :: jids(200), jgdt(200), jpdt(200) integer :: kgds(200) integer :: kpds(200) integer :: message_num integer :: numbytes integer :: numpts integer :: isgrib logical :: unpack real, allocatable :: dummy(:,:) real :: dum type(gribfield) :: gfld use_nesdis = .true. if ( len_trim(nesdis_snow_file) == 0 ) then print*,"- WILL NOT USE NESDIS/IMS DATA." use_nesdis = .false. return end if print*,"- OPEN AND READ NESDIS/IMS SNOW FILE ", trim(nesdis_snow_file) call grib_check(nesdis_snow_file, isgrib) if (isgrib==0) then print*,'- FATAL ERROR: IMS FILE MUST BE GRIB 1 OR GRIB2 FORMAT' call w3tage('SNOW2MDL') call errexit(41) end if call baopenr (iunit, nesdis_snow_file, iret) if (iret /= 0) then print*,'- FATAL ERROR: BAD OPEN OF FILE, IRET IS ', iret call w3tage('SNOW2MDL') call errexit(73) end if if (isgrib==1) then ! grib 1 format !----------------------------------------------------------------------- ! tell degribber to look for requested data. !----------------------------------------------------------------------- lskip = -1 jpds = -1 jgds = -1 jpds(5) = 91 ! ice cover kpds = jpds kgds = jgds print*,"- GET GRIB HEADER" call getgbh(iunit, lugi, lskip, jpds, jgds, numbytes, & numpts, message_num, kpds, kgds, iret) if (iret /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF HEADER. IRET IS ", iret call w3tage('SNOW2MDL') call errexit(72) end if kgds_nesdis = kgds inesdis = kgds(2) jnesdis = kgds(3) mesh_nesdis = inesdis / 64 nesdis_res = 381. / float(mesh_nesdis) ! in km print*,"- DATA VALID AT (YYMMDDHH): ", kpds(8:11) allocate (dummy(inesdis,jnesdis)) allocate (sea_ice_nesdis(inesdis,jnesdis)) allocate (bitmap_nesdis(inesdis,jnesdis)) print*,"- DEGRIB SEA ICE." call getgb(iunit, lugi, (inesdis*jnesdis), lskip, jpds, jgds, & numpts, lskip, kpds, kgds, bitmap_nesdis, dummy, iret) if (iret /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ", iret call w3tage('SNOW2MDL') call errexit(71) end if sea_ice_nesdis = nint(dummy) ! only needed as yes/no flag deallocate (dummy) lskip = -1 jpds = -1 jgds = -1 jpds(5) = 238 ! snow cover kpds = jpds kgds = jgds allocate (snow_cvr_nesdis(inesdis,jnesdis)) print*,"- DEGRIB SNOW COVER." call getgb(iunit, lugi, (inesdis*jnesdis), lskip, jpds, jgds, & numpts, lskip, kpds, kgds, bitmap_nesdis, snow_cvr_nesdis, iret) if (iret /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ", iret call w3tage('SNOW2MDL') call errexit(70) end if elseif (isgrib==2) then ! grib 2 format print*,"- DEGRIB SNOW COVER." j = 0 ! search at beginning of file jdisc = 0 ! search for discipline; 0 - meteorological products jpdtn = 0 ! search for product definition template number; 0 - analysis at one level jgdtn = 20 ! search for grid definition template number; 20 - polar stereographic grid jids = -9999 ! array of values in identification section, set to wildcard jgdt = -9999 ! array of values in grid definition template 3.m jpdt = -9999 ! array of values in product definition template 4.n jpdt(1) = 1 ! search for parameter category - moisture jpdt(2) = 201 ! search for parameter number - snow cover in percent. unpack = .true. ! unpack data call grib2_null(gfld) call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, iret) if (iret /=0) then print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret call w3tage('SNOW2MDL') call errexit(70) endif print*,"- DATA VALID AT (YYYYMMDDHH): ", gfld%idsect(6),gfld%idsect(7), & gfld%idsect(8),gfld%idsect(9) !----------------------------------------------------------------------- ! set the grib1 kgds array from the g2 grid definition template array. ! the kgds array is used by ipolates. !----------------------------------------------------------------------- call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds_nesdis, & inesdis, jnesdis, dum) mesh_nesdis = inesdis / 64 nesdis_res = 381. / float(mesh_nesdis) ! in km if (mesh_nesdis==16) kgds_nesdis(6)=136 ! the ims 16th mesh grib2 data ! is gribbed with an elliptical ! earth. that is wrong. hardwire ! a fix here. allocate (snow_cvr_nesdis(inesdis,jnesdis)) allocate (sea_ice_nesdis(inesdis,jnesdis)) allocate (bitmap_nesdis(inesdis,jnesdis)) bitmap_nesdis = reshape (gfld%bmap , (/inesdis,jnesdis/) ) snow_cvr_nesdis = reshape (gfld%fld , (/inesdis,jnesdis/) ) call grib2_free(gfld) print*,"- DEGRIB SEA ICE." j = 0 ! search at beginning of file jdisc = 10 ! search for discipline; 10 - ocean products jpdtn = 0 ! search for product definition template number; 0 - analysis at one level jgdtn = 20 ! search for grid definition template number; 20 - polar stereographic grid jids = -9999 ! array of values in identification section, set to wildcard jgdt = -9999 ! array of values in grid definition template 3.m jpdt = -9999 ! array of values in product definition template 4.n jpdt(1) = 2 ! search for parameter category - ice jpdt(2) = 0 ! search for parameter number - ice cover in percent. unpack = .true. ! unpack data call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, iret) if (iret /=0) then print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret call w3tage('SNOW2MDL') call errexit(71) endif sea_ice_nesdis = reshape (gfld%fld , (/inesdis,jnesdis/) ) call grib2_free(gfld) end if call baclose(iunit,iret) !----------------------------------------------------------------------- ! the 16th mesh nesdis/ims grib data does not have a proper ! bitmap section. therefore, need to read in the mask ! from file. but the 96th mesh data has a proper bitmap, so use it. !----------------------------------------------------------------------- if (mesh_nesdis == 16) then print*,"- OPEN NESDIS/IMS 16TH MESH LAND MASK: ", trim(nesdis_lsmask_file) open(iunit2, file=trim(nesdis_lsmask_file), form="formatted", & iostat = iret) if (iret /= 0) then print*,"- FATAL ERROR OPENING NESDIS/IMS LAND MASK FILE. ISTAT IS: ", iret call errexit(87) end if print*,"- READ NESDIS/IMS 16TH MESH LAND MASK." allocate (dummy4(inesdis,jnesdis)) do j = 1, 1024 read(iunit2, 123, iostat=iret) (dummy4(i,j),i=1,1024) if (iret /= 0) then print*,"- FATAL ERROR READING NESDIS/IMS LAND MASK FILE. ISTAT IS: ", iret call errexit(88) end if enddo close (iunit2) !----------------------------------------------------------------------- ! the file has 0-sea, 1-land, 9-off hemi. this code expects ! 0-non-land (or don't use data), 1-land (use data). !----------------------------------------------------------------------- bitmap_nesdis=.false. do j = 1, 1024 do i = 1, 1024 if (dummy4(i,j) == 1) bitmap_nesdis(i,j) = .true. enddo enddo deallocate(dummy4) 123 FORMAT(80I1) endif ! is nesdis/ims data 16th mesh? bad_nesdis=.false. call nh_climo_check(kgds_nesdis,snow_cvr_nesdis,bitmap_nesdis,inesdis,jnesdis,2,bad_nesdis) !----------------------------------------------------------------------- ! for the 2009 nmm-b implementation, it was decided to not run with ! afwa only. so even if afwa data is current and not corrupt, ! but the ims is bad, then abort program. exception, if ims is very old ! (there is a catastropic outage) then program will run with afwa ! only. this is done by setting the nesdis_snow_file variable to ! a zero length string (i.e., ims data not selected). this variable ! setting is accomplished in the run script. !----------------------------------------------------------------------- if (bad_nesdis) then print*,'- FATAL ERROR: NESDIS/IMS DATA BAD, DO NOT USE.' print*,'- DONT RUN PROGRAM.' use_nesdis=.false. call w3tage('SNOW2MDL') call errexit(53) stop endif return end subroutine readnesdis !> Read snow depth data and masks. !! !! @note Read afwa snow depth data and !! land sea mask. !! !! program history log: !! !! 2005-dec-16 gayno - initial version !! 2007-nov-28 gayno - read 16th mesh afwa data in grib format !! !! files: !! input: !! - global afwa data in grib 2 (if selected) !! - nh afwa data in grib 1 (if selected) !! - sh afwa data in grib 1 (if selected) !! !! condition codes: !! 60 - bad open afwa file !! 61 - bad degrib of afwa file !! !! @author George Gayno org: w/np2 @date 2005-Dec-16 subroutine readafwa use grib_mod implicit none integer, parameter :: iunit=17 integer :: jgds(200), jpds(200), kgds(200), kpds(200) integer :: istat, isgrib integer :: lugi, lskip, numbytes, numpts, message_num integer :: j, k, jdisc, jpdtn, jgdtn integer :: jpdt(200), jgdt(200), jids(200) logical :: unpack type(gribfield) :: gfld bad_afwa_nh=.false. bad_afwa_sh=.false. bad_afwa_global=.false. use_global_afwa=.true. use_nh_afwa = .true. use_sh_afwa = .true. if (len_trim(afwa_snow_nh_file) == 0 .and. & len_trim(afwa_snow_sh_file) == 0 .and. & len_trim(afwa_snow_global_file) == 0) then print*,"- WILL NOT USE AFWA DATA." use_nh_afwa = .false. use_sh_afwa = .false. use_global_afwa = .false. return end if !----------------------------------------------------------------------- ! If chosen, read global AFWA GRIB2 file. !----------------------------------------------------------------------- if ( len_trim(afwa_snow_global_file) > 0 ) then print*,"- OPEN AND READ global AFWA SNOW FILE ", trim(afwa_snow_global_file) call baopenr (iunit, afwa_snow_global_file, istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD OPEN OF FILE, ISTAT IS ', istat call w3tage('SNOW2MDL') call errexit(60) end if call grib2_null(gfld) jdisc = 0 ! Search for discipline; 0 - meteorological products j = 0 ! Search at beginning of file. lugi = 0 ! No grib index file. jids = -9999 ! Identification section, set to wildcard. jgdt = -9999 ! Grid definition template, set to wildcard. jgdtn = -1 ! Grid definition template number, set to wildcard. jpdtn = 0 ! Search for product definition template number 0 - analysis or forecast jpdt = -9999 ! Product definition template (Sec 4), initialize to wildcard. jpdt(1) = 1 ! Search for parameter category 1 (Sec 4 oct 10) - ! moisture. jpdt(2) = 11 ! Search for parameter 11 (Sec 4 oct 11) - snow depth. unpack = .true. ! Unpack data. call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, istat) if (istat /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF GLOBAL DATA. ISTAT IS ", istat call w3tage('SNOW2MDL') call errexit(61) end if print*,"- DATA VALID AT (YYMMDDHH): ", gfld%idsect(6:9) print*,"- DEGRIB SNOW DEPTH." call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds_afwa_global, & iafwa, jafwa, afwa_res) allocate(bitmap_afwa_global(iafwa,jafwa)) allocate(snow_dep_afwa_global(iafwa,jafwa)) snow_dep_afwa_global = reshape(gfld%fld, (/iafwa,jafwa/)) bitmap_afwa_global = reshape(gfld%bmap, (/iafwa,jafwa/)) call baclose(iunit, istat) call nh_climo_check(kgds_afwa_global,snow_dep_afwa_global,bitmap_afwa_global,iafwa,jafwa,1,bad_afwa_global) if (bad_afwa_global) then print*,'- WARNING: AFWA DATA BAD, DO NOT USE.' use_global_afwa = .false. endif use_nh_afwa=.false. ! Use global or hemispheric files. not both. use_sh_afwa=.false. return ! Use global or hemispheric files. not both. else use_global_afwa=.false. endif if ( len_trim(afwa_snow_nh_file) > 0 ) then ! afwa nh data selected call grib_check(afwa_snow_nh_file, isgrib) if (isgrib==0) then ! old ncep binary format iafwa = 512 jafwa = 512 afwa_res = 47.625 ! in kilometers kgds_afwa_nh = kgds_afwa_nh_8th allocate (snow_dep_afwa_nh(iafwa,jafwa)) call read_afwa_binary(afwa_snow_nh_file, snow_dep_afwa_nh) allocate (bitmap_afwa_nh(iafwa,jafwa)) call read_afwa_mask(afwa_lsmask_nh_file, bitmap_afwa_nh) else ! afwa data is grib print*,"- OPEN AND READ AFWA SNOW FILE ", trim(afwa_snow_nh_file) call baopenr (iunit, afwa_snow_nh_file, istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD OPEN OF FILE, ISTAT IS ', istat call w3tage('SNOW2MDL') call errexit(60) end if !----------------------------------------------------------------------- ! tell degribber to look for requested data. !----------------------------------------------------------------------- lugi = 0 lskip = -1 jpds = -1 jgds = -1 jpds(5) = 66 ! snow depth kpds = jpds kgds = jgds print*,"- GET GRIB HEADER" call getgbh(iunit, lugi, lskip, jpds, jgds, numbytes, & numpts, message_num, kpds, kgds, istat) if (istat /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF HEADER. ISTAT IS ", istat call w3tage('SNOW2MDL') call errexit(61) end if iafwa = kgds(2) jafwa = kgds(3) afwa_res = float(kgds(8))*0.001 ! in km. print*,"- DATA VALID AT (YYMMDDHH): ", kpds(8:11) print*,"- DEGRIB SNOW DEPTH." allocate(bitmap_afwa_nh(iafwa,jafwa)) allocate(snow_dep_afwa_nh(iafwa,jafwa)) call getgb(iunit, lugi, (iafwa*jafwa), lskip, jpds, jgds, & numpts, lskip, kpds, kgds, bitmap_afwa_nh, snow_dep_afwa_nh, istat) if (istat /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF DATA. ISTAT IS ", istat call w3tage('SNOW2MDL') call errexit(61) end if kgds_afwa_nh = kgds kgds_afwa_nh(7) = -80000 ! ipolates definition of orientation angle is ! 180 degrees off from grib standard. call baclose(iunit, istat) endif ! is nh afwa data grib? call nh_climo_check(kgds_afwa_nh,snow_dep_afwa_nh,bitmap_afwa_nh,iafwa,jafwa,1,bad_afwa_nh) else use_nh_afwa=.false. endif !----------------------------------------------------------------------- ! now, read southern hemisphere data. !----------------------------------------------------------------------- if ( len_trim(afwa_snow_sh_file) > 0 ) then call grib_check(afwa_snow_sh_file, isgrib) if (isgrib==0) then ! old ncep binary format iafwa = 512 jafwa = 512 afwa_res = 47.625 kgds_afwa_sh = kgds_afwa_sh_8th allocate (snow_dep_afwa_sh(iafwa,jafwa)) call read_afwa_binary(afwa_snow_sh_file, snow_dep_afwa_sh) allocate (bitmap_afwa_sh(iafwa,jafwa)) call read_afwa_mask(afwa_lsmask_sh_file, bitmap_afwa_sh) else ! sh afwa data is grib print*,"- OPEN AND READ AFWA SNOW FILE ", trim(afwa_snow_sh_file) call baopenr (iunit, afwa_snow_sh_file, istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD OPEN OF FILE, ISTAT IS ', istat call w3tage('SNOW2MDL') call errexit(60) end if !----------------------------------------------------------------------- ! tell degribber to look for requested data. !----------------------------------------------------------------------- lugi = 0 lskip = -1 jpds = -1 jgds = -1 jpds(5) = 66 ! snow cover kpds = jpds kgds = jgds print*,"- GET GRIB HEADER" call getgbh(iunit, lugi, lskip, jpds, jgds, numbytes, & numpts, message_num, kpds, kgds, istat) if (istat /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF HEADER. ISTAT IS ", istat call w3tage('SNOW2MDL') call errexit(61) end if iafwa = kgds(2) jafwa = kgds(3) afwa_res = float(kgds(8))*0.001 ! in km. print*,"- DATA VALID AT (YYMMDDHH): ", kpds(8:11) print*,"- DEGRIB SNOW DEPTH." allocate(bitmap_afwa_sh(iafwa,jafwa)) allocate(snow_dep_afwa_sh(iafwa,jafwa)) call getgb(iunit, lugi, (iafwa*jafwa), lskip, jpds, jgds, & numpts, lskip, kpds, kgds, bitmap_afwa_sh, snow_dep_afwa_sh, istat) if (istat /= 0) then print*,"- FATAL ERROR: BAD DEGRIB OF DATA. ISTAT IS ", istat call w3tage('SNOW2MDL') call errexit(61) end if kgds_afwa_sh = kgds kgds_afwa_sh(7) = -80000 ! ipolates definition of orientation angle is ! 180 degrees off from grib standard. call baclose(iunit, istat) endif ! is sh afwa data grib or not? call afwa_check(2) else use_sh_afwa = .false. endif !------------------------------------------------------------------- !if either hemisphere is bad, don't trust all hemispheres !------------------------------------------------------------------- if (bad_afwa_nh .or. bad_afwa_sh) then print*,'- WARNING: AFWA DATA BAD, DO NOT USE.' use_nh_afwa = .false. use_sh_afwa = .false. endif return end subroutine readafwa !> Check for corrupt nh snow cover data. !! !! @note Check for corrupt nh data by comparing it !! to climatology. !! !! program history log: !! 2009-jun-3 gayno - initial version !! 2011-apr-26 gayno - Perform gross check first, !! then check against climo. !! 2014-sep-30 gayno - Weekly climo file converted !! to grib 2. !! !! @param[in] kgds_data Grib 1 grid description sect of data to be qcd. !! @param[in] snow_data Snow cover to be qcd. !! @param[in] bitmap_data bitmap of data to be qcd. !! @param[in] idata I dimension of data to be qcd. !! @param[in] jdata J dimension of data to be qcd. !! @param[in] isrc Flag indicating data source; 1- afwa depth, 2-ims cover. !! @param[out] bad When true, data failed check. !! !! files: !! input: !! - NH weekly climatological snow cover file (grib 2). !! !! @author George Gayno org: w/np2 @date 2009-Jun-3 subroutine nh_climo_check(kgds_data,snow_data,bitmap_data,idata,jdata,isrc,bad) use gdswzd_mod use program_setup, only : climo_qc_file, & grib_year, grib_month, grib_day, & grib_century use grib_mod ! for grib2 library implicit none ! describes the climo data grid. integer, parameter :: iclim = 1080 integer, parameter :: jclim = 270 real, parameter :: lat11_clim = 90.0 real, parameter :: lon11_clim = -180.0 real, parameter :: dx_clim = 1./3. real, parameter :: dy_clim = 1./3. integer, intent(in) :: idata, jdata, kgds_data(200), isrc logical*1, intent(in) :: bitmap_data(idata,jdata) logical, intent(out) :: bad real, intent(in) :: snow_data(idata,jdata) ! local variables integer :: idat(8), jdow, jdoy, jday integer :: century, year, week, iret, lugb, i, j, ii, jj integer :: lugi, jdisc, jpdtn, jgdtn, k, nret integer :: jids(200), jgdt(200), jpdt(200) integer :: count_nosnow_climo, count_nosnow_data integer :: count_snow_climo, count_snow_data, count_grosschk_data logical*1, allocatable :: bitmap_clim(:,:) logical :: unpack real, allocatable :: climo(:,:) real :: fill, percent, x, y real, allocatable :: xpts(:,:),ypts(:,:),rlon_data(:,:),rlat_data(:,:) real :: thresh_gross, thresh type(gribfield) :: gfld bad=.false. if (len_trim(climo_qc_file)==0) return print*,"- QC SNOW DATA IN NH." if (isrc==1) then thresh_gross=50.0 ! afwa data is depth in meters elseif (isrc==2) then thresh_gross=100.0 ! nesdis/ims data is coverage in percent endif fill=999. allocate(xpts(idata,jdata)) allocate(ypts(idata,jdata)) allocate(rlon_data(idata,jdata)) allocate(rlat_data(idata,jdata)) do j=1,jdata do i=1,idata xpts(i,j)=i ypts(i,j)=j enddo enddo print*,"- CALC LAT/LONS OF SOURCE POINTS." call gdswzd(kgds_data,1,(idata*jdata),fill,xpts,ypts,rlon_data,rlat_data,nret) deallocate(xpts,ypts) if (nret /= (idata*jdata)) then print*,"- WARNING: CALC FAILED. WILL NOT PERFORM QC." deallocate (rlon_data,rlat_data) return endif count_grosschk_data=0 do j=1,jdata do i=1,idata if (rlat_data(i,j)>0.0 .and. bitmap_data(i,j)) then if (snow_data(i,j) < 0.0 .or. snow_data(i,j) > thresh_gross) then count_grosschk_data=count_grosschk_data+1 endif endif enddo enddo if (count_grosschk_data > 1) then print*,'- NUMBER OF DATA POINTS THAT FAIL GROSS CHECK ',count_grosschk_data deallocate (rlon_data,rlat_data) bad=.true. return endif print*,"- QC DATA SOURCE AGAINST CLIMO." print*,"- OPEN CLIMO SNOW COVER FILE ",trim(climo_qc_file) lugb=11 call baopenr(lugb,climo_qc_file,iret) if (iret /= 0) then print*,"- WARNING: BAD OPEN, WILL NOT PERFORM QC ", iret deallocate (rlon_data,rlat_data) return endif !--------------------------------------------------------------- ! climo file is weekly. so calculate the current week ! then read that record from the climo file. !--------------------------------------------------------------- if (grib_year == 100) then century = grib_century else century = grib_century-1 endif year = century*100 + grib_year idat=0 idat(1)=year idat(2)=grib_month idat(3)=grib_day call w3doxdat(idat,jdow,jdoy,jday) ! the climo file date is the beginning of the 7 day period week = nint((jdoy+3.)/7.) if (week==0) week=52 if (week==53) week=1 print*,"- READ CLIMO FOR WEEK ",week call grib2_null(gfld) j = week-1 ! search for specific week (# records to skip) lugi = 0 ! no grib index file jdisc = 0 ! search for discipline; 0 - meteorological products jpdtn = 8 ! search for product definition template number; 8 - average jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid jids = -9999 ! array of values in identification section, set to wildcard jgdt = -9999 ! array of values in grid definition template 3.m jgdt(8) = iclim ! search for assumed grid specs - i/j dimensions and corner ! point lat/lons must match. jgdt(9) = jclim jgdt(12) = nint(lat11_clim * 1e6) jgdt(13) = nint(abs(lon11_clim) * 1e6) jpdt = -9999 ! array of values in product definition template 4.n jpdt(1) = 1 ! search for parameter category - moisture jpdt(2) = 201 ! search for parameter number - snow cover in percent. unpack = .true. ! unpack data call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, iret) if (iret /= 0) then print*,"- WARNING: PROBLEM READING GRIB FILE ", iret print*,"- WILL NOT PERFORM QC." deallocate(rlon_data,rlat_data) call baclose(lugb,iret) return endif call baclose(lugb,iret) allocate(climo(iclim,jclim)) climo = reshape (gfld%fld , (/iclim,jclim/) ) allocate(bitmap_clim(iclim,jclim)) bitmap_clim = reshape (gfld%bmap , (/iclim,jclim/) ) call grib2_free(gfld) !--------------------------------------------------------------- ! loop over all data points in nh. gross check data. ! afwa is a depth in meters, ims is % coverage. there should be ! no neg values or very large values. if point passes gross check, ! then check against climatology. find the ! nearest point on the climo snow cover grid. if ! climo indicates snow is likely (100% coverage), then ! check if afwa/ims has snow. if climo indicates snow is ! impossible (0% coverage), then check if afwa/ims has no snow. if ! afwa/ims differs from climo too much, then afwa/ims is ! considered suspect and will not be used. !--------------------------------------------------------------- count_nosnow_climo=0 count_nosnow_data=0 count_snow_data=0 count_snow_climo=0 if (isrc==1) then thresh=.005 elseif (isrc==2) then thresh=50.0 endif do j=1,jdata do i=1,idata if (rlat_data(i,j)>0.0 .and. bitmap_data(i,j)) then y = (lat11_clim-rlat_data(i,j))/dy_clim + 1.0 if (rlon_data(i,j)>180.0) rlon_data(i,j)=rlon_data(i,j)-360.0 x = (rlon_data(i,j)-lon11_clim)/dx_clim + 1.0 jj=nint(y) if (jj<1) jj=1 if (jj>jclim) jj=jclim ii=nint(x) if (ii<1) ii=ii+iclim if (ii>iclim) ii=ii-iclim if (bitmap_clim(ii,jj)) then ! climo point is land if (climo(ii,jj) <1.0) then ! climo point is snow impossible count_nosnow_climo=count_nosnow_climo+1 if (snow_data(i,j) == 0.0) then count_nosnow_data=count_nosnow_data+1 endif endif if (climo(ii,jj) > 99.) then ! climo point is snow likely count_snow_climo=count_snow_climo+1 if (snow_data(i,j) >thresh) then count_snow_data=count_snow_data+1 endif endif endif endif enddo enddo percent = float(count_snow_climo-count_snow_data) / float(count_snow_climo) percent = percent*100. write(6,200) '- NUMBER OF DATA POINTS THAT SHOULD HAVE SNOW',count_snow_climo write(6,201) '- NUMBER OF THESE POINTS THAT ARE BARE GROUND',(count_snow_climo-count_snow_data), & 'OR', percent, '%' 200 format(1x,a45,1x,i10) 201 format(1x,a45,1x,i10,1x,a2,1x,f6.2,a1) if (percent>50.0) then print*,"- WARNING: PERCENTAGE OF BARE GROUND POINTS EXCEEDS ACCEPTABLE LEVEL." print*,"- WILL NOT USE SOURCE DATA." bad=.true. endif percent = float(count_nosnow_climo-count_nosnow_data) / float(count_nosnow_climo) percent = percent*100. write(6,202) '- NUMBER OF DATA POINTS THAT SHOULD *NOT* HAVE SNOW',count_nosnow_climo write(6,203) '- NUMBER OF THESE POINTS WITH SNOW',(count_nosnow_climo-count_nosnow_data), & 'OR', percent, '%' 202 format(1x,a51,1x,i10) 203 format(1x,a34,1x,i10,1x,a2,1x,f6.2,a1) if (percent>20.0) then print*,"- WARNING: PERCENTAGE OF POINTS WITH SNOW EXCEEDS ACCEPTABLE LEVEL." print*,"- WILL NOT USE SOURCE DATA." bad=.true. endif if (allocated(rlat_data)) deallocate (rlat_data) if (allocated(rlon_data)) deallocate (rlon_data) if (allocated(climo)) deallocate (climo) if (allocated(bitmap_clim)) deallocate (bitmap_clim) return end subroutine nh_climo_check !> Check for corrupt afwa data. !! !! @param[in] hemi (1-nh, 2-sh) !! @author George Gayno org: w/np2 @date 2009-Jun-3 subroutine afwa_check(hemi) use gdswzd_mod implicit none integer, intent(in) :: hemi integer :: kgds(200), nret integer, parameter :: npts=1 real :: fill, xpts(npts), ypts(npts) real :: rlon(npts), rlat(npts) kgds=0 fill=9999. if (hemi==1) then print*,'- QC DATA IN NH.' kgds=kgds_afwa_nh rlat=75.0 rlon=-40. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_nh(nint(xpts(1)),nint(ypts(1))) < 0.001) then print*,'- WARNING: NO SNOW IN GREENLAND: ',snow_dep_afwa_nh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_nh=.true. endif rlat=3.0 rlon=-60. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_nh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN S AMERICA: ',snow_dep_afwa_nh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_nh=.true. endif rlat=23.0 rlon=10. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_nh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN SAHARA: ',snow_dep_afwa_nh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_nh=.true. endif rlat=15.0 rlon=10. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_nh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN S INDIA: ',snow_dep_afwa_nh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_nh=.true. endif endif if (hemi==2) then print*,'- QC DATA IN SH.' kgds=kgds_afwa_sh rlat=-88.0 rlon=0. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_sh(nint(xpts(1)),nint(ypts(1))) < 0.001) then print*,'- WARNING: NO SNOW IN ANTARCTICA: ',snow_dep_afwa_sh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_sh=.true. endif rlat=-10. rlon=-45. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_sh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN SOUTH AMERICA: ',snow_dep_afwa_sh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_sh=.true. endif rlat=-20.0 rlon=130. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_sh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN AUSTRALIA: ',snow_dep_afwa_sh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_sh=.true. endif rlat=-9.0 rlon=25. call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) if (snow_dep_afwa_sh(nint(xpts(1)),nint(ypts(1))) > 0.0) then print*,'- WARNING: SNOW IN AFRICA: ',snow_dep_afwa_sh(nint(xpts),nint(ypts)) print*,'- DONT USE AFWA DATA.' bad_afwa_sh=.true. endif endif end subroutine afwa_check !> Read afwa binary snow depth file. !! !! @param[in] file_name file name !! @param[out] snow_dep_afwa snow depth in meters !! !! files: !! input: !! - nh/sh afwa data in simple binary format !! !! condition codes: all fatal !! 60 - bad open of afwa file !! 61 - bad read of afwa file !! !! @note Read logic for binary data is taken from hua-lu's code, !! /nwprod/sorc/grib_snowgrib.fd. !! !! @author George Gayno org: w/np2 @date 2007-Nov-28 subroutine read_afwa_binary(file_name, snow_dep_afwa) implicit none character*8 :: afwa_file_info(2) character*(*), intent(in) :: file_name integer*2, allocatable :: dummy(:,:) integer :: i,j, istat integer, parameter :: iafwa = 512 integer, parameter :: jafwa = 512 integer, parameter :: iunit=11 ! input afwa data file real, intent(out) :: snow_dep_afwa(iafwa,jafwa) print*,"- OPEN AFWA BINARY FILE ", trim(file_name) open (iunit, file=trim(file_name), access="direct", recl=iafwa*2, iostat=istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD OPEN. ISTAT IS ',istat call w3tage('SNOW2MDL') call errexit(60) end if print*,"- READ AFWA BINARY FILE ", trim(file_name) read(iunit, rec=2, iostat = istat) afwa_file_info if (istat /= 0) then print*,'- FATAL ERROR: BAD READ. ISTAT IS ',istat call w3tage('SNOW2MDL') call errexit(61) end if print*,"- AFWA DATA IS ", afwa_file_info(1), " AT TIME ", afwa_file_info(2)(2:7) allocate(dummy(iafwa,jafwa)) do j = 1, jafwa read(iunit, rec=j+2, iostat=istat) (dummy(i,j),i=1,iafwa) if (istat /= 0) then print*,'- FATAL ERROR: BAD READ. ISTAT IS ',istat call w3tage('SNOW2MDL') call errexit(61) end if enddo close(iunit) !----------------------------------------------------------------------- ! "4090" is the sea ice flag. we don't use the afwa sea ice. !----------------------------------------------------------------------- where (dummy == 4090) dummy = 0 snow_dep_afwa = float(dummy) !--------------------------------------------------------------------- ! afwa data is a snow depth in units of tenths of inches. ! convert this to meters. !--------------------------------------------------------------------- snow_dep_afwa = snow_dep_afwa * 2.54 / 1000.0 deallocate (dummy) return end subroutine read_afwa_binary !> Read afwa land mask file to get a bitmap. !! !! @param[in] file_name land mask file name !! @param[out] bitmap_afwa .true. if land !! !! files: !! input: !! - afwa landmask in simple binary format !! !! condition codes: all fatal !! 62 - bad open of afwa landmask file !! 63 - bad read of afwa landmask file !! !! @author George Gayno org: w/np2 @date 2007-Nov-28 subroutine read_afwa_mask(file_name, bitmap_afwa) implicit none character*(*), intent(in) :: file_name integer, parameter :: iunit=11 ! input mask file integer, parameter :: iafwa = 512 integer, parameter :: jafwa = 512 integer :: i, j, istat integer*4, allocatable :: dummy4(:,:) logical*1, intent(out) :: bitmap_afwa(iafwa,jafwa) allocate (dummy4(iafwa,jafwa)) print*,'- OPEN AFWA MASK FILE ', trim(file_name) open(iunit, file=trim(file_name), access='direct', & recl=iafwa*jafwa*4, iostat=istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD OPEN. ISTAT IS ', istat call w3tage('SNOW2MDL') call errexit(62) end if print*,'- READ AFWA MASK FILE ', trim(file_name) read(iunit, rec=1, iostat=istat) dummy4 if (istat /= 0) then print*,'- FATAL ERROR: BAD READ. ISTAT IS ', istat call w3tage('SNOW2MDL') call errexit(63) end if close(iunit) !----------------------------------------------------------------------- ! here -1-offhemi, 1-ocean, 2-land, 4-coast. !----------------------------------------------------------------------- bitmap_afwa = .false. do j = 1, jafwa do i = 1, iafwa if (dummy4(i,j) > 1) then bitmap_afwa(i,j) = .true. endif enddo enddo deallocate (dummy4) end subroutine read_afwa_mask end module snowdat