write(unit=69,rec=record) char(int(127.0d0+scalingValue*ez(i,j,k)))                  record = record + 1               end do            end do         end do       	 !! close the output file for this iteration:         close(unit=69) 	 !! write the dimensions and name of the output file for this 	 !!     iteration to the viz file:         write(unit=33,fmt=334) nx+1,ny+1,nz,filename334      format(I4.4,"x",I4.4,"x",I4.4," ",A14)	 !! write identification of the corners of the mesh and the max	 !!     and min values for this iteration to the viz file:         write(unit=33,fmt=335) dx*nx,dy*ny,dz*nz,min,max335      format("bbox: 0.0 0.0 0.0 ",F7.4,F7.4,F7.4," ",F8.4," ",F8.4)      end if !! end bob output section      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      !! Compute the stimulus: a plane wave emanates from the x=0 face:      !!     The length of the guide lies in the x-direction, the width of the       !!     guide lies in the y-direction, and the height of the guide lies      !!     in the z-direction.  So the guide is sourced by all the ez       !!     components on the stimulus face.        stimulus = sin(omega*currentSimulatedTime)      do i=0,0,1         do j=0,ny,1            do k=0,(nz-1),1               ez(i,j,k) = stimulus            end do         end do      end do          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      !! Update the interior of the mesh:      !!    all vector components except those on the faces of the mesh.        !!      !! Update all the H field vector components within the mesh:      !!     Since all H vectors are internal, all H values are updated here.      !!     Note that the normal H vectors on the faces of the mesh are not       !!     computed here, and in fact were never allocated -- the normal       !!     H components on the faces of the mesh are never used to update       !!     any other value, so they are left out of the memory allocation        !!     entirely.             !! Update the hx values:      do i=0,(nx-2),1         do j=0,(ny-1),1            do k=0,(nz-1),1               hx(i,j,k) = hx(i,j,k) + dtmudz*(ey(i+1,j,k+1) - ey(i+1,j,k)) &                                   & - dtmudy*(ez(i+1,j+1,k) - ez(i+1,j,k))            end do         end do      end do          !! Update the hy values:      do i=0,(nx-1),1         do j=0,(ny-2),1            do k=0,(nz-1),1               hy(i,j,k) = hy(i,j,k) + dtmudx*(ez(i+1,j+1,k) - ez(i,j+1,k)) &                                   & - dtmudz*(ex(i,j+1,k+1) - ex(i,j+1,k))            end do         end do      end do          !! Update the hz values:      do i=0,(nx-1),1         do j=0,(ny-1),1            do k=0,(nz-2),1               hz(i,j,k) = hz(i,j,k) + dtmudy*(ex(i,j+1,k+1) - ex(i,j,k+1)) &                                   & - dtmudx*(ey(i+1,j,k+1) - ey(i,j,k+1))            end do         end do      end do          !! Update the E field vector components.        !! The values on the faces of the mesh are not updated here; they       !!      are handled by the boundary condition computation       !!      (and stimulus computation).            !! Update the ex values:      do i=0,(nx-1),1         do j=1,(ny-1),1            do k=1,(nz-1),1               ex(i,j,k) = ex(i,j,k) + dtepsdy*(hz(i,j,k-1) - hz(i,j-1,k-1)) &                                   & - dtepsdz*(hy(i,j-1,k) - hy(i,j-1,k-1))            end do         end do      end do          !! Update the ey values:      do i=1,(nx-1),1         do j=0,(ny-1),1            do k=1,(nz-1),1               ey(i,j,k) = ey(i,j,k) + dtepsdz*(hx(i-1,j,k) - hx(i-1,j,k-1)) &                                   & - dtepsdx*(hz(i,j,k-1) - hz(i-1,j,k-1))            end do         end do      end do          !! Update the ez values:      do i=1,(nx-1),1         do j=1,(ny-1),1            do k=0,(nz-1),1               ez(i,j,k) = ez(i,j,k) + dtepsdx*(hy(i,j-1,k) - hy(i-1,j-1,k)) &                                   & - dtepsdy*(hx(i-1,j,k) - hx(i-1,j-1,k))            end do         end do      end do          print*," "          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      !! Compute the boundary conditions:           !! OK, so I'm yanking your chain on this one.  The PEC condition is       !! enforced by setting the tangential E field components on all the      !! faces of the mesh to zero every timestep (except the stimulus       !! face).  But the lazy/efficient way out is to initialize those       !! vectors to zero and never compute them again, which is exactly       !! what happens in this code.    end do !! end mainloop  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! Output section:    !! The output routine is repeated one last time to write out    !! the last data computed.       !! time in simulated seconds that the simulation has progressed:   currentSimulatedTime = dt*iteration                         !! print to standard output the iteration number   !!     and current simulated time:  write(*,372) iteration, currentSimulatedTime  !! 3D data output for the last timestep:   !! create the filename for this iteration,   !!     which includes the iteration number:  !! Construct the filename:  !! The numerical part is obtained by successively  !!      dividing the current iteration by ten, and the  !!      remainders of these divisions are the digits.  !!      For example, 12345 would be converted as follows:  !!      12345/10 = 1234+5/10 (5 is the first digit in the number)  !!      1234/10 = 123 + 4/10 (4 is the next digit)  !!      ....and so on..this is performed six times to create  !!      a filename out of any 6-digit number.  temp = iteration  filename(1:4) = "f90_"  filename(10:10) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(9:9) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(8:8) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(7:7) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(6:6) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(5:5) = numbers(mod(temp,10)+1:mod(temp,10)+1)  temp = int(temp/10)  filename(11:14) = ".bob"    !! open a new data file for this iteration:  open(unit=69,file=filename,access="direct",iostat=ios,recl=1)  !! if the file fails to open, print an error message to   !!     standard output:  if(ios /= 0) then    print*,"Difficulty opening ",filename    STOP  end if    !!! find the min and max values for this iteration  min = FLT_MAX  max = -FLT_MAX  do i=0,(nx),1     do j=0,(ny),1        do k=0,(nz-1),1           if(ez(i,j,k) < min) then              min = ez(i,j,k)           end if           if(ez(i,j,k) > max) then              max = ez(i,j,k)           end if        end do     end do  end do      !!  update the tracking variables for minimum and maximum   !!  values for the entire simulation:  if(min < simulationMin) then     simulationMin = min  end if  if(max > simulationMax) then     simulationMax = max  end if    !! set norm to be max or min, whichever is greater in magnitude:  if(abs(max) > abs(min)) then     norm = abs(max)  else     norm = abs(min)  end if    !! if everything is zero, give norm a tiny value   !!     to avoid division by zero:  if(norm == 0.0) then     norm = DBL_EPSILON  end if  scalingValue = 127.0d0/norm  !! write to standard output the minimum and maximum values   !!     from this iteration and the minimum and maximum values  !!     that will be written to the bob file this iteration:  write(*, 373) min, int(127.0d0 + scalingValue * min), max, int(127.0d0 + scalingValue * max)    !! scale each ez value in the mesh to the range of integers   !!     from zero through 254 and write them to the output   !!     file for this iteration:  record = 1  do k=0,(nz-1),1     do j=0,(ny),1        do i=0,(nx),1           write(unit=69,rec=record) char(int(127.0d0+scalingValue*ez(i,j,k)))           record = record + 1        end do     end do  end do    !! close the output file for this iteration:  close(unit=69)   !! write the dimensions and name of the output file for this   !!     iteration to the viz file:  write(unit=33,fmt=334) nx+1,ny+1,nz,filename  !! write identification of the corners of the mesh and the max  !!     and min values for this iteration to the viz file:  write(unit=33,fmt=335) dx*nx,dy*ny,dz*nz,min,max  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! close the viz file for this simulation:  close(unit=33)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! write some progress notes to standard output:  print*," "  print*," "  !! print out how much memory has been allocated:   print*,"Dynamically allocated ", allocatedBytes, " bytes"  !! print out some simulation parameters:   print*,"PLOT_MODULUS = ", PLOT_MODULUS  print*,"rectangular waveguide"  print*,"Stimulus = ", FREQUENCY, " Hertz continuous plane wave"  print*," "  print*,"Meshing parameters:"  print '(I4,"x",I3,"x",I3," cells")', nx, ny, nz  print '(" dx=",F11.8," dy=",F11.8," dz=",F11.8," meters")', dx, dy, dz  print '(F11.8, " x ", F11.8, " x ", F11.8, " meter^3 simulation region")', &       & GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda  print*," "  print '(" Time simulated was ",E12.5E2," seconds, ",I5.5," timesteps")', &       & totalSimulatedTime, MAXIMUM_ITERATION  print*," "  print*,"3D output scaling parameters:"  print*,"Autoscaling every timestep"  print*," "  print*," "  !! print out simulationMin and simulationMax:   print'(" Minimum output value was : ", f10.5)', simulationMin  print'(" Maximum output value was : ", f10.5)', simulationMax  print*," "  print*," "  !! print out a bob command line, including the dimensions  !!      of the output files:  print'(" bob -cmap chengGbry.cmap -s ",I5.5,"x",I5.5,"x",I5.5," *.bob")',nx+1,ny+1,nz  print*," "  !! print out a viz command line:  print*,"viz ToyFDTD1f90.viz"  print*," "  print*," "  print*," "  print*," "end program ToyFDTD1 !! end main