HOMEWORK  for  FORTRAN week 2.

You will continue the top-down approach of FORTRAN code writing that you started last week.  We will continue to write code to find wind power.  The numbering of these assignments is a continuation of the numbering from Fortran week 1.   

To allow you to all start on the same footing, and to give you all a head start, I will give you the fortran code for my version 12 that you can use as your starting point this week.  

We will start the first few assignments below together in lab, with you following along at first.  Then you can finish the remaining assignments as your homework.

Assignment Intro

• If you want, you may use the copy of "wp12s.f95"  as a starting point for fortran2.   (see link under Assignment Intro, above.) 
• Copy the contents of this whole file to the clipboard, and then in emacs paste it into the "wp12s.f95" that you already have open, and save it (inside the "fortran2" directory). 
• Then, minize your emacs window and your browser window.
• Follow along with the instructor to create alias on the eidolon terminal window by typing:
   alias g12s='gfortran wp12s.f95 -o runwp12s'  
• Compile this code, by typing g12s  
• Run the code, by typing  ./runwp12s  
• When it asks for a file name, see what happens when you type a wrong name.  Then type a good file name, such as "porthardy.txt", and see what happens.
• Run the code again by typing  ./runwp12s  , and respond with filename "darwin.txt" .
  ...
  resume lecture on I/O for multiple slides, until it says "exercise 13".

• In emacs, start with the "wp12s.f95" file (or with your own working version 12) that you already have, and Save As  "wp13.f95" .
  
• Create a new alias on eidolon by typing:
   alias g13='gfortran wp13.f95 -o runwp13'  
