Hello to all,
This is not in response to any specific question. However,
I have received several questions regarding using fortran
codes from NCL. I thought I would post several emails about
this subject. The emails will appear several days apart.
The first post uses a fortran 77 subroutines to:
(i) read data from an ascii file, (2) do some (trivial) computations
and (3) then print out the data/results.
The examples are simple but illustrate the basic concepts.
The following URL has additional details:
http://www.cgd.ucar.edu/csm/support/Document/pdf/fortran.pdf
For convenience:
(1) I will present the various steps in a C-shell script (see below).
The pedantic advantage of this approach is that
the fortran, WRAPIT and NCL are contained in one script.
(2) I will assume that users have a relatively recent
version of NCL. FYI: as of 10 April 2002 the release is
4.2.0.a023
When might a user want to use fortran from an NCL script?
NCL (like IDL, Matlab, ferret, yorick) is an interpreted
language. As such, use of 'do loops' should be minimized
because interpreted languages do not optimize loops. Thus,
when multi-level do loops are required and time is
important using fortran is the best option.
Reading 'complicated' ascii files may be best done via
fortran formatted read statements.
Formatting text via NCL is possible but can be tedious.
This approach would involve using NCL's built-in
functions "sprinti" and "sprintf". However, creating nice
tabular output via fortran is often easier.
Sometimes you may have an existing fortran subroutine that
performs a calculation and you do not want to rewrite it in NCL.
etc
-----------------
Example 1: Reading an ascii file that has imbedded characters.
Consider the following ascii file:
YY/MM/DD HH:MM HH.MM Kd Kn Max Ratio
01/08/12 00:00 0.00 547.5 439.3 528.9 1.04
01/08/12 00:05 0.08 414.2 321.6 517.5 0.80
01/08/12 00:10 0.17 562.6 445.5 494.5 1.14
01/08/12 00:15 0.25 426.3 338.7 482.7 0.88
01/08/12 00:20 0.33 486.2 382.8 470.9 1.03
01/08/12 00:25 0.42 480.9 383.2 459.0 1.05
01/08/12 00:30 0.50 339.6 262.6 434.8 0.78
01/08/12 00:35 0.58 382.8 296.0 422.6 0.91
01/08/12 00:40 0.67 360.7 284.4 410.3 0.88
01/08/12 00:45 0.75 308.9 237.7 397.9 0.78
01/08/12 00:50 0.83 407.4 314.5 372.8 1.09
01/08/12 00:55 0.92 357.3 275.9 360.2 0.99
01/08/12 01:00 1.00 249.2 187.5 347.5 0.72
01/08/12 01:05 1.08 242.9 182.8 334.7 0.73
Yes, you can read this via NCL but for illustration we
will read it via fortran.
The following comments are associated with the script below:
(1) Note the delimiters 'C NCLFORTSTART' and 'C NCLEND'
in the beginning of the subroutines. These
enclose all declarations pertaining to the arguments
being passed. They should appear only in subroutines
being called from NCL. To the fortran compiler,
they are merely comment lines.
(2) After 'C NCLEND' the code may include other declarations
(eg: data statement, common blocks, additional local memory, etc)
(3) You may have any number of fortran subroutines.
(4) WRAPIT is a script which is distributed with more recent versions
of NCL. It performs the various steps necessary
to create the shared object (.so). In addition, WRAPIT
removes any temporary files. Incidentally, the suffix
'.so' is merely a common convention.
By default the shared object prefix will be the same prefix
as the 1st file (here "sub"). You can use the -n option
to specify a different name.
WRAPIT runs on a number of different systems.
Mark Stevens (CGD) provides many examples within
the WRAPIT function. Feel free to look at it.
It is located in $NCARG_ROOT/bin/WRAPIT
(5) Use "new" to preallocate memory for the
variables being returned from the fortran subroutine.
The variable type (eg: float, double) should agree between
NCL and fortran: integer-to-integer, float-to-real,
double-to-double precision.
(6) the external line in the NCL script provides an *arbitrary*
identifier to the shared object. By convention it is capitalized
but this is not required. This identifier can have any name.
#!/usr/bin/csh
# ===========================FORTRAN===================================
cat >! sub.f << "END_FORTRAN"
C NCLFORTSTART
subroutine rddata(filnam,nn ,yy ,mo ,dd ,hh ,mm
+ ,hhmm, kd, kn, zmax, rat)
implicit none
character*(*) filnam
integer nn
integer yy(nn), mo(nn), dd(nn),hh(nn), mm(nn)
real hhmm(nn), kd(nn), kn(nn), zmax(nn), rat(nn)
C NCLEND
integer k
character*80 skip
open (unit=10,file=filnam,form="formatted",access="sequential")
c read/ignore the first record. This is just ascii header stuff
read (10,"(a)") skip
c read the desired variables
do k=1,nn
read(10,"(1x,5(i2,1x),f5.2,f7.2,3(1x,f7.2) )")
+ yy(k), mo(k), dd(k), hh(k), mm(k)
+ ,hhmm(k), kd(k), kn(k), zmax(k), rat(k)
end do
return
end
C NCLFORTSTART
subroutine compute(nn,kd, zmax, kd2, ratio)
implicit none
integer nn
real kd(nn), zmax(nn)
real kd2(nn), ratio(nn)
C NCLEND
integer k
C trivial computations
do k=1,nn
kd2(k) = 2.*kd(k)
ratio(k) = zmax(k)/kd2(k)
end do
return
end
C NCLFORTSTART
subroutine prdata(head ,nn ,yy ,mo ,dd ,hh ,mm
+ ,hhmm, kd, kn, zmax, rat)
implicit none
character*(*) head
integer nn
integer yy(nn), mo(nn), dd(nn),hh(nn), mm(nn)
real hhmm(nn), kd(nn), kn(nn), zmax(nn), rat(nn)
C NCLEND
integer k
c space down 3 lines then print out header
write (*,"(///,10x,a)") head
c print column headers
write(*,"(3x,'YY/MM/DD',1x,'HH:MM',1x,'hh.mm',5x,'Kd*2'
+ 3x,'sqrt(Kn)',5x,'Max',2x,'NewRatio')")
c write the desired variables
do k=1,nn
write(*,"(3x,5(i2,1x),f5.2,2x,f7.2,2x,f8.2,2(2x,f7.2) )")
+ yy(k), mo(k), dd(k), hh(k), mm(k)
+ ,hhmm(k), kd(k), kn(k), zmax(k), rat(k)
end do
return
end
"END_FORTRAN"
# =========================== WRAPIT===================================
# $NCARG_ROOT/bin/WRAPIT sub.f # create sub.so
WRAPIT sub.f
# ==============================NCL====================================
cat >! main.ncl << "END_NCL"
; commented out cuz not used here
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
external SUB "./sub.so" ; fortran shared object produced by WRAPIT
begin
diri = "./"
fili = "ascii.data"
fname= diri+fili
; access the file, determine # records
a = asciiread(fname, -1, "string")
nRow = dimsizes(a) ; # rows (records)
nrow = nRow-1 ; remember: skip the 1st rec
delete(a)
; preallocate to hold data
; returned from fortran
yy = new (nrow, "integer")
mo = new (nrow, "integer")
dd = new (nrow, "integer")
hh = new (nrow, "integer")
mm = new (nrow, "integer")
hhmm = new (nrow, "float")
Kd = new (nrow, "float")
Kn = new (nrow, "float")
zmax = new (nrow, "float")
rat = new (nrow, "float")
SUB::rddata(fname, nrow ,yy ,mo ,dd ,hh ,mm \
,hhmm, Kd, Kn, zmax, rat)
print ("======> NCL print <=====")
print (yy+" "+mo+" "+dd+" "+hh+" "+mm+" " \
+hhmm+" "+Kd+" "+Kn+" "+zmax+" "+rat)
Kn = sqrt(Kn) ; overwrite Kn with sqrt (array syntax)
Kd2 = new (dimsizes(Kd), typeof(Kd) ) ; Kd2 = 2.*Kd
RAT = new (dimsizes(Kd), typeof(Kd) ) ; RAT = zmax/Kd2
SUB::compute(nrow,Kd,zmax,Kd2,RAT)
title = "+++ Fortran print of Bogus New Data +++"
SUB::prdata(title, nrow ,yy ,mo ,dd ,hh ,mm \
,hhmm, Kd2, Kn, zmax, RAT)
end
"END_NCL"
# ============================EXECUTE==================================
ncl < main.ncl
; ncl < main.ncl >&! z.out # send ncl output to a specific file
# ====================CLEAN UP TEMPORARY FILES=========================
/usr/bin/rm sub.f
/usr/bin/rm sub.so
/usr/bin/rm main.ncl
exit
++++++++++++++++++++++++++++++++++++++++++++=
Here is the output produced by the above script.
++++++++++++++++++++++++++++++++++++++++++++=
See http://ngwww.ucar.edu/ncl/ for more details.
(0) ======> NCL print <=====
(0) 1 8 12 0 0 0 547.5 439.3 528.9 1.04
(1) 1 8 12 0 5 0.08 414.2 321.6 517.5 0.8
(2) 1 8 12 0 10 0.17 562.6 445.5 494.5 1.14
(3) 1 8 12 0 15 0.25 426.3 338.7 482.7 0.88
(4) 1 8 12 0 20 0.33 486.2 382.8 470.9 1.03
(5) 1 8 12 0 25 0.42 480.9 383.2 459 1.05
(6) 1 8 12 0 30 0.5 339.6 262.6 434.8 0.78
(7) 1 8 12 0 35 0.58 382.8 296 422.6 0.91
(8) 1 8 12 0 40 0.67 360.7 284.4 410.3 0.88
(9) 1 8 12 0 45 0.75 308.9 237.7 397.9 0.78
(10) 1 8 12 0 50 0.83 407.4 314.5 372.8 1.09
(11) 1 8 12 0 55 0.92 357.3 275.9 360.2 0.99
(12) 1 8 12 1 0 1 249.2 187.5 347.5 0.72
(13) 1 8 12 1 5 1.08 242.9 182.8 334.7 0.73
(0) ========
+++ Fortran print of Bogus New Data +++
YY/MM/DD HH:MM hh.mm Kd*2 sqrt(Kn) Max NewRatio
1 8 12 0 0 0.00 1095.00 439.30 528.90 0.48
1 8 12 0 5 0.08 828.40 321.60 517.50 0.62
1 8 12 0 10 0.17 1125.20 445.50 494.50 0.44
1 8 12 0 20 0.33 972.40 382.80 470.90 0.48
1 8 12 0 25 0.42 961.80 383.20 459.00 0.48
1 8 12 0 30 0.50 679.20 262.60 434.80 0.64
1 8 12 0 35 0.58 765.60 296.00 422.60 0.55
1 8 12 0 40 0.67 721.40 284.40 410.30 0.57
1 8 12 0 45 0.75 617.80 237.70 397.90 0.64
1 8 12 0 50 0.83 814.80 314.50 372.80 0.46
1 8 12 0 55 0.92 714.60 275.90 360.20 0.50
1 8 12 1 0 1.00 498.40 187.50 347.50 0.70
1 8 12 1 5 1.08 485.80 182.80 334.70 0.69
This archive was generated by hypermail 2b29 : Tue Apr 16 2002 - 14:06:53 MDT