out2_nc.f.html

out2_nc.f.html
Source file: out2_nc.f
Directory: /home/rjl/git/rjleveque/clawpack-4.x/clawpack/2d/lib
Converted: Sun May 15 2011 at 19:15:43 using clawcode2html
This documentation file will not reflect any later changes in the source file.
 
c
c
c =========================================================
      subroutine out2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
     &                 dx,dy,q,t,iframe,aux,maux)
c =========================================================
c
c     # Output the results for a general system of conservation laws
c     # in 2 dimensions
c
c     # Write the results to the file fort.q.nc</font>
<font color="blue">c     # Use format required by matlab script  plotclaw2.m</font>
<font color="blue">c     # The same format is used by the amrclaw package.</font>
<font color="blue">c     # Here it's adapted to output just the single grid.</font>
<font color="blue">c     # set outaux = .true. to also output the aux arrays to fort.a<iframe></font>
<font color="blue">c</font>

<font color="blue">c</font>
<font color="blue">c -----------------------------------------------------</font>
<font color="blue">c     Routine to write netcdf files in the classic format</font>
<font color="blue">!        #jj-2011.03.29</font>
<font color="blue">!        # Each file written by the fortran code has </font>
<font color="blue">!        # Dimensions:</font>
<font color="blue">!        #           timedimension : UNLIMITED</font>
<font color="blue">!        #           meqn          : The number of equations</font>
<font color="blue">!        #           dimx_<gridno> : X dimension for grid number <gridno></font>
<font color="blue">!        #           dimy_<gridno> : Y dimension for grid number <gridno></font>
<font color="blue">!        # Variables:</font>
<font color="blue">!        #           timedimension : Stores the time of the frame</font>
<font color="blue">!        #           ngrids        : Number of grids in this frame</font>
<font color="blue">!        #           naux          : Number of Auxilary Variables</font>
<font color="blue">!        #           ndim          : Number of Dimensions in the frame</font>
<font color="blue">!        #           grid_<gridno> : A grid of (dimx,dimy,meqn)</font>
<font color="blue">!        # Attributes:</font>
<font color="blue">!        # (grid_<no>) gridno      : The number of this grid <grid_no></font>
<font color="blue">!        #           level         : The AMR level</font>
<font color="blue">!        #           dim_names     : a list of dimensions [dimx,dimy]</font>
<font color="blue">!        #           dim<x,y>.low  : The lowest dimension value </font>
<font color="blue">!        #           dim<x,y>.d    : The distance between grid points </font>
<font color="blue">c -----------------------------------------------------</font>

<font color="blue">c</font>
      implicit double precision (a-h,o-z)
      include 'netcdf.inc'
      
      dimension   q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      dimension aux(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, maux)
      character*10 fname1, fname2, fname3
      logical outaux,do_ascii
            
      real(kind=8) time
      integer ncid,rcode
      integer timeid,tVarID,meqnID,ngridsVarID,nauxVarID,ndimVarID
      integer dimxid,dimyid,xlowid,ylowid,dxid,dyid
      integer gridid
      integer ntimes
      character*2 gridstr
      character*40 dim_names
      REAL(kind=8), ALLOCATABLE  ::grid(:,:,:)
      integer nx, ny
      real xlow,ylow
      

      nx=mx
      ny=my
      outaux = .false.
      do_ascii = .false.
<font color="blue">c</font>
<font color="blue">c     # first create the file name and open file</font>
<font color="blue">c</font>
         fname1 = 'fort.qxxxx'
         fname2 = 'fort.txxxx'
         fname3 = 'fort.axxxx'
         nstp = iframe
         do 55 ipos = 10, 7, -1
            idigit = mod(nstp,10)
            fname1(ipos:ipos) = char(ichar('0') + idigit)
            fname2(ipos:ipos) = char(ichar('0') + idigit)
            fname3(ipos:ipos) = char(ichar('0') + idigit)
            nstp = nstp / 10
 55      continue

         if (do_ascii) then  open(unit=50,file=fname1,status='unknown',form='formatted')
         if (do_ascii) then  open(unit=60,file=fname2,status='unknown',form='formatted')

<font color="blue">c</font>
<font color="blue">c     # the following parameters are used in amrclaw where there are</font>
<font color="blue">c     # multiple grids.  Here they are all set to 1:</font>
      ngrids = 1
      mptr = 1
      level = 1

      ntimes=1