• .
• Follow along with the instructor.  In emacs, edit version 13 as follows.  In getsounding: 
• from inside  the loop where you read past the first 5 header lines, save the first line of the sounding input as a new character variable called title.    Don't forget to declare title as a character variable of length 100 characters.  Then, just after that loop, print the title one the screen.    save, compile (g13), run (./runwp13).
• copy from the main program the array declarations for maxlines parameter, and for height zmsl and wind speed, and paste these into getsounding.   save, compile (g13), run.
• after skipping past the header lines (which you've already done), read the first 10 lines of the sounding to get height, and speed.  To do this first declare a new character string called fmt to hold the format for the read statement.  Define fmt to use the tab (T)  format to get to column 8 to read the height, to read the height as format F7.0, to then use the tab to get to column 50, and finally to read the wind speed as F7.0.  
  Next, create a counting do loop to read the first 20 lines of data values from unit 1, using the fmt that you just created.  Also, inside that loop, convert each speed from knots to m/s (where knots times 0.5144 = m/s), and  echo these values to the screen, using the skip (X) format to nicely space the columns of output.   
• Also just before that counting do loop, write to the screen some column headers with the heights and speeds (with their units).   save, compile (g13), run.
• Fix any bugs, resave,  recompile and test run until working successfully. 
• Finally, minimize all except emacs.
  ...
  resume lecture for one slide on modules, until it says "exercise 14".
  .
• 
    

• In emacs, start with the "wp13.f95" file that you already have, and Save As  "wp14.f95" .
  .
• Follow along with the instructor.  In emacs, edit version 14 as follows:
• Add header comment lines at the very beginning of your code, describing this program, and put your name  and date, and your copyright.  save, compile, run.
• In your fortran code, just after the header comment lines and just before the main program, add a module called "soundmod".  Copy the type declarations statements for filename, title, maxlines, zmsl, speed from subroutine "getsounding", and paste them into the new module .   save, compile, run.
• In both the main program and in "getsounding", just before "implicit none", add the line "use soundmod".  Also, from those two routines, remove the type declarations that are already included in the module.
• Save, compile, run,
• Fix any bugs, resave,  recompile and test run until working successfully.
  .
  .
  
• On your own, just after the soundmod module, create a new module called turbinemod that holds the turbine specs ( zhub and r), and use that module in the main program, in getturbinespecs and in findpower .
• Save, compile, run.  
• Debug as needed, and save final working result.  Minimize all except emacs.
  ....
  resume lecture for 2 slides on intrinsic functions and user-defined functions, until it says "exercise 15".

• In emacs, start with the "wp14.f95" file that you already have, and Save As  "wp15.f95" .  Update the version number and date in the header, if needed.
  .
  .
• First create a loop in findpower to give 20 equally-spaced height values between the bottom and the top of the turbine disk and display these heights to the screen so you can see if they are reasonable.  (follow along with instructor).    For those who couldn't attend class, here is the code for the updated subroutine findpower.
   
  !=======================================
  subroutine findpower
    use turbinemod            !module holding turbine specs
    implicit none             !enforce strong typing
    integer :: i              ! dummy index
    real :: delz              ! thickness (m) of each height strip across turbine
    real :: zref              ! height (m) AGL of center of strip below bottom strip
    real :: zz                !an arbitrary height (m) AGL
    real :: M                 !the function that finds the wind speed
    real :: speedz            !the wind speed (m/s) at height z (m)
   
    write(*,*)
    write(*,*) "findpower:  Calculate the wind power."
  
  !Divide the diameter of the turbine disk into 20 height strips of equal thickness 
    delz = 2.0 * r / 20.0         !vertical thickness of each height strip
    zref = zhub - r - 0.5*delz    !height (m AGL) of center of strip below bottom strip
    do i = 20,1,-1                !for each of 20 height strips across turbine, from top down
      zz = zref + real(i) * delz  !find the height (m AGL) of the center of the strip
      speedz = M(zz)              !use M function to find interpolated with speed
      write(*,*) "zz(mAGL), M(m/s) = ",zz, speedz    !display the height and wind speed
    enddo
   
  end subroutine findpower
  
  
  
• Notice that this code calls a User-defined Function called M to find the interpolated wind speed.  
• Create this user-defined function as follows, at the very end of your whole program.  This function uses the info in the soundmod module.  For those who couldn't attend class, here is the code.
  
  
  !=======================================
  real function M (zr)                !interpolate wind speed M for any height zr (m AGL)
    use soundmod                      !access the sounding data module
    implicit none                     !enforce strong typing
    real, intent(in) :: zr            !height (m) above ground level AGL
    real :: zrmsl                     !height (m) above mean sea level MSL
    integer :: i                      !dummy height index
    integer :: itop                   !index of the sounding height that is just above zr
    real :: ratio                     !relative dist of zr between bracketing sounding heights
   
    M = 0.0                           !initialize
    zrmsl = zr + zmsl(1)              !convert height zr to above mean sea level (MSL)
   
  !find array indices that bracket the desired height
    do i = 1,20                       !start at the bottom of the sounding, and work up
      itop = i                        !save the index value, in case it is for the 1st sounding ht above zrmsl
      if (zmsl(i) .ge. zrmsl) exit    !yes, found the 1st sounding height above, zrmsl, so stop looping
    enddo
  
  !do the interpolation using sounding data to find the wind speed at height zrmsl
    if (itop .gt. 1) then              !when zrmsl is bracketed by 2 sounding levels having wind speed data
      ratio = (zrmsl - zmsl(itop-1) ) / (zmsl(itop) - zmsl(itop-1) )
      M = speed(itop-1) + ratio*(speed(itop) - speed(itop-1))
    else                               !if z is below lowest sounding height
      ratio = zrmsl / zmsl(itop)       !assume z bottom = 0 ,
      M = ratio * speed(itop)          !assume M at (zbottom) = 0 m/s
    endif
     
  end function M
  
  
  !======end
  
• Save, compile, run.  Debug and resave if needed.
  
• ...
• Move on to HW 16 below.   (no PowerPoint slides at this point)

• In emacs, start with the "wp15.f95" file that you already have, and Save As  "wp16.f95" .  Update the version number and date in the header, if needed.
  .
• On your own, create a new user-defined function density that finds the air density rho (kg/m3) at any height zr AGL, based on the approximation  density = denref*exp(-zr/H), where the reference density at sea level is denref = 1.225 kg/m3, and the scale height for density is roughly H = 8550 m.
• Then access that function from   findpower for each height inside the height loop (i.e., inside the loop was given in Version 15 above), to find the value of rho at that height.
• Also, from inside that loop, add to the existing write statement to also display the air density at each height.  (don't forget to update the format for the write statement if needed.)
• Save, compile, run, fix program if needed and save
  .
• On your own, still within that loop in findpower, calculate the length  of the chord of the circle, and include the chord length  in the display to the screen.  Recall from the PowerPoint slides from week 6 that:  length = 2*sqrt(r**2 - B**2), where B = zhub - zz , where these two heights are both AGL..   Save, compile, run, fix program if needed and save.
  .
• On your own, still within that loop in findpower for, calculate the amount of wind power associated with that one height strip.  Name that wind power variable delpower, and also display it.  (Recall from the PowerPoint slides from Fortran week 6,  wind power P = 0.5 * rho * area * (speedz**3) ).  The portion of the wind-turbine area associated with that one height strip is  delz * length  .
• Save, compile, run, fix program if needed and save.
  ...
• Move on to HW 17 below.   (no powerpoint slides at this point)

• In emacs, start with the "wp16.f95" file that you already have, and Save As  "wp17.f95" .  Update the version number and date in the header, if needed.
  .
  On Your Own:
  
• In the main program remove the statement that said power = 0.0, and remove the line that wrote that power to the screen.  (we will let the subroutines calcualte and display the power from now on.  
• Inside find power, declare a new variable for total power, and initialize it to zero.  Then, inside the existing counting do loop,  add each delpower value to the variable power,  to gradually accumulate the total power as we add the contribution from each horizontal strip of wind turbine disk.   
• After that end of that loop, Convert the power from W to kW.  Display the total power in kW on the screen.
• Save, compile, run, fix program if needed and save.  Test for both darwin.txt and porthardy.txt
  .
• Return to 3 powerpoint slides on string manipulation and file handling.

• In emacs, start with the "wp17.f95" file that you already have, and Save As  "wp18.f95" .  Update the version number and date in the header, if needed.
  .
•  Follow along with the instructor (or see my code below).   Modify subroutine saveresults as follows: use the file name of the input file, and use substring methods to extract the root part of it (namely, don't include the suffix), and display that root file name to the screen.   Save, compile, run, fix program if needed and save.  
  
• Then use concatenation to add the characters "out.txt" to the root name to create a new output file name.  Display that output file name to the screen.   Save, compile, run, fix program if needed and save.
• Open a new file with that output file name, and write the location title line from the input file back into this output file.  Also write a line for the future column headers for height, wind speed, chord length, density, and delpower.  Then close the file.   
  (Because we didn't get to this in class, I give a copy of this code below:)
  
  !=======================================
  subroutine saveresults
    use soundmod                      !access data from sounding module
    use turbinemod                    !access data from turbine specs module
    implicit none                     !enforce strong typing
    character (len=70) :: outfile     !output file name
    character (len=50) :: rootfile    !file name without suffixes
    character (len=1) :: tab=achar(9) !ascii tab character
    integer :: strlength              !character length of rootfile name
    integer :: ero                    !error flag for opening a file
   
    write(*,*)
    write(*,*) "saveresults:  Write to disk and screen the wind power potential"
    write(*,*)
   
  !Create an output file name from the root of the input filename.  The root is just 
  !the original input file name without the last 4 characters ".txt"
    strlength = len_trim(filename) - 4     !the number of characters in the root name
    rootfile = filename(1:strlength)       !extract substring
    write(*,*) "  rootfile = ", rootfile   !debugging output
    write(*,*)
   
    outfile = trim(rootfile) // "out.txt"        !concatenate suffix to create output file name
    write(*,*) "  Results are saved in file: ", outfile     !debugging output
   
  !Open a file to hold output results
    open(1, file=outfile, status="replace", action="write", iostat=ero)    !open file holding the sounding
  
  !Save debugging output
    write(1,*) title                    !write the location of the sounding
    write(1,*)
    write(1,*) "Hub height agl(m)=", zhub, "    Turbine radius (m)=", r
    write(1,*)
  
       !part (a) of your code for version 19 will go here. <===========
  
    write(1,*)
    write(1,*) "  zz(m)",tab,"M(m/s)",tab,"L(m)",tab,"rho(kg/m3)",tab,"delpower(W) = "  !write column headers
  
      !parts (b) and (c) of your code for version 19 will go here. <===========
  
    write(1,*)
    close(1)
   
    write(*,*)
    write(*,*) "Finished Wind Power Program.  Bye."
    write(*,*)
  end subroutine saveresults
  
  
  
  
• Save, compile, run, fix program if needed and save.   To see if it worked, use your linux terminal window to list the files (ls).  Hopefully this new output file is there.  Then use cat to view the contents of that file on the linux screen.  Debug, fix, and recompile if needed.   
  .
  
  

• In emacs, start with the "wp18.f95" file that you already have, and Save As  "wp19.f95" .  Update the version number and date in the header, if needed.
  .
• Modify saveresults to also write to the output file the following:
  (a)  the sounding data (z and speed) for the first sounding 10 heights, identical to what you displayed  on the screen from the getsounding subroutine.  save, compile, run.  Then in linux use cat to display the results of the output file that you just created.  If the output file doesn't contain what it needs to , then fix the code and save, recompile, run, and then again check the output file.
  .
  (b) the table from the findpower subroutine that gives the height, wind speed, chord length, air density, and delpower for each of the 20 height strips across the wind turbine, idential to what you displayed on the screen from the findpower subroutine.   [You will need to figure out how to get these values (which you already calculated in findpower) from findpower into saveresults. ]   Save, compile, run, debug, and check the output file, and debug if needed.
  .
  (c) also write the total power in kilowatts to this file.   Save, compile, run, debug, and check the output file, and debug if needed.
  .
• Your program is finished.  Be sure to test it for both Port Hardy and Darwin.

Marking:

The bottom line is that we will look for the following files in the "fortran2" directory under your user directory on eidolon:

1. Seven ascii text files holding your fortran code versions:
• wp13.f95
• wp14.f95
• wp15.f95
• wp16.f95
• wp17.f95
• wp18.f95
• wp19.f95
2. ... and Seven executable files:
• runwp13
• runwp14
• runwp15
• runwp16
• runwp17
• runwp18
• runwp19
3. ... and the two output files:
• porthardyout.txt
• darwinout.txt

The marker will grade your work by looking at EACH of your fortran versions, and compiling them himself to be sure they run.  Also, marks will be deducted if your code is not adequately documented with in-line comments.

The above files are due (i.e., must be inside your /fortran2 directory on eidolon) by the deadline discussed in class.

Hints:

• here are some tips on how to read compiler error messages .
• look at the sample large FORTRAN program.  
• refer to the the tutorials listed below (Reading Assignment) .
• check the full, online fortran users reference manual listed in the Resources link of the course web page.  Although this is for a different compiler than the one we are using, its FORTRAN 95 description should be close enough.
• Check the Fortran code segments the PDF copies of the lectures for both weeks.

Reading Assignment  

• Skim chapters 6 & 8 of the following FORTRAN tutorials.
  (courtesy of Stephen Brooks, Univ. of St. Andrews, Scotland. http://www-solar.mcs.st-and.ac.uk/%7Esteveb/course/course.html )
  2. Basics: program structure, data types, programming style, conditional statements, loops, arrays .
  3. Arithmetic expressions, intrinsic functions, simple input and output (I/O) .
  4. Characters, strings, print .
  5. Internal subroutines and functions .
  6. Modules . 
  7. More about arrays, dynamic allocation of arrays, key words . 
  8. More about I/O . 
  9. Case statement .

• Read Chapter 5 on Debugging, from Kernighan and Pike "The Practice of Programming" in our course material pack.

-end-