PROGRAM Zmosaic2arps,19
!
!
!-----------------------------------------------------------------------
!
! Program for generating a reflectivity file in ARPS grid from
! NSSL Mosaica reflectivity
!
! Author: Ming Hu, CAPS. University of Oklahma.
! First written: 04/06/2006.
!
! History:
!
! 5/15/2007 Fanyou Kong
! Rewrite to remove mosaicfile, and to include 2D field read
! (e.g., cref, 1h_rad_hsr, etc)
! Also remove NAMELIST jobname (not used)
! Add NetCDF file existence check, and uncompress if applicable
!
! O5/16/2007 Y. Wang
! Added check for "/" at the end of mosaicPath.
! Changed wrtvar1 to wrtvar2.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
! ARPS grid
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'mp.inc'
INCLUDE 'netcdf.inc'
!--------------------------------------------------------------
! ARPS grid
!--------------------------------------------------------------
REAL, allocatable :: lath(:,:) ! Latitude of each terrain point
REAL, allocatable :: lonh(:,:) ! Longitude of each terrain point
REAL, allocatable :: zp(:,:,:) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL, allocatable :: ref3d(:,:,:) !3D reflectivity in mosaic vertical grid
REAL, allocatable :: ref_arps_3d(:,:,:) ! 3D reflectivity in arps grid
REAL, allocatable :: var2d(:,:) ! 2D variable
REAL, allocatable :: compref(:,:) ! composite reflectivity
REAL, allocatable :: hsr_1hr(:,:) ! 1-h accumulate precipitation
!--------------------------------------------------------------
! Namelist
!--------------------------------------------------------------
!
! grid
!
INTEGER :: nx, ny, nz
! NAMELIST /jobname/ runname
NAMELIST /grid/ nx,ny, nz, dx,dy,dz, &
strhopt,dzmin,zrefsfc,dlayer1,dlayer2,strhtune,zflat, &
ctrlat,ctrlon
NAMELIST /projection/ mapproj,trulat1,trulat2,trulon,sclfct
NAMELIST /terrain/ ternopt,terndta,ternfmt
!
! For reflectiivty mosaic
!
CHARACTER*256 mosaicPath
CHARACTER*13 mosaictime(50)
CHARACTER*256 mosaicTile(14)
INTEGER :: tversion,mosaic_opt,ifcompositeRef,numvolume,ntiles
NAMELIST /mosaicpara/mosaicPath,tversion,mosaic_opt, &
ifcompositeRef,numvolume, &
mosaictime,ntiles,mosaicTile
INTEGER :: dmpfmt
NAMELIST /output/ dirname, dmpfmt
!--------------------------------------------------------------
! Grid of NSSL mosaic
!--------------------------------------------------------------
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 :: lonMin,latMin,lonMax,latMax,dlon,dlat
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
CHARACTER*256 tempfile
INTEGER :: NCID,iSTATUS
LOGICAL :: fexist
INTEGER :: maxlvl
REAL :: swx,swy,ctrx,ctry ! Temporary variables.
INTEGER :: i,j,k,ifl,im,ifn
REAL :: rlon,rlat
INTEGER :: ip,jp
REAL :: rip,rjp
REAL :: dip,djp
REAL :: w1,w2,w3,w4
REAL :: ref1,ref2,ref3,ref4
INTEGER :: iyr,imon,iday,ihr,imin,isec
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! CALL mpinit_var ! calling of writtrn needs some variables
!
!-----------------------------------------------------------------------
!
! Set up the default values for all the variables to be read in
! using the namelist method. In case the user does not specify a
! particular value, this value will be used.
!
!-----------------------------------------------------------------------
!
nx = 67
ny = 67
nz = 53
! runname = 'may20'
dx = 1000.0
dy = 1000.0
dz = 400.0
ctrlat = 35.0
ctrlon = -100.0
mapproj = 0
trulat1 = 30.0
trulat2 = 60.0
trulon = -100.0
sclfct = 1.0
!
!-----------------------------------------------------------------------
!
! Read and Print of Namelist variables and parameters.
!
!-----------------------------------------------------------------------
!
! READ (5,jobname,ERR=980,END=990)
! WRITE(6,'(/a,a)') 'The name of this run is: ', runname
READ (5,mosaicpara,ERR=980,END=990)
WRITE(6,*) 'NetCDF data type: ', mosaic_opt
WRITE(6,*) 'Number of times: ', numvolume
WRITE(6,*) 'Number of tiles: ', ntiles
READ(5,grid,ERR=980,END=990)
READ(5,projection,ERR=980,END=990)
READ(5,terrain,ERR=980,END=990)
READ(5,output,ERR=980,END=990)
lfnkey = LEN_TRIM(mosaicpath)
IF (mosaicpath(lfnkey:lfnkey) /= '/') mosaicpath = TRIM(mosaicpath) // '/'
!
! CALL gtlfnkey( runname, lfnkey )
!
! get the latitude and longitude of ARPS domain
!
allocate(lath(0:nx,0:ny))
allocate(lonh(0:nx,0:ny))
allocate(compref(nx,ny))
allocate(hsr_1hr(nx,ny))
CALL getcoordinate
(nx,ny,lath,lonh)
CALL a3dmax0(lath,0,nx,0,nx,0,ny,0,ny,1,1,1,1,latmax,latmin)
CALL a3dmax0(lonh,0,nx,0,nx,0,ny,0,ny,1,1,1,1,lonmax,lonmin)
! write(*,*) lonmax,lonmin
! write(*,*) latmax,latmin
IF(mosaic_opt == 1) THEN
allocate(zp(nx,ny,nz))
allocate(ref_arps_3d(nx,ny,nz))
CALL getVertcoordinate(nx,ny,nz,zp)
END IF
!
! begin to read in mosaic file
!
nv: DO ifl=1, numvolume
!
! get dimension of mosaic field and allocate arrays
!
DO im=1, ntiles
tempfile=trim(mosaicPath)//trim(mosaictime(ifl))//trim(mosaicTile(im))
write(*,*) 'process file No.',im,'with name:',trim(tempfile)
INQUIRE(FILE=trim(tempfile), EXIST = fexist )
IF( .NOT.fexist) THEN
INQUIRE(FILE=trim(tempfile)//'.gz', EXIST = fexist )
IF(fexist) THEN
CALL uncmprs(trim(tempfile)//'.gz')
ELSE
PRINT *,trim(tempfile),' or its .gz file does not exist, SKIP'
CYCLE nv
END IF
END IF
IF( tversion == 14 ) then
call GET_DIM_ATT_Mosaic(tempfile,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat)
ELSEIF( tversion == 8 ) then
call GET_DIM_ATT_Mosaic8(tempfile,mscNlon,mscNlat,mscNlev, &
lonMin,latMin,lonMax,latMax,dlon,dlat,mosaic_opt)
ELSE
write(*,*) ' unknown tile version !!!'
stop 123
ENDIF
allocate(msclon(mscNlon))
allocate(msclat(mscNlat))
allocate(msclev(mscNlev))
allocate(mscValue(mscNlon,mscNlat))
if(im == 1) then
maxlvl=mscNlev
allocate(ref3d(nx,ny,maxlvl))
ref3d=-99999.0
else
if(maxlvl .ne. mscNlev) then
write(*,*) ' The vertical dimensions are different in two tiles'
stop
endif
endif
msclon(1)=lonMin
DO i=1,mscNlon-1
msclon(i+1)=msclon(i)+dlon
ENDDO
msclat(1)=latMin
DO i=1,mscNlat-1
msclat(i+1)=msclat(i)+dlat
ENDDO
!
! ingest mosaic file and interpolation
!
call OPEN_Mosaic(tempfile, 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,mosaic_opt)
endif
DO k=1, mscNlev
IF(mosaic_opt == 1) THEN
call GET_Mosaic_sngl_Mosaic(NCID,mscNlon,mscNlat,k,mscValue)
ELSE IF(mosaic_opt == 2) THEN
call GET_Mosaic_cref_Mosaic(NCID,mscNlon,mscNlat,mscValue)
ELSE IF(mosaic_opt == 3) THEN
call GET_Mosaic_hsr_Mosaic(NCID,mscNlon,mscNlat,mscValue)
ELSE
print *,'mosaic_opt=',mosaic_opt,' is invalid, STOP!!!'
stop
END IF
! DO j=0,ny
! DO i=0,nx
DO j=1,ny
DO i=1,nx
rlat=lath(i,j)
rlon=lonh(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 ) then
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
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(ip+1,jp)
ref3=mscValue(ip+1,jp+1)
ref4=mscValue(ip,jp+1)
if(ref1 > -500.0 .and. ref2 > -500.0 .and. &
ref3 > -500.0 .and. ref4 > -500.0 ) then
ref3d(i,j,k)=(ref1*w1+ref2*w2+ref3*w3+ref4*w4)/10.0
endif
endif
endif
ENDDO
ENDDO
ENDDO ! mscNlev
!
call CLOSE_Mosaic(NCID)
deallocate(msclon)
deallocate(msclat)
deallocate(msclev)
deallocate(mscValue)
ENDDO ! ntiles
write(*,*)
write(*,*) 'successfully process Mosaic data at ', &
trim(mosaictime(ifl))
! write(*,*)
!
IF(mosaic_opt == 1) THEN
call vert_interp_ref(nx,ny,nz,mscNlev,ref3d,ref_arps_3d,zp)
! write(*,*)
write(*,*) 'successfully get reflectivity in ARPS grid at ', &
trim(mosaictime(ifl))
! write(*,*)
!
! dump out colume data for data analysis
!
read(mosaictime(ifl),'(i4,i2,i2,1x,i2,i2)') iyr,imon,iday,ihr,imin
isec=0
call wtradcol_Mosaic(nx,ny,nz, &
mosaictime(ifl),iyr,imon,iday,ihr,imin,isec, &
lath,lonh,zp,ref_arps_3d)
ref_arps_3d=-99999.0
call readadcol_Mosaic(nx,ny,nz,mosaictime(ifl),ref_arps_3d)
!
! get composite reflectivity
!
if(ifcompositeRef==1) then
compref=-99999.0
DO k=2, nz-1
DO j=1,ny
DO i=1,nx
compref(i,j)=max(compref(i,j),ref_arps_3d(i,j,k))
ENDDO
ENDDO
ENDDO
!
! dump out composite reflectivity
!
DO j=1,ny
DO i=1,nx
if( compref(i,j) < 0.0 ) compref(i,j)=0
if( compref(i,j) > 100.0 ) write(*,*) i,j,compref(i,j)
ENDDO
ENDDO
call wrtvar2(nx,ny,1,compref, 'compst','composite reflectivity', &
'dBZ',0,mosaictime(ifl),dirname,dmpfmt,0,0,iSTATUS)
endif ! ifcompositeRef=1
ELSE IF(mosaic_opt == 2) THEN
DO j=1,ny
DO i=1,nx
compref(i,j) = max(ref3d(i,j,1), 0.0)
ENDDO
ENDDO
call wrtvar2(nx,ny,1,compref, 'compst','composite reflectivity', &
'dBZ',0,mosaictime(ifl),dirname,dmpfmt,0,0,iSTATUS)
ELSE IF(mosaic_opt == 3) THEN
DO j=1,ny
DO i=1,nx
hsr_1hr(i,j) = max(ref3d(i,j,1)/25.4, 0.0) ! convert to inch
ENDDO
ENDDO
call wrtvar2(nx,ny,1,hsr_1hr, 'hsr1hr','1-h precipitation', &
'in',0,mosaictime(ifl),dirname,dmpfmt,0,0,iSTATUS)
ELSE
print *,'mosaic_opt=',mosaic_opt,' is invalid, STOP!!!'
stop
END IF
deallocate(ref3d)
222 CONTINUE
ENDDO nv !numvolume
! close(ifn)
deallocate(lath)
deallocate(lonh)
deallocate(compref)
deallocate(hsr_1hr)
IF(mosaic_opt == 1) THEN
deallocate(zp)
deallocate(ref_arps_3d)
END IF
CALL exit(0)
980 CONTINUE
WRITE(6,'(/1x,a,a)') &
'Error occured when reading namelist input file. Program stopped.'
CALL exit(980)
990 CONTINUE
WRITE(6,'(/1x,a,a)') &
'End of file reached when reading namelist input file. ', &
'Program stopped.'
CALL exit(990)
981 CONTINUE
WRITE(6,'(/1x,a,a)') &
'Error occured when reading mosaic file name. Program stopped.'
CALL exit(981)
991 CONTINUE
WRITE(6,'(/1x,a,a)') &
'End of file reached when reading mosaic file path. ', &
'Program stopped.'
CALL exit(991)
STOP
END PROGRAM Zmosaic2arps