format output from successive do loops to write in separate columns

[DustBunnie]

I am new to this forum and need help with formatting output. I have written a program in Fortran 90 that calculates wind erosion. It uses 24 hourly wind data sets and computes the dust emission for each model grid cell for each hour. Currently the output files write these data in one file in a single column. I would like to change it so each hourly iteration of the program writes to a column and therefore I would have 24 columns of data in each the output file. In previous applications of this model I have had a low enough number of grid cells that I could import the output files into excel and use a macro to reformat the file. Now I have 257664 grid cells so I need to reformat the model's output file.
Example code:

! --------- temporal and spatial wind fields ustar computation ------------------------
!* This subroutine is the engine of the program, it controls the cycle through input days and hours *
!* it also opens all the external files for data output, calculates the ustar from the input wind *
!* files and then calls the subroutine that calculate the horizontal flux *
!* it USE the arrays module which defines the dimensions of the the arrays shared between program *
!* units, the arrays module also includes the erodpara.inc file that defines all the parameters used*
!* by program units. *
!****************************************************************************************************
subroutine beginBC

use arraysBC ! module containing arrays and include file erodparaBC.inc
implicit none
integer:: ij, ih, lat, long, i, j
real:: van, ustarns, maxwind, minsoilthres, lnval, tempz0, v_z0
! open all external output files
open (70, file='vflux.dat',status='unknown', form='formatted')
open (75, file='F_tot.dat',status='unknown', form='formatted')
do 10 ij=1,nday ! do loop based on days modelled, defined parameter in erodpara.inc
write (6,'(a5,i3)') 'day=',ij ! written to screen so user can know programs status
do 20 ih=1,nhour ! do loop based on hours modelled, defined parameter in erodpara.inc
write (6,'(a5,i3)') 'hr=',ih ! written to screen so user can know programs status
write(70,'(a5,1X,i3,1X,a5,1X,i3)')'day=',ij, 'hour=',ih
write(75,'(a5,1X,i3,1X,a5,1X,i3)')'day=',ij, 'hour=',ih
do 30 lat=1,nlat !these are rows
do 40 long=1,nlong !these are columns
van=wind(long,lat,ij,ih)
do 50 i=1,ntyp !this is the number of soiltypes&vegtypes in a grid cell


tempz0=(z/v_z0) ! z0 fixed for this study at met tower derivied 0.0002m; z = 10m
lnval = log(tempz0) ! this compiler not recognizing LN function
Ustarns=(VK*van)/ lnval

Ustar (long,lat,i)=Ustarns ! fill array of Ustar values based on 10m wind field generated by Calmet

50 continue
40 continue
30 continue
!filter used to check if the winds strong enough for erosion on bare soil, if not move to next time step
maxwind=Maxval(ustar)

if (maxwind .lt. umin)then
!fill the output files with 0.0s

vert_flux=0.0
write(70,'(E13.3)') (vert_flux)! try to write to output file of vertical flux by hour
tot_flux=0.0
write(75,'(E13.3)') (tot_flux)
percent=0

UthR=0.0


Goto 20 ! go to next hour time step

else

call Calc_UthBC !


endif

20 continue
10 continue

end subroutine beginBC



Subroutine Calc_UthBC

use arraysBC
implicit none
real(8), Dimension (numgas):: soilarray
integer:: lat, long, i, r, j, k, count
real:: UthRs, qwhites, cal_q, ustarpass, rep_uths, tmp_percent, rcount, rnumgas,v_flux
real:: cal_vfluxsol2, cal_vfluxsol3, cal_vfluxsol4, cal_vfluxsol5


do 30 lat=1,nlat !rows
do 40 long=1,nlong !columns
do 50 i=1,ntyp !soiltypes in a cell

count=0 !reset counter

if ((soiltyp(long,lat,i)).EQ.99.)then !off the grid
UthRs =-1.d0
UthR(long,lat,i)=UthRs !for array of partitioned shear velocity
percent(long,lat,i)=-1 ! for array of percent emission
v_flux=-1
vert_flux(long,lat,i)=v_flux



elseif ((soiltyp(long,lat,i)).EQ.2.)then

UthRs= min2 ! min2 is min threshold shear velocity for soiltype2;

IF (Ustar(long,lat,i).gt.UthRs) count=count+1

rcount=real(count)
rnumgas=real(numgas)
tmp_percent=(rcount/rnumgas)*100.0
percent(long,lat,i)=nint(tmp_percent) !converted to an integer value

UthRs=0.0 ! reset
soilarray=0.0 !reset
rep_uths=0.0 !reset
soilarray=soil2 !define array to be passed to subroutine
call ran_index(soilarray,rep_uths)

UthRs= rep_uths ! *(dsqrt((1-m*sigma(i)*lambda(i))*(1+m*betaR*lambda(i)))) ! rep_uths is a random value drawn from threshold shear velocity for soiltype array
IF (Ustar(long,lat,i).lt.UthRs) then !no erosion based on sampled threshold value
v_flux = 0.d0
vert_flux(long,lat,i)=v_flux
UthR(long,lat,i)=UthRs
ELSE !erosion, because wind force is greater than threshold shear velocity
!keep track of number of times this occurs for probability purposes
ustarpass=ustar(long,lat,i)

v_flux=cal_vfluxsol2(ustarpass)! for BC hydro changing cal_flux to use power functions for soil type
UthR(long,lat,i)=UthRs
vert_flux(long,lat,i)=v_flux


ENDIF
! note there are many more soil classes but I've only shown 2 in this posted example

end if

50 continue


write(80,'(i3)') (percent(long,lat,i),i=1,ntyp) !Integer number output
write(70,'(E13.3)') (vert_flux(long,lat,i),i=1,ntyp)! Scientific number output ntyp=1
write(90,'(F10.2)') (UthR(long,lat,i),i=1,ntyp) !Real number output



40 continue
30 continue

return
end subroutine Calc_UthBC

Thank you for you time

 

[xwb]

Which output file are you having problems with: 70, 75, 80 or 90?

 

[DustBunnie]

Hi xwb,
output for both 70 and 75.

 

[xwb]

Let's look at 70 first. In beginBC, it prints

write(70,'(a5,1X,i3,1X,a5,1X,i3)')'day=',ij, 'hour=',ih

for every hour of every day and

write(70,'(E13.3)') (vert_flux)! try to write to output file of vertical flux by hour

at the end.

For calcUTHBC, it prints

write(70,'(E13.3)') (vert_flux(long,lat,i),i=1,ntyp)! Scientific number output ntyp=1

for every lat and lon.

One way to reduce the amount of output is to buffer them up and prefix the format with the maximum number of possible values. Eg if ntyp is 5 and nlong is 6 then make it 30E13.3. Also move the wri(te statement after the lon loop. Also, you don't really need the labels: that is just for compatibility with Fortran IV (66). So you have something like

do lat=1,nlat
   do long=1,nlong
      do i = 1,ntyp
         ...
      end do
      ! Old write was here
   end do
   write (70, '(30E13.3)') ((vert_flux(long, lat, i), i = 1, ntyp), long = 1, nlong)
end do

For 75 it is similar: just move it out by one loop and put the number of times in front of the format

do ij = 1, nday
   do ih = 1, nhour
       ...
   end do
   write (75, '("day = ", 25I4)') ij, (ih = 1, nhour)
end do