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 !