PROGRAM ARPSCALCTRAJC ,105
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSTRAJC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! (4/08/2004)
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL, allocatable :: x (:),y (:),z (:)
REAL, allocatable :: xs (:), ys (:), zs1d(:) ! x,y coord for scalar points
REAL, allocatable :: zp(:,:,:), ptbar(:,:,:), pbar(:,:,:),qvbar(:,:,:),rhobar(:,:,:)
REAL, allocatable :: tem1 (:,:,:), tem2(:,:,:), tem3(:,:,:), tem4(:,:,:)
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
CHARACTER (LEN=256) :: grdbasfn,hisfile
INTEGER :: hdmpinopt, nf
CHARACTER (LEN=256) :: hdmpftrailer
CHARACTER (LEN=256) :: hdmpfheader
INTEGER :: ntimes
INTEGER, PARAMETER :: nmax_times=20
REAL :: tzero, tend, tstart_calc, tend_calc, tinc_calc, reftime(nmax_times)
CHARACTER(LEN=256) :: trajc_fn_in(nmax_times), trajc_fn_out
NAMELIST /input/ hdmpfheader,hdmpftrailer,tstart_calc,tend_calc,tinc_calc,reftime,trajc_fn_in,ntimes
INTEGER :: ntrajc0,ntrajcs, npoints, npoints_in
REAL, allocatable :: xtrajc(:,:,:),ytrajc(:,:,:),ztrajc(:,:,:),ttrajc(:)
REAL, allocatable :: xtrajc1(:,:),ytrajc1(:,:),ztrajc1(:,:)
REAL, allocatable :: utrajc(:,:),vtrajc(:,:),wtrajc(:,:)
REAL, allocatable :: pt_trajc(:,:),p_trajc(:,:),qv_trajc(:,:),q_trajc(:,:),ptprt_trajc(:,:)
REAL, allocatable :: xweight(:,:),yweight(:,:),zweight(:,:)
INTEGER, allocatable :: itrajc(:,:),jtrajc(:,:),ktrajc(:,:)
REAL, allocatable :: vortx_trajc(:,:),vorty_trajc(:,:),vortz_trajc(:,:), &
vorts_trajc(:,:),vortc_trajc(:,:)
REAL, allocatable :: buoy_trajc(:,:),buoyq_trajc(:,:),frcs_trajc(:,:),pprt_trajc(:,:)
REAL, allocatable :: upgrad_trajc(:,:),vpgrad_trajc(:,:),wpgrad_trajc(:,:)
REAL, allocatable :: vorts_gen_trajc(:,:),vortc_gen_trajc(:,:)
REAL, allocatable :: vortz_tilt_trajc(:,:),vortz_stch_trajc(:,:)
INTEGER :: istatus, cur_level
REAL :: xlow, xhigh, ylow, yhigh, zlow, zhigh, pttem
INTEGER :: nstart, nend, npnt, ntrj, ninc
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.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, do_diagnostics
REAL :: ttrajc0_in
INTEGER :: ntrajcs_in, ntrajcs_each_j
CHARACTER (LEN=80) :: timsnd
INTEGER :: tmstrln, nunit(nmax_times),ltrajc_fn,istat,tmstrln3
INTEGER :: kl, counter
REAL missing_value, pi
INTEGER :: tmstrln0,tmstrln1,tmstrln2, idot
CHARACTER (LEN=6) :: timsnd0,timsnd1,timsnd2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
missing_value = -99999.0
hinfmt = 3
reftime = 0.0
READ(5,input,ERR=100)
WRITE(6,'(a)')'Namelist input was successfully read.'
write(6,input)
DO k=1,ntimes
CALL getunit
(nunit(k))
!-----------------------------------------------------------------------
! Read trajectory data
!-----------------------------------------------------------------------
OPEN(UNIT=nunit(k),FILE=trim(trajc_fn_in(k)),STATUS='old', &
FORM='formatted',IOSTAT= istat )
READ(nunit(k),'(a)') runname
READ(nunit(k),'(6e17.6)') xlow, xhigh, ylow, yhigh, zlow, zhigh
write(6,'(6e17.6)') xlow, xhigh, ylow, yhigh, zlow, zhigh
READ(nunit(k),'(3e17.6)') dx, dy, dz
write(6,'(3e17.6)') dx, dy, dz
READ(nunit(k),'(3e17.6)') tstart, tzero, tend
READ(nunit(k),'(i10)') npoints
READ(nunit(k),'(i10)') ntrajcs
ENDDO
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(ttrajc(npoints),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ttrajc")
allocate(xtrajc1(ntrajcs,ntimes),stat=istatus)
allocate(ytrajc1(ntrajcs,ntimes),stat=istatus)
allocate(ztrajc1(ntrajcs,ntimes),stat=istatus)
DO k=1,ntimes
DO j=1,npoints
READ(nunit(k),'(4e17.6)',err=115,end=115) ttrajc(j)
READ(nunit(k),'(i10)',err=115,end=115) ntrajcs_in
IF( ntrajcs_in /= ntrajcs ) then
print*,'ntrajcs read in .ne. ntrajcs in program.'
print*,'Job stopped'
STOP
ENDIF
READ(nunit(k),'(6e17.6)',err=115,end=115) ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
ENDDO
CLOSE(UNIT=nunit(k))
CALL retunit(nunit(k))
ENDDO
npoints_in = npoints
GOTO 125
115 continue
npoints_in = max(1,j-1)
125 continue
print*,'done reading trajectory data'
! IF( tstart_calc /= tend_calc ) then
nstart = npoints_in
DO j=1,npoints_in
IF( ttrajc(j) >= tstart_calc ) then
nstart = j
exit
ENDIF
ENDDO
nend = 1
DO j=npoints_in,1,-1
IF( ttrajc(j) <= tend_calc ) then
nend = j
exit
ENDIF
ENDDO
! ELSE
! nstart = 1
! nend = npoints
! ENDIF
print*,'ttrajc(nstart),ttrajc(nend), tinc_calc=', ttrajc(nstart),ttrajc(nend), tinc_calc
ninc = max(1, nint( tinc_calc/(ttrajc(nstart+1)-ttrajc(nstart)) ))
print*,'tstart_calc, tend_calc, nstart, nend =', tstart_calc, tend_calc, nstart, nend
grdbasfn = trim(hdmpfheader)//'.hdfgrdbas'//trim(hdmpftrailer)
CALL get_dims_from_data(hinfmt,trim(grdbasfn), &
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(zs1d(nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:zs")
zs1d = 0.0
allocate(zp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:zp")
zp=0.0
allocate(ptbar(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ptbar")
ptbar=0.0
allocate(pbar(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:pbar")
pbar=0.0
allocate(qvbar(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:qvbar")
qvbar=0.0
allocate(rhobar(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:rhobar")
rhobar=0.0
allocate(tem1 (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:tem1 ")
tem1 =0.0
allocate(tem2 (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:tem2 ")
tem2 =0.0
allocate(tem3 (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:tem3 ")
tem3 =0.0
allocate(tem4 (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:tem4 ")
tem4 =0.0
allocate(itmp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:itmp")
itmp=0
allocate(vortx_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortx_trajc")
vortx_trajc = missing_value
allocate(vorty_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vorty_trajc")
vorty_trajc = missing_value
allocate(vortz_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortz_trajc")
vortz_trajc = missing_value
allocate(vorts_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vorts_trajc")
vorts_trajc= missing_value
allocate(vortc_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortc_trajc")
vortc_trajc= missing_value
allocate(buoy_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:buoy_trajc")
buoy_trajc = missing_value
allocate(buoyq_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:buoya_trajc")
buoyq_trajc= missing_value
allocate(vorts_gen_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vorts_gen_trajc")
vorts_gen_trajc= missing_value
allocate(vortc_gen_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortc_gen_trajc")
vortc_gen_trajc= missing_value
allocate(vortz_tilt_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortz_tilt_trajc")
vortz_tilt_trajc= missing_value
allocate(vortz_stch_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vortz_stch_trajc")
vortz_stch_trajc= missing_value
allocate(frcs_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:frcs_trajc")
frcs_trajc = missing_value
allocate(pprt_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:pprt_trajc")
pprt_trajc = missing_value
allocate(ptprt_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:ptprt_trajc")
ptprt_trajc = missing_value
allocate(pt_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:pt_trajc")
pt_trajc = missing_value
allocate(qv_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:qv_trajc")
qv_trajc = missing_value
allocate(p_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:p_trajc")
p_trajc = missing_value
allocate(q_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:q_trajc")
q_trajc = missing_value
allocate(upgrad_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:upgrad_trajc")
upgrad_trajc = missing_value
allocate(vpgrad_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vpgrad_trajc")
vpgrad_trajc = missing_value
allocate(wpgrad_trajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:wpgrad_trajc")
wpgrad_trajc = missing_value
allocate(utrajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:utrajc")
allocate(vtrajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:vtrajc")
allocate(wtrajc(ntrajcs,ntimes),stat=istatus)
CALL check_alloc_status(istatus, "arpstrajc:wtrajc")
! allocate(xweight(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:xweight")
! allocate(yweight(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:yweight")
! allocate(zweight(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:zweight")
! allocate(itrajc(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:itrajc")
! allocate(jtrajc(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:jtrajc")
! allocate(ktrajc(ntrajcs,ntimes),stat=istatus)
! CALL check_alloc_status(istatus, "arpstrajc:ktrajc")
!-----------------------------------------------------------------------
! Read x,y,z
!-----------------------------------------------------------------------
CALL hdfreadxyz(nx,ny,nz,x,y,z, trim(grdbasfn), time)
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
!-----------------------------------------------------------------------
! First re-read data at starting (reference) time if necessary.
!-----------------------------------------------------------------------
CALL hdfreadvar(nx,ny,nz,trim(grdbasfn),time,'zp', tem1, itmp )
! print*,'mark 1'
zp(:,:,:) = tem1(:,:,:)
CALL hdfreadvar(nx,ny,nz,trim(grdbasfn),time,'ptbar', tem1, itmp )
ptbar(:,:,:) = tem1(:,:,:)
! print*,'mark 2'
CALL hdfreadvar(nx,ny,nz,trim(grdbasfn),time,'pbar', tem1, itmp )
pbar(:,:,:) = tem1(:,:,:)
CALL hdfreadvar(nx,ny,nz,trim(grdbasfn),time,'qvbar', tem1, itmp )
qvbar(:,:,:) = tem1(:,:,:)
! print*,'mark 3'
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
! DO k=1,nz-1
! zs1d(k) = 0.5*(zp1d(k)+zp1d(k+1))
! ENDDO
! print*,'mark 4'
print*,'nstart, nend, ninc =', nstart, nend, ninc
!-----------------------------------------------------------------------
! Write out data along trajectories
!-----------------------------------------------------------------------
DO k=1,ntimes
idot = index( trim(trajc_fn_in(k)) , '.')
CALL cvttsnd( ttrajc(nstart) ,timsnd0, tmstrln0 )
CALL cvttsnd( ttrajc(nend) ,timsnd1, tmstrln1 )
CALL cvttsnd( reftime(k), timsnd2, tmstrln2 )
trajc_fn_out = trajc_fn_in(k)(1:idot-1)//'.trajc_'//trim(timsnd0)//'-'//trim(timsnd1)//'_'//trim(timsnd2)//'.data'
CALL getunit(nunit(k))
ltrajc_fn = len_trim(trajc_fn_out)
CALL fnversn( trajc_fn_out, ltrajc_fn )
OPEN(UNIT=nunit(k),FILE=trim(trajc_fn_out),STATUS='unknown', &
FORM='formatted',IOSTAT= istat )
counter = 0
DO j = nstart, nend, ninc
counter = counter + 1
ENDDO
WRITE(nunit(k),'(a)') trim(runname)
WRITE(nunit(k),'(6e17.6)') xs(1),xs(nx-1),ys(1),ys(ny-1),zp(1,1,1),zp(1,1,nz-1)
WRITE(nunit(k),'(3e17.6)') dx, dy, z(2)-z(1)
WRITE(nunit(k),'(3e17.6)') ttrajc(nstart),ttrajc(nstart),ttrajc(nend)
WRITE(nunit(k),'(i10)') counter
WRITE(nunit(k),'(i10)') ntrajcs
write(nunit(k),'(a,a,a,a,a,a)') &
' t, x, y, z, ', &
'u, v, w, pt, ', &
'p, qv, ptprt, pprt, ', &
'vortx, vorty, vortz, vorts, vortc, ', &
'buoy, buoyq, frcs, upgrad, vpgrad, ', &
' wpgrad, vorts_gen, vortc_gen, vortz_tilt, vortz_stch '
ENDDO
!-----------------------------------------------------------------------
! Start calculations along trajectories
!-----------------------------------------------------------------------
DO npnt = nstart, nend, ninc
DO k=1,ntimes
DO j=1,ntrajcs
xtrajc1(j,k)=xtrajc(npnt,j,k)
ytrajc1(j,k)=ytrajc(npnt,j,k)
ztrajc1(j,k)=ztrajc(npnt,j,k)
ENDDO
ENDDO
time=ttrajc(npnt)
! print*,'Inside npnt loop, npnt = ', npnt
! May have ninc and hdmpfheader, hdmpftrailer depends on time
CALL cvttsnd( ttrajc(npnt), timsnd, tmstrln )
hisfile = trim(hdmpfheader)//'.hdf'//timsnd(1:tmstrln)//trim(hdmpftrailer)
WRITE(6,'(/a,a,a)') ' Data set ', trim(hisfile),' to be read.'
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'u', tem1, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'v', tem2, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'w', tem3, itmp )
CALL a3dmax0(tem1,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(tem2,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(tem3,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
!-----------------------------------------------------------------------
! Do diagnostics
!-----------------------------------------------------------------------
do_diagnostics = 10
IF( do_diagnostics == 1 .or. do_diagnostics == 10 ) then
!CALL cal_xvort(tem2,tem3,y,zp,nx,ny,nz, tem4)
CALL cal_xvort_flat(tem2,tem3,y,zp,nx,ny,nz, tem4)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortx_trajc )
! Tilting of vortx into vertical
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-2
tem4(i,j,k) = tem4(i,j,k)*(tem3(i+1,j,k)+tem3(i+1,j,k+1) &
-tem3(i-1,j,k)-tem3(i-1,j,k+1))/(4*dx)
ENDDO
ENDDO
ENDDO
CALL edgfill(tem4,1,nx,2,nx-2, 1,ny,1,ny-1, 1,nz,1,nz-1)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortz_trajc ) ! vortz_trajc is used as work array
!CALL cal_yvort(tem1,tem3,x,zp,nx,ny,nz, tem4)
CALL cal_yvort_flat(tem1,tem3,x,zp,nx,ny,nz, tem4)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vorty_trajc )
! Tilting of vorty into vertical
DO k=1,nz-1
DO j=2,ny-2
DO i=1,nx-1
tem4(i,j,k) = tem4(i,j,k)*(tem3(i,j+1,k)+tem3(i,j+1,k+1) &
-tem3(i,j-1,k)-tem3(i,j-1,k+1))/(4*dy)
ENDDO
ENDDO
ENDDO
CALL edgfill(tem4,1,nx,1,nx-1, 1,ny,2,ny-2, 1,nz,1,nz-1)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortz_tilt_trajc )
vortz_tilt_trajc(:,:) = vortz_tilt_trajc(:,:) + vortz_trajc(:,:)
! vertical vorticity
!CALL cal_zvort(tem1,tem2,x,y ,nx,ny,nz, tem4,zp)
CALL cal_zvort_flat(tem1,tem2,x,y ,nx,ny,nz, tem4)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortz_trajc )
! Vertical stretching for vertical vorticity
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k) = tem4(i,j,k)*(tem3(i,j,k+1)-tem3(i,j,k))/(zp(i,j,k+1)-zp(i,j,k))
ENDDO
ENDDO
ENDDO
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortz_stch_trajc )
CALL avgx(tem1,0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem4 )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, utrajc )
tem1 = tem4
CALL avgy(tem2,0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem4 )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vtrajc )
tem2 = tem4
CALL avgz(tem3,0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem4 )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, wtrajc )
tem3 = tem4
IF( .false. ) then
! angle/direction of horizontal wind
pi = 4.0 * atan(1.0)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
IF (tem1(i,j,k).eq.0.0 .and. tem2(i,j,k).gt.0.0) THEN
tem4 (i,j,k) = (pi/2.0)
ELSE IF (tem1(i,j,k).eq.0.0 .and. tem2(i,j,k).lt.0.0) THEN
tem4 (i,j,k) = 3.0*pi/2.0
ELSE IF (tem1(i,j,k).gt.0.0 .and. tem2(i,j,k).eq.0.0) THEN
tem4 (i,j,k) = 0.0
ELSE IF (tem1(i,j,k).lt.0.0 .and. tem2(i,j,k).eq.0.0) THEN
tem4 (i,j,k) = pi
ELSE IF (tem1(i,j,k).gt.0.0 .and. tem2(i,j,k).gt.0.0) THEN
tem4 (i,j,k) = atan(tem2(i,j,k)/tem1(i,j,k))
ELSE IF (tem1(i,j,k).gt.0.0 .and. tem2(i,j,k).lt.0.0) THEN
tem4 (i,j,k) = atan(tem2(i,j,k)/tem1(i,j,k)) + (2.0*pi)
ELSE IF (tem1(i,j,k).lt.0.0 .and. tem2(i,j,k).lt.0.0) THEN
tem4 (i,j,k) = atan(tem2(i,j,k)/tem1(i,j,k)) + (pi)
ELSE IF (tem1(i,j,k).lt.0.0 .and. tem2(i,j,k).gt.0.0) THEN
tem4 (i,j,k) = atan(tem2(i,j,k)/tem1(i,j,k)) + (pi)
ENDIF
ENDDO
ENDDO
ENDDO
!
! horizontal wind speed
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k) = sqrt( (tem1(i,j,k))**2 + (tem2(i,j,k))**2 )
ENDDO
ENDDO
ENDDO
ENDIF
!-----------------------------------------------------------------------
! Horizontal streamwise and cross-stream vorticity
!-----------------------------------------------------------------------
DO k=1,ntimes
DO ntrj=1,ntrajcs
utem = utrajc(ntrj,k)
vtem = vtrajc(ntrj,k)
wtem = wtrajc(ntrj,k)
vorts_trajc(ntrj,k)=( vortx_trajc(ntrj,k)*utem &
+vorty_trajc(ntrj,k)*vtem) &
/(1.0e-10+sqrt(utem**2+vtem**2))
vortc_trajc(ntrj,k)=(-vortx_trajc(ntrj,k)*vtem &
+vorty_trajc(ntrj,k)*utem) &
/(1.0e-10+sqrt(utem**2+vtem**2))
ENDDO
ENDDO
ENDIF
IF( do_diagnostics == 2 .or. do_diagnostics == 10 ) then
!-----------------------------------------------------------------------
! Buoyancy term
!-----------------------------------------------------------------------
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'qc', tem1, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'qr', tem2, itmp )
tem4 = tem1 + tem2
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k) = -g*tem4(i,j,k)/(1+qvbar(i,j,k))
END DO
END DO
END DO
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, buoyq_trajc )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'pt', tem1, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'p', tem2, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'qv', tem3, itmp )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem1, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, pt_trajc )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem2, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, p_trajc )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem3, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, qv_trajc )
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k) = tem1(i,j,k)-ptbar(i,j,k)
tem2(i,j,k) = tem2(i,j,k)- pbar(i,j,k)
tem3(i,j,k) = tem3(i,j,k)-qvbar(i,j,k)
ENDDO
ENDDO
ENDDO
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem1, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, ptprt_trajc )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem2, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, pprt_trajc )
!
!-----------------------------------------------------------------------
! wbuoy = g ( ptprt/ptbar-pprt/(rddrv*pbar)+
! qvprt/(rddrv+qvbar)-(qvprt+qc+qr+qs+qi+qh)/(1+qvbar)
! -(ptprt*ptprt)/(ptbar*ptbar) !2nd-order
! +0.5*(ptprt*pprt)/(cpdcv*ptbar*pbar)) !2nd-order
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pttem = tem1(i,j,k)/ptbar(i,j,k)
tem1(i,j,k) = g*( &
pttem - tem2(i,j,k)/(rddrv*pbar(i,j,k)) &
-pttem*pttem+0.5*pttem*tem2(i,j,k)/(cpdcv*pbar(i,j,k)) &
+ tem3(i,j,k)/(rddrv+qvbar(i,j,k))) &
+ tem4(i,j,k)
END DO
END DO
END DO
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem1, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, buoy_trajc )
!-----------------------------------------------------------------------
! Baroclinic vorticity generation in streamwise direction
!-----------------------------------------------------------------------
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'u', tem2, itmp )
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'v', tem3, itmp )
do k=1,nz-1
do j=2,ny-2
do i=2,nx-2
utem = 0.5*(tem2(i+1,j,k)+tem2(i,j,k))
vtem = 0.5*(tem3(i,j+1,k)+tem3(i,j,k))
tem4(i,j,k)=( -vtem*(tem1(i+1,j,k)-tem1(i-1,j,k))/(x(i+1)-x(i)) &
+utem*(tem1(i,j+1,k)-tem1(i,j-1,k))/(y(j+1)-y(j)) )/ &
( 2.0 * (1.0e-10+sqrt( utem*utem + vtem*vtem)) )
enddo
enddo
enddo
CALL edgfill(tem4,1,nx,2,nx-2, 1,ny,2,ny-2, 1,nz,1,nz-1)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vorts_gen_trajc )
!-----------------------------------------------------------------------
! Baroclinic vorticity generation in cross-wise direction
!-----------------------------------------------------------------------
do k=1,nz-1
do j=2,ny-2
do i=2,nx-2
utem = 0.5*(tem2(i+1,j,k)+tem2(i,j,k))
vtem = 0.5*(tem3(i,j+1,k)+tem3(i,j,k))
tem4(i,j,k)=( utem*(tem1(i+1,j,k)-tem1(i-1,j,k))/(x(i+1)-x(i)) &
+vtem*(tem1(i,j+1,k)-tem1(i,j-1,k))/(y(j+1)-y(j)) )/ &
( 2.0 * (1.0e-10+sqrt( utem*utem + vtem*vtem)) )
enddo
enddo
enddo
CALL edgfill(tem4,1,nx,2,nx-2, 1,ny,2,ny-2, 1,nz,1,nz-1)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vortc_gen_trajc )
ENDIF
IF( do_diagnostics == 3 .or. do_diagnostics == 10 ) then
!-----------------------------------------------------------------------
! Vertical pressure gradient force
!-----------------------------------------------------------------------
CALL hdfreadvar(nx,ny,nz,trim(hisfile),time,'p', tem2, itmp )
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem2(i,j,k) = tem2(i,j,k)- pbar(i,j,k)
ENDDO
ENDDO
ENDDO
DO k=2,nz-2
if( k == 2 ) then
kl = 2
factor = 1.0
else
kl = k-1
factor = 2.0
endif
DO j=1,ny-1
DO i=1,nx-1
rhobar(i,j,k) = pbar(i,j,k)/(rd*ptbar(i,j,k)*(pbar(i,j,k)/p0)**rddcp)
tem1(i,j,k)=-(tem2(i,j,k+1)-tem2(i,j,kl))/ &
(factor*rhobar(i,j,k)*(zp(i,j,k+1)-zp(i,j,kl)))
END DO
END DO
DO j=1,ny-1
DO i=2,nx-2
tem3(i,j,k)=-(tem2(i+1,j,k)-tem2(i-1,j,k))/(2*dx*rhobar(i,j,k))
END DO
END DO
DO j=2,ny-2
DO i=1,nx-1
tem4(i,j,k)=-(tem2(i,j,k+1)-tem2(i,j,kl))/(2*dy*rhobar(i,j,k))
END DO
END DO
END DO
CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,2,nz-2)
CALL edgfill(tem3,1,nx,2,nx-2, 1,ny,1,ny-1, 1,nz,2,nz-2)
CALL edgfill(tem4,1,nx,1,nx-1, 1,ny,2,ny-2, 1,nz,2,nz-2)
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem1, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, wpgrad_trajc )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem3, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, upgrad_trajc )
CALL intrp_trajc(nx,ny,nz,xs,ys,zp, tem4, &
xtrajc1,ytrajc1,ztrajc1, ntrajcs, ntimes, vpgrad_trajc )
DO k=1,ntimes
DO ntrj=1,ntrajcs
utem = utrajc(ntrj,k)
vtem = vtrajc(ntrj,k)
wtem = wtrajc(ntrj,k)
frcs_trajc(ntrj,k)=(upgrad_trajc(ntrj,k)*utem &
+vpgrad_trajc(ntrj,k)*vtem+ &
(wpgrad_trajc(ntrj,k)+buoy_trajc(ntrj,k))*wtem) &
/(1.0e-10+sqrt(utem**2+vtem**2+wtem**2))
ENDDO
ENDDO
ENDIF
!-----------------------------------------------------------------------
! Write out data long trajectories
!-----------------------------------------------------------------------
j = npnt
DO k=1,ntimes
write(nunit(k),'(f8.1)') ttrajc(j)
DO i = 1,ntrajcs
write(nunit(k),'(f8.1,6f10.2,32e14.6)') ttrajc(j), &
xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k), &
utrajc(i,k),vtrajc(i,k),wtrajc(i,k), &
pt_trajc(i,k),p_trajc(i,k),qv_trajc(i,k), &
ptprt_trajc(i,k),pprt_trajc(i,k), &
vortx_trajc(i,k),vorty_trajc(i,k),vortz_trajc(i,k), &
vorts_trajc(i,k), vortc_trajc(i,k), &
buoy_trajc(i,k),buoyq_trajc(i,k),frcs_trajc(i,k), &
upgrad_trajc(i,k),vpgrad_trajc(i,k),wpgrad_trajc(i,k), &
vorts_gen_trajc(i,k),vortc_gen_trajc(i,k) , &
vortz_tilt_trajc(i,k),vortz_stch_trajc(i,k)
ENDDO
ENDDO
ENDDO
!-----------------------------------------------------------------------
! Close the trajectory file
!-----------------------------------------------------------------------
DO k=1,ntimes
CALL flush(nunit(k))
CLOSE(UNIT=nunit(k))
CALL retunit( nunit(k))
ENDDO
print*,'done writing trajectory data'
STOP
100 WRITE(6,'(a)') 'Error reading NAMELIST file. Program ARPSTRAJC stopped.'
STOP
END PROGRAM ARPSCALCTRAJC
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DTAREADWIND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dtareadwind(nx,ny,nz,x,y,z,zp1d, grdbasfn,datafn, time,u,v,w, itmp, tem1) 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 :: zp1d(nz)
REAL :: tem1(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.
!
!-----------------------------------------------------------------------
!
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 DTAREAD.'
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,zp1d,grdbas,trim(grdbasfn), &
time,u,v,w, itmp,tem1)
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.'
CALL arpsstop('arpsstop called from dtareadwind during base read-2',1)
100 CONTINUE
CALL hdfreadwind(nx,ny,nz,x,y,z,zp1d,grdbas,trim(datafn), &
time,u,v,w, itmp,tem1)
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,zp1d,grdbas,filename, & 4,122
time,u,v,w, itmp, tem1)
!-----------------------------------------------------------------------
! 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 :: zp1d(nz)
REAL :: tem1(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,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
zp1d(:)=tem1(2,2,:)
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'
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)
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.'
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
INTEGER FUNCTION get_ktrajc(z, zs, nz)
IMPLICIT NONE
INTEGER :: nz, k
REAL :: z, zs(nz)
IF( z < zs(1) ) then
get_ktrajc = 1
RETURN
ENDIF
IF( z >= zs(nz) ) then
get_ktrajc = nz-1
RETURN
ENDIF
DO k=1,nz-1
IF( z >= zs(k) .and. z < zs(k+1) ) then
get_ktrajc = k
EXIT
endif
ENDDO
RETURN
END FUNCTION get_ktrajc
SUBROUTINE intrp_trajc(nx,ny,nz,xs,ys,zp, u, & 21
xtrajc,ytrajc,ztrajc, ntrajcs,ntimes, utrajc )
!
!-----------------------------------------------------------------------
! 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) ! input - field to be interpolated to trajectory points
INTEGER :: ntrajcs, ntimes
REAL :: xtrajc(ntrajcs,ntimes),ytrajc(ntrajcs,ntimes),ztrajc(ntrajcs,ntimes)
REAL :: utrajc(ntrajcs,ntimes) ! output - u interpolated to trajectory points
INTEGER :: get_ktrajc
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 :: ntrajc,itrajc1,jtrajc1,ktrajc1
REAL :: xweight1,yweight1,zweight1
!
REAL :: zs1d(nz),zy1(nz),zy2(nz) ! automatic work array
INTEGER :: i,k,l
dx = xs(2)-xs(1)
dy = ys(2)-ys(1)
missing_value = -99999.0
DO k=1,ntimes
DO i= 1,ntrajcs
IF( xtrajc(i,k) == missing_value .or. &
ytrajc(i,k) == missing_value .or. &
ztrajc(i,k) == missing_value ) then
utrajc(i,k) = missing_value
ELSE
itrajc1 = MAX(1, MIN(nx-2, INT((xtrajc(i,k)-xs(1))/dx)+1 ))
jtrajc1 = MAX(1, MIN(ny-2, INT((ytrajc(i,k)-ys(1))/dy)+1 ))
! ktrajc1 = get_ktrajc(ztrajc(i,k),zs,nz-1)
xweight1 = (xtrajc(i,k)-xs(itrajc1))/dx
yweight1 = (ytrajc(i,k)-ys(jtrajc1))/dy
DO l=1,nz-1
zy1(l) = (1.0-xweight1)*(zp(itrajc1 ,jtrajc1 ,l)+zp(itrajc1 ,jtrajc1 ,l+1))*0.5 &
+xweight1 *(zp(itrajc1+1,jtrajc1 ,l)+zp(itrajc1+1,jtrajc1 ,l+1))*0.5
zy2(l) = (1.0-xweight1)*(zp(itrajc1 ,jtrajc1+1,l)+zp(itrajc1 ,jtrajc1+1,l+1))*0.5 &
+xweight1 *(zp(itrajc1+1,jtrajc1+1,l)+zp(itrajc1+1,jtrajc1+1,l+1))*0.5
zs1d(l) = ( 1.0-yweight1 )*zy1(l) + yweight1*zy2(l)
END DO
ktrajc1 = get_ktrajc(ztrajc(i,k),zs1d,nz-1)
uy1z1 = (1.0-xweight1)*u(itrajc1,jtrajc1 ,ktrajc1 )+xweight1*u(itrajc1+1,jtrajc1 ,ktrajc1 )
uy2z1 = (1.0-xweight1)*u(itrajc1,jtrajc1+1,ktrajc1 )+xweight1*u(itrajc1+1,jtrajc1+1,ktrajc1 )
uy1z2 = (1.0-xweight1)*u(itrajc1,jtrajc1 ,ktrajc1+1)+xweight1*u(itrajc1+1,jtrajc1 ,ktrajc1+1)
uy2z2 = (1.0-xweight1)*u(itrajc1,jtrajc1+1,ktrajc1+1)+xweight1*u(itrajc1+1,jtrajc1+1,ktrajc1+1)
uz1 = ( 1.0-yweight1 )*uy1z1 + yweight1*uy2z1
uz2 = ( 1.0-yweight1 )*uy1z2 + yweight1*uy2z2
zweight1 = (ztrajc(i,k)-zs1d(ktrajc1))/(zs1d(ktrajc1+1)-zs1d(ktrajc1))
utrajc(i,k) = (1.0-zweight1)*uz1+zweight1*uz2
ENDIF
ENDDO
ENDDO
RETURN
END SUBROUTINE intrp_trajc
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreadvar(nx,ny,nz,filename, time, varname, var, itmp ) 15,8
!-----------------------------------------------------------------------
! 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
CHARACTER (LEN=*) :: varname
CHARACTER (LEN=*) :: filename
REAL :: time
REAL :: var(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmax(nz), hmin(nz) ! Temporary array
INTEGER :: ireturn
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin
INTEGER :: istat, sd_id
INTEGER :: varflg, istatus
REAL :: timein
!-----------------------------------------------------------------------
!
! 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
CALL hdfrdr(sd_id,"time",timein,istat)
IF( timein /= time ) then
print*,'Warning: time in data does not match time passed into READHDFVAR.'
Print*,'time in program =', time, ', time in data =', timein
print*,'time in program reset to ', timein
time = timein
ENDIF
!-----------------------------------------------------------------------
! 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 total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,varname,nx,ny,nz,var,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
!-----------------------------------------------------------------------
! Friendly exit message
!-----------------------------------------------------------------------
930 CONTINUE
WRITE(6,'(/a,a,a,F8.1,a/)') &
' Variable ', varname,' 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 HDFREADVAR'
STOP
ireturn=1
130 CONTINUE
CALL hdfclose(sd_id,istat)
IF (ireturn == 0) THEN
IF (istat == 0) THEN
WRITE(6,'(/a/a)') &
"HDFREADVAR: Successfully read ", trim(filename)
ELSE
WRITE(6,'(/a,i3,a/,a)') &
"HDFREADVAR: ERROR (status=", istat, ") closing ", trim(filename)
END IF
END IF
RETURN
END SUBROUTINE hdfreadvar
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREADXYZ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreadxyz(nx,ny,nz,x,y,z, filename, time) 1,10
!-----------------------------------------------------------------------
! 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.
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
CHARACTER (LEN=*) :: filename
REAL :: x(nx) ! x coord.
REAL :: y(ny) ! y coord.
REAL :: z(nz) ! z coord.
REAL :: time
INTEGER :: ireturn
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 :: 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
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).
!-----------------------------------------------------------------------
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
!-----------------------------------------------------------------------
! 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'
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)
END IF
END IF
RETURN
END SUBROUTINE hdfreadxyz
SUBROUTINE cal_xvort(v,w,y,zp,nx,ny,nz, xvort)
!
! Rewritten (DTD 10-10-05) to include terrain-following case
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: xvort(nx,ny,nz)
REAL :: v(nx,ny,nz), w(nx,ny,nz)
REAL :: y(ny),zp(nx,ny,nz)
REAL :: zweightl,zweightr,zs(nx,ny,nz)
REAL :: ws(nx,ny,nz),wr,wl
REAL :: dwdy, dvdz
INTEGER :: i,j,k,kl,k1,k2
REAL :: factor
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j=2,ny-2
DO i=1,nx-1
DO k=2,nz-2
! Determine value of z and w at center scalar point
zs(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
ws(i,j,k)=0.5*(w(i,j,k)+w(i,j,k+1))
END DO
END DO
END DO
DO j=2,ny-2
DO i=1,nx-1
DO k=2,nz-2
k1=k
k2=k
! Determine k levels for vertical interpolation on either side of center point
IF(zs(i,j,k) < zp(i,j-1,k)) THEN
DO WHILE (zs(i,j,k) < zp(i,j-1,k1))
k1=k1-1
END DO
ELSE IF(zs(i,j,k) > zp(i,j-1,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i,j-1,k1+1))
k1=k1+1
END DO
END IF
IF(zs(i,j,k) < zp(i,j+1,k)) THEN
DO WHILE (zs(i,j,k) < zp(i,j+1,k2))
k2=k2-1
END DO
ELSE IF(zs(i,j,k) > zp(i,j+1,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i,j+1,k2+1))
k2=k2+1
END DO
END IF
! Precalculate left-side and right-side interpolated w values
IF(k1 >= 2 .and. k1 <= nz-2) THEN
zweightl=(zs(i,j,k)-zp(i,j-1,k1))/(zp(i,j-1,k1+1)-zp(i,j-1,k1))
wl=(1.0-zweightl)*zp(i,j-1,k1)+zweightl*zp(i,j-1,k1+1)
END IF
IF(k2 >= 2 .and. k2 <= nz-2) THEN
zweightr=(zs(i,j,k)-zp(i,j+1,k2))/(zp(i,j+1,k2+1)-zp(i,j+1,k2))
wr=(1.0-zweightr)*zp(i,j+1,k2)+zweightr*zp(i,j+1,k2+1)
END IF
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
! Make sure the k levels are not below ground or above the model top
IF((k1 < 2 .and. k2 < 2) .or. (k1 > nz-2 .and. k2 > nz-2)) THEN
! No calculation can be performed! Set vorticity to zero at this point.
xvort(i,j,k) = 0.0
ELSE IF(k1 < 2) THEN ! Do a right-sided difference
dwdy=(wr-ws(i,j,k))/(y(j+1)-y(j))
ELSE IF(k2 < 2) THEN ! Do a left-sided difference
dwdy=(ws(i,j,k)-wl)/(y(j)-y(j-1))
ELSE ! Do a normal centered difference
dwdy=(wr-wl)/(y(j+1)-y(j-1))
END IF
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
! Calculate dvdz (there will be some error here due to assumption of uniform grid in the vertical)
dvdz=0.25*factor*(v(i,j,k+1)+v(i,j+1,k+1)-v(i,j,kl )-v(i,j+1,kl ))/(zs(i,j,k+1)-zs(i,j,k))
! Calculate xvort
xvort(i,j,k)=dwdy-dvdz
END DO
END DO
END DO
DO k=2,nz-2
DO i=1,nx-1
xvort(i, 1,k)=xvort(i, 2,k)
xvort(i,ny-1,k)=xvort(i,ny-2,k)
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
xvort(i,j, 1)=xvort(i,j, 2)
xvort(i,j,nz-1)=xvort(i,j,nz-2)
END DO
END DO
RETURN
END SUBROUTINE cal_xvort
SUBROUTINE cal_yvort(u,w,x,zp,nx,ny,nz, yvort)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: yvort(nx,ny,nz)
REAL :: zp(nx,ny,nz)
REAL :: u(nx,ny,nz), w(nx,ny,nz)
REAL :: x(ny),zs(nx,ny,nz)
REAL :: wl,wr,ws(nx,ny,nz)
REAL :: dwdx,dudz,zweightl,zweightr
INTEGER :: k1,k2
INTEGER :: i,j,k, kl
REAL :: factor
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j=1,ny-1
DO i=2,nx-2
DO k=2,nz-2
! Determine value of z and w at center scalar point
zs(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
ws(i,j,k)=0.5*(w(i,j,k)+w(i,j,k+1))
END DO
END DO
END DO
DO j=1,ny-1
DO i=2,nx-2
DO k=2,nz-2
k1=k
k2=k
! Determine k levels for vertical interpolation on either side of center point
IF(zs(i,j,k) < zp(i-1,j,k)) THEN
DO WHILE (zs(i,j,k) < zp(i-1,j,k1))
k1=k1-1
END DO
ELSE IF(zs(i,j,k) > zp(i-1,j,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i-1,j,k1+1))
k1=k1+1
END DO
END IF
IF(zs(i,j,k) < zp(i+1,j,k)) THEN
DO WHILE (zs(i,j,k) < zp(i+1,j,k2))
k2=k2-1
END DO
ELSE IF(zs(i,j,k) > zp(i+1,j,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i+1,j,k2+1))
k2=k2+1
END DO
END IF
! Precalculate left-side and right-side interpolated w values
IF(k1 >= 2 .and. k1 <= nz-2) THEN
zweightl=(zs(i,j,k)-zp(i-1,j,k1))/(zp(i-1,j,k1+1)-zp(i-1,j,k1))
wl=(1.0-zweightl)*zp(i-1,j,k1)+zweightl*zp(i-1,j,k1+1)
END IF
IF(k2 >= 2 .and. k2 <= nz-2) THEN
zweightr=(zs(i,j,k)-zp(i+1,j,k2))/(zp(i+1,j,k2+1)-zp(i+1,j,k2))
wr=(1.0-zweightr)*zp(i+1,j,k2)+zweightr*zp(i+1,j,k2+1)
END IF
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
! Make sure the k levels are not below ground or above the model top
IF((k1 < 2 .and. k2 < 2) .or. (k1 > nz-2 .and. k2 > nz-2)) THEN
! No calculation can be performed! Set vorticity to zero at this point.
yvort(i,j,k) = 0.0
ELSE IF(k1 < 2) THEN ! Do a right-sided difference
dwdx=(wr-ws(i,j,k))/(x(i+1)-x(i))
ELSE IF(k2 < 2) THEN ! Do a left-sided difference
dwdx=(ws(i,j,k)-wl)/(x(i)-x(i-1))
ELSE ! Do a normal centered difference
dwdx=(wr-wl)/(x(i+1)-x(i-1))
END IF
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
! Calculate dudz
dudz=0.25*factor*(u(i,j,k+1)+u(i+1,j,k+1)-u(i,j,kl)-u(i+1,j,kl))/(zs(i,j,k+1)-zs(i,j,k))
! Calculate xvort
yvort(i,j,k)=dudz-dwdx
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
yvort( 1,j,k)=yvort( 2,j,k)
yvort(nx-1,j,k)=yvort(nx-2,j,k)
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
yvort(i,j, 1)=yvort(i,j, 2)
yvort(i,j,nz-1)=yvort(i,j,nz-2)
END DO
END DO
RETURN
END SUBROUTINE cal_yvort
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINES cal_zvort ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Calculate Vort*10^5 (1/s) value.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Min Zou
! 3/2/98
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
SUBROUTINE cal_zvort(u,v,x,y,nx,ny,nz,zvort,zp)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: zvort(nx,ny,nz),zp(nx,ny,nz)
REAL :: u(nx,ny,nz), v(nx,ny,nz)
REAL :: x(nx), y(ny)
INTEGER :: kw, ke, ks, kn
REAL :: zs(nx,ny,nz)
REAL :: us(nx,ny,nz),vs(nx,ny,nz)
REAL :: viw,vie,uin,uis,zweightw,zweighte,zweightn,zweights
REAL :: dudy,dvdx
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Precalculate z, u, and v, at the scalar points
DO j=2,ny-2
DO i=2,nx-2
DO k=2,nz-2
zs(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
us(i,j,k)=0.5*(u(i+1,j,k)+u(i,j,k))
vs(i,j,k)=0.5*(v(i,j+1,k)+v(i,j,k))
END DO
END DO
END DO
DO j=2,ny-2
DO i=2,nx-2
DO k=2,nz-2
kw=k
ke=k
ks=k
kn=k
! Determine k levels for vertical interpolation of u and v on either side of center point
IF(zs(i,j,k) < zp(i-1,j,k)) THEN
DO WHILE (zs(i,j,k) < zp(i-1,j,kw))
kw=kw-1
END DO
ELSE IF(zs(i,j,k) > zp(i-1,j,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i-1,j,kw+1))
kw=kw+1
END DO
END IF
IF(zs(i,j,k) < zp(i+1,j,k)) THEN
DO WHILE (zs(i,j,k) < zp(i+1,j,ke))
ke=ke-1
END DO
ELSE IF(zs(i,j,k) > zp(i+1,j,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i+1,j,ke+1))
ke=ke+1
END DO
END IF
IF(zs(i,j,k) < zp(i,j-1,k)) THEN
DO WHILE (zs(i,j,k) < zp(i,j-1,ks))
ks=ks-1
END DO
ELSE IF(zs(i,j,k) > zp(i,j-1,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i,j-1,ks+1))
ks=ks+1
END DO
END IF
IF(zs(i,j,k) < zp(i,j+1,k)) THEN
DO WHILE (zs(i,j,k) < zp(i,j+1,kn))
kn=kn-1
END DO
ELSE IF(zs(i,j,k) > zp(i,j+1,k+1)) THEN
DO WHILE (zs(i,j,k) > zp(i,j+1,kn+1))
kn=kn+1
END DO
END IF
! Precalculate vertically-interpolated u and v values
IF(kw >= 2 .and. kw <= nz-2) THEN
zweightw=(zs(i,j,k)-zp(i-1,j,kw))/(zp(i-1,j,kw+1)-zp(i-1,j,kw))
viw=(1.0-zweightw)*vs(i-1,j,kw)+zweightw*vs(i-1,j,kw+1)
END IF
IF(ke >= 2 .and. ke <= nz-2) THEN
zweighte=(zs(i,j,k)-zp(i+1,j,ke))/(zp(i+1,j,ke+1)-zp(i+1,j,ke))
vie=(1.0-zweighte)*vs(i+1,j,ke)+zweighte*vs(i,j,ke+1)
END IF
IF(ks >= 2 .and. ks <= nz-2) THEN
zweights=(zs(i,j,k)-zp(i,j-1,ks))/(zp(i,j-1,ks+1)-zp(i,j-1,ks))
uis=(1.0-zweights)*us(i,j-1,ks)+zweights*us(i,j-1,ks+1)
END IF
IF(kn >= 2 .and. kn <= nz-2) THEN
zweightn=(zs(i,j,k)-zp(i,j+1,kn))/(zp(i,j+1,kn+1)-zp(i,j+1,kn))
uin=(1.0-zweightn)*us(i,j+1,kn)+zweightn*us(i,j+1,kn+1)
END IF
! Make sure the k levels are not below ground or above the model top and calculate dudy and dvdx
IF((ks < 2 .and. kn < 2) .or. (ks > nz-2 .and. kn > nz-2)) THEN
dudy = 0.0
ELSE IF(ks < 2) THEN ! Do a right-sided difference
dudy=(uin-us(i,j,k))/(y(j+1)-y(j))
ELSE IF(kn < 2) THEN ! Do a left-sided difference
dudy=(us(i,j,k)-uis)/(y(j)-y(j-1))
ELSE ! Do a normal centered difference
dudy=(uin-uis)/(y(j+1)-y(j-1))
END IF
IF((kw < 2 .and. ke < 2) .or. (kw > nz-2 .and. ke > nz-2)) THEN
dvdx = 0.0
ELSE IF(kw < 2) THEN ! Do a right-sided difference
dvdx=(vie-vs(i,j,k))/(x(i+1)-x(i))
ELSE IF(ke < 2) THEN ! Do a left-sided difference
dvdx=(vs(i,j,k)-viw)/(x(i)-x(i-1))
ELSE ! Do a normal centered difference
dvdx=(vie-viw)/(x(i+1)-x(i-1))
END IF
! Finally, calculate the vorticity
zvort(i,j,k)=dvdx-dudy
END DO
END DO
END DO
DO j=2,ny-2
DO i=2,nx-2
zvort(i,j, 1)=zvort(i,j, 2)
zvort(i,j,nz-1)=zvort(i,j,nz-2)
END DO
END DO
DO k=1,nz-1
DO j=2,ny-2
zvort( 1,j,k)=zvort( 2,j,k)
zvort(nx-1,j,k)=zvort(nx-2,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
zvort(i, 1,k)=zvort(i, 2,k)
zvort(i,ny-1,k)=zvort(i,ny-2,k)
END DO
END DO
RETURN
END SUBROUTINE cal_zvort
SUBROUTINE cal_xvort_flat(v,w,y,zp,nx,ny,nz, xvort) 1
!
! Only valid for Cartesian grid, flat terrain case
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: xvort(nx,ny,nz)
REAL :: v(nx,ny,nz), w(nx,ny,nz)
REAL :: y(ny),zp(nx,ny,nz)
INTEGER :: i,j,k,kl
REAL :: factor
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=2,nz-2
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
DO j=2,ny-2
DO i=1,nx-1
xvort(i,j,k)=0.25*( &
(w(i,j+1,k)+w(i,j+1,k+1)-w(i,j-1,k)-w(i,j-1,k+1))/(y(j+1)-y(j)) &
-factor*(v(i,j,k+1)+v(i,j+1,k+1)-v(i,j,kl )-v(i,j+1,kl ))/(zp(i,j,k+1)-zp(i,j,k)) )
END DO
END DO
END DO
DO k=2,nz-2
DO i=1,nx-1
xvort(i, 1,k)=xvort(i, 2,k)
xvort(i,ny-1,k)=xvort(i,ny-2,k)
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
xvort(i,j, 1)=xvort(i,j, 2)
xvort(i,j,nz-1)=xvort(i,j,nz-2)
END DO
END DO
RETURN
END SUBROUTINE cal_xvort_flat
SUBROUTINE cal_yvort_flat(u,w,x,zp,nx,ny,nz, yvort) 1
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: yvort(nx,ny,nz)
REAL :: u(nx,ny,nz), w(nx,ny,nz)
REAL :: x(ny),zp(nx,ny,nz)
INTEGER :: i,j,k, kl
REAL :: factor
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=2,nz-2
if( k == 2 ) then
kl = 2
factor = 2.0
else
kl = k-1
factor = 1.0
endif
DO j=1,ny-1
DO i=2,nx-2
yvort(i,j,k)=0.25*( &
-(w(i+1,j,k)+w(i+1,j,k+1)-w(i-1,j,k)-w(i-1,j,k+1))/(x(i+1)-x(i)) &
+factor*(u(i,j,k+1)+u(i+1,j,k+1)-u(i,j,kl )-u(i+1,j,kl )) &
/(zp(i,j,k+1)-zp(i,j,k)) )
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
yvort( 1,j,k)=yvort( 2,j,k)
yvort(nx-1,j,k)=yvort(nx-2,j,k)
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
yvort(i,j, 1)=yvort(i,j, 2)
yvort(i,j,nz-1)=yvort(i,j,nz-2)
END DO
END DO
RETURN
END SUBROUTINE cal_yvort_flat
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINES cal_zvort_flat ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Calculate Vort*10^5 (1/s) value.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Min Zou
! 3/2/98
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
SUBROUTINE cal_zvort_flat(u,v,x,y,nx,ny,nz,zvort) 1
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: zvort(nx,ny,nz)
REAL :: u(nx,ny,nz), v(nx,ny,nz)
REAL :: x(nx), y(ny)
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
zvort(i,j,k)= ( &
(v(i+1,j,k)-v(i-1,j,k)+v(i+1,j+1,k)-v(i-1,j+1,k))/ &
(4*(x(i+1)-x(i)))- &
(u(i,j+1,k)-u(i,j-1,k)+u(i+1,j+1,k)-u(i+1,j-1,k))/ &
(4*(y(j+1)-y(j))) )
END DO
END DO
END DO
DO j=2,ny-2
DO i=2,nx-2
zvort(i,j, 1)=zvort(i,j, 2)
zvort(i,j,nz-1)=zvort(i,j,nz-2)
END DO
END DO
DO k=1,nz-1
DO j=2,ny-2
zvort( 1,j,k)=zvort( 2,j,k)
zvort(nx-1,j,k)=zvort(nx-2,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
zvort(i, 1,k)=zvort(i, 2,k)
zvort(i,ny-1,k)=zvort(i,ny-2,k)
END DO
END DO
RETURN
END SUBROUTINE cal_zvort_flat
!