module io_ascii ! ! Reads/writes old ascii or binary restart files for gryffin ! (io_netcdf is preferred, as it is more portable, its as fast as binary, ! and it allows variables to be added without affecting the file format.) ! ! caution: io_binary.f90 is generated from io_ascii.f90 using io_binary.sh, ! plus a little hand editing. ! ! Warning: the ascii *.res files don't print out quite enough digits ! to preserve the full precision of the Cray 64-bit floating point numbers. ! Someday could upgrade the format statements (though then they wouldn't ! be backwards compatible, and we are trying to use the more portable ! io_netcdf routines anyway). ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! implicit none real, private, allocatable, dimension(:,:) :: ri1 real, private, allocatable, dimension(:,:,:) :: ri2 real, private, allocatable, dimension(:,:,:,:) :: ri3 real, private, allocatable, dimension(:,:,:,:,:) :: ri4 contains subroutine wread(filename, ntim, ntimr, tim, dt1, nsteps, nfile) ! ! Write/read results to/from the *.res file. ! The format of the *.res file is essentially given by the ! read/writes in this subroutine, and must be compatible between ! output subroutine wpunch.f and the input subroutine wread.f ! ! (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, ! and Gregory W. Hammett. ALL RIGHTS RESERVED. ! use itg_data use gryffin_layouts use convert use mp, only: proc0, broadcast use gryffin_grid use fields use diagnostics_arrays ! subroutine arguments character(*), intent(in) :: filename ! name of file to read integer ntim ! on entry, the index of last time point previously used. ! on exit, the index of the last time point now used. integer ntimr ! the same as for ntim, but for a copy of phi(t) at ! one point with potentially higher time resolution ! for frequency calculation. integer nsteps ! (# of timesteps/nfreq) * (# of runs) ! for runs that have been restarted. integer nfile ! number of restart files ! output subroutine arguments: real tim ! final time when dt is the same for all ky's. real dt1 ! time step when dt is the same for all ky's. ! local variables: integer iou ! fortran i/o unit number of the *.res or *.resb file integer i,id,l,m,n,n1,nt1,nt2,ntf1,ntf2 integer :: mloc, nloc, ploc integer, dimension(:), allocatable :: m_min, m_max logical unit_exist,file_open real xxx real, dimension(:), allocatable :: timeall, dtall real, dimension(:,:), allocatable :: rdiffall ! For backwards compatibility with old cosine/sine coding ! (where the amplitude of the sine part is the negative of the ! amplitude of the imaginary part), reverse the sign of the imaginary ! parts of some of the fields. ! find a free fortran unit number to use for i/o: do iou=20,99 inquire(unit=iou,exist=unit_exist,opened=file_open) if(unit_exist .and. .not. file_open) goto 33 enddo write(*,*) 'ERROR in io_ascii.f90 or io_binary.f90, in subroutine wread' write(*,*) 'unable to find a free fortran i/o unit number' 33 continue open(unit=iou,file=filename,status='old',form='formatted',err=98) ! read first few lines of data to determine array sizes: read(iou,1,err=96)ld,ldb,kd,nmin,nmax,nd,lin,kdpass read(iou,1,err=96)nstp,ne,nfreq,ntotal,md,nspecies,ikx,iperiod read(iou,1)nparmom,nperpmom,nemom,iodd,iflr,iphi00,ifilter,iexp,igradon allocate(m_min(nd), m_max(nd)) read(iou,1)(m_min(n),n=1,nd) read(iou,1)(m_max(n),n=1,nd) ! need to read beta_e before allocating arrays (so we know whether ! to allocate electromagnetic arrays): read(iou,2)tim,shr,qsf,epsn,eps,epse,nueeff,alpha, & etai,nuii,rmu1,xxx,beta_e,tiovte,vy,vz,rmime,etae read(iou,2)dt0,dt1,x0,y0,z0,xp ! ! If mmin is not allocated, 'allocate' has not been called. ! if(.not.allocated(mmin)) then call init_gryffin_layouts (md, nd) call alloc_grid ! see if the above init_gryffin_layouts is enough for use in n2cres, ! without needing to call init_grid: ! Initialize grid quantities: ! call init_grid (epse, iperiod, lin, linlay) call alloc_itg ! Allocate arrays (mostly diagnostics) call alloc_fields ! Allocate main 3-D arrays nsteps = max(nsteps, ntotal) ! ! This could be done more reliably by storing two new variables. ! One would count the total number of time steps taken to that point, ! and would replace nsteps, and the other would count the total number ! of times prstep was called, through all the continuation runs. This ! would replace nfile*ne here. ! call alloc_diagnostics(nfile*ne, nsteps) allocate (rdiff(m_low:m_alloc, n_low:n_alloc)) endif mmin = m_min mmax = m_max deallocate(m_min, m_max) ! now read the rest of the file n1=nspecies nt1=ntim+1 nt2=ntim+ne ntim=ntim+ne ntf1=ntimr+1 ntf2=ntimr+ntotal ntimr=ntimr+ntotal mloc = m_high - m_low + 1 nloc = n_high - n_low + 1 ! ploc > 0 only if both mloc and nloc are ploc = min(mloc, nloc) ! eventually, put the following arrays in a more logical order: ! For now, keep backward compability. Put new arrays at end: allocate(ri4(2,ld,md,nd,n1)) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(density, ri4(:,:,m_low:m_high,n_low:n_high,:)) allocate (ri3(2,ld,md,nd)) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) ri3(2,:,:,:) = -ri3(2,:,:,:) if(ploc > 0) call r2c (potential, ri3(:,:,m_low:m_high,n_low:n_high)) deallocate (ri3) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) !t_par read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(t_par, ri4(:,:,m_low:m_high,n_low:n_high,:)) ! redundant write, just to keep other arrays in the same position ! in the file for backwards compatability, for now: read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(u_par, ri4(:,:,m_low:m_high,n_low:n_high,:)) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(t_perp, ri4(:,:,m_low:m_high,n_low:n_high,:)) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(q_perp, ri4(:,:,m_low:m_high,n_low:n_high,:)) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) ! if(ploc > 0) call r2c(q_par, ri4(:,:,m_low:m_high,n_low:n_high,:)) ! more redundant writes, for backwards compatibility for now: read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) ! call r2c(q_par, ri4) read(iou,2)((((ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) read(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ri4(2,:,:,:,:) = -ri4(2,:,:,:,:) if(ploc > 0) call r2c(q_par, ri4(:,:,m_low:m_high,n_low:n_high,:)) deallocate (ri4) ! The above replace 3rd and 4th perp moments, not yet used toroidally: ! read(iou,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! read(iou,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! read(iou,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! read(iou,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) if(epse > 0.) then allocate (ri3(2,kdpass, md, nd)) read(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) if(ploc > 0) call r2c(e_density, ri3(:,:,m_low:m_high,n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) if(ploc > 0) call r2c(e_p, ri3(:,:,m_low:m_high,n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) if(ploc > 0) call r2c(e_r, ri3(:,:,m_low:m_high,n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) if(ploc > 0) call r2c(e_t, ri3(:,:,m_low:m_high,n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) if(ploc > 0) call r2c(phi_ba, ri3(:,:,m_low:m_high,n_low:n_high)) deallocate (ri3) allocate (ri3(2,ld,md,nd)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c(phi_bk, ri3(:,:,m_low:m_high,n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c(e_denk, ri3(:,:,m_low:m_high,n_low:n_high)) deallocate (ri3) allocate (ri2(2,md,nd)) read(iou,2)((ri2(1,m,n),m=1,md),n=1,nd) read(iou,2)((ri2(2,m,n),m=1,md),n=1,nd) if(ploc > 0) call r2c (phiav_old, ri2(:,m_low:m_high,n_low:n_high)) deallocate (ri2) endif read(iou,2)(timo(i),i=nt1,nt2) !time history plot array read(iou,2)((wpfx(i,n),i=nt1,nt2),n=1,n1) !ion thermal flux read(iou,2)((fluxi(i,n),i=nt1,nt2),n=1,n1) !ion particle flux read(iou,2)(fluxe(i),i=nt1,nt2) !electron particle flux read(iou,2)(qfluxe(i),i=nt1,nt2) !electron thermal flux read(iou,2)(gamx(i),i=nt1,nt2) !total growth rate? allocate(ri2(md, nd, nt1:nt2)) read(iou,2)(((ri2(m,n,i),m=1,md),n=1,nd),i=nt1,nt2) !gamma vs. mode if(ploc > 0) mgamx(:,:,nt1:nt2) = ri2(m_low:m_high,n_low:n_high,nt1:nt2) deallocate (ri2) ! read(iou,2)(wrhx(i),i=nt1,nt2) ! read(iou,2)(wrhy(i),i=nt1,nt2) ! read(iou,2)(wdux(i),i=nt1,nt2) ! read(iou,2)(wu2x(i),i=nt1,nt2) ! read(iou,2)(wp2x(i),i=nt1,nt2) ! read(iou,2)(wv2x(i),i=nt1,nt2) ! read(iou,2)(wlux(i),i=nt1,nt2) ! read(iou,2)(wdpx(i),i=nt1,nt2) ! read(iou,2)(wdvx(i),i=nt1,nt2) read(iou,2)(wenx(i),i=nt1,nt2) ! total energy ! read(iou,2)(wdfx(i),i=nt1,nt2) read(iou,2)(pcerr(i),i=nt1,nt2) ! energy conservation allocate(ri2(nt1:nt2, md, nd)) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !energy vs. mode if(ploc > 0) wtif(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !kin.en. vs. mode if(ploc > 0) wkif(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !p^2 vs mode if(ploc > 0) wpif(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !phi*phi vs mode if(ploc > 0) psp(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !energy vs. mode if(ploc > 0) wenrk(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !eta drive v mode if(ploc > 0) dketa(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !par damp. v mode if(ploc > 0) dkpar(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor drive v mode if(ploc > 0) dktor(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor damp. v mode if(ploc > 0) dktdp(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !wupsi v mode if(ploc > 0) wkups(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tot. drive v mode if(ploc > 0) dktot(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)((cx(i,n),i=nt1,nt2),n=1,4*nd) !phi(x) not aved. in y read(iou,2)((cy(i,m),i=nt1,nt2),m=1,4*md) !phi(y) not aved. in x read(iou,2)((cz(i,n),i=nt1,nt2),n=1,ld) !phi(z) not aved. in x,y read(iou,2)((czc(i,n),i=nt1,nt2),n=1,ld) !phi(0)phi(z) aved. in x,y read(iou,2)((czcn(i,n),i=nt1,nt2),n=1,ld) !phi(z)^2 aved. in x,y read(iou,2)((czn(i,n),i=nt1,nt2),n=1,ld) !ne(0)ne(z) aved. in x,y read(iou,2)((cznn(i,n),i=nt1,nt2),n=1,ld) !ne(z)^2 aved. in x,y read(iou,2)(ct(i),i=nt1,nt2) !phi(0)phi(t) aved. in x,y read(iou,1) ninterv read(iou,2)((phi00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. potential read(iou,2)((den00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. density read(iou,2)((upar00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. velocity read(iou,2)((tpar00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. t_par read(iou,2)((tperp00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. t_perp allocate (ri1(2,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(phih0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(uparh0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(denh0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(tparh0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(tperph0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(qparh0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(qperph0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(veh0, ri1(:,nt1:nt2)) read(iou,2)(ri1(1,i),i=nt1,nt2) read(iou,2)(ri1(2,i),i=nt1,nt2) call r2c(uparph0, ri1(:,nt1:nt2)) deallocate (ri1) read(iou,2)((drt(id,i),id=1,14),i=nt1,nt2) !drive terms read(iou,2)(phirms(i),i=nt1,nt2) ! rms phi ! read(iou,2)(((wvif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) ! read(iou,2)(((wplx(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) read(iou,2)(wakx(i),i=nt1,nt2) ! mode width in x read(iou,2)(waky(i),i=nt1,nt2) ! mode width in y read(iou,2)(wxsp(i),i=nt1,nt2) ! mode spread read(iou,2)(utor(i),i=nt1,nt2) ! toroidal flow read(iou,2)(upol(i),i=nt1,nt2) ! poloidal flow read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) ! amplitude if(ploc > 0) grtmx(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) deallocate(ri2) allocate (ri3(2, ntf1:ntf2, md, nd)) read(iou,2)(((ri3(1,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd) ! Phi(t) read(iou,2)(((ri3(2,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd) if(ploc > 0) call r2c (utim(ntf1:ntf2,:,:), ri3(:,ntf1:ntf2,m_low:m_high,n_low:n_high)) deallocate (ri3) read(iou,2)(timf(i),i=ntf1,ntf2) ! time base for freq(t) and phi(t) read(iou,2)avgflag read(iou,3)date read(iou,4)note ! ! Write multi-species stuff at the end for backwards compatibility: ! read(iou,2)(rmass(n),n=1,nspecies) read(iou,2)(charge(n),n=1,nspecies) read(iou,2)(n_I(n),n=1,nspecies) read(iou,2)(Ln(n),n=1,nspecies) read(iou,2)(eta(n),n=1,nspecies) read(iou,2)(eta_par(n),n=1,nspecies) read(iou,2)(tau(n),n=1,nspecies) allocate(ri2(nt1:nt2, md, nd)) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) if(ploc > 0) phisq(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) allocate(ri3(nt1:nt2, md, nd, nspecies)) read(iou,2)((((ri3(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),l=1,nspecies) if(ploc > 0) fluximn(nt1:nt2,:,:,:) = ri3(nt1:nt2,m_low:m_high,n_low:n_high,:) read(iou,2)((((ri3(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),l=1,nspecies) if(ploc > 0) qfluximn(nt1:nt2,:,:,:) = ri3(nt1:nt2,m_low:m_high,n_low:n_high,:) deallocate(ri3) if(epse > 0. .or. beta_e > 0.) then read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) if(ploc > 0) fluxemn(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) read(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) if(ploc > 0) qfluxemn(nt1:nt2,:,:) = ri2(nt1:nt2,m_low:m_high,n_low:n_high) endif deallocate(ri2) ! more new variables tacked on at the end for upward compabitility: read(iou,1) icrit,ntimf allocate(timeall(md),dtall(md)) read(iou,2)(timeall(m),m=1,md) read(iou,2)(dtall(m),m=1,md) if(mloc > 0) then time(:) = timeall(m_low:m_high) dt(:) = dtall(m_low:m_high) endif deallocate(timeall,dtall) allocate(ri2(ld, md, nd)) read(iou,2)(((ri2(l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) rkperp2 = ri2(:,m_low:m_high,n_low:n_high) deallocate(ri2) allocate(rdiffall(md, nd)) read(iou,2) ((rdiffall(m,n),m=1,md),n=1,nd) if(ploc > 0) rdiff(:,:) = rdiffall(m_low:m_high,n_low:n_high) deallocate(rdiffall) read(iou,2) chi_int read(iou,1) malias, nalias if(beta_e > 0.) then read(iou,1) inuei, semi_imp read(iou,2) beta_e, meovmi, nuei allocate (ri3(2,ld,md,nd)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c (n_e, ri3(:,:,m_low:m_high, n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c (u_e, ri3(:,:,m_low:m_high, n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c (tpar_e, ri3(:,:,m_low:m_high, n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c (tperp_e, ri3(:,:,m_low:m_high, n_low:n_high)) read(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) read(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) if(ploc > 0) call r2c (apar, ri3(:,:,m_low:m_high, n_low:n_high)) deallocate (ri3) endif close(unit=iou) return 96 write(*,*) 'ERROR in io_ascii.f90 or io_binary.f90, in subroutine wread' write(*,*) 'error while reading file ', filename stop 98 write(*,*) 'ERROR in io_ascii.f90 or io_binary.f90, in subroutine wread' write(*,*) 'unable to open input file ', filename stop 1 format(1x,3i10) 2 format(1x,3e22.14) 3 format(1x,a30) 4 format(1x,a80) end subroutine wread subroutine wpunch(filename,ntim,ntimr,tim,dt1) use itg_data use gryffin_layouts use convert use mp, only: proc0, broadcast use gryffin_grid use fields use diagnostics_arrays use funnel ! ! Write/read results to/from the *.res file. ! The format of the *.res file is essentially given by the ! read/writes in this subroutine, and must be compatible between ! output subroutine wpunch.f and the input subroutine wread.f ! ! (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, ! and Gregory W. Hammett. ALL RIGHTS RESERVED. ! ! subroutine arguments character(*), intent(in) :: filename ! name of file to read integer ntim ! on entry, the index of last time point previously used. ! on exit, the index of the last time point now used. integer ntimr ! the same as for ntim, but for a copy of phi(t) at ! one point with potentially higher time resolution ! for frequency calculation. ! output subroutine arguments: real tim ! final time when dt is the same for all ky's. real dt1 ! time step when dt is the same for all ky's. ! local variables: integer iou ! fortran i/o unit number of the *.res or *.resb file integer i,id,l,m,n,n1,nt1,nt2,ntf1,ntf2 logical unit_exist,file_open real xxx real, dimension(:,:,:), allocatable :: rall3 real, dimension(:,:), allocatable :: rall2 real, dimension(:), allocatable :: rall1 ! defaults values for backwards compatibility: xxx=0.0 if(idx_local(m_lo, md, 1)) then tim=time(md) dt1=dt(md) endif call broadcast(tim, proc_id(m_lo, md, 1)) call broadcast(dt1, proc_id(m_lo, md, 1)) ! find a free fortran unit number to use for i/o: do iou=20,99 inquire(unit=iou,exist=unit_exist,opened=file_open) if(unit_exist .and. .not. file_open) goto 33 enddo write(*,*) 'ERROR in io_ascii.f90 or io_binary.f90, in subroutine wpunch' write(*,*) 'unable to find a free fortran i/o unit number' 33 continue ! start writing to the *.res or *.resb file: open(unit=iou,file=filename,status='unknown',form='formatted',err=98) if(proc0) then write(iou,1)ld,ldb,kd,nmin,nmax,nd,lin,kdpass write(iou,1)nstp,ne,nfreq,ntotal,md,nspecies,ikx,iperiod write(iou,1)nparmom,nperpmom,nemom,iodd,iflr,iphi00,ifilter,iexp,igradon write(iou,1)(mmin(n),n=1,nd) write(iou,1)(mmax(n),n=1,nd) write(iou,2)tim,shr,qsf,epsn,eps,epse,nueeff,alpha, & etai,nuii,rmu1,xxx,beta_e,tiovte,vy,vz,rmime,etae ! Starting Feb. 18, 1999, beta_e is written here, early in ! the file so we know whether to allocate the electromagnetic ! arrays. beta_e is written again later in the file. write(iou,2)dt0,dt1,x0,y0,z0,xp endif n1=nspecies nt1=ntim+1 nt2=ntim+ne ntim=ntim+ne ntf1=ntimr+1 ntf2=ntimr+ntotal ntimr=ntimr+ntotal ! eventually, put the following arrays in a more logical order: ! For now, keep backward compability. Put new arrays at end: if(proc0) allocate(ri4(2,ld,md,nd,n1)) ! density call funnel20(density, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif if(proc0) allocate(ri3(2,ld,md,nd)) ! potential call funnel20(potential, ri3, 1, ld) if(proc0) then write(iou,2)(((-ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) deallocate (ri3) endif ! t_par call funnel20(t_par, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif ! redundant write, just to keep other arrays in the same position ! in the file for backwards compatability, for now: ! u_par call funnel20(u_par, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif ! t_perp call funnel20(t_perp, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif ! q_perp call funnel20(q_perp, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif ! q_par call funnel20(q_par, ri4) if(proc0) then write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! more redundant writes, for backwards compatibility for now: write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((-ri4(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) write(iou,2)((((ri4(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) endif if(proc0) deallocate (ri4) ! The above replace 3rd and 4th perp moments, not yet used toroidally: ! write(iou,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! write(iou,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! write(iou,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) ! write(iou,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1) if(epse > 0.) then if(proc0) allocate (ri3(2, kdpass, md, nd)) ! e_density call funnel20(e_density, ri3, 1, kdpass) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) endif ! e_p call funnel20(e_p, ri3, 1, kdpass) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) endif ! e_r call funnel20(e_r, ri3, 1, kdpass) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) endif ! e_t call funnel20(e_t, ri3, 1, kdpass) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) endif ! phi_ba call funnel20(phi_ba, ri3, 1, kdpass) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd) endif if(proc0) deallocate (ri3) if(proc0) allocate (ri3(2,ld,md,nd)) ! phi_bk call funnel20(phi_bk, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif ! e_denk call funnel20(e_denk, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif if(proc0) deallocate (ri3) if(proc0) allocate (ri2(2,md, nd)) ! phiav_old call funnel20(phiav_old, ri2) if(proc0) then write(iou,2)((ri2(1,m,n),m=1,md),n=1,nd) write(iou,2)((ri2(2,m,n),m=1,md),n=1,nd) endif if(proc0) deallocate (ri2) endif if(proc0) then write(iou,2)(timo(i),i=nt1,nt2) !time history plot array write(iou,2)((wpfx(i,n),i=nt1,nt2),n=1,n1) !ion thermal flux write(iou,2)((fluxi(i,n),i=nt1,nt2),n=1,n1) !ion particle flux write(iou,2)(fluxe(i),i=nt1,nt2) !electron particle flux write(iou,2)(qfluxe(i),i=nt1,nt2) !electron thermal flux write(iou,2)(gamx(i),i=nt1,nt2) !total growth rate? endif if(proc0) allocate(rall3(md, nd, nt1:nt2)) call mfunnel3(mgamx, rall3, nt1, nt2) if(proc0) then write(iou,2)(((rall3(m,n,i),m=1,md),n=1,nd),i=nt1,nt2) !gamma vs. mode endif if(proc0) deallocate (rall3) if(proc0) then ! write(iou,2)(wrhx(i),i=nt1,nt2) ! write(iou,2)(wrhy(i),i=nt1,nt2) ! write(iou,2)(wdux(i),i=nt1,nt2) ! write(iou,2)(wu2x(i),i=nt1,nt2) ! write(iou,2)(wp2x(i),i=nt1,nt2) ! write(iou,2)(wv2x(i),i=nt1,nt2) ! write(iou,2)(wlux(i),i=nt1,nt2) ! write(iou,2)(wdpx(i),i=nt1,nt2) ! write(iou,2)(wdvx(i),i=nt1,nt2) write(iou,2)(wenx(i),i=nt1,nt2) ! total energy ! write(iou,2)(wdfx(i),i=nt1,nt2) write(iou,2)(pcerr(i),i=nt1,nt2) ! energy conservation endif if(proc0) allocate(ri2(nt1:nt2, md, nd)) call funnel20(wtif, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !energy vs. mode call funnel20(wkif, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !kin.en. vs. mode call funnel20(wpif, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !p^2 vs mode call funnel20(psp, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !phi*phi vs mode call funnel20(wenrk, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !energy vs. mode call funnel20(dketa, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !eta drive v mode call funnel20(dkpar, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !par damp. v mode call funnel20(dktor, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor drive v mode call funnel20(dktdp, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor damp. v mode call funnel20(wkups, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !wupsi v mode call funnel20(dktot, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tot. drive v mode ! it would be better to send only to proc0, but for now this is ok... do n = m_lo%nl_world, m_lo%nu_world call broadcast(phi00(:,n), proc_id(m_lo, 1, n)) call broadcast(den00(:,n), proc_id(m_lo, 1, n)) call broadcast(upar00(:,n), proc_id(m_lo, 1, n)) call broadcast(tpar00(:,n), proc_id(m_lo, 1, n)) call broadcast(tperp00(:,n), proc_id(m_lo, 1, n)) enddo if(proc0) then write(iou,2)((cx(i,n),i=nt1,nt2),n=1,4*nd) !phi(x) not aved. in y write(iou,2)((cy(i,m),i=nt1,nt2),m=1,4*md) !phi(y) not aved. in x write(iou,2)((cz(i,n),i=nt1,nt2),n=1,ld) !phi(z) not aved. in x,y write(iou,2)((czc(i,n),i=nt1,nt2),n=1,ld) !phi(0)phi(z) aved. in x,y write(iou,2)((czcn(i,n),i=nt1,nt2),n=1,ld) !phi(z)^2 aved. in x,y write(iou,2)((czn(i,n),i=nt1,nt2),n=1,ld) !ne(0)ne(z) aved. in x,y write(iou,2)((cznn(i,n),i=nt1,nt2),n=1,ld) !ne(z)^2 aved. in x,y write(iou,2)(ct(i),i=nt1,nt2) !phi(0)phi(t) aved. in x,y write(iou,1) ninterv write(iou,2)((phi00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. potential write(iou,2)((den00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. density write(iou,2)((upar00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. velocity write(iou,2)((tpar00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. t_par write(iou,2)((tperp00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. t_perp allocate (ri1(2,nt1:nt2)) call c2r(phih0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(uparh0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(denh0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(tparh0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(tperph0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(qparh0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(qperph0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(veh0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) call c2r(uparph0(nt1:nt2), ri1) write(iou,2)(ri1(1,i),i=nt1,nt2) write(iou,2)(ri1(2,i),i=nt1,nt2) deallocate (ri1) write(iou,2)((drt(id,i),id=1,14),i=nt1,nt2) !drive terms write(iou,2)(phirms(i),i=nt1,nt2) ! rms phi endif ! write(iou,2)(((wvif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) ! write(iou,2)(((wplx(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) if(proc0) then write(iou,2)(wakx(i),i=nt1,nt2) ! mode width in x write(iou,2)(waky(i),i=nt1,nt2) ! mode width in y write(iou,2)(wxsp(i),i=nt1,nt2) ! mode spread write(iou,2)(utor(i),i=nt1,nt2) ! toroidal flow write(iou,2)(upol(i),i=nt1,nt2) ! poloidal flow endif call funnel20(grtmx, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) ! amplitude if(proc0) deallocate(ri2) if(proc0) allocate (ri3(2,ntf1:ntf2,md,nd)) call funnel20(utim, ri3, ntf1, ntf2) if(proc0) then write(iou,2)(((ri3(1,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd) ! Phi(t) write(iou,2)(((ri3(2,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd) endif if(proc0) deallocate (ri3) if(proc0) then write(iou,2)(timf(i),i=ntf1,ntf2) ! time base for freq(t) and phi(t) write(iou,2)avgflag write(iou,3)date write(iou,4)note ! ! Write multi-species stuff at the end for backwards compatibility: ! write(iou,2)(rmass(n),n=1,nspecies) write(iou,2)(charge(n),n=1,nspecies) write(iou,2)(n_I(n),n=1,nspecies) write(iou,2)(Ln(n),n=1,nspecies) write(iou,2)(eta(n),n=1,nspecies) write(iou,2)(eta_par(n),n=1,nspecies) write(iou,2)(tau(n),n=1,nspecies) endif if(proc0) allocate(ri2(nt1:nt2, md, nd)) call funnel20(phisq, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) if(proc0) deallocate(ri2) if(proc0) allocate(ri3(nt1:nt2, md, nd, nspecies)) call funnel20(fluximn, ri3, nt1, nt2) if(proc0) write(iou,2)((((ri3(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),l=1,nspecies) call funnel20(qfluximn, ri3, nt1, nt2) if(proc0) write(iou,2)((((ri3(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),l=1,nspecies) if(proc0) deallocate(ri3) if(proc0) allocate(ri2(nt1:nt2, md, nd)) if(epse > 0. .or. beta_e > 0.) then call funnel20(fluxemn, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) call funnel20(qfluxemn, ri2, nt1, nt2) if(proc0) write(iou,2)(((ri2(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) endif if(proc0) deallocate(ri2) ! more new variables tacked on at the end for upward compabitility: if(proc0) write(iou,1) icrit,ntimf if(proc0) allocate(rall1(md)) call funnel20(time, rall1) if(proc0) write(iou,2)(rall1(m),m=1,md) call funnel20(dt, rall1) if(proc0) write(iou,2)(rall1(m),m=1,md) if(proc0) deallocate(rall1) if(proc0) allocate(rall3(ld,md,nd)) call funnel20(rkperp2, rall3, 1, ld) if(proc0) write(iou,2)(((rall3(l,m,n),l=1,ld),m=1,md),n=1,nd) if(proc0) deallocate(rall3) if(proc0) allocate(rall2(md, nd)) call funnel20(rdiff, rall2) if(proc0) write(iou,2) ((rall2(m,n),m=1,md),n=1,nd) if(proc0) deallocate(rall2) if(proc0) then write(iou,2) chi_int write(iou,1) malias, nalias endif if (beta_e > 0.) then if(proc0) then write(iou,1) inuei, semi_imp write(iou,2) beta_e, meovmi, nuei endif if(proc0) allocate(ri3(2,ld,md,nd)) ! n_e call funnel20(n_e, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif ! u_e call funnel20(u_e, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif ! tpar_e call funnel20(tpar_e, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif ! tperp_e call funnel20(tperp_e, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif ! apar call funnel20(apar, ri3, 1, ld) if(proc0) then write(iou,2)(((ri3(1,l,m,n),l=1,ld),m=1,md),n=1,nd) write(iou,2)(((ri3(2,l,m,n),l=1,ld),m=1,md),n=1,nd) endif if(proc0) deallocate (ri3) endif close(unit=iou) return 98 write(*,*) 'ERROR in io_ascii.f90 or io_binary.f90, in subroutine wpunch' write(*,*) 'unable to open input file ', filename stop 1 format(1x,3i10) 2 format(1x,3e22.14) 3 format(1x,a30) 4 format(1x,a80) end subroutine wpunch end module io_ascii