program process_NSSL_mosaic
!
! PRGMMR: Ming Hu ORG: GSD DATE: 2007-12-17
!
! ABSTRACT:
! This routine read in NSSL reflectiivty mosaic fiels and
! interpolate them into GSI mass grid
!
! tversion=8 : NSSL 8 tiles netcdf
! tversion=81 : NCEP 8 tiles binary
! tversion=4 : NSSL 4 tiles binary
! tversion=1 : NSSL 1 tile grib2
!
! PROGRAM HISTORY LOG:
!
! variable list
!
! USAGE:
! INPUT FILES: mosaic_files
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90 + EXTENSIONS
! MACHINE: wJET
!
!$$$
!
!_____________________________________________________________________
!
! use constants, only: zero,one_tenth,one,deg2rad,rad2deg
! use gridmod, only: regional,nlon,nlat,nsig, &
! tll2xy,txy2ll, &
! regional_time,nhr_assimilation, &
! regional_fhr, &
! ylat,xlon
use mpi
use kinds, only: r_kind,i_kind
implicit none
!
INCLUDE 'netcdf.inc'
!
! MPI variables
integer :: npe, mype, mypeLocal,ierror
real :: rad2deg = 180.0/3.1415926
!
character*256 output_file
!
! grid
integer(i_kind) :: nlon,nlat
real,allocatable:: xlon(:,:) !
real,allocatable:: ylat(:,:) !
REAL, allocatable :: ref3d(:,:,:) ! 3D reflectivity
REAL, allocatable :: ref0(:,:,:) ! 3D reflectivity
REAL, allocatable :: maxref(:,:) ! composite reflectivity
integer , allocatable :: imaxref(:,:) ! composite reflectivity
REAL(r_kind), allocatable :: ref3d_column(:,:) ! 3D reflectivity in column
CHARACTER*180 geofile
!
! For reflectiivty mosaic
!
CHARACTER*256 mosaicfile
CHARACTER*256 filenameall(200)
INTEGER :: mscNlon ! number of longitude of mosaic data
INTEGER :: mscNlat ! number of latitude of mosaic data
INTEGER :: mscNlev ! number of vertical levels of mosaic data
REAL, allocatable :: msclon(:) ! longitude of mosaic data
REAL, allocatable :: msclat(:) ! latitude of mosaic data
REAL, allocatable :: msclev(:) ! level of mosaic data
REAL, allocatable :: mscValue(:,:) ! reflectivity
REAL, allocatable :: mscValue3d(:,:,:) ! reflectivity
REAL :: lonMin,latMin,lonMax,latMax
REAL*8 :: dlon,dlat
integer :: height
integer,parameter :: maxMosaiclvl=33
integer :: levelheight(maxMosaiclvl)
data levelheight /500, 750, 1000, 1250, 1500, 1750, 2000, 2250, 2500,&
2750,3000,3500, 4000,4500,5000,5500,6000,6500,7000,&
7500,8000,8500,9000,10000,11000,12000,13000,14000, &
15000,16000,17000,18000,19000/
!
! 4 Tile binary format
!
integer :: ntot, ntot2d, mt
integer*4 :: nx,ny,nz
integer*4 :: yr, mo, da, hr, mn, sc
real*8 :: rdx,rdy
real*4 :: rdx4,rdy4
real :: rlatmax,rlonmin
integer*4 :: var_scale
integer*2, dimension(:), allocatable :: var
!
!
! namelist files
!
INTEGER(i_kind) :: tversion
character*10 :: analysis_time
CHARACTER*180 dataPath
namelist/setup/ tversion,analysis_time,dataPath
integer(i_kind) :: idate
!
!
! ** misc
real ::rix ! x-grid coordinate (grid units)
real ::riy ! y-grid coordinate (grid units)
logical ::outside ! .false., then point is inside x-y domain
! .true., then point is outside x-y domain
logical :: fileexist
integer i,j,k,itype,iymdh,ier,jret,ifn,kk
integer iz,n,nlv,isao,nflag,np,ilen,iflag,iostat
integer :: NCID
REAL :: rlat,rlon
INTEGER :: ip,jp,ipp1,jpp1
REAL :: rip,rjp
REAL :: dip,djp
REAL :: w1,w2,w3,w4
REAL :: ref1,ref2,ref3,ref4,refl_ltng
INTEGER(i_kind) :: maxlvl, nlevel,worklevel
INTEGER(i_kind) :: numlvl,numref
integer :: status
REAL :: rthresh_ref,rthresh_miss
integer :: maxcores
!**********************************************************************
!
! 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) write(*,*) mype, 'deal with mosaic'
open(15, file='mosaic.namelist')
read(15,setup)
close(15)
!
! safty check for cores used in this run
!
read(analysis_time,'(I10)') idate
if(mype==0) write(6,*) 'cycle time is :', idate
if( tversion == 8 .or. tversion == 14) then
maxcores=8
elseif( tversion == 81 ) then
maxcores=8
elseif( tversion == 4 ) then
maxcores=4
elseif( tversion == 1 ) then
maxcores=33
else
write(*,*) 'unknow tversion !'
stop 1234
endif
write(6,*) 'total cores for this run is ',npe
if(npe < maxcores) then
write(6,*) 'ERROR, this run must use ',maxcores,' or more cores !!!'
call MPI_FINALIZE(ierror)
stop 1234
endif
!
!
if( tversion == 8 .or. tversion == 14) then
maxlvl = 31
rthresh_ref=-500.0
rthresh_miss=-5000.0
elseif( tversion == 81 ) then
maxlvl = 31
rthresh_ref=-90.0
rthresh_miss=-900.0
elseif( tversion == 4 ) then
maxlvl = 33
rthresh_ref=-500.0
rthresh_miss=-5000.0
elseif( tversion == 1 ) then
maxlvl = 33
rthresh_ref=-90.0
rthresh_miss=-900.0
if(npe < maxlvl) then
write(*,*) 'npe must larger or equal to maxlvl'
stop 1234
endif
else
write(*,*) 'unknow tversion !'
stop 1234
endif
!
! set geogrid fle name
!
write(geofile,'(a,a)') './', 'geo_em.d01.nc'
if(mype==0) write(*,*) 'geofile', trim(geofile)
call GET_DIM_ATT_geo(geofile,NLON,NLAT)
if(mype==0) write(*,*) 'NLON,NLAT',NLON,NLAT
!
! get GSI horizontal grid in latitude and longitude
!
allocate(xlon(nlon,nlat))
allocate(ylat(nlon,nlat))
call OPEN_geo(geofile, NCID)
call GET_geo_sngl_geo(NCID,Nlon,Nlat,ylat,xlon)
call CLOSE_geo(NCID)
!
mypeLocal=mype+1
if ( tversion == 1 ) then
open(10,file='filelist_mrms',form='formatted',err=300)
do n=1,200
read(10,'(a)',err=200,end=400) filenameall(n)
enddo
300 write(6,*) 'read_grib2 open filelist_mrms failed ',mype
stop(555)
200 write(6,*) 'read_grib2 read msmr file failed ',n,mype
stop(555)
400 nlevel=n-1
close(10)
if(nlevel .gt. maxlvl) then
write(*,*) 'vertical level is too large:',nlevel,maxlvl
stop 666
endif
if(mypeLocal <= npe) then
mosaicfile=trim(filenameall(mypeLocal))
write(*,*) 'process level:',mypeLocal,trim(mosaicfile)
else
mosaicfile=''
endif
if(npe < nlevel ) then
write(6,*) 'Error, need more cores to run',npe
stop 234
endif
elseif(tversion==81) then
if(mypeLocal <= 8) then
write(mosaicfile,'(a,a,I1)') trim(dataPath), 'mosaic_t',mypeLocal
write(*,*) 'process tile:',trim(mosaicfile)
else
mosaicfile=''
endif
if(npe <8 ) then
write(6,*) 'Error, need more cores to run',npe
stop 234
endif
else
if(mypeLocal <= tversion) then
write(mosaicfile,'(a,a,I1)') trim(dataPath), 'mosaic_t',mypeLocal
write(*,*) 'process tile:',trim(mosaicfile)
else
mosaicfile=''
endif
if(npe < tversion ) then
write(6,*) 'Error, need more cores to run',npe
stop 234
endif
endif
!
! deal with certain tile
!
fileexist=.false.
if( tversion == 8 .or. tversion == 14) then
call ifexist_file(mosaicfile,STATUS)
fileexist=STATUS .EQ. NF_NOERR
elseif( tversion == 81) then
inquire(file=trim(mosaicfile),exist=fileexist)
if(mypeLocal > 8) fileexist=.false.
elseif(tversion == 4) then
open(99,file=trim(mosaicfile),form='unformatted',access='direct',&
recl=6*4,status='old',err=225)
rewind(99)
read(99,rec=1,err=225) yr, mo, da, hr, mn, sc
fileexist=.true.
225 continue
close(99)
elseif(tversion == 1) then
inquire(file=trim(mosaicfile),exist=fileexist)
if(mypeLocal > nlevel) fileexist=.false.
endif
if(fileexist) then
IF( tversion == 14 ) then
call GET_DIM_ATT_Mosaic(mosaicfile,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat)
var_scale=10
ELSEIF( tversion == 8 ) then
call GET_DIM_ATT_Mosaic8(mosaicfile,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat)
var_scale=10
ELSEIF( tversion == 81 ) then
call read_ncep_binary_head(mype,mosaicfile,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,rdx4,rdy4)
nx=mscNlon
ny=mscNlat
nz=mscNlev
dlon=rdx4
dlat=rdy4
var_scale=1
ELSEIF( tversion == 4 ) then
call read_head_Mosaic4(mosaicfile,nx,ny,nz,rlonmin,rlatmax,&
rdx,rdy,var_scale)
mscNlon=nx
mscNlat=ny
mscNlev=nz
dlon=rdx
dlat=rdy
lonMin=rlonmin
lonMax=lonMin+dlon*(mscNlon-1)
latMax=rlatmax
latMin=latMax-dlat*(mscNlat-1)
ELSEIF( tversion == 1 ) then
call read_grib2_head(mosaicfile,nx,ny,nz,rlonmin,rlatmax,&
rdx,rdy)
var_scale=1
mscNlon=nx
mscNlat=ny
mscNlev=nz
dlon=rdx
dlat=rdy
lonMin=rlonmin
lonMax=lonMin+dlon*(mscNlon-1)
latMax=rlatmax
latMin=latMax-dlat*(mscNlat-1)
ELSE
write(*,*) ' unknown tile version !!!'
stop 123
ENDIF
if(mype==0) then
write(*,*) 'mscNlon,mscNlat,mscNlev'
write(*,*) mscNlon,mscNlat,mscNlev
write(*,*) 'dlon,dlat,lonMin,lonMax,latMax,latMin'
write(*,*) dlon,dlat,lonMin,lonMax,latMax,latMin
endif
if( maxlvl == mscNlev .or. maxlvl >= nlevel ) then
allocate(ref3d(nlon,nlat,maxlvl))
else
write(*,*) 'Wrong vertical layers:', maxlvl, mscNlev,nlevel
stop 1234
endif
ref3d=-999.0
allocate(msclon(mscNlon))
allocate(msclat(mscNlat))
allocate(msclev(mscNlev))
allocate(mscValue(mscNlon,mscNlat))
DO i=1,mscNlon
msclon(i)=lonMin+(i-1)*dlon
ENDDO
DO i=1,mscNlat
msclat(i)=latMin+(i-1)*dlat
ENDDO
!
! ingest mosaic file and interpolation
!
if( tversion == 8 .or. tversion == 14) then
call OPEN_Mosaic(mosaicfile, NCID)
if(tversion == 14 ) then
call Check_DIM_ATT_Mosaic(NCID,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat)
elseif(tversion == 8 ) then
call Check_DIM_ATT_Mosaic8(NCID,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat)
endif
write(*,*) mscNlon,mscNlat,mscNlev
write(*,*) 'Area of tile=',lonMin,latMin,lonMax,latMax,dlon,dlat
elseif(tversion == 81) then
allocate(mscValue3d(mscNlon,mscNlat,mscNlev))
call read_ncep_binary_value(mype,mosaicfile,mscNlon,mscNlat,mscNlev, &
mscValue3d)
elseif(tversion == 4) then
ntot = nx*ny*nz
allocate(var(ntot))
call read_data_Mosaic4(mosaicfile,ntot,var)
elseif(tversion == 1) then
ntot = nx*ny*nz
call read_grib2_sngle(mosaicfile,ntot,height,mscValue)
! write(*,*) 'height,max,min',height,maxval(mscValue),minval(mscValue)
if(levelheight(mypeLocal) .eq. height) then
worklevel=mypeLocal
else
worklevel=0
do k=1,maxMosaiclvl
if (levelheight(k) .eq. height) worklevel=k
enddo
if(worklevel==0) then
write(6,*) 'Error, cannot find working level', &
mypeLocal,levelheight(mypeLocal), height
stop 12345
else
write(6,*) 'Find new level for core ',mypeLocal,' new level is',worklevel
endif
endif
else
write(*,*) 'unknown type'
stop 1234
endif
endif
call mpi_barrier(MPI_COMM_WORLD,ierror)
! stop 999
if(fileexist) then
!
if(tversion == 1) mscNlev=1
DO k=1, mscNlev
! if(tversion > 1) write(*,*) mype, 'deal with level:', k,mscNlon,mscNlat
if( tversion == 8 .or. tversion == 14) then
call GET_Mosaic_sngl_Mosaic(NCID,mscNlon,mscNlat,k,mscValue)
elseif(tversion == 81) then
do j=1,ny
do i=1,nx
mscValue(i,j) = mscValue3d(i,j,k)
enddo
enddo
elseif(tversion == 4) then
ntot2d=nx*ny*(k-1)
do j=1,ny
do i=1,nx
mscValue(i,j) = var(ntot2d+(j-1)*nx+i)
enddo
enddo
elseif(tversion == 1) then
write(*,*) 'level max min height',mypeLocal,maxval(mscValue),minval(mscValue),height
endif
DO j=1,nlat
DO i=1,nlon
rlat=ylat(i,j)
rlon=xlon(i,j)
if(tversion == 14 ) then
rip=(rlon-lonMin)/dlon+1
rjp=(rlat-latMin)/dlat+1
ip=int(rip)
jp=int(rjp)
dip=rip-ip
djp=rjp-jp
elseif(tversion == 8 .or. tversion == 81) then
rip=(rlon-lonMin)/dlon+1
rjp=(latMax-rlat)/dlat+1
ip=int(rip)
jp=int(rjp)
dip=rip-ip
djp=rjp-jp
elseif(tversion == 4 ) then
rip=(rlon-lonMin)/dlon+1
rjp=(rlat-latMin)/dlat+1
ip=int(rip)
jp=int(rjp)
dip=rip-ip
djp=rjp-jp
elseif(tversion == 1 ) then
if(rlon<0.0) rlon=360.0+rlon
rip=(rlon-lonMin)/dlon+1
rjp=(latMax-rlat)/dlat+1
ip=int(rip)
jp=int(rjp)
dip=rip-ip
djp=rjp-jp
else
write(*,*) ' Unknown Mosaic format !!'
stop 123
endif
if( ip >= 1 .and. ip <= mscNlon ) then
if( jp >= 1 .and. jp <= mscNlat ) then
! inside mosaic domain
ipp1=min(ip+1,mscNlon)
jpp1=min(jp+1,mscNlat)
w1=(1.0-dip)*(1.0-djp)
w2=dip*(1.0-djp)
w3=dip*djp
w4=(1.0-dip)*djp
ref1=mscValue(ip,jp)
ref2=mscValue(ipp1,jp)
ref3=mscValue(ipp1,jpp1)
ref4=mscValue(ip,jpp1)
kk=k
if(tversion == 1) kk=worklevel
if(ref1 > rthresh_ref .and. ref2 > rthresh_ref .and. &
ref3 > rthresh_ref .and. ref4 > rthresh_ref ) then
ref3d(i,j,kk)=(ref1*w1+ref2*w2+ref3*w3+ref4*w4)/float(var_scale)
elseif(ref1 > rthresh_miss .and. ref2 > rthresh_miss .and. &
ref3 > rthresh_miss .and. ref4 > rthresh_miss ) then
ref3d(i,j,kk)=-99.0 ! clear
else
ref3d(i,j,kk)=-999.0 ! no observation
endif
endif
endif
ENDDO
ENDDO
ENDDO ! mscNlev
if( tversion == 8 .or. tversion == 14) then
call CLOSE_Mosaic(NCID)
endif
if(tversion == 4) deallocate(var)
if(tversion == 81) deallocate(mscValue3d)
deallocate(msclon)
deallocate(msclat)
deallocate(msclev)
deallocate(mscValue)
else
allocate(ref3d(nlon,nlat,maxlvl))
ref3d=-999.0
write(*,*) trim(mosaicfile), ' does not exist!!!'
ENDIF
call mpi_barrier(MPI_COMM_WORLD,ierror)
!
! collect data from all processes to root (0)
!
if(mype==0) then
allocate( ref0(nlon,nlat,maxlvl) )
allocate( maxref(nlon,nlat) )
allocate( imaxref(nlon,nlat) )
endif
call MPI_REDUCE(ref3d, ref0, nlon*nlat*maxlvl, MPI_REAL, MPI_MAX, 0, &
MPI_COMM_WORLD, ierror)
deallocate(ref3d)
!
if(mype==0) then
OPEN(10,file='./'//'RefInGSI3D.dat',form='unformatted')
write(10) maxlvl,nlon,nlat
write(10) ref0
close(10)
DO k=1,maxlvl
write(*,*) k,maxval(ref0(:,:,k)),minval(ref0(:,:,k))
ENDDO
endif
if(mype==0 .and. 1==1) then
!
allocate(ref3d_column(maxlvl+2,nlon*nlat))
ref3d_column=-999.0
numref=0
! DO j=1,nlat
! DO i=1,nlon
DO j=2,nlat-1
DO i=2,nlon-1
numlvl=0
DO k=1,maxlvl
if(abs(ref0(i,j,k)) < 888.0 ) numlvl=numlvl+1
ENDDO
if(numlvl > 0 ) then
numref=numref+1
ref3d_column(1,numref)=float(i)
ref3d_column(2,numref)=float(j)
DO k=1,maxlvl
ref3d_column(2+k,numref)=ref0(i,j,k)
ENDDO
endif
ENDDO
ENDDO
write(*,*) 'Dump out results', numref, 'out of', nlon*nlat
OPEN(10,file='./'//'RefInGSI.dat',form='unformatted')
write(10) maxlvl,nlon,nlat,numref,1,2
write(10) ((ref3d_column(k,i),k=1,maxlvl+2),i=1,numref)
close(10)
write(*,*) 'Start write_bufr_nsslref'
call write_bufr_nsslref(maxlvl,nlon,nlat,numref,ref3d_column,idate)
endif
if(mype==0) write(6,*) "=== RAPHRRR PREPROCCESS SUCCESS ==="
call MPI_FINALIZE(ierror)
!
end program process_NSSL_mosaic