<font color="blue">c        write(50,1001) mptr,level,mx,my</font>
<font color="blue">c        write(50,1002) xlower,ylower,dx,dy</font>
        <font color="blue">!!!!Define netcdf file</font>
         rcode=NF_CREATE(fname1//'.nc',NF_NOCLOBBER,ncid)
         if(rcode.ne.NF_NOERR) print *,'ERROR OPENING NETCDF FILE'
         rcode=NF_DEF_DIM(ncid,'timedimension',NF_UNLIMITED,timeid)
         rcode=NF_DEF_VAR(ncid,'timedimension',NF_DOUBLE,1,timeid,
     &   tVarID)
         rcode=NF_DEF_DIM(ncid,'meqn',meqn,meqnid)
         rcode=NF_DEF_VAR(ncid,'ngrids',NF_INT,0,0,ngridsVarID)
         rcode=NF_DEF_VAR(ncid,'naux',NF_INT,0,0,nauxVarID)
         rcode=NF_DEF_VAR(ncid,'ndim',NF_INT,0,0,ndimVarID)
        write(gridstr,67) mptr
              
67            format(I2.2)              
              rcode=NF_DEF_DIM(ncid,'dimx_'//trim(gridstr),nx,dimxid)
              if(rcode.ne.NF_NOERR) print *,'ERROR  DEFINE DIMS'
              rcode=NF_DEF_DIM(ncid,'dimy_'//trim(gridstr),ny,dimyid)
              if(rcode.ne.NF_NOERR) print *,'ERROR  DEFINE DIMS'
              

              rcode=NF_DEF_Var(ncid,'grid_'//trim(gridstr),NF_DOUBLE,4,
     &              (/dimxid,dimyid,meqnid,timeid/),gridid)
               if(rcode.ne.NF_NOERR) print *,'ERROR  DEFINE VAR'
               
              rcode=NF_PUT_ATT_INT(ncid,gridid,'gridno',NF_INT,1,
     &              mptr)
     
              rcode=NF_PUT_ATT_INT(ncid,gridid,'level',NF_INT,1,level)
              
              dim_names="['dimx','dimy']"
              rcode=NF_PUT_ATT_TEXT(ncid,gridid,'dim_names',
     &         LEN_TRIM(dim_names),TRIM(dim_names))
     
              rcode=NF_PUT_ATT_DOUBLE(ncid,gridid,'dimx.lower',NF_FLOAT,
     &              1,xlower)     
              rcode=NF_PUT_ATT_DOUBLE(ncid,gridid,'dimy.lower',NF_FLOAT,
     &              1,ylower)
     
              rcode=NF_PUT_ATT_DOUBLE(ncid,gridid,'dimx.d',NF_FLOAT,1,
     &          dx)
              rcode=NF_PUT_ATT_DOUBLE(ncid,gridid,'dimy.d',NF_FLOAT,1,
     &          dy) 
     
              rcode=NF_ENDDEF(ncid)
            write(0,*) "xlower, ylower:",xlower,ylower,sizeof(xlower)        
      allocate(grid(nx,ny,meqn))
      time=t      
<font color="blue">c</font>
      do j=1,my
        do i=1,mx
          do m=1,meqn
<font color="blue">c            # exponents with more than 2 digits cause problems reading</font>
<font color="blue">c            # into matlab... reset tiny values to zero:</font>
             if (dabs(q(i,j,m)) .lt. 1d-99) q(i,j,m) = 0.d0
             grid(i,j,m)=q(i,j,m)
             enddo
              if (do_ascii) then write(50,1005) (grid(i,j,m), m=1,meqn)
         enddo
         if (do_ascii) then write(50,*) ' '
      enddo
      if (do_ascii) then write(50,*) ' '
      
      rcode=NF_PUT_VARA_DOUBLE(ncid,gridid,(/1,1,1,1/),
     & (/nx,ny,meqn,1/),grid)
     
      deallocate(grid)
      
      rcode=NF_PUT_VAR_DOUBLE(ncid,tVarID,t)
      if(rcode.ne.NF_NOERR) print *,'ERROR  Write Time'
      rcode=NF_PUT_VAR_INT(ncid,ngridsVarID,int(ngrids))
      if(rcode.ne.NF_NOERR) print *,'ERROR  Write GridNo'
      rcode=NF_PUT_VAR_INT(ncid,nauxVarID,maux)
      rcode=NF_PUT_VAR_INT(ncid,ndimVarID,2)
      rcode=NF_CLOSE(ncid)

 1001 format(i5,'                 grid_number',/,
     &       i5,'                 AMR_level',/,
     &       i5,'                 mx',/,
     &       i5,'                 my')

 1002 format(e26.16,'    xlow', /,
     &       e26.16,'    ylow', /,
     &       e26.16,'    dx', /,
     &       e26.16,'    dy',/)
 1005     format(4e26.16)

            
      if (outaux) then 
<font color="blue">c     # also output the aux arrays:</font>
      open(unit=70,file=fname3,status='unknown',form='formatted')
      write(70,1001) mptr,level,mx,my
      write(70,1002) xlower,ylower,dx,dy
      do 120 j=1,my
         do 110 i=1,mx
            do m=1,maux
<font color="blue">c              # exponents with more than 2 digits cause problems reading</font>
<font color="blue">c              # into matlab... reset tiny values to zero:</font>
               if (dabs(aux(i,j,m)) .lt. 1d-99) aux(i,j,m) = 0.d0
            enddo
<font color="blue">c</font>
            write(70,1005) (aux(i,j,m), m=1,maux)
<font color="blue">c</font>
  110       continue
         write(70,*) ' '
  120    continue
      write(70,*) ' '
      close(unit=70)
      endif

      if (do_ascii) then write(60,1000) t,meqn,ngrids,maux,2

 1000 format(e26.16,'    time', /,
     &       i5,'                 meqn'/,
     &       i5,'                 ngrids'/,
     &       i5,'                 maux'/,
     &       i5,'                 ndim'/,/)
<font color="blue">!c</font>
      if (do_ascii) then close(unit=50)
      if (do_ascii) then close(unit=60)

      return
      end
</pre></html>