efficient writing command

[LEILABLYTHE]

Hello everyone,
You may find these questions stupid but I prefer to ask and not to walk as a blind person in FORTRAN! would be really kind of you if you can help me.
I am using Fortran 90 to simulate the plasma chemistry. Simulation consists of 6 species and 32 reactions. I would like to store the information about a matrix in time which reveal the info about participation of each reactions in production or consumption of each species. I am using the following commands:

double precision :: SourceTerm(species_max,reactions_max)
! SourceTerm Matrix will be produced in a time loop

open(13,file=species_name(1)//"SourceTerm1")
open(14,file=species_name(1)//"SourceTerm2")
open(15,file=species_name(1)//"SourceTerm3")
open(16,file=species_name(1)//"SourceTerm4")

write(13,'(99(A16))',advance='yes')(trim(reaction_sign(i)),i=1, 8)
write(14,'(99(A16))',advance='yes')(trim(reaction_sign(i)),i=9,16)
write(15,'(99(A16))',advance='yes')(trim(reaction_sign(i)),i=17,24 )
write(16,'(99(A18))',advance='yes')(trim(reaction_sign(i)),i=25,32 )

write(13,'(99(1pe16.6))',advance='yes') SourceTerm(1,1:8)
write(14,'(99(1pe16.6))',advance='yes') SourceTerm(1,9:16)
write(15,'(99(1pe16.6))',advance='yes') SourceTerm(1,17:24)
write(16,'(99(1pe19.6))',advance='yes') SourceTerm(1,25:32)

close(13)
close(14)
close(15)
close(16)

There are some problems
1- These really basic commands and if I want to change the program to another would take to much time to change the writing part only! e.g. by species_name(1) I will get the info about species 1, but since I can not write info about 32 reactions in 1 file I divide it to 4 files but I want to have the produce the name of the files also in a loop but I do not know how to write numbers (1-4) after trim in the file name.
2- The time step are 1.0d-15 but I want this info only at e certain times, e.g. 1.0d-9, 1.0d-6, 1.0d-3, etc. How can I do it?
3- I copy and paste the format of the output files so I do not understand exactly what '(99(A16))',advance='yes' and '(99(1pe16.6))' means! Is there any other way to store this info? The numbers that I got are in the order of 1.0d-12 to 1.0d12.

Thanks in advance.

Leila

 

[xwb]

Is reaction_sign a character string and SourceTerm a real?

 

[LEILABLYTHE]

yes, source term Integer, and reaction_sign a character!

 

[xwb]

Minor point - you have an array called SourceTerm and a subroutine called SourceTerm. Not a good idea. It is confusing for someone else who reads the program. If it builds it may give weird results.

First the formats. Advance='yes' is on by default.

99(A16) means print up to 99 items of 16characters(max) on each line. If one of the items was ABCDEFGHIJKLMNOPQ, you would only get ABCDEFGHIJKLMNOP

99(1PE16.6) means print one decimal place. You could try the following to see what happens. Play with the 16 and the .6. See what you can get away with. Note that this is also used for spacing out the data (eg instead of 2X, I5, you could use I7)

program main
   real:: x
   x = 123.456789
   
   ! No places before decimal
   write(*, '(E16.6)') x
   ! 1 place before decimal
   write(*, '(1PE16.6)') x
   ! 2 places before decimal
   write(*, '(2PE16.6)') x
   stop
end program

Another way of storing the data: it really depends on what you want. What are you going to do with the data? Is another program going to read it in to do some more processing, is it going into a database or it is being used to plot graphs or is it going to be the results of a published paper?

I can't see a time step in your program so I can't really advise you on how to go about it. Normally it would be a parameter to the routine.

Looping the code: this isn't necessarily more efficient in machine terms. It is just more efficient in human terms. If the last E19.6 was actually meant to be E16.6, then it would be something like

   subroutine PrintSourceTerm (species_max, reaction_step, reaction_max)
      integer, intent(in):: species_max, reaction_step, reaction_max
      integer, parameter:: CHAN_START = 12
      integer:: chan, species, b, block_max, r, rbeg, rend
      character*80 filename
      
      block_max = reaction_max / reaction_step
      do species = 1, species_max
         rend = 0
         do b = 1, block_max
            ! Create a filename
            write (filename, '(A,A,I1)') species_name(species), 'SourceTerm', b
            chan = CHAN_START + b
            open (chan, file=filename)
            ! Find out start and end
            rbeg = rend + 1  ! if you want a formula use (b - 1) * reaction_step + 1
            rend = b * reaction_step
            ! Write out reaction sign, allow block size of up  to 99
            ! Trim isn't necessary as A16 truncates
            write (chan, '(99A16)') (reaction_sign(r), r = rbeg, rend)
            ! Print with shifted decimal
            write (chan, '(99(1PE16.6))') SourceTerm(1,rbeg:rend)
            close (chan)
         end do
      end do
      
   end subroutine PrintSourceTerm

 

[xwb]

Sorry, just read the bit that SourceTerm is an integer. In that case you do not need the shifted decimal. Use something like

write (chan, '(99I8)') SourceTerm(1, rbeg:rend)

 

[LEILABLYTHE]

Dear xwb,

Thank you very much for your suggestion. I need to present some of the data as an input and the rest I need to plot them. Here I am trying to combine this subroutine with the time loop but apparently I am not using an efficient way because the time loop is repeating more than once!!! Also, there are some conflicts in the name of the files.

program test
implicit none

integer :: i,icount = 0
DOUBLE PRECISION :: time, dtime
DOUBLE PRECISION :: density_old(species_max)
DOUBLE PRECISION :: species_source(species_max)
DOUBLE PRECISION :: source_terms(reactions_max), varstep
DOUBLE PRECISION :: SourceTerm(species_max,reactions_max)

DOUBLE PRECISION, PARAMETER :: density_min = 1.0d5
DOUBLE PRECISION, PARAMETER :: absvar = 0.1, relvar = 0.1
DOUBLE PRECISION, PARAMETER :: dtime_min = 1.0d-19
DOUBLE PRECISION, PARAMETER :: time_end = 1.0d4

INTEGER :: SPECIES_MAX, REACTION_STEP, REACTION_MAX
INTEGER :: CHANR, CHANS, SPECIES, B, BLOCK_REACTION, BLOCK_MAX, INTEGER :: RBEG, REND, R
INTEGER, PARAMETER :: CHAN_START = 12, icount_save = 100
character*80 filenameR
character*80 filenameS
CALL INPUTFILENAME(SPECIES_STEP, SPECIES_MAX, REACTION_STEP, REACTION_MAX)

BLOCK_MAX = CEILING(REACTION_MAX / DBLE(REACTION_STEP))

DO SPECIES = 1, SPECIES_MAX
rend = 0
do b = 1, BLOCK_MAX
write (filenameR, '(A,A,I3)') species_name(species),'Source', b
CHANR = CHAN_START + B
open (CHANR, file=filename)

RBEG = REND + 1
REND = B * REACTION_STEP
write (CHANR, '(99A16)') (REACTION_SIGN(R), R = RBEG, REND)

write (filenameS, '(99A4)') 'DENSITY',(SPECIES_NAME(R),R = RBEG, REND)
CHANS = 10*CHAN_START + B
open (CHANS, file=filename)
write (CHANS, '(99A16)') (SPECIES_NAME(R),R = RBEG, REND)

time = 0.0d0
dtime = dtime_min

do while(time .lt. time_end)
density_old) = density)
call timestep(time,dtime)

time = time + dtime

Call rates(SOURCE_TERMS=species_source),REACTION_RATES=source_terms))
call rates(SOURCE_TERMS_MATRIX=SourceTerm,)

