PROGRAM process_NESDIS_imssnow ! ! PRGMMR: Ming Hu ORG: GSD DATE: 2009-04-15 ! ! ABSTRACT: ! This routine read in NESDIS NESDIS SNOW/ICE data from a grib file and ! map them into RR mass grid ! ! ! PROGRAM HISTORY LOG: ! ! variable list ! ! USAGE: ! INPUT FILES: imssnow ! ! OUTPUT FILES: RRimssnow ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 + EXTENSIONS ! MACHINE: wJET ! !$$$ ! !_____________________________________________________________________ use mpi use grib2_read_mod, only : read_grib2_head_dim,read_grib2_head_time use grib2_read_mod, only : read_grib2_sngle implicit none ! MPI variables integer :: npe, mype, mypeLocal,ierror ! ! ! for grib2 ! real :: rlatmin,rlonmin real*8 :: rdx,rdy integer :: nxobs,nyobs logical :: fileexist integer :: ntot integer :: kgds(200) character(len=100) :: snowfile ! RR grid integer :: nlonRR,nlatRR ! parameter (nlonRR=1799,nlatRR=1059) real,allocatable :: xlonRR(:,:) ! real,allocatable :: ylatRR(:,:) ! real,allocatable :: xlandRR(:,:) ! real,allocatable :: luseRR(:,:) ! real,allocatable :: xlandIMS(:,:) ! real, allocatable :: snowiceRR(:,:) ! snow/ice in RR ! integer :: i,j ! ! sea ice and snow, and mask REAL, allocatable :: ICEC(:,:) REAL, allocatable :: SNOWICEC(:,:) Logical*1, allocatable :: lmask(:,:) REAL, allocatable :: maskims(:,:) ! integer, allocatable :: imaskims(:,:) integer :: iyear, imonth, iday, ihr, imm !********************************************************************** ! ! END OF DECLARATIONS....start of program ! MPI setup call MPI_INIT(ierror) call MPI_COMM_SIZE(mpi_comm_world,npe,ierror) call MPI_COMM_RANK(mpi_comm_world,mype,ierror) ! if(mype==0) then snowfile='imssnow2' inquire(file=trim(snowfile),exist=fileexist) if(fileexist) then call read_grib2_head_dim(snowfile,nxobs,nyobs,rlonmin,rlatmin,rdx,rdy,kgds) write(*,*) 'grib2 file grid =',nxobs,nyobs,rlonmin,rlatmin,rdx,rdy write(*,*) 'kgds=',kgds(1:11) ! ! if grib1 file ! nxobs=6144 ! nyobs=6144 ! allocate(ICEC(nxobs,nyobs)) allocate(SNOWICEC(nxobs,nyobs)) allocate(lmask(nxobs,nyobs)) allocate(maskims(nxobs,nyobs)) ICEC=0 SNOWICEC=0 call read_grib2_head_time(snowfile,iyear,imonth,iday,ihr,imm) write(*,'(a,5I5)') 'date: iyear,imonth,iday,ihr,imm=',& iyear,imonth,iday,ihr,imm ! read field ntot=nxobs*nyobs call read_grib2_sngle(snowfile,ntot,ICEC,SNOWICEC,lmask) write(*,*) 'read in ice =',maxval(icec),minval(icec) write(*,*) 'read in snow =',maxval(SNOWICEC),minval(SNOWICEC) ! ! read from grib1 file ! call read_issnow(icec,snowicec,lmask,kgds,iyear,imonth,iday,ihr) ! write(*,*) 'NESDIS SNOW and ICE in ', iyear,imonth,iday,ihr ! ! ! Map NESDIS SNOW/ICE to RR grid: Based on the following RUC routine: ! /whome/rucdev/code/13km/hybpre_code/snohires.f ! DO J = 1, nxobs DO I = 1, nyobs IF (icec(I,J) > 0. .OR. snowicec(I,J) > 0.) THEN snowicec(I,J) = 1. ELSE snowicec(I,J) = 0. ENDIF !tgs Logical lmask --> .true. for land, .false. for water IF (lmask(I,J)) THEN maskims(i,j) = 1. ! land ELSE maskims(i,j) = 0. ! water ENDIF ENDDO ENDDO deallocate(icec) deallocate(lmask) ! ! call GET_DIM_ATT_geo('./geo_em.d01.nc',nlonRR,nlatRR) write(*,*) 'grid dimension =',nlonRR,nlatRR allocate(xlonRR(nlonRR,nlatRR)) allocate(ylatRR(nlonRR,nlatRR)) allocate(xlandRR(nlonRR,nlatRR)) allocate(luseRR(nlonRR,nlatRR)) allocate(xlandIMS(nlonRR,nlatRR)) call GET_RR_GRID(xlonRR,ylatRR,nlonRR,nlatRR,xlandRR,luseRR) ! map to RR grid ! allocate(snowiceRR(nlonRR,nlatRR)) snowiceRR=0 ! pecentage call map2RR(snowicec,maskims,kgds,xlandRR,nlonRR,nlatRR,xlonRR,ylatRR,snowiceRR,xlandIMS) deallocate(snowicec) deallocate(maskims) deallocate(xlonRR,ylatRR) ! ! ! trim snow cover field based on NESDIS snow cover data ! call update_SNOWICE_netcdf_mass(snowiceRR, xlandRR, luseRR, nlonRR, nlatRR,xlandIMS) ! else write(*,*) 'file is not exist: ',trim(snowfile) write(*,*) 'stop update snow/ice' endif write(6,*) "=== RAPHRRR PREPROCCESS SUCCESS ===" endif ! mype==0 call MPI_FINALIZE(ierror) ! END PROGRAM process_NESDIS_imssnow SUBROUTINE map2RR(f,maskims,kgds_src,xlandRR,nlonRR,nlatRR,xlonRR,ylatRR,snowiceRR,xlandIMS) ! ! PRGMMR: Ming Hu ORG: GSD DATE: 2009-04-15 ! ! ABSTRACT: ! This routine map NESDIS SNOW/ICE data to RR grid ! ! PROGRAM HISTORY LOG: ! ! variable list ! ! USAGE: ! INPUT FILES: imssnow ! ! OUTPUT FILES: RRimssnow ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 + EXTENSIONS ! MACHINE: wJET ! !$$$ ! !_____________________________________________________________________ use gdswzd_mod, only : gdswzd implicit none INTEGER JPDS(200),JGDS(200),KPDS_SRC(200),KGDS_SRC(200) REAL, intent(in):: F(6144,6144) REAL, intent(in):: maskims(6144,6144) ! grid integer, intent(in) :: nlonRR,nlatRR real, intent(in):: xlandRR(nlonRR,nlatRR) real, intent(in):: xlonRR(nlonRR,nlatRR) ! real, intent(in):: ylatRR(nlonRR,nlatRR) ! ! real, intent(out):: snowiceRR(nlonRR,nlatRR) ! ! real :: xlandIMS(nlonRR,nlatRR) REAL :: DTR REAL :: XPNMC8,YPNMC8,ENNMC8,ALNMC8,ORIENT8 REAL :: XPNMCAF,YPNMCAF,ENNMCAF,ALNMCAF,ORIENTAF real :: YYLAT(1),XLONG(1), RM, RAD, XX(1),YY(1),X, Y integer :: IS,IP1,JS,JP1, IPOINT, JPOINT ! integer :: iland,KOUNT real :: XRATIO,YRATIO,AREA11,AREA21,AREA12,AREA22,AREA integer :: i,j,k,ifound,LL,JPE,JPB,IPE,IPB,NK,MK ! ! integer :: nret,nout real :: dum real, parameter :: undefined_value = -1.0 ! print *,'NESDIS grid info: kgds_src', kgds_src(1:11) nout=0 DO j=1,nlatRR DO i=1,nlonRR ! call gdswiz(kgds_src,-1,1,undefined_value,X,Y, & ! xlonRR(i,j),ylatRR(i,j),nret,0,dum,dum) YYLAT(1)=ylatRR(i,j) XLONG(1)=xlonRR(i,j) call gdswzd(kgds_src,-1,1,undefined_value,XX,YY, & XLONG,YYLAT,nret) X=XX(1) Y=YY(1) if (nret /= 1) then nout=nout+1 ! snowiceRR(i,j) = 0. cycle else ! check if X and Y are defined if(X == undefined_value .or. Y == undefined_value) then ! snowiceRR(i,j) = 0.0 cycle endif if(X < 1.00001 .or. Y < 1.00001) then cycle endif IS = NINT(X) IP1 = IS + 1 JS = NINT(Y) JP1 = JS + 1 ipoint = nint(X) jpoint = nint(Y) ! if (f(ipoint,jpoint) == 1) then ! snowiceRR(i,j) = 1.0 ! else ! snowiceRR(i,j) = 0.0 ! end if end if ! AND ONLY SEA POINTS ARE INTERPOLATED TO SEA POINTS (FOR ICE) ! (NESDIS/IMS LAND MASK: SEA=0,LAND=1, WHILE THE RR MASK in XLAND IS: SEA=2, ! LAND=1). ! ILAND = 1 IF( int(xlandRR(i,j)) == 0 ) ILAND=0 !tgs - buggy IF( int(maskims(i,j)) == ILAND ) THEN IF( int(maskims(is,js)) == ILAND ) THEN snowiceRR(i,j) = F(IPOINT,JPOINT) ELSE ! ! NEAREST NEIGHBOR NOT SAME SFC TYPE, SO USE ALL 4 SURROUNDING POINTS ! KOUNT = 0 ! XRATIO = X - REAL(IS) YRATIO = Y - REAL(JS) ! AREA11 = (1.0E0 - XRATIO) * (1.0E0 - YRATIO) AREA21 = XRATIO * (1.0E0 - YRATIO) AREA12 = (1.0E0 - XRATIO) * YRATIO AREA22 = XRATIO * YRATIO ! ! IF( int(maskims(IS, JS)) .EQ. ILAND) THEN KOUNT = KOUNT + 1 AREA = AREA11 IPOINT = IS JPOINT = JS END IF ! IF( int(maskims(IS, JP1)) .EQ. ILAND ) THEN KOUNT = KOUNT +1 IF (KOUNT .EQ. 1) THEN IPOINT = IS JPOINT = JP1 ELSEIF (AREA12 .GT. AREA) THEN AREA = AREA12 IPOINT = IS JPOINT = JP1 END IF END IF ! IF( int(maskims(IP1, JS)) .EQ. ILAND ) THEN KOUNT = KOUNT + 1 IF (KOUNT .EQ. 1) THEN AREA = AREA21 IPOINT = IP1 JPOINT = JS ELSEIF (AREA21 .GT. AREA) THEN AREA = AREA21 IPOINT = IP1 JPOINT = JS END IF END IF ! ! IF( int(maskims(IP1, JP1)) .EQ. ILAND ) THEN KOUNT = KOUNT + 1 IF (KOUNT .EQ. 1) THEN AREA = AREA22 IPOINT = IP1 JPOINT = JP1 ELSEIF (AREA22 .GT. AREA) THEN AREA = AREA22 IPOINT = IP1 JPOINT = JP1 END IF END IF ! ! DETERMINE SNO/ICE USING NEAREST NEIGHBOR WITH SAME SFC TYPE ! IF(KOUNT .GT. 0) THEN snowiceRR(i,j) = F(IPOINT,JPOINT) ELSE ! ! NO IMMEDIATELY SURROUNDING POINTS IN THE 6144 X 6144 FIELD OF ! SNOW/ICE HAVE THE SAME LAND-SEA TYPE AS THE model POINT. THE ! model POINT MAY BE SMALL ISLAND OR LAKE OR SMALL BAY OR PENNIN. ! (INVARIABLY A SMALL LAKE IN ETA GRID) ! SO EXPAND SEARCH RADIUS AND TAKE FIRST SFC TYPE MATCH ! IPOINT = NINT(X) JPOINT = NINT(Y) ! ! Define the frame (no. of grid points) over which to search for ! a matching land/water type from IMS data for the model gridpoint. ifound=0 DO LL=1,16 JPE = MIN (6144, JPOINT+LL) JPB = MAX (1 , JPOINT-LL) IPE = MIN (6144, IPOINT+LL) IPB = MAX (1 , IPOINT-LL) ! DO NK=IPB,IPE DO MK=JPB,JPE IF ( int(maskims(nk,mk)) == ILAND .and. ifound ==0 ) THEN snowiceRR(i,j) = F(NK,MK) ifound=1 ENDIF ENDDO ! MK ENDDO ! NK ENDDO ! LL ! ! NO LAND/SEA MASK MATCHES FOUND, SO ! A) NORTH OF 55N, WE ASSIGN SNOW/ICE IRRESPECTIVE OF SFC TYPE ! B) SOUTH OF 55N, WE KEEP A PRIORI ZERO DEFAULT ! (THE "B" OPTION BEST FOR WARMER LATS OF U.S., WHERE THIS CONDITION ! IS VIRTUALLY ALWAYS A SMALL ETA LAKE WITH NO COUNTERPART WATER ! NEARBY IN THE NESDIS/IMS GRID, E.G., SALT LAKE, WHERE WE MUST ! AVOID GETTING SEA-ICE OWING TO SURROUNDING SNOW COVER) ! ! IF (YYLAT .GE. 55.0 .and. ifound==0) THEN IF ( ifound==0 ) THEN snowiceRR(i,j) = F(IPOINT,JPOINT) ENDIF ENDIF ! KOUNT .GT. 0 ! ENDIF ! !tgs - save NESDIS 4-km land/water mask xlandIMS(i,j)=maskims(IPOINT,JPOINT) ! ENDDO ! nlon ENDDO ! nlat print*,"- WARNING: NUM of MODEL POINT OUTSIDE NESDIS GRID.",nout end subroutine map2RR SUBROUTINE read_issnow(icec,snowc,lb,KGDS,iyear,imonth,iday,ihr) ! ! PRGMMR: Ming Hu ORG: GSD DATE: 2009-04-15 ! ! ABSTRACT: ! This routine read in NESDIS NESDIS SNOW/ICE data from a grib file ! ! PROGRAM HISTORY LOG: ! ! variable list ! ! USAGE: ! INPUT ! OUTPUT ! icec : ice cover ! snowc: snow cover ! iyear: Year ! imonth: month ! iday: day ! ihr: hour ! INPUT FILES: imssnow ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 + EXTENSIONS ! MACHINE: wJET ! !$$$ ! !_____________________________________________________________________ implicit none INTEGER JF, MBUF PARAMETER (JF=6144*6144) ! PARAMETER (JF=1024*1024) INTEGER JPDS(200),JGDS(200),KPDS(200),KGDS(200) LOGICAL*1 LB(JF) !water/land mask REAL, intent(out):: ICEC(JF) REAL, intent(out):: SNOWC(JF) ! PARAMETER(MBUF=256*1024) PARAMETER(MBUF=256*6144) CHARACTER CBUF(MBUF) INTEGER :: LUGB,LUGI, J, IRET, K, KF character(200):: FNAME1 INTEGER :: I, iyear, imonth, iday, ihr !----------------------------------------------------------------------- LUGB = 11 LUGI = 0 J = 0 ! FNAME1='/public/data/grids/ncep/snow/grib/latest.SNOW_IMS' FNAME1='./imssnow' call baopen(LUGB,trim(fname1),IRET) if ( IRET .ne. 0) then write(6,*) 'bad baopen!!! ', IRET STOP endif JPDS=-1 JGDS=-1 JPDS(5) = 91 ! ice cover call GETGB(LUGB,LUGI,JF,J,JPDS,JGDS, & KF,K,KPDS,KGDS,LB,ICEC,IRET) if (IRET .ne. 0) then write(6,*) 'bad getgb ', IRET STOP endif JPDS(5) = 238 ! ice cover call GETGB(LUGB,LUGI,JF,J,JPDS,JGDS, & KF,K,KPDS,KGDS,LB,SNOWC,IRET) if (IRET .ne. 0) then write(6,*) 'bad getgb ', IRET STOP endif iyear=2000 + KPDS(8) imonth=KPDS(9) iday=KPDS(10) ihr=KPDS(11) ! POLAR STEREOGRAPHIC GRIDS ! KGDS(i) ! (2) - N(I) NR POINTS ALONG LAT CIRCLE ! (3) - N(J) NR POINTS ALONG LON CIRCLE ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) ! (7) - LOV GRID ORIENTATION ! (8) - DX - X DIRECTION INCREMENT ! (9) - DY - Y DIRECTION INCREMENT ! (10) - PROJECTION CENTER FLAG ! (11) - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28) ! do j=1,6144,100 ! write(*,*) (lb(((j-1)*6144 + i)+0.2),i=1,6144,50) ! enddo KGDS=0 KGDS(1)=5 KGDS(2) = 6144 KGDS(3) = 6144 KGDS(4) = -21492 KGDS(5) = 235000 KGDS(6) = 72 KGDS(7) = 280000 KGDS(8) = 4000 KGDS(9) = 4000 KGDS(10) = 0 KGDS(11) = 64 ! ! end subroutine SUBROUTINE GET_RR_GRID(xlon,ylat,nlon,nlat,xland,luse) ! ! PRGMMR: Ming Hu ORG: GSD DATE: 2009-04-15 ! ! ABSTRACT: ! This routine read in Rapid Refresh grid and land mask ! ! PROGRAM HISTORY LOG: ! ! variable list ! ! USAGE: ! INPUT: ! ! OUTPUT: ! xlon: longitude in each grid ! ylat: latitude in each grid ! xland: land mask ! ! INPUT FILES: imssnow ! ! OUTPUT FILES: RRimssnow ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 + EXTENSIONS ! MACHINE: wJET ! !$$$ ! !_____________________________________________________________________ ! implicit none ! INCLUDE 'netcdf.inc' ! grid integer, intent(in) :: nlon,nlat real, intent(out):: xlon(nlon,nlat) ! real, intent(out):: ylat(nlon,nlat) ! real, intent(out):: xland(nlon,nlat) ! real, intent(out):: luse(nlon,nlat) ! ! integer :: tnlon,tnlat ! integer :: NCID CHARACTER*180 geofile CHARACTER*180 workPath ! ! set geogrid fle name ! workPath='./' write(geofile,'(a,a)') trim(workPath), 'geo_em.d01.nc' write(*,*) 'geofile', trim(geofile) call GET_DIM_ATT_geo(geofile,TNLON,TNLAT) if( (TNLON.ne.NLON) .or. (TNLAT.ne.NLAT)) then write(6,*) ' unmatched dimension' write(*,*) 'NLON,NLAT',NLON,NLAT, 'TNLON,TNLAT',TNLON,TNLAT stop 123 endif ! ! get GSI horizontal grid in latitude and longitude ! call OPEN_geo(geofile, NCID) call GET_geo_sngl_geo(NCID,Nlon,Nlat,ylat,xlon,xland,luse) call CLOSE_geo(NCID) END SUBROUTINE GET_RR_GRID