PROGRAM ARPSTRAJC ,46
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSTRAJC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! (4/08/2004)
!
! MODIFICATION HISTORY:
!
! 9/28/2005
! Modified to handle terrain-following grid. 3D zp is defined.
! Removed used dignositc codes and arrays
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL, allocatable :: x(:),y(:),z(:)
REAL, allocatable :: xs (:), ys (:) ! x,y coord for scalar points
REAL, allocatable :: u(:,:,:)
REAL, allocatable :: v(:,:,:)
REAL, allocatable :: w(:,:,:)
REAL, allocatable :: zp(:,:,:)
! REAL, allocatable :: ptbar(:),pbar(:),rhobar(:),qvbar(:)
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
INTEGER :: nx,ny,nz,nzsoil ! Grid dimensions.
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf
PARAMETER (nhisfile_max=2000)
CHARACTER (LEN=256) :: grdbasfn,hisfile(nhisfile_max)
INTEGER :: hdmpinopt, nf
CHARACTER (LEN=256) :: hdmpftrailer
CHARACTER (LEN=256) :: hdmpfheader
REAL :: tintv_dmpin, tbgn_dmpin, tend_dmpin
NAMELIST /history_data/ hinfmt, nhisfile, grdbasfn, hisfile, &
hdmpinopt,hdmpfheader,hdmpftrailer, &
tintv_dmpin, tbgn_dmpin, tend_dmpin
INTEGER, PARAMETER :: nmax_trajcs=1000, nmax_times=20
REAL :: xtrajc0(nmax_trajcs,nmax_times), &
ytrajc0(nmax_trajcs,nmax_times), &
ztrajc0(nmax_trajcs,nmax_times),ttrajc0
INTEGER :: ntrajc0,ntrajcs, npoints
!
! initrajc = 1 ! specified locations
! = 2 ! circle of radius radius and athimuthal angle increments of theta_inc
! = 3 ! square of radius radius and athimuthal angle increments of theta_inc
! = 4 ! read in from trajectory data file
!
CHARACTER(LEN=256) :: trajc_fn_in(nmax_times)
REAL :: radius,xctr,yctr,zctr(100),theta_inc
INTEGER :: initrajc,nzctr, ntimes
CHARACTER (LEN=256) :: trajcfn_header
REAL :: reftime(nmax_times)
NAMELIST /trajectories/ trajcfn_header,ntimes,reftime, initrajc, ntrajcs, &
xtrajc0,ytrajc0,ztrajc0,ttrajc0,ntrajc0, &
radius,xctr,yctr,zctr,nzctr,theta_inc, trajc_fn_in
REAL, allocatable :: xtrajc(:,:,:),ytrajc(:,:,:),ztrajc(:,:,:),ttrajc(:)
REAL, allocatable :: utrajc(:,:,:),vtrajc(:,:,:),wtrajc(:,:,:)
! REAL, allocatable :: xweight(:,:),yweight(:,:),zweight(:,:)
! INTEGER, allocatable :: itrajc(:,:),jtrajc(:,:),ktrajc(:,:)
INTEGER :: istatus, cur_level
REAL :: tinc
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INTEGER :: i,j,k,ireturn
REAL :: time, time1, time2, amin, amax, degree, degree2radian, theta
LOGICAL :: iexist
REAL :: u1_changed, utem,vtem,wtem, factor
INTEGER :: nfile
REAL :: ttrajc0_in
INTEGER :: ntrajcs_in, ntrajcs_each_j
CHARACTER (LEN=80) :: timsnd1,timsnd2,timsnd0,timsnd3
INTEGER :: tmstrln1,tmstrln2,tmstrln0, nunit(nmax_times),ltrajc_fn,istat,tmstrln3
INTEGER :: kl, notinteg, nunit1(nmax_times)
CHARACTER (LEN=256) :: trajc_fn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(/9(/5x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to ARPSTRAJC ###', &
'### ###', &
'###############################################################', &
'###############################################################'
!
!-----------------------------------------------------------------------
! Get the names of the input data files.
!-----------------------------------------------------------------------
!
READ(5,history_data,ERR=100)
WRITE(6,'(a,i3)') 'Input hdmpinopt=', hdmpinopt
WRITE(6,'(a,i3)') 'Input hinfmt =', hinfmt
WRITE(6,'(a,i3)') 'Input nhisfile =', nhisfile
IF( hinfmt /= 3 ) then
Print*,'Only HDF format (3) is supported.'
STOP
ENDIF
WRITE(6,'(a)')'Namelist history_data was successfully read.'
IF( hdmpinopt == 1 ) THEN
CALL gthinfns
(hdmpfheader,hdmpftrailer,hinfmt, &
tintv_dmpin, tbgn_dmpin, tend_dmpin, &
grdbasfn,hisfile,nhisfile)
END IF
IF( nhisfile > nhisfile_max ) THEN
WRITE(6,'(/a,a,I5,/a/)') &
'The number of history files to be processed', &
' exceeded the maximum number ',nhisfile_max, &
'Please increase the size of array hisfile.'
END IF
lengbf =len_trim( grdbasfn )
WRITE(6,'(/a,a)')'The grid base-state file name is ', &
grdbasfn(1:lengbf)
DO nfile=1,nhisfile
WRITE(6,'(a,i3,a,a)') &
'History file No.',nfile,' is ',trim(hisfile(nfile))
END DO
initrajc = 1
reftime = 0.0
trajcfn_header = ' '
nzctr = 1
ntimes = 1
ttrajc0=0.0
ntrajc0=1
xtrajc0(:,:)=0.0
ytrajc0(:,:)=0.0
ztrajc0(:,:)=0.0
READ(5,trajectories, END=100)
WRITE(6,'(/a,a/)') 'NAMELIST block trajectories successfully read.'
trajc_fn_in(:)='default_trajecotry_filename'
write(6,*) 'Print out namelist block trajectories parameters read in.'
write(6,trajectories)
IF( ntimes > nmax_times ) then
print*,'ntimes in input file exceeded max allowed value of ', nmax_times
STOP
ENDIF
IF( initrajc == 1 .and. ntrajcs > nmax_trajcs ) then
Print*,'ntrajcs larger than maximum allowed. Increase nmax_trajcs in program to ',ntrajcs
STOP
endif
degree2radian = 4.0*atan(1.0)/180.0
IF( initrajc == 2 ) then ! circle
ntrajcs_each_j = nint( 360.0/theta_inc ) + 1
ntrajcs = nzctr * ntrajcs_each_j
IF( ntrajcs > nmax_trajcs ) then
Print*,'ntrajcs larger than maximum allowed. Increase nmax_trajcs in program to ',ntrajcs
STOP
endif
print*,'ntrajcs,nzctr,initrajc=',ntrajcs,nzctr,initrajc
do k=1,ntimes
do j=1,nzctr
xtrajc0((j-1)*(ntrajcs_each_j)+1,k)= xctr
ytrajc0((j-1)*(ntrajcs_each_j)+1,k)= yctr
ztrajc0((j-1)*(ntrajcs_each_j)+1,k)= zctr(j)
do i=1,ntrajcs
theta = (i-1)*theta_inc * degree2radian
xtrajc0((j-1)*(ntrajcs_each_j)+(i+1),k)= xctr+radius*sin( theta )
ytrajc0((j-1)*(ntrajcs_each_j)+(i+1),k)= yctr+radius*cos( theta )
ztrajc0((j-1)*(ntrajcs_each_j)+(i+1),k)= zctr(j)
enddo
enddo
enddo
endif
IF( initrajc == 3 ) then ! box
ntrajcs = nint( 360.0/theta_inc )
IF( ntrajcs > nmax_trajcs ) then
Print*,'ntrajcs larger than maximum allowed. Increase nmax_trajcs in program to ',ntrajcs
STOP
endif
do k=1,ntimes
do i=1,ntrajcs
degree = (i-1)*theta_inc
theta = degree * degree2radian
if( degree <= 45.0 ) then
xtrajc0(i,k)= xctr+radius
ytrajc0(i,k)= yctr+radius*tan( theta )
elseif( degree <= 135.0 ) then
xtrajc0(i,k)= xctr+radius*cos( theta )
ytrajc0(i,k)= yctr+radius
elseif( degree <= 225.0 ) then
xtrajc0(i,k)= xctr-radius
ytrajc0(i,k)= yctr+radius*tan( theta )
elseif( degree <= 315.0 ) then
xtrajc0(i,k)= xctr+radius*cos( theta )
ytrajc0(i,k)= yctr-radius
else
xtrajc0(i,k)= xctr+radius
ytrajc0(i,k)= yctr+radius*tan( theta )
endif
ztrajc0(i,k)=zctr(1)
enddo
enddo
endif
print*,'I am here. initrajc, ntrajcs=', initrajc, ntrajcs
IF( initrajc == 4 ) then ! read from file
DO k=1,ntimes
CALL getunit(nunit(k))
OPEN(UNIT=nunit(k),FILE=trim(trajc_fn_in(k)),STATUS='old',FORM='formatted')
READ(nunit(k),*)
READ(nunit(k),*)
READ(nunit(k),*)
READ(nunit(k),*)
READ(nunit(k),*)
READ(nunit(k),*)
READ(nunit(k),'(4e17.6)') ttrajc0_in
READ(nunit(k),'(i10)') ntrajcs_in
print*,'ntrajcs_in=', ntrajcs_in
IF( ntrajcs_in /= ntrajcs ) then
print*,'ntrajcs read in .ne. ntrajcs in program. Value in data used.'
ntrajcs = ntrajcs_in
ENDIF
READ(nunit(k),'(6e17.6)') ((xtrajc0(i,k),ytrajc0(i,k),ztrajc0(i,k)),i=1,ntrajcs)
CLOSE(nunit(k))
call retunit(nunit(k))
ENDDO
ENDIF
DO k=1,ntimes
print*,'I am here now. ntrajcs=',ntrajcs,' No. of times=', k
write(6,'(6e17.6)') ((xtrajc0(i,k),ytrajc0(i,k),ztrajc0(i,k)),i=1,ntrajcs)
ENDDO
IF( hdmpinopt == 1 ) THEN
CALL gthinfns(hdmpfheader,hdmpftrailer,hinfmt, &
tintv_dmpin, tbgn_dmpin, tend_dmpin, &
grdbasfn,hisfile,nhisfile)
END IF
IF( nhisfile > nhisfile_max ) THEN
WRITE(6,'(/a,a,I5,/a/)') &
'The number of history files to be processed', &
' exceeded the maximum number ',nhisfile_max, &
'Please increase the size of array hisfile.'
END IF
lengbf =len_trim( grdbasfn )
WRITE(6,'(/a,a)')'The grid base-state file name is ', &
grdbasfn(1:lengbf)
DO nf=1,nhisfile
WRITE(6,'(a,i3,a,a)') &
'History file No.',nf,' is ',trim(hisfile(nf))
END DO
npoints = nhisfile
lengbf = len_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
Print*,'nx,ny,nz of input data were ', nx,ny,nz
allocate(x(nx ),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:x")
x = 0.0
allocate(y(ny ),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:y")
y = 0.0
allocate(z(nz ),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:z")
z = 0.0
allocate(xs(nx),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:xs")
xs = 0.0
allocate(ys(ny),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ys")
ys = 0.0
allocate(u(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:u")
u=0.0
allocate(v(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:v")
v=0.0
allocate(w(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:w")
w=0.0
allocate(zp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:zp")
zp=0.0
allocate(itmp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:itmp")
itmp=0
print*,'npoints,ntrajcs=', npoints,ntrajcs
allocate(xtrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:xtrajc")
allocate(ytrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ytrajc")
allocate(ztrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ztrajc")
allocate(utrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:utrajc")
allocate(vtrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vtrajc")
allocate(wtrajc(npoints,ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:wtrajc")
! allocate(xweight(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:xweight")
! allocate(yweight(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:yweight")
! allocate(zweight(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:zweight")
! allocate(itrajc(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:itrajc")
! allocate(jtrajc(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:jtrajc")
! allocate(ktrajc(npoints,ntrajcs),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:ktrajc")
allocate(ttrajc(npoints),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ttrajc")
if( hdmpinopt == 1 ) then ! determine ntrajc0 from ttrajc0
! ttrajc0 is known
ntrajc0 = nint( (ttrajc0 - tbgn_dmpin)/tintv_dmpin ) + 1
else
! ntrajc0 is known from namelist input
! Get ttrajc0 from NO. ntrajc0 history data
endif
print*,'ntrajc0 =', ntrajc0
DO k=1,ntimes
DO j=1,ntrajcs
xtrajc(ntrajc0,j,k)=xtrajc0(j,k)
ytrajc(ntrajc0,j,k)=ytrajc0(j,k)
ztrajc(ntrajc0,j,k)=ztrajc0(j,k)
ENDDO
ENDDO
CALL get_gridinfo_from_hdf(grdbasfn(1:lengbf),nx,ny,nz,x,y,z,ireturn)
Print*,'x,y,z of input data read in.'
print*,'x(1 )=',x(1)
print*,'x(nx)=',x(nx)
print*,'y(1 )=',y(1)
print*,'y(ny)=',y(ny)
dx = x(2) - x(1)
dy = y(2) - y(1)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz ,'nzsoil=',nzsoil
!-----------------------------------------------------------------------
! Find start and end times for trajectories from data
! for hdmpinopt =2.
!-----------------------------------------------------------------------
IF( hdmpinopt == 1 ) then
ttrajc(1) = tbgn_dmpin
ttrajc(npoints) = tend_dmpin
ELSE
nfile = nhisfile
CALL dtareadwind(nx,ny,nz,x,y,z,zp, trim(grdbasfn),trim(hisfile(nfile)), &
time,u,v,w, itmp)
ttrajc(npoints) = time
nfile = 1
print*,'calling dtareadwind to read file ', trim(hisfile(nfile))
CALL dtareadwind(nx,ny,nz,x,y,z,zp, trim(grdbasfn),trim(hisfile(nfile)), &
time,u,v,w, itmp)
ttrajc(1) = time
ENDIF
!-----------------------------------------------------------------------
! Read history data at the starting (reference) time of trajectories.
!-----------------------------------------------------------------------
nfile = ntrajc0
WRITE(6,'(/a,a,a)') ' Data set ', trim(hisfile(nfile)),' to be read.'
CALL dtareadwind(nx,ny,nz,x,y,z,zp, trim(grdbasfn),trim(hisfile(nfile)), &
time,u,v,w, itmp)
ttrajc0 = time
time1 = time
ttrajc(nfile)=time
print*,'ntrajc0, ttrajc0=', ntrajc0, ttrajc0
print*,'At time=', ttrajc0, 'xtrajc(1),ytrajc(1),ztrajc(1) =', &
xtrajc(nfile,1,1),ytrajc(nfile,1,1),ztrajc(nfile,1,1)
print*,'Data at Time =', time, ' read.'
write(0,*) 'before a3dmax0u'
CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'umin = ', amin,', umax =',amax
write(0,*) 'before a3dmax0v'
CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vmin = ', amin,', vmax =',amax
write(0,*) 'before a3dmax0w'
CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
print*,'Data at Time =', time, ' read.'
write(6,*) 'here'
DO i=1,nx-1
xs(i) = 0.5*(x(i)+x(i+1))
ENDDO
DO j=1,ny-1
ys(j) = 0.5*(y(j)+y(j+1))
ENDDO
write(6,*) 'ntimes = ',ntimes
!-----------------------------------------------------------------------
! Open trajectory data file and write header information
!-----------------------------------------------------------------------
IF( trajcfn_header == ' ' ) trajcfn_header = runname
DO k=1,ntimes
CALL cvttsnd( ttrajc0 , timsnd0, tmstrln0 )
CALL cvttsnd( ttrajc(1) , timsnd1, tmstrln1 )
CALL cvttsnd( ttrajc(npoints), timsnd2, tmstrln2 )
CALL cvttsnd( reftime(k), timsnd3, tmstrln3 )
trajc_fn = trim(trajcfn_header)//'.trajc_'//trim(timsnd1)//'-'//trim(timsnd2)//'_'//trim(timsnd3)
CALL getunit( nunit(k))
ltrajc_fn = len_trim(trajc_fn)
CALL fnversn( trajc_fn, ltrajc_fn )
OPEN(UNIT=nunit(k),FILE=trim(trajc_fn),STATUS='unknown', &
FORM='formatted',IOSTAT= istat )
CALL getunit( nunit1(k))
OPEN(UNIT=nunit1(k),FILE=trim(trajc_fn)//'.backward',STATUS='unknown', &
FORM='formatted',IOSTAT= istat )
WRITE(nunit(k),'(a)') trim(trajcfn_header)
WRITE(nunit(k),'(6e17.6)') x(2),x(nx-1),y(2),y(ny-1),z(2),z(nz-1)
WRITE(nunit(k),'(3e17.6)') dx, dy, z(2)-z(1)
WRITE(nunit(k),'(3e17.6)') ttrajc(1),ttrajc0,ttrajc(npoints)
WRITE(nunit(k),'(i10)') npoints
WRITE(nunit(k),'(i10)') ntrajcs
WRITE(nunit1(k),'(a)') trim(trajcfn_header)
WRITE(nunit1(k),'(6e17.6)') x(2),x(nx-1),y(2),y(ny-1),z(2),z(nz-1)
WRITE(nunit1(k),'(3e17.6)') dx, dy, z(2)-z(1)
WRITE(nunit1(k),'(3e17.6)') ttrajc(1),ttrajc0,ttrajc(npoints)
WRITE(nunit1(k),'(i10)') npoints
WRITE(nunit1(k),'(i10)') ntrajcs
notinteg = 1
tinc = 1.0 ! should not be used because notinteg = 1
cur_level = ntrajc0
CALL calc_trajc(nx,ny,nz,xs,ys,zp,u,v,w,tinc,cur_level, &
xtrajc(1,1,k),ytrajc(1,1,k),ztrajc(1,1,k), &
utrajc(1,1,k),vtrajc(1,1,k),wtrajc(1,1,k), &
! itrajc,jtrajc,ktrajc,xweight,yweight,zweight, &
npoints,ntrajcs, notinteg)
j=ntrajc0
write(nunit1(k),'(4e17.6)') ttrajc(j)
write(nunit1(k),'(i10)') ntrajcs
write(nunit1(k),'(6e17.6)') ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
ENDDO
!-----------------------------------------------------------------------
! Calculate backward trajectories.
!-----------------------------------------------------------------------
time1 = ttrajc(ntrajc0)
DO nfile = ntrajc0-1,1,-1
WRITE(6,'(/a,a,a)') ' Data set ', trim(hisfile(nfile)),' to be read.'
CALL dtareadwind(nx,ny,nz,x,y,z,zp, trim(grdbasfn),trim(hisfile(nfile)), &
time,u,v,w, itmp)
time2 = time
ttrajc(nfile)=time
print*,'Data at Time =', time, ' read.'
CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'umin = ', amin,', umax =',amax
CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vmin = ', amin,', vmax =',amax
CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
write(0,*) 'here'
!-----------------------------------------------------------------------
! Calculate trajectories
!-----------------------------------------------------------------------
write(0,*) 'calling calc_trajc ',time2, time1
tinc = time2-time1
cur_level = nfile+1
notinteg = 0
DO k=1,ntimes
write(0,*) 'calling calc_trajc'
CALL calc_trajc(nx,ny,nz,xs,ys,zp,u,v,w,tinc,cur_level, &
xtrajc(1,1,k),ytrajc(1,1,k),ztrajc(1,1,k), &
utrajc(1,1,k),vtrajc(1,1,k),wtrajc(1,1,k), &
! itrajc,jtrajc,ktrajc,xweight,yweight,zweight, &
npoints,ntrajcs, notinteg)
print*,'At time=', time2, 'xtrajc(1),ytrajc(1) =', xtrajc(nfile,1,1),ytrajc(nfile,1,1)
time1 = time2
j = nfile
write(nunit1(k),'(4e17.6)') ttrajc(j)
write(nunit1(k),'(i10)') ntrajcs
write(nunit1(k),'(6e17.6)') ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
ENDDO
ENDDO
DO k=1,ntimes
!-----------------------------------------------------------------------
! Write out the backward trajectory points into the combined trajectory file
!-----------------------------------------------------------------------
do j=1,ntrajc0
write(nunit(k),'(4e17.6)') ttrajc(j)
write(nunit(k),'(i10)') ntrajcs
write(nunit(k),'(6e17.6)') ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
enddo
CLOSE(nunit1(k),status='delete') ! No need for this file anymore
CALL retunit(nunit1(k))
ENDDO
!-----------------------------------------------------------------------
! Calculate forward trajectories
! First re-read data at starting (reference) time if necessary.
!-----------------------------------------------------------------------
time1 = ttrajc(ntrajc0)
DO nfile = ntrajc0+1,npoints,+1
WRITE(6,'(/a,a,a)') ' Data set ', trim(hisfile(nfile)),' to be read.'
CALL dtareadwind(nx,ny,nz,x,y,z,zp, trim(grdbasfn),trim(hisfile(nfile)), &
time,u,v,w, itmp)
time2 = time
ttrajc(nfile)=time
print*,'Data at Time =', time, ' read.'
CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'umin = ', amin,', umax =',amax
CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vmin = ', amin,', vmax =',amax
CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
write(0,*) 'here'
!-----------------------------------------------------------------------
! Calculate trajectories
!-----------------------------------------------------------------------
write(0,*) 'calling calc_trajc ',time2, time1
tinc = time2-time1
cur_level = nfile-1
notinteg = 0
! print*,'cur_level =', cur_level
! print*,'xtrajc(cur_level,:) =', xtrajc(cur_level,1), xtrajc(cur_level,2)
DO k=1,ntimes
write(0,*) 'calling calc_trajc'
CALL calc_trajc(nx,ny,nz,xs,ys,zp,u,v,w,tinc,cur_level, &
xtrajc(1,1,k),ytrajc(1,1,k),ztrajc(1,1,k), &
utrajc(1,1,k),vtrajc(1,1,k),wtrajc(1,1,k), &
! itrajc,jtrajc,ktrajc,xweight,yweight,zweight, &
npoints,ntrajcs, notinteg)
time1 = time2
print*,'At time=', time2, 'xtrajc(1),ytrajc(1),ztrajc(1) =', &
xtrajc(nfile,1,k),ytrajc(nfile,1,k),ztrajc(nfile,1,k)
!-----------------------------------------------------------------------
! Write out the forward trajectory points
!-----------------------------------------------------------------------
j = nfile
write(nunit(k),'(4e17.6)') ttrajc(j)
write(nunit(k),'(i10)') ntrajcs
write(nunit(k),'(6e17.6)') ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
ENDDO
ENDDO
!-----------------------------------------------------------------------
! Close the trajectory file
!-----------------------------------------------------------------------
DO k=1,ntimes
CLOSE(UNIT=nunit(k))
CALL retunit( nunit(k))
ENDDO
STOP
100 WRITE(6,'(a)') 'Error reading NAMELIST file. Program ARPSTRAJC stopped.'
STOP
END PROGRAM ARPSTRAJC
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DTAREADWIND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dtareadwind(nx,ny,nz,x,y,z,zp, grdbasfn,datafn, time,u,v,w, itmp) 5,16
!
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: hinfmt ! The format of the history data dump
CHARACTER(LEN=*) :: grdbasfn ! Name of the grid/base state array file
CHARACTER(LEN=*) :: datafn ! Name of the other time dependent data file
REAL :: x(nx),y(ny),z(nz)
REAL :: time
REAL :: u(nx,ny,nz)
REAL :: v(nx,ny,nz)
REAL :: w(nx,ny,nz)
REAL :: zp(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
INTEGER :: grdread,iread
SAVE grdread
INTEGER :: istat
INTEGER :: ireturn ! Return status indicator
INTEGER :: grdbas ! Wether this is a grid/base state
! array dump
INTEGER :: i,j,k
LOGICAL :: fexist
INTEGER :: is
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'mp.inc' ! mpi parameters.
INCLUDE 'exbc.inc'
INCLUDE 'phycst.inc'
DATA grdread /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ireturn = 0
hinfmt = 3
!
!-----------------------------------------------------------------------
!
! Open and read grid and base state data file depending on the
! values of parameters grdin and basin, which are read in from the
! time dependent data set. If grdin or basin is zero, the grid and
! base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
Print*,'In dtareadwind, grdbasfn =', trim(grdbasfn)
Print*,'In dtareadwind, datafn =', trim(datafn)
IF( grdread == 0 ) THEN
! print*,'grdread inside if block=', grdread
grdbas = 1
INQUIRE(FILE=grdbasfn, EXIST = fexist )
IF( fexist ) GO TO 200
INQUIRE(FILE=grdbasfn//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn//'.Z' )
GO TO 200
END IF
INQUIRE(FILE=grdbasfn//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn//'.gz' )
GO TO 200
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//grdbasfn// &
' or its compressed version not found.', &
'Program stopped in DTAREADWIND location 1.'
CALL arpsstop('arpsstop called from dtareadwind during base state read',1)
200 CONTINUE
!
!-----------------------------------------------------------------------
! Read grid and base state fields.
!-----------------------------------------------------------------------
!
CALL hdfreadwind(nx,ny,nz,x,y,z,zp,grdbas,trim(grdbasfn), &
time,u,v,w, itmp)
grdread = 1
END IF
!
!-----------------------------------------------------------------------
! Read data fields.
!-----------------------------------------------------------------------
!
grdbas = 0
INQUIRE(FILE=trim(datafn), EXIST = fexist )
IF( fexist ) GO TO 100
INQUIRE(FILE=trim(datafn)//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn)//'.Z' )
GO TO 100
END IF
INQUIRE(FILE=trim(datafn)//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn)//'.gz' )
GO TO 100
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//trim(datafn) &
//' or its compressed version not found.', &
'Program stopped in DTAREADWIND location 2.'
CALL arpsstop('arpsstop called from dtareadwind during base read-2',1)
100 CONTINUE
CALL hdfreadwind(nx,ny,nz,x,y,z,zp,grdbas,trim(datafn), &
time,u,v,w, itmp)
RETURN
END SUBROUTINE dtareadwind
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreadwind(nx,ny,nz,x,y,z,zp,grdbas,filename, & 4,122
time,u,v,w, itmp)
!-----------------------------------------------------------------------
! PURPOSE:
! Read in history data in the NCSA HDF4 format.
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 2000/04/15
!
! MODIFICATION HISTORY:
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! grdbas Data read flag.
! =1, only grid and base state arrays will be read
! =0, all arrays will be read based on data
! parameter setting.
! filename Character variable nhming the input HDF file
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: grdbas
CHARACTER (LEN=*) :: filename
REAL :: x (nx) ! x coord.
REAL :: y (ny) ! y coord.
REAL :: z (nz) ! z coord.
REAL :: time
REAL :: u(nx,ny,nz)
REAL :: v(nx,ny,nz)
REAL :: w(nx,ny,nz)
REAL :: zp(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmax(nz), hmin(nz) ! Temporary array
INTEGER :: ireturn
!-----------------------------------------------------------------------
! Parameters describing routine that wrote the gridded data
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version.
! intver??: an integer to allow faster comparison than fmtver??,
! which are strings.
!
! Verion 5.00: significant change in soil variables since version 4.10.
!
!-----------------------------------------------------------------------
CHARACTER (LEN=40) :: fmtver410,fmtver500
INTEGER :: intver,intver410,intver500
PARAMETER (fmtver410='004.10 HDF4 Coded Data',intver410=410)
PARAMETER (fmtver500='005.00 HDF4 Coded Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
!-----------------------------------------------------------------------
! Misc. local variables
!-----------------------------------------------------------------------
INTEGER :: lchanl
PARAMETER (lchanl=6) ! Channel number for formatted printing.
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin
INTEGER :: bgrdin,bbasin,bvarin,bicein,btrbin,btkein
INTEGER :: istat, sd_id
INTEGER :: varflg, istatus
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'indtflg.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Beginning of executable code...
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(*,*) 'HDFREAD: Reading HDF file: ', trim(filename)
!-----------------------------------------------------------------------
! Open file for reading
!-----------------------------------------------------------------------
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "HDFREAD: ERROR opening ", trim(filename)," for reading."
GO TO 110
END IF
fmtverin = fmtver500
WRITE(6,'(/1x,a,a/)') 'Incoming data format, fmtverin=',fmtverin
CALL hdfrdc(sd_id,40,"runname",runname,istat)
CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
IF( nocmnt > 0 ) THEN
CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
END IF
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') trim(runname)
WRITE (6,*) "Comments:"
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE (6,*) " "
CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
CALL hdfrdr(sd_id,"time",time,istat)
!-----------------------------------------------------------------------
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdi(sd_id,"nz",nzin,istat)
IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') ' Dimensions in HDFREAD inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz
WRITE(6,'(1x,a)') ' Program aborted in HDFREAD.'
CALL arpsstop('arpsstop called from HDFREAD due to nxin...',1)
END IF
!-----------------------------------------------------------------------
! Read in x,y and z at grid cell centers (scalar points).
!-----------------------------------------------------------------------
IF( grdin == 1 .OR. grdbas == 1 ) THEN
CALL hdfrd1d(sd_id,"x",nx,x,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"y",ny,y,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"z",nz,z,istat)
IF (istat /= 0) GO TO 110
END IF ! grdin
!-----------------------------------------------------------------------
! Read in flags for different data groups
!-----------------------------------------------------------------------
IF ( grdbas == 1 ) THEN ! Read grid and base state arrays
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)') &
'To read grid and base state data at time ', time, &
' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",bgrdin,istat)
CALL hdfrdi(sd_id,"basflg",bbasin,istat)
CALL hdfrdi(sd_id,"varflg",bvarin,istat)
CALL hdfrdi(sd_id,"mstflg",mstin,istat)
CALL hdfrdi(sd_id,"iceflg",bicein,istat)
CALL hdfrdi(sd_id,"trbflg",btrbin,istat)
CALL hdfrdi(sd_id,"landflg",landin,istat)
CALL hdfrdi(sd_id,"totflg",totin,istat)
CALL hdfrdi(sd_id,"tkeflg",btkein,istat)
ELSE ! Normal data reading
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', &
time,' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",grdin,istat)
CALL hdfrdi(sd_id,"basflg",basin,istat)
CALL hdfrdi(sd_id,"varflg",varin,istat)
CALL hdfrdi(sd_id,"mstflg",mstin,istat)
CALL hdfrdi(sd_id,"iceflg",icein,istat)
CALL hdfrdi(sd_id,"trbflg",trbin,istat)
CALL hdfrdi(sd_id,"sfcflg",sfcin,istat)
CALL hdfrdi(sd_id,"rainflg",rainin,istat)
CALL hdfrdi(sd_id,"totflg",totin,istat)
CALL hdfrdi(sd_id,"tkeflg",tkein,istat)
print*,'done reading parameters'
END IF
CALL hdfrdi(sd_id,"prcflg",prcin,istat)
CALL hdfrdi(sd_id,"radflg",radin,istat)
CALL hdfrdi(sd_id,"flxflg",flxin,istat)
CALL hdfrdi(sd_id,"snowflg",snowin,istat)
CALL hdfrdi(sd_id,"month",month,istat)
CALL hdfrdi(sd_id,"day",day,istat)
CALL hdfrdi(sd_id,"year",year,istat)
CALL hdfrdi(sd_id,"hour",hour,istat)
CALL hdfrdi(sd_id,"minute",minute,istat)
CALL hdfrdi(sd_id,"second",second,istat)
CALL hdfrdr(sd_id,"umove",umove,istat)
CALL hdfrdr(sd_id,"vmove",vmove,istat)
CALL hdfrdr(sd_id,"xgrdorg",xgrdorg,istat)
CALL hdfrdr(sd_id,"ygrdorg",ygrdorg,istat)
CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
CALL hdfrdr(sd_id,"trulon",trulon,istat)
CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
CALL hdfrdr(sd_id,"tstop",tstop,istat)
CALL hdfrdr(sd_id,"thisdmp",thisdmp,istat)
CALL hdfrdr(sd_id,"latitud",latitud,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)
IF( grdin == 1 .OR. grdbas == 1 ) THEN
CALL hdfrd3d(sd_id,"zp",nx,ny,nz,zp,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
ENDIF
!-----------------------------------------------------------------------
! Read in base state fields
!-----------------------------------------------------------------------
! Print*,'start doing 3d arrays'
IF( (basin == 1 .OR. grdbas == 1) .and. hdmpfmt /= 0) THEN
! CALL hdfrd3d(sd_id,"ubar",nx,ny,nz,ubar,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"vbar",nx,ny,nz,vbar,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"ptbar",nx,ny,nz,ptbar,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"pbar",nx,ny,nz,pbar,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! IF( mstin == 1 ) THEN
! CALL hdfrd3d(sd_id,"qvbar",nx,ny,nz,qvbar,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! ELSE
! qvbar= 0.0
! END IF
ENDIF
landout = 0 ! skipping landout for know
IF( grdbas == 1 ) GO TO 930
IF( varin == 1 ) then
!-----------------------------------------------------------------------
! Read in total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,"u",nx,ny,nz,u,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
CALL hdfrd3d(sd_id,"v",nx,ny,nz,v,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
CALL hdfrd3d(sd_id,"w",nx,ny,nz,w,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"pt",nx,ny,nz,pt,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"p",nx,ny,nz,p,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! IF( mstin == 1 ) THEN
! CALL hdfrd3d(sd_id,"qv",nx,ny,nz,qv,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
!
! CALL hdfrd3d(sd_id,"qc",nx,ny,nz,qc,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
!
! CALL hdfrd3d(sd_id,"qr",nx,ny,nz,qr,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! IF( icein == 1 ) THEN
! CALL hdfrd3d(sd_id,"qi",nx,ny,nz,qi,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
!
! CALL hdfrd3d(sd_id,"qs",nx,ny,nz,qs,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
!
! CALL hdfrd3d(sd_id,"qh",nx,ny,nz,qh,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! END IF
! END IF
END IF
! IF( tkein == 1 ) THEN
! CALL hdfrd3d(sd_id,"tke",nx,ny,nz,tke,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! END IF
! IF( trbin == 1 ) THEN
! CALL hdfrd3d(sd_id,"kmh",nx,ny,nz,kmh,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"kmv",nx,ny,nz,,istat,itmp,hmax,hmin)
! IF (istat /= 0) GO TO 110
! END IF
!
! Surface arrays skipped
!
!-----------------------------------------------------------------------
!
! Friendly exit message
!
!-----------------------------------------------------------------------
930 CONTINUE
WRITE(6,'(/a,F8.1,a/)') &
' Data at time=', time/60,' (min) were successfully read.'
ireturn = 0
GO TO 130
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in HDFREAD'
STOP
ireturn=1
130 CONTINUE
CALL hdfclose(sd_id,istat)
IF (ireturn == 0) THEN
IF (istat == 0) THEN
WRITE(6,'(/a/a)') &
"HDFREADWIND: Successfully read ", trim(filename)
ELSE
WRITE(6,'(/a,i3,a/,a)') &
"HDFREADWIND: ERROR (status=", istat, ") closing ", trim(filename)
STOP
END IF
END IF
RETURN
END SUBROUTINE hdfreadwind
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GET_GRIDINFO_FROM_HDF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE get_gridinfo_from_hdf(filename,nx,ny,nz,x,y,z,ireturn) 2,15
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! filename Channel number for binary reading.
!
! OUTPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: stat, sd_id
CHARACTER (LEN=*) :: filename
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: x(nx),y(ny),z(nz)
INTEGER :: ireturn ! Return status indicator
INTEGER istat
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "get_gridinfo_from_hdf: ERROR opening ", &
trim(filename)," for reading."
GO TO 110
ELSE
WRITE(6,*) 'File ',filename,' openned.'
END IF
! print*,'sd_id, nx =', sd_id, nx
CALL hdfrd1d(sd_id,"x",nx,x,istat)
! print*,'istat after reading x =', istat
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"y",ny,y,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"z",nz,z,istat)
IF (istat /= 0) GO TO 110
ireturn = 0
GO TO 130
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in GET_GRIDINFO_FROM_HDF.'
print*,'Program stopped in GET_GRIDINFO_FROM_HDF.'
STOP
ireturn=1
130 CONTINUE
!tmp stat = sfendacc(sd_id) ! is this necessary?
CALL hdfclose(sd_id,stat)
RETURN
END SUBROUTINE get_gridinfo_from_hdf
SUBROUTINE copyarray(arrayin,nx,ny,nz, arrayout, nx1,ny1,nz1, & 32
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
INTEGER :: nx,ny,nz, nx1,ny1,nz1,nxbgn,nxend,nybgn,nyend,nzbgn,nzend
REAL :: arrayin(nx,ny,nz)
REAL :: arrayout(nx1,ny1,nz1)
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
arrayout(i,j,k)=arrayin(nxbgn+i-1,nybgn+j-1,nzbgn+k-1)
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE copyarray
SUBROUTINE intrpx3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgtx,ix, & 3
aout,nx1,is1,ie1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the first dimension
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nx1,is1,ie1
REAL :: ain (nx ,ny,nz)
REAL :: aout(nx1,ny,nz)
REAL :: wgtx(nx1)
INTEGER :: ix(nx1)
INTEGER :: i,j,k
DO k=ks ,ke
DO j=js ,je
DO i=is1,ie1
aout(i,j,k)= wgtx(i) *ain(ix(i) ,j,k) &
+(1.0-wgtx(i))*ain(ix(i)+1,j,k)
END DO
END DO
END DO
RETURN
END SUBROUTINE intrpx3d
SUBROUTINE intrpy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgty,jy, & 3
aout,ny1,js1,je1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the second dimension
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: ny1,js1,je1
REAL :: ain (nx,ny ,nz)
REAL :: aout(nx,ny1,nz)
REAL :: wgty(ny1)
INTEGER :: jy(ny1)
INTEGER :: i,j,k
DO k=ks ,ke
DO j=js1,je1
DO i=is ,ie
aout(i,j,k)= wgty(j) *ain(i,jy(j) ,k) &
+(1.0-wgty(j))*ain(i,jy(j)+1,k)
END DO
END DO
END DO
RETURN
END SUBROUTINE intrpy3d
SUBROUTINE intrpxy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, & 53,8
wgtx,ix,wgty,jy, &
aout,nx1,is1,ie1, ny1,js1,je1, &
temx1yz)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the first and second dimensions
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nx1,is1,ie1, ny1,js1,je1
REAL :: ain (nx ,ny ,nz)
REAL :: aout(nx1,ny1,ks:ke)
REAL :: wgtx(nx1),wgty(ny1)
INTEGER :: ix(nx1),jy(ny1)
INTEGER :: k
REAL :: temx1yz(nx1,ny)
DO k=ks,ke
CALL intrpx3d(ain(1,1,k),nx,is,ie, ny,js,je, 1,1,1, &
wgtx,ix,temx1yz,nx1,is1,ie1)
CALL intrpy3d(temx1yz,nx1,is1,ie1, ny,js,je, 1,1,1, &
wgty,jy,aout(1,1,k),ny1,js1,je1)
ENDDO
CALL edgfill(aout,1,nx1,is1,ie1,1,ny1,js1,je1,ks,ke,ks,ke)
RETURN
END SUBROUTINE intrpxy3d
SUBROUTINE calc_trajc(nx,ny,nz,xs,ys,zp,u,v,w,tinc,clevel, & 3
xtrajc,ytrajc,ztrajc,utrajc,vtrajc,wtrajc, &
! itrajc,jtrajc,ktrajc,xweight,yweight,zweight, &
ntrajc_points,ntrajcs, notinteg)
!
!-----------------------------------------------------------------------
! PURPOSE:
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! (4/08/2004)
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
! Variable Declarations:
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: xs(nx),ys(ny),zp(nx,ny,nz)
REAL :: u(nx,ny,nz)
REAL :: v(nx,ny,nz)
REAL :: w(nx,ny,nz)
INTEGER :: ntrajc_points,ntrajcs
REAL :: xtrajc(ntrajc_points,ntrajcs),ytrajc(ntrajc_points,ntrajcs)
REAL :: ztrajc(ntrajc_points,ntrajcs)
REAL :: utrajc(ntrajc_points,ntrajcs),vtrajc(ntrajc_points,ntrajcs),wtrajc(ntrajc_points,ntrajcs)
! REAL :: xweight(ntrajc_points,ntrajcs),yweight(ntrajc_points,ntrajcs),zweight(ntrajc_points,ntrajcs)
! INTEGER :: itrajc(ntrajc_points,ntrajcs),jtrajc(ntrajc_points,ntrajcs),ktrajc(ntrajc_points,ntrajcs)
REAL :: tinc
INTEGER :: clevel, flevel,get_ktrajc
INTEGER :: niteration,niteration_max,ntrajc,level_inc
REAL :: dx, dy, utrajc1,vtrajc1,wtrajc1
REAL :: uy1z1,uy2z1,uy1z2,uy2z2,vy1z1,vy2z1,vy1z2,vy2z2,wy1z1,wy2z1,wy1z2,wy2z2
REAL :: uz1,uz2,vz1,vz2,wz1,wz2
REAL :: missing_value
INTEGER :: notinteg ! for the first time step where xtrajc,ytrajc,ztrajc
! are known and no time integration of trajectory is performed.
INTEGER :: niterat,itrajc1,jtrajc1,ktrajc1
REAL :: xweight1,yweight1,zweight1
REAL :: zs1d(nz),zy1(nz),zy2(nz) ! automatic work array
INTEGER :: i,j,k, vlevel
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Begin of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
write(0,*) 'entering calc_trajc'
dx = xs(2)-xs(1)
dy = ys(2)-ys(1)
IF( tinc == 0 ) then
Print*,'tinc cannot be equal 0, job stopped in calc_trajc.'
STOP
ENDIF
! print*,'mark 1'
IF( tinc > 0 ) level_inc = 1
IF( tinc < 0 ) level_inc = -1
flevel = clevel + level_inc
if( notinteg == 1 ) flevel = clevel
niteration_max = 7
IF( notinteg == 1 ) THEN
niterat = 1
ELSE
niterat = niteration_max
ENDIF
! print*,'ntrajc_points,ntrajcs=', ntrajc_points,ntrajcs
! print*,'mark 2, flevel =', flevel, clevel, level_inc
xtrajc(flevel,:) = xtrajc(clevel,:)
ytrajc(flevel,:) = ytrajc(clevel,:)
ztrajc(flevel,:) = ztrajc(clevel,:)
! print*,'xtrajc(clevel,:), xtrajc(flevel,:)=', xtrajc(clevel,:), xtrajc(flevel,:)
! print*,'ytrajc(clevel,:), ytrajc(flevel,:)=', ytrajc(clevel,:), ytrajc(flevel,:)
! print*,'ztrajc(clevel,:), ztrajc(flevel,:)=', ztrajc(clevel,:), ztrajc(flevel,:)
! print*,' I am here'
missing_value = -99999.0
write(0,*) 'inside calc_trajc ',ntrajcs
DO ntrajc = 1,ntrajcs
write(0,*) 'ntrajc = ',ntrajc
IF( xtrajc(clevel,ntrajc) == missing_value .or. &
ytrajc(clevel,ntrajc) == missing_value .or. &
ztrajc(clevel,ntrajc) == missing_value ) cycle
DO niteration = 1, niterat
write(0,*) 'niteration = ',niteration, niterat
IF( niteration == 1 .and. notinteg /= 1 ) THEN
utrajc1 = utrajc(clevel,ntrajc)
vtrajc1 = vtrajc(clevel,ntrajc)
wtrajc1 = wtrajc(clevel,ntrajc)
ELSE
itrajc1 = MAX(1, MIN(nx-2, INT((xtrajc(flevel,ntrajc)-xs(1))/dx)+1 ))
jtrajc1 = MAX(1, MIN(ny-2, INT((ytrajc(flevel,ntrajc)-ys(1))/dy)+1 ))
write(0,*) 'xtrajc = ',flevel,xtrajc(flevel,ntrajc)
xweight1 = (xtrajc(flevel,ntrajc)-xs(itrajc1))/dx
yweight1 = (ytrajc(flevel,ntrajc)-ys(jtrajc1))/dy
! DTD test initial trajectory height
IF(notinteg == 1) THEN
print*,'Before: itrajc1,jtrajc1,ztrajc,zp(k=1)',itrajc1,jtrajc1,ztrajc(clevel,ntrajc),zp(itrajc1,jtrajc1,1)
! adjust initial trajectory height to be consistent with zp array, which is in ASL (user input should
! be interpreted as being AGL)
ztrajc(clevel,ntrajc)=ztrajc(clevel,ntrajc)+zp(itrajc1,jtrajc1,1)
print*,'After: itrajc1,jtrajc1,ztrajc,zp(k=1)',itrajc1,jtrajc1,ztrajc(clevel,ntrajc),zp(itrajc1,jtrajc1,1)
END IF
DO k=1,nz-1
zy1(k) = (1.0-xweight1)*(zp(itrajc1 ,jtrajc1 ,k)+zp(itrajc1 ,jtrajc1 ,k+1))*0.5 &
+xweight1 *(zp(itrajc1+1,jtrajc1 ,k)+zp(itrajc1+1,jtrajc1 ,k+1))*0.5
zy2(k) = (1.0-xweight1)*(zp(itrajc1 ,jtrajc1+1,k)+zp(itrajc1 ,jtrajc1+1,k+1))*0.5 &
+xweight1 *(zp(itrajc1+1,jtrajc1+1,k)+zp(itrajc1+1,jtrajc1+1,k+1))*0.5
zs1d(k) = ( 1.0-yweight1 )*zy1(k) + yweight1*zy2(k)
END DO
write(0,*) ' here 2'
ktrajc1 = get_ktrajc(ztrajc(flevel,ntrajc),zs1d,nz-1)
! print*,'niteration =', niteration
! print*,'xtrajc, ytrajc, ztrajc=', xtrajc(flevel,ntrajc),ytrajc(flevel,ntrajc),ztrajc(flevel,ntrajc)
! print*,'itrajc,jtrajc ,ktrajc =',itrajc1,jtrajc1 ,ktrajc1
uy1z1 = (1.0-xweight1)*(u(itrajc1 ,jtrajc1 ,ktrajc1 )+u(itrajc1+1,jtrajc1 ,ktrajc1 ))*0.5 &
+xweight1*(u(itrajc1+1,jtrajc1 ,ktrajc1 )+u(itrajc1+2,jtrajc1 ,ktrajc1 ))*0.5
uy2z1 = (1.0-xweight1)*(u(itrajc1 ,jtrajc1+1,ktrajc1 )+u(itrajc1+1,jtrajc1+1,ktrajc1 ))*0.5 &
+xweight1 *(u(itrajc1+1,jtrajc1+1,ktrajc1 )+u(itrajc1+2,jtrajc1+1,ktrajc1 ))*0.5
uy1z2 = (1.0-xweight1)*(u(itrajc1 ,jtrajc1 ,ktrajc1+1)+u(itrajc1+1,jtrajc1 ,ktrajc1+1))*0.5 &
+xweight1* (u(itrajc1+1,jtrajc1 ,ktrajc1+1)+u(itrajc1+2,jtrajc1 ,ktrajc1+1))*0.5
uy2z2 = (1.0-xweight1)*(u(itrajc1 ,jtrajc1+1,ktrajc1+1)+u(itrajc1+1,jtrajc1+1,ktrajc1+1))*0.5 &
+xweight1* (u(itrajc1+1,jtrajc1+1,ktrajc1+1)+u(itrajc1+2,jtrajc1+1,ktrajc1+1))*0.5
vy1z1 = (1.0-xweight1)*(v(itrajc1 ,jtrajc1 ,ktrajc1 )+v(itrajc1 ,jtrajc1+1,ktrajc1 ))*0.5 &
+xweight1* (v(itrajc1+1,jtrajc1 ,ktrajc1 )+v(itrajc1+1,jtrajc1+1,ktrajc1 ))*0.5
vy2z1 = (1.0-xweight1)*(v(itrajc1 ,jtrajc1+1,ktrajc1 )+v(itrajc1 ,jtrajc1+2,ktrajc1 ))*0.5 &
+xweight1* (v(itrajc1+1,jtrajc1+1,ktrajc1 )+v(itrajc1+1,jtrajc1+2,ktrajc1 ))*0.5
vy1z2 = (1.0-xweight1)*(v(itrajc1 ,jtrajc1 ,ktrajc1+1)+v(itrajc1 ,jtrajc1+1,ktrajc1+1))*0.5 &
+xweight1* (v(itrajc1+1,jtrajc1 ,ktrajc1+1)+v(itrajc1+1,jtrajc1+1,ktrajc1+1))*0.5
vy2z2 = (1.0-xweight1)*(v(itrajc1 ,jtrajc1+1,ktrajc1+1)+v(itrajc1 ,jtrajc1+2,ktrajc1+1))*0.5 &
+xweight1* (v(itrajc1+1,jtrajc1+1,ktrajc1+1)+v(itrajc1+1,jtrajc1+2,ktrajc1+1))*0.5
wy1z1 = (1.0-xweight1)*(w(itrajc1 ,jtrajc1 ,ktrajc1 )+w(itrajc1 ,jtrajc1 ,ktrajc1+1))*0.5 &
+xweight1* (w(itrajc1+1,jtrajc1 ,ktrajc1 )+w(itrajc1+1,jtrajc1 ,ktrajc1+1))*0.5
wy2z1 = (1.0-xweight1)*(w(itrajc1 ,jtrajc1+1,ktrajc1 )+w(itrajc1 ,jtrajc1+1,ktrajc1+1))*0.5 &
+xweight1* (w(itrajc1+1,jtrajc1+1,ktrajc1 )+w(itrajc1+1,jtrajc1+1,ktrajc1+1))*0.5
wy1z2 = (1.0-xweight1)*(w(itrajc1 ,jtrajc1 ,ktrajc1+1)+w(itrajc1 ,jtrajc1 ,ktrajc1+2))*0.5 &
+xweight1* (w(itrajc1+1,jtrajc1 ,ktrajc1+1)+w(itrajc1+1,jtrajc1 ,ktrajc1+2))*0.5
wy2z2 = (1.0-xweight1)*(w(itrajc1 ,jtrajc1+1,ktrajc1+1)+w(itrajc1 ,jtrajc1+1,ktrajc1+2))*0.5 &
+xweight1* (w(itrajc1+1,jtrajc1+1,ktrajc1+1)+w(itrajc1+1,jtrajc1+1,ktrajc1+2))*0.5
write(0,*) 'wylz2,wy2z2',wy1z2,wy2z2, xweight1
write(0,*) 'here 3 ',itrajc1, jtrajc1, ktrajc1, w(itrajc1 ,jtrajc1 ,ktrajc1+1)
uz1 = ( 1.0-yweight1 )*uy1z1 + yweight1*uy2z1
uz2 = ( 1.0-yweight1 )*uy1z2 + yweight1*uy2z2
vz1 = ( 1.0-yweight1 )*vy1z1 + yweight1*vy2z1
vz2 = ( 1.0-yweight1 )*vy1z2 + yweight1*vy2z2
wz1 = ( 1.0-yweight1 )*wy1z1 + yweight1*wy2z1
wz2 = ( 1.0-yweight1 )*wy1z2 + yweight1*wy2z2
write(0,*) 'here 4 ',ktrajc1,zs1d(ktrajc1+1), zs1d(ktrajc1)
zweight1 = (ztrajc(flevel,ntrajc)-zs1d(ktrajc1))/(zs1d(ktrajc1+1)-zs1d(ktrajc1))
write(0,*) 'here 4.1'
utrajc1 = (1.0-zweight1)*uz1+zweight1*uz2
write(0,*) 'here 4.2'
vtrajc1 = (1.0-zweight1)*vz1+zweight1*vz2
write(0,*) 'here 4.3', wz1, wz2
wtrajc1 = (1.0-zweight1)*wz1+zweight1*wz2
write(0,*) 'here 5'
ENDIF
! print*,'xweight,yweight,zweight=', xweight1,yweight1,zweight1
! print*,'tinc =', tinc
write(0,*) 'notinteg = ',notinteg
IF( notinteg /= 1 ) THEN
IF( xtrajc(clevel,ntrajc) .ne. missing_value ) then
xtrajc(flevel,ntrajc) = xtrajc(clevel,ntrajc) + tinc * 0.5*(utrajc(clevel,ntrajc)+utrajc1)
ELSE
xtrajc(flevel,ntrajc) = missing_value
ENDIF
IF( ytrajc(clevel,ntrajc) .ne. missing_value ) then
ytrajc(flevel,ntrajc) = ytrajc(clevel,ntrajc) + tinc * 0.5*(vtrajc(clevel,ntrajc)+vtrajc1)
ELSE
ytrajc(flevel,ntrajc) = missing_value
ENDIF
ztrajc(flevel,ntrajc) = ztrajc(clevel,ntrajc) + tinc * 0.5*(wtrajc(clevel,ntrajc)+wtrajc1)
! print*,'xtrajc, ytrajc, ztrajc=', xtrajc(flevel,ntrajc),ytrajc(flevel,ntrajc),ztrajc(flevel,ntrajc)
ENDIF
! print*,'clevel,flevel,niteration,ntrajc=',clevel,flevel,niteration,ntrajc
! print*,'utrajc0,vtrajc0,wtrajc0=', utrajc(clevel,ntrajc),vtrajc(clevel,ntrajc),wtrajc(clevel,ntrajc)
! print*,'utrajc1,vtrajc1,wtrajc1=', utrajc1,vtrajc1,wtrajc1
! print*,'xtrajc, ytrajc, ztrajc=', xtrajc(flevel,ntrajc),ytrajc(flevel,ntrajc),ztrajc(flevel,ntrajc)
ENDDO ! niteration
write(0,*) 'here 1'
utrajc(flevel,ntrajc) = utrajc1
vtrajc(flevel,ntrajc) = vtrajc1
wtrajc(flevel,ntrajc) = wtrajc1
! xweight(flevel,ntrajc) = xweight1 ! save the calculated weights in array for reuse
! yweight(flevel,ntrajc) = yweight1 ! save the calculated weights in array for reuse
! zweight(flevel,ntrajc) = zweight1 ! save the calculated weights in array for reuse
! itrajc(flevel,ntrajc) = itrajc1 ! save the calculated weights in array for reuse
! jtrajc(flevel,ntrajc) = jtrajc1 ! save the calculated weights in array for reuse
! ktrajc(flevel,ntrajc) = ktrajc1 ! save the calculated weights in array for reuse
IF( notinteg /= 1 ) THEN
ztrajc(flevel,ntrajc) = max(ztrajc(flevel,ntrajc), (zs1d(1)+zs1d(2))*0.5+0.1*(zs1d(2)-zs1d(1)))
ztrajc(flevel,ntrajc) = min(ztrajc(flevel,ntrajc), (zs1d(nz-1)+zs1d(nz-2))*0.5 )
if( xtrajc(flevel,ntrajc) .lt. xs(1) .or. xtrajc(flevel,ntrajc).gt. xs(nx-1) ) then
xtrajc(flevel,ntrajc) = missing_value
endif
if( ytrajc(flevel,ntrajc) .lt. ys(1) .or. ytrajc(flevel,ntrajc).gt. ys(ny-1) ) then
ytrajc(flevel,ntrajc) = missing_value
endif
ENDIF
ENDDO
RETURN
END SUBROUTINE calc_trajc
INTEGER FUNCTION get_ktrajc(z, zs1d, nz)
IMPLICIT NONE
INTEGER :: nz, k
REAL :: z, zs1d(nz)
IF( z < zs1d(1) ) then
get_ktrajc = 1
RETURN
ENDIF
IF( z >= zs1d(nz) ) then
get_ktrajc = nz-1
RETURN
ENDIF
DO k=1,nz-1
IF( z >= zs1d(k) .and. z < zs1d(k+1) ) then
get_ktrajc = k
EXIT
endif
ENDDO
RETURN
END FUNCTION get_ktrajc