!output
if(mod(icount,icount_save) .eq. 0) then
write (CHANR, '(A,99(I8))') time, SourceTerm(1,RBEG:REND)
write (CHANS, '(A,99(I8))') time, SPECIES_DENSITY(1,RBEG:REND)
endif
!------------------------------
!new timestep
!------------------------------
varstep = maxval(abs(density)-density_old))/(0.5d0*(density)+density_old))+ density_min))

dtime = dtime/(varstep/relvar + absvar)

dtime = max(dtime,dtime_min)
icount = icount + 1
enddo
write(*,*) time,density(species_electrons)

CLOSE(CHANR)
CLOSE(CHANS)
end do
end do
write(*,'(A,$)') 'PRESS ENTER TO EXIT ...'
read(*,*)

end program test



Could you please have a look and advice me how to do it?

Best,
Leila

 

[xwb]

1) Declare your constants before you use them
2) Array limits have to be declared.
3) Check that all your variables are declared (where is SPECIES_DENSITY defined)

If you wish to use optional arguments for rates, there is an example in

http://www.netlib.org/ccm/page/api/optional.html.

I'm not sure what you are doing with rates. Is it meant to be

call rates ( &
SOURCE_TERMS=species_source), &
REACTION_RATES=source_terms), &
SOURCE_TERMS_MATRIX=SourceTerm,)