fortran - 悲催的科学匠人 - 冷水's blog
[zz]FORTRAN版非递归快速排序
在二叉树离散网格单元时，最初使用递归的，但是网格单元到百万量级时，排序程序就出错了。怀疑是递归耗尽系统堆栈的问题，搜到一个非递归版本： http://www.fortran.com/fortran/quick_sort2.f 百万量级也不是问题了。
PROGRAM test_nrqsort
REAL,ALLOCATABLE :: array(:)
INTEGER,ALLOCATABLE :: idx(:)
INTEGER N,i
INTEGER,PARAMETER :: seed = 86456
!READ(*,*) N
N = 5000000
ALLOCATE(array(N),idx(N))
CALL srand(seed)
DO i=1,N;
array(i) = rand()
ENDDO
WRITE(*,*) 'sorting'
CALL SORTRX(N,array,idx)
!DO i=1,N
! WRITE(*,*) array(idx(i))
!ENDDO
DO i=2,N
IF(array(idx(i)) .lt. array(idx(i-1)) ) THEN
WRITE(*,*) 'error', array(idx(i)), array(idx(i-1))
ENDIF
ENDDO
DEALLOCATE(array,idx)
END
C From Leonard J. Moss of SLAC:
C Here's a hybrid QuickSort I wrote a number of years ago. It's
C based on suggestions in Knuth, Volume 3, and performs much better
C than a pure QuickSort on short or partially ordered input arrays.
SUBROUTINE SORTRX(N,DATA,INDEX)
C===================================================================
C
C SORTRX -- SORT, Real input, indeX output
C
C
C Input: N INTEGER
C DATA REAL
C
C Output: INDEX INTEGER (DIMENSION N)
C
C This routine performs an in-memory sort of the first N elements of
C array DATA, returning into array INDEX the indices of elements of
C DATA arranged in ascending order. Thus,
C
C DATA(INDEX(1)) will be the smallest number in array DATA;
C DATA(INDEX(N)) will be the largest number in DATA.
C
C The original data is not physically rearranged. The original order
C of equal input values is not necessarily preserved.
C
C===================================================================
C
C SORTRX uses a hybrid QuickSort algorithm, based on several
C suggestions in Knuth, Volume 3, Section 5.2.2. In particular, the
C "pivot key" [my term] for dividing each subsequence is chosen to be
C the median of the first, last, and middle values of the subsequence;
C and the QuickSort is cut off when a subsequence has 9 or fewer
C elements, and a straight insertion sort of the entire array is done
C at the end. The result is comparable to a pure insertion sort for
C very short arrays, and very fast for very large arrays (of order 12
C micro-sec/element on the 3081K for arrays of 10K elements). It is
C also not subject to the poor performance of the pure QuickSort on
C partially ordered data.
C
C Created: 15 Jul 1986 Len Moss
C
C===================================================================
INTEGER N,INDEX(N)
REAL DATA(N)
INTEGER LSTK(31),RSTK(31),ISTK
INTEGER L,R,I,J,P,INDEXP,INDEXT
REAL DATAP
C QuickSort Cutoff
C
C Quit QuickSort-ing when a subsequence contains M or fewer
C elements and finish off at end with straight insertion sort.
C According to Knuth, V.3, the optimum value of M is around 9.
INTEGER M
PARAMETER (M=9)
C===================================================================
C
C Make initial guess for INDEX
DO 50 I=1,N
INDEX(I)=I
50 CONTINUE
C If array is short, skip QuickSort and go directly to
C the straight insertion sort.
IF (N.LE.M) GOTO 900
C===================================================================
C
C QuickSort
C
C The "Qn:"s correspond roughly to steps in Algorithm Q,
C Knuth, V.3, PP.116-117, modified to select the median
C of the first, last, and middle elements as the "pivot
C key" (in Knuth's notation, "K"). Also modified to leave
C data in place and produce an INDEX array. To simplify
C comments, let DATA[I]=DATA(INDEX(I)).
C Q1: Initialize
ISTK=0
L=1
R=N
200 CONTINUE
C Q2: Sort the subsequence DATA[L]..DATA[R].
C
C At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L,
C r > R, and L <= m <= R. (First time through, there is no
C DATA for l < L or r > R.)
I=L
J=R
C Q2.5: Select pivot key
C
C Let the pivot, P, be the midpoint of this subsequence,
C P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R)
C so the corresponding DATA values are in increasing order.
C The pivot key, DATAP, is then DATA[P].
P=(L+R)/2
INDEXP=INDEX(P)
DATAP=DATA(INDEXP)
IF (DATA(INDEX(L)) .GT. DATAP) THEN
INDEX(P)=INDEX(L)
INDEX(L)=INDEXP
INDEXP=INDEX(P)
DATAP=DATA(INDEXP)
ENDIF
IF (DATAP .GT. DATA(INDEX(R))) THEN
IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN
INDEX(P)=INDEX(L)
INDEX(L)=INDEX(R)
ELSE
INDEX(P)=INDEX(R)
ENDIF
INDEX(R)=INDEXP
INDEXP=INDEX(P)
DATAP=DATA(INDEXP)
ENDIF
C Now we swap values between the right and left sides and/or
C move DATAP until all smaller values are on the left and all
C larger values are on the right. Neither the left or right
C side will be internally ordered yet; however, DATAP will be
C in its final position.
300 CONTINUE
C Q3: Search for datum on left >= DATAP
C
C At this point, DATA[L] <= DATAP. We can therefore start scanning
C up from L, looking for a value >= DATAP (this scan is guaranteed
C to terminate since we initially placed DATAP near the middle of
C the subsequence).
I=I+1
IF (DATA(INDEX(I)).LT.DATAP) GOTO 300
400 CONTINUE
C Q4: Search for datum on right <= DATAP
C
C At this point, DATA[R] >= DATAP. We can therefore start scanning
C down from R, looking for a value <= DATAP (this scan is guaranteed
C to terminate since we initially placed DATAP near the middle of
C the subsequence).
J=J-1
IF (DATA(INDEX(J)).GT.DATAP) GOTO 400
C Q5: Have the two scans collided?
IF (I.LT.J) THEN
C Q6: No, interchange DATA[I] <--> DATA[J] and continue
INDEXT=INDEX(I)
INDEX(I)=INDEX(J)
INDEX(J)=INDEXT
GOTO 300
ELSE
C Q7: Yes, select next subsequence to sort
C
C At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r],
C for all L <= l < I and J < r <= R. If both subsequences are
C more than M elements long, push the longer one on the stack and
C go back to QuickSort the shorter; if only one is more than M
C elements long, go back and QuickSort it; otherwise, pop a
C subsequence off the stack and QuickSort it.
IF (R-J .GE. I-L .AND. I-L .GT. M) THEN
ISTK=ISTK+1
LSTK(ISTK)=J+1
RSTK(ISTK)=R
R=I-1
ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN
ISTK=ISTK+1
LSTK(ISTK)=L
RSTK(ISTK)=I-1
L=J+1
ELSE IF (R-J .GT. M) THEN
L=J+1
ELSE IF (I-L .GT. M) THEN
R=I-1
ELSE
C Q8: Pop the stack, or terminate QuickSort if empty
IF (ISTK.LT.1) GOTO 900
L=LSTK(ISTK)
R=RSTK(ISTK)
ISTK=ISTK-1
ENDIF
GOTO 200
ENDIF
900 CONTINUE
C===================================================================
C
C Q9: Straight Insertion sort
DO 950 I=2,N
IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN
INDEXP=INDEX(I)
DATAP=DATA(INDEXP)
P=I-1
920 CONTINUE
INDEX(P+1) = INDEX(P)
P=P-1
IF (P.GT.0) THEN
IF (DATA(INDEX(P)).GT.DATAP) GOTO 920
ENDIF
INDEX(P+1) = INDEXP
ENDIF
950 CONTINUE
C===================================================================
C
C All done
END
PLOT3D文件格式[ZZ]
http://hi.baidu.com/necrohan/blog/item/4997711f294102fde0fe0b8e.html
其实早已熟悉，只是借来存档一下。此外，以前不知iblank的unformatted格式写法，这个博文居然有。
PLOT3D文件格式
FLOW3D可以导出的数据格式有限，说明中说支持tecplot，但是是需要tecplot插件的，而且要360版本才行，找了2个插件可是不能用，而且插件的功能好像很有限。
PLOT3D格式的数据可以在后处理时导出，为此查了一下PLOT3D格式的数据，实在不行再转为tecplot格式的数据。
PLOT3D 数据格式
PLOT3D数据格式源于NASA，广泛用于规则网格的CFD数据文件。PLOT3D文件可以是ASCII, 也可是Fortran unformatted 或 C binary形式。
PLOT3D文件分为网格文件(XYZ 文件), 空气动力学结果文件 (Q 文件)和通用结果文件(函数文件 + 函数名称文件)。网格文件中可加入所谓的IBlank参数。
IBlank参数的定义：
IBlank是在每一网格点上给出的一个正数值，定义如下：
0 - 计算域之外，非流体点
1 - 正常点
2 - 固面边界点
负数 - 分块网格界面点，其数值为相邻网格块的块号
以下为各文件使用FORTRAN读入的语句，所有文件均为无格式文件。
网格文件（XYZ文件）：
XYZ 文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) (((X(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Y(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Z(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX)
XYZ 文件, 单块（single-block）, 加 IBlank:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) (((X(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Y(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Z(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((IBLANK(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX)
XYZ 文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) ((X(I,J),I=1,IMAX),J=1,JMAX),&
((Y(I,J),I=1,IMAX),J=1,JMAX)
XYZ 文件, 多块（multi-block）
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) (((X(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Y(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Z(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N))
ENDDO
XYZ 文件, 多块（multi-block）, 加 IBlank:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) (((X(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Y(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Z(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((IBLANK(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N))
ENDDO
XYZ 文件, 二维, 多块（multi-block）
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) ((X(I,J),I=1,IMAX(N)),J=1,JMAX(N)),&
((Y(I,J),I=1,IMAX(N)),J=1,JMAX(N))
ENDDO
空气动力学结果文件 (Q 文件)(Q 文件定义过窄，现已很少使用。):
Q 文件专为外流空气动力学设计，对三维流动，数组变量如下：
Q1 - 无量纲 密度
Q2 - 无量纲 X-动量
Q3 - 无量纲 Y-动量
Q4 - 无量纲 Z-动量
Q5 - 无量纲 总能
另加4个参数：
FSMACH - 自由流马赫数
ALPHA - 攻角
RE - 雷诺数
TIME - 时间
Q 文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT) ((((Q(I,J,K,M),I=1,IMAX),J=1,JMAX),K=1,KMAX),M=1,5)
Q 文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT) (((Q(I,J,M),I=1,IMAX),J=1,JMAX),M=1,4)
Q 文件, 多块（multi-block）:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT)((((Q(I,J,K,M),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),M=1,5)
ENDDO
Q 文件定义过窄，现已很少使用。大多CFD工作者使用通用结果文件，即
所谓的函数文件（function file）。函数文件中可定义任意数量的数组变量，
其名称在另一函数名称文件（function name file）中给出。
函数文件:
函数文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) ((((F(I,J,K,M),I=1,IMAX),J=1,JMAX),K=1,KMAX),M=1,NFUN)
函数文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) (((F(I,J,M),I=1,IMAX),J=1,JMAX),M=1,NFUN)
函数文件, 多块（multi-block）:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) ((((F(I,J,K,M),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),M=1,NFUN)
ENDDO
函数名称文件:
函数名称文件是一ASCII文件，列有函数文件中数组变量的对应名称。以下为一例子：
density
pressure
u;velocity vector
v
w
temperature
turbulence energy
...
...
注意事项：
1 名称的数量和排列次序与函数文件中数组变量相同。
2 分割号 ";" 表示向量行的开始。";"之右为向量名，";"之左为X分量。其下两行为Y分量和Z分量。
http://www.grc.nasa.gov/WWW/wind/valid/plot3d.html
下面例子假设网格文件句柄为7，结果文件句柄为8。
2D, Whole, Formatted, Single-Block Grid and Solution
parameter ( imax = 100 )
parameter ( jmax = 100 )
integer i
integer j
integer m
integer n
integer ni
integer nj
real mach ! freestream Mach number
real alpha ! freestream angle-of-attack
real reyn ! freestream Reynolds number
real time ! time
real x(imax,jmax)
real y(imax,jmax)
real q(imax,jmax,4)
open ( unit=7, form='formatted', file='2D.x' )
open ( unit=8, form='formatted', file='2D.q' )
read(7,*) ni, nj
read(7,*)
& (( x(i,j), i=1,ni), j=1,nj),
& (( y(i,j), i=1,ni), j=1,nj)
read(8,*) ni, nj
read(8,*) mach, alpha, reyn, time
read(8,*) ((( q(i,j,n), i=1,ni), j=1,nj), n=1,4)
3D, Whole, Unformatted, Multi-Block Grid and Solution
parameter ( imax = 100 )
parameter ( jmax = 100 )
parameter ( kmax = 100 )
parameter ( nbmax = 10 )
integer i
integer j
integer m
integer n
integer nblocks
integer ni (nbmax)
integer nj (nbmax)
integer nk (nbmax)
real mach ! freestream Mach number
real alpha ! freestream angle-of-attack
real reyn ! freestream Reynolds number
real time ! time
real x(imax,jmax,kmax,nbmax)
real y(imax,jmax,kmax,nbmax)
real z(imax,jmax,kmax,nbmax)
real q(imax,jmax,kmax,nbmax,5)
open ( unit=7, form='unformatted', file='3D.x' )
open ( unit=8, form='unformatted', file='3D.q' )
read(7) nblocks
read(7) ( ni(m), nj(m), nk(m), m = 1, nblocks )
do m = 1, nblocks
read(7)
& ((( x(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)),
& ((( y(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)),
& ((( z(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m))
enddo
read(8) nblocks
read(8) ( ni(m), nj(m), nk(m), m = 1, nblocks )
do m = 1, nblocks
read(8) mach, alpha, reyn, time
read(8)
& (((( q(i,j,k,m,n), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)), n=1,5)
enddo
TECPLOT 的 PREPLOT 可以把 PLOT3D 文件转为 TECPLOT 格式。（没有测试）
20090526 add
由于后来发现可以用flow3d自带的命令行方式后处理提取数据，所以这个方法放弃。
PLOT3D格式的数据可以在后处理时导出，为此查了一下PLOT3D格式的数据，实在不行再转为tecplot格式的数据。
PLOT3D 数据格式
PLOT3D数据格式源于NASA，广泛用于规则网格的CFD数据文件。PLOT3D文件可以是ASCII, 也可是Fortran unformatted 或 C binary形式。
PLOT3D文件分为网格文件(XYZ 文件), 空气动力学结果文件 (Q 文件)和通用结果文件(函数文件 + 函数名称文件)。网格文件中可加入所谓的IBlank参数。
IBlank参数的定义：
IBlank是在每一网格点上给出的一个正数值，定义如下：
0 - 计算域之外，非流体点
1 - 正常点
2 - 固面边界点
负数 - 分块网格界面点，其数值为相邻网格块的块号
以下为各文件使用FORTRAN读入的语句，所有文件均为无格式文件。
网格文件（XYZ文件）：
XYZ 文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) (((X(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Y(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Z(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX)
XYZ 文件, 单块（single-block）, 加 IBlank:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) (((X(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Y(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((Z(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX),&
(((IBLANK(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX)
XYZ 文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) ((X(I,J),I=1,IMAX),J=1,JMAX),&
((Y(I,J),I=1,IMAX),J=1,JMAX)
XYZ 文件, 多块（multi-block）
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) (((X(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Y(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Z(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N))
ENDDO
XYZ 文件, 多块（multi-block）, 加 IBlank:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) (((X(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Y(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((Z(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),&
(((IBLANK(I,J,K),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N))
ENDDO
XYZ 文件, 二维, 多块（multi-block）
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) ((X(I,J),I=1,IMAX(N)),J=1,JMAX(N)),&
((Y(I,J),I=1,IMAX(N)),J=1,JMAX(N))
ENDDO
空气动力学结果文件 (Q 文件)(Q 文件定义过窄，现已很少使用。):
Q 文件专为外流空气动力学设计，对三维流动，数组变量如下：
Q1 - 无量纲 密度
Q2 - 无量纲 X-动量
Q3 - 无量纲 Y-动量
Q4 - 无量纲 Z-动量
Q5 - 无量纲 总能
另加4个参数：
FSMACH - 自由流马赫数
ALPHA - 攻角
RE - 雷诺数
TIME - 时间
Q 文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT) ((((Q(I,J,K,M),I=1,IMAX),J=1,JMAX),K=1,KMAX),M=1,5)
Q 文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT) (((Q(I,J,M),I=1,IMAX),J=1,JMAX),M=1,4)
Q 文件, 多块（multi-block）:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) FSMACH,ALPHA,RE,TIME
READ(IUNIT)((((Q(I,J,K,M),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),M=1,5)
ENDDO
Q 文件定义过窄，现已很少使用。大多CFD工作者使用通用结果文件，即
所谓的函数文件（function file）。函数文件中可定义任意数量的数组变量，
其名称在另一函数名称文件（function name file）中给出。
函数文件:
函数文件, 单块（single-block）:
READ(IUNIT) IMAX,JMAX,KMAX
READ(IUNIT) ((((F(I,J,K,M),I=1,IMAX),J=1,JMAX),K=1,KMAX),M=1,NFUN)
函数文件, 二维, 单块（single-block）:
READ(IUNIT) IMAX,JMAX
READ(IUNIT) (((F(I,J,M),I=1,IMAX),J=1,JMAX),M=1,NFUN)
函数文件, 多块（multi-block）:
READ(IUNIT) NBLOCK
READ(IUNIT) (IMAX(N),JMAX(N),KMAX(N),N=1,NBLOCK)
DO N=1,NBLOCK
READ(IUNIT) ((((F(I,J,K,M),I=1,IMAX(N)),J=1,JMAX(N)),K=1,KMAX(N)),M=1,NFUN)
ENDDO
函数名称文件:
函数名称文件是一ASCII文件，列有函数文件中数组变量的对应名称。以下为一例子：
density
pressure
u;velocity vector
v
w
temperature
turbulence energy
...
...
注意事项：
1 名称的数量和排列次序与函数文件中数组变量相同。
2 分割号 ";" 表示向量行的开始。";"之右为向量名，";"之左为X分量。其下两行为Y分量和Z分量。
http://www.grc.nasa.gov/WWW/wind/valid/plot3d.html
下面例子假设网格文件句柄为7，结果文件句柄为8。
2D, Whole, Formatted, Single-Block Grid and Solution
parameter ( imax = 100 )
parameter ( jmax = 100 )
integer i
integer j
integer m
integer n
integer ni
integer nj
real mach ! freestream Mach number
real alpha ! freestream angle-of-attack
real reyn ! freestream Reynolds number
real time ! time
real x(imax,jmax)
real y(imax,jmax)
real q(imax,jmax,4)
open ( unit=7, form='formatted', file='2D.x' )
open ( unit=8, form='formatted', file='2D.q' )
read(7,*) ni, nj
read(7,*)
& (( x(i,j), i=1,ni), j=1,nj),
& (( y(i,j), i=1,ni), j=1,nj)
read(8,*) ni, nj
read(8,*) mach, alpha, reyn, time
read(8,*) ((( q(i,j,n), i=1,ni), j=1,nj), n=1,4)
3D, Whole, Unformatted, Multi-Block Grid and Solution
parameter ( imax = 100 )
parameter ( jmax = 100 )
parameter ( kmax = 100 )
parameter ( nbmax = 10 )
integer i
integer j
integer m
integer n
integer nblocks
integer ni (nbmax)
integer nj (nbmax)
integer nk (nbmax)
real mach ! freestream Mach number
real alpha ! freestream angle-of-attack
real reyn ! freestream Reynolds number
real time ! time
real x(imax,jmax,kmax,nbmax)
real y(imax,jmax,kmax,nbmax)
real z(imax,jmax,kmax,nbmax)
real q(imax,jmax,kmax,nbmax,5)
open ( unit=7, form='unformatted', file='3D.x' )
open ( unit=8, form='unformatted', file='3D.q' )
read(7) nblocks
read(7) ( ni(m), nj(m), nk(m), m = 1, nblocks )
do m = 1, nblocks
read(7)
& ((( x(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)),
& ((( y(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)),
& ((( z(i,j,k,m), i=1,ni(m)), j=1,nj(m)), k=1,nk(m))
enddo
read(8) nblocks
read(8) ( ni(m), nj(m), nk(m), m = 1, nblocks )
do m = 1, nblocks
read(8) mach, alpha, reyn, time
read(8)
& (((( q(i,j,k,m,n), i=1,ni(m)), j=1,nj(m)), k=1,nk(m)), n=1,5)
enddo
TECPLOT 的 PREPLOT 可以把 PLOT3D 文件转为 TECPLOT 格式。（没有测试）
20090526 add
由于后来发现可以用flow3d自带的命令行方式后处理提取数据，所以这个方法放弃。
一个简单的动态数据容器
能够动态的添加成员数据，目前只支持integer(4)和real(8)
! dynamic data type
! Qu Kun, April 2012
! 一个简单的动态数据结构容器，能够容纳INTEGER(4)和REAL(8)数据，维数范围为0到5
! 多个容器可链接成为链表
! 在链表中创建插入容器 dyd_new，指定数据类型、维数、指标范围和名称
! 获取容器指针 dyd_findnode 在链表中搜索指定名称的容器
! 释放容器内的数据空间 dyd_empty
! 开辟容器数据空间 dyd_construct
! 获取数据指针 dyd_getdata， 通过名称得到数据指针
!
! 一般使用最常用的是 dyd_new和dyd_getdata
! dyd_new的使用方法，建议采用关键字方式指定参数
! 1 变量指定数组范围，range11,range12,range21,range22,...,range51,range52
! CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
! range11=0,range12=9)
! 2 字符串指定数组范围
! CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
! 3 数组方式指定数组范围
! CALL dyd_new(head=list, atype=int4, adim=1,arange=irange, aname='X')
! 比较方便的是前2种
!
! 获取数据指针的方法为
! INTEGER(4),DIMENSION(:),POINTER :: X
! ...
! CALL getdata(list,'X', X)
! DO i=0,9; X(i) = i; ENDDO
MODULE qkdynadata
IMPLICIT NONE
INTEGER,PARAMETER,PRIVATE :: dyd_maxdim=5,dyd_namelen=32
INTEGER,PARAMETER,PRIVATE :: int4=1,real8=2
TYPE dynadata
INTEGER :: mdim=-1
INTEGER :: mrange(2,dyd_maxdim)=-1
INTEGER :: metype=-1 ! data type of the element, integer(4), real(8),real(4)
CHARACTER(LEN=dyd_namelen) mname
INTEGER(4),POINTER :: mint4d0p=>NULL()
INTEGER(4),POINTER,DIMENSION(:) :: mint4d1p=>NULL()
INTEGER(4),POINTER,DIMENSION(:,:) :: mint4d2p=>NULL()
INTEGER(4),POINTER,DIMENSION(:,:,:) :: mint4d3p=>NULL()
INTEGER(4),POINTER,DIMENSION(:,:,:,:) :: mint4d4p=>NULL()
INTEGER(4),POINTER,DIMENSION(:,:,:,:,:) :: mint4d5p=>NULL()
REAL(8),POINTER :: mreal8d0p=>NULL()
REAL(8),POINTER,DIMENSION(:) :: mreal8d1p=>NULL()
REAL(8),POINTER,DIMENSION(:,:) :: mreal8d2p=>NULL()
REAL(8),POINTER,DIMENSION(:,:,:) :: mreal8d3p=>NULL()
REAL(8),POINTER,DIMENSION(:,:,:,:) :: mreal8d4p=>NULL()
REAL(8),POINTER,DIMENSION(:,:,:,:,:) :: mreal8d5p=>NULL()
!pointer for inner dynamic data member
TYPE(dynadata),POINTER :: mdynadata
! pointer for linklist
TYPE(dynadata),POINTER :: mprev,mnext,mhead
END TYPE
INTERFACE getdata
MODULE PROCEDURE getdata_int4p0
MODULE PROCEDURE getdata_int4p1
MODULE PROCEDURE getdata_int4p2
MODULE PROCEDURE getdata_int4p3
MODULE PROCEDURE getdata_int4p4
MODULE PROCEDURE getdata_int4p5
MODULE PROCEDURE getdata_real8p0
MODULE PROCEDURE getdata_real8p1
MODULE PROCEDURE getdata_real8p2
MODULE PROCEDURE getdata_real8p3
MODULE PROCEDURE getdata_real8p4
MODULE PROCEDURE getdata_real8p5
END INTERFACE
CONTAINS
SUBROUTINE getdata_int4p0(head,aname, dp)
INTEGER(4),POINTER :: dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
IF(ASSOCIATED(iter)) dp => iter%mint4d0p
END SUBROUTINE
SUBROUTINE getdata_int4p1(head,aname,dp)
INTEGER(4),DIMENSION(:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
! WRITE(*,*) iter%mname
dp => iter%mint4d1p
END SUBROUTINE
SUBROUTINE getdata_int4p2(head,aname,dp)
INTEGER(4),DIMENSION(:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mint4d2p
END SUBROUTINE
SUBROUTINE getdata_int4p3(head,aname,dp)
INTEGER(4),DIMENSION(:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mint4d3p
END SUBROUTINE
SUBROUTINE getdata_int4p4(head,aname,dp)
INTEGER(4),DIMENSION(:,:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mint4d4p
END SUBROUTINE
SUBROUTINE getdata_int4p5(head,aname,dp)
INTEGER(4),DIMENSION(:,:,:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mint4d5p
END SUBROUTINE
SUBROUTINE getdata_real8p0(head,aname, dp)
REAL(8),POINTER :: dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d0p
END SUBROUTINE
SUBROUTINE getdata_real8p1(head,aname,dp)
REAL(8),DIMENSION(:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d1p
END SUBROUTINE
SUBROUTINE getdata_real8p2(head,aname,dp)
REAL(8),DIMENSION(:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d2p
END SUBROUTINE
SUBROUTINE getdata_real8p3(head,aname,dp)
REAL(8),DIMENSION(:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d3p
END SUBROUTINE
SUBROUTINE getdata_real8p4(head,aname,dp)
REAL(8),DIMENSION(:,:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d4p
END SUBROUTINE
SUBROUTINE getdata_real8p5(head,aname,dp)
REAL(8),DIMENSION(:,:,:,:,:),POINTER ::dp
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dp)
iter => dyd_findnode(head,aname)
dp => iter%mreal8d5p
END SUBROUTINE
SUBROUTINE dyd_empty(this)
TYPE(dynadata),INTENT(INOUT) :: this
! local
LOGICAL fk
INTEGER ierr
this%mdim=-1; this%mrange=-1; this%metype=-1; this%mname=''
NULLIFY(this%mdynadata)
NULLIFY(this%mprev)
NULLIFY(this%mnext)
NULLIFY(this%mhead)
IF(ASSOCIATED(this%mint4d0p)) THEN
DEALLOCATE(this%mint4d0p,STAT=ierr); NULLIFY(this%mint4d0p)
ENDIF
IF(ASSOCIATED(this%mint4d1p)) THEN
DEALLOCATE(this%mint4d1p,STAT=ierr)
NULLIFY(this%mint4d1p)
ENDIF
IF(ASSOCIATED(this%mint4d2p)) THEN
DEALLOCATE(this%mint4d2p,STAT=ierr)
NULLIFY(this%mint4d2p)
ENDIF
IF(ASSOCIATED(this%mint4d3p)) THEN
DEALLOCATE(this%mint4d3p,STAT=ierr)
NULLIFY(this%mint4d3p)
ENDIF
IF(ASSOCIATED(this%mint4d4p)) THEN
DEALLOCATE(this%mint4d4p,STAT=ierr)
NULLIFY(this%mint4d4p)
ENDIF
IF(ASSOCIATED(this%mint4d5p)) THEN
DEALLOCATE(this%mint4d5p,STAT=ierr)
NULLIFY(this%mint4d5p)
ENDIF
IF(ASSOCIATED(this%mreal8d0p)) THEN
DEALLOCATE(this%mreal8d0p,STAT=ierr); NULLIFY(this%mreal8d0p)
ENDIF
IF(ASSOCIATED(this%mreal8d1p)) THEN
DEALLOCATE(this%mreal8d1p,STAT=ierr); NULLIFY(this%mreal8d1p)
ENDIF
IF(ASSOCIATED(this%mreal8d2p)) THEN
DEALLOCATE(this%mreal8d2p,STAT=ierr); NULLIFY(this%mreal8d2p)
ENDIF
IF(ASSOCIATED(this%mreal8d3p)) THEN
DEALLOCATE(this%mreal8d3p,STAT=ierr); NULLIFY(this%mreal8d3p)
ENDIF
IF(ASSOCIATED(this%mreal8d4p)) THEN
DEALLOCATE(this%mreal8d4p,STAT=ierr); NULLIFY(this%mreal8d4p)
ENDIF
IF(ASSOCIATED(this%mreal8d5p)) THEN
DEALLOCATE(this%mreal8d5p,STAT=ierr); NULLIFY(this%mreal8d5p)
ENDIF
END SUBROUTINE
SUBROUTINE dyd_new(head,elementtype, dataname, ierr, &
range11, range12,range21, range22,range31, range32, &
range41, range42,range51, range52)
TYPE(dynadata),POINTER,INTENT(INOUT) :: head
CHARACTER(LEN=*),INTENT(IN) :: elementtype, dataname
INTEGER,INTENT(IN),OPTIONAL :: range11, range12,range21, range22, &
range31, range32, range41, range42,range51, range52
INTEGER,INTENT(OUT) :: ierr
! local
INTEGER :: atype,adim
INTEGER :: arange(2,10)
TYPE(dynadata),POINTER :: pdyd
CHARACTER(LEN=32) :: tstr
ierr=0
IF(.NOT.dyd_isvalidname(head,dataname)) THEN
ierr=1
RETURN
ENDIF
! get type
tstr = TRIM(elementtype)
CALL upper_case(tstr)
IF(tstr=='INTEGER(4)') THEN; atype=int4
ELSEIF(tstr=='REAL(8)') THEN; atype=real8
ELSE
ierr=2
RETURN
ENDIF
! get
adim=0
IF( PRESENT(range11) .AND. PRESENT(range12) ) THEN
IF(range11>=range12) THEN
WRITE(*,*) 'Invalid range'; ierr=3
RETURN
ENDIF
adim=1; arange(1,1)=range11; arange(2,1)=range12
! DIM 2
IF( PRESENT(range21) .AND. PRESENT(range22) ) THEN
IF(range21>=range22) THEN
WRITE(*,*) 'Invalid range'; ierr=3
RETURN
ENDIF
adim=2; arange(1,2)=range21; arange(2,2)=range22
! DIM 3
IF( PRESENT(range31) .AND. PRESENT(range32) ) THEN
IF(range31>=range32) THEN
WRITE(*,*) 'Invalid range'; ierr=3
RETURN
ENDIF
adim=3; arange(1,3)=range31; arange(2,3)=range32
! DIM 4
IF( PRESENT(range41) .AND. PRESENT(range42) ) THEN
IF(range41>=range42) THEN
WRITE(*,*) 'Invalid range'; ierr=3
RETURN
ENDIF
adim=4; arange(1,4)=range41; arange(2,4)=range42
! DIM 5
IF( PRESENT(range51) .AND. PRESENT(range52) ) THEN
IF(range51>=range52) THEN
WRITE(*,*) 'Invalid range'; ierr=3
RETURN
ENDIF
adim=5; arange(1,5)=range51; arange(2,5)=range52
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
!WRITE(*,*) 'dyd_new ',atype,adim,arange(:,1:adim),dataname
CALL dyd_new_0(head,atype,adim,arange,dataname)
END SUBROUTINE
SUBROUTINE dyd_new_1(head,elementtype, arrayrange, dataname, ierr)
TYPE(dynadata),POINTER,INTENT(INOUT) :: head
CHARACTER(LEN=*),INTENT(IN) :: elementtype, arrayrange,dataname
INTEGER,INTENT(OUT) :: ierr
! local
INTEGER :: atype,offset,suboffset,i,adim
INTEGER :: arange(2,10), loc(20),com(10)
TYPE(dynadata),POINTER :: pdyd
CHARACTER(LEN=32) :: tstr
!WRITE(*,*) 'new1 ',dataname
ierr=0
IF(.NOT.dyd_isvalidname(head,dataname)) THEN
ierr=1
RETURN
ENDIF
! get type
tstr = TRIM(elementtype)
CALL upper_case(tstr)
IF(tstr=='INTEGER(4)') THEN; atype=int4
ELSEIF(tstr=='REAL(8)') THEN; atype=real8
ELSE
ierr=2
RETURN
ENDIF
! get
adim=0; offset=1
DO
suboffset = INDEX(arrayrange(offset:),':')
IF(suboffset==0) EXIT
adim = adim+1
IF(adim>5) STOP
offset = offset+suboffset
ENDDO
!WRITE(*,*) 'adim=',adim,' ',arrayrange
IF(adim>0) THEN
loc(1) = INDEX(arrayrange,'(')
loc(adim+1) = INDEX(arrayrange,')')
com(1) = INDEX(arrayrange,':')
IF(adim>1) THEN
DO i=2,adim
loc(i) = loc(i-1)+INDEX(arrayrange(loc(i-1)+1:),',')
!write(*,*) loc(i),' ',arrayrange(loc(i-1)+1:)
com(i) = loc(i)+INDEX(arrayrange(loc(i)+1:),':')
!write(*,*) com(i),' ',arrayrange(loc(i)+1:)
ENDDO
ENDIF
write(*,*) loc(1:adim+1)
write(*,*) com(1:adim)
DO i=1,adim
READ( arrayrange(loc(i)+1 : com(i)-1 ), *) arange(1,i)
READ( arrayrange(com(i)+1 : loc(i+1)-1), *) arange(2,i)
ENDDO
ENDIF
CALL dyd_new_0(head,atype,adim,arange,dataname)
!NULLIFY(pdyd)
!ALLOCATE(pdyd)
!CALL dyd_construct(pdyd, atype,adim,arange,dataname)
!CALL dyd_insert(pdyd,head)
END SUBROUTINE
SUBROUTINE dyd_new_0(head,atype,adim,arange,aname)
TYPE(dynadata),POINTER,INTENT(INOUT) :: head
INTEGER,INTENT(IN) :: atype,adim
INTEGER,INTENT(IN),OPTIONAL :: arange(2,adim)
CHARACTER(LEN=*),INTENT(IN) :: aname
! local
TYPE(dynadata),POINTER :: pdyd
NULLIFY(pdyd)
IF(dyd_isvalidname(head,aname)) THEN
ALLOCATE(pdyd)
CALL dyd_construct(pdyd, atype,adim,arange,aname)
CALL dyd_insert(pdyd,head)
ENDIF
END SUBROUTINE
SUBROUTINE dyd_construct(this,atype,adim,arange,aname)
TYPE(dynadata),INTENT(INOUT) :: this
INTEGER,INTENT(IN) :: atype,adim
INTEGER,INTENT(IN),OPTIONAL :: arange(2,adim)
CHARACTER(LEN=*),INTENT(IN) :: aname
!
INTEGER i
CALL dyd_empty(this)
IF(atype<=0 .OR. atype>2) THEN
WRITE(*,*) 'Invalid type: ',atype
STOP
ENDIF
IF(adim<0) THEN
WRITE(*,*) 'Invalid dim: ',adim
STOP
ENDIF
IF(adim>0) THEN; DO i=1,adim
IF(arange(1,i)>arange(2,i)) THEN
WRITE(*,*) 'Invalid range: ',arange(1:2,i), ' in DIM ',i
STOP
ENDIF
ENDDO; ENDIF
this%mname = aname
this%metype=atype
this%mdim=adim
IF(adim>0) this%mrange(1:2,1:adim) = arange(1:2,1:adim)
SELECT CASE(this%metype)
CASE (int4)
IF(this%mdim==0) THEN
ALLOCATE(this%mint4d0p)
ENDIF
IF(this%mdim==1) THEN
ALLOCATE(this%mint4d1p(this%mrange(1,1):this%mrange(2,1) ) )
ENDIF
IF(this%mdim==2) THEN
ALLOCATE(this%mint4d2p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2) ) )
ENDIF
IF(this%mdim==3) THEN
ALLOCATE(this%mint4d3p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3) ) )
ENDIF
IF(this%mdim==4) THEN
ALLOCATE(this%mint4d4p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3), &
this%mrange(1,4):this%mrange(2,4) ) )
ENDIF
IF(this%mdim==5) THEN
ALLOCATE(this%mint4d5p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3), &
this%mrange(1,4):this%mrange(2,4), &
this%mrange(1,5):this%mrange(2,5) ) )
ENDIF
CASE (real8)
IF(this%mdim==0) THEN
ALLOCATE(this%mreal8d0p)
ENDIF
IF(this%mdim==1) THEN
ALLOCATE(this%mreal8d1p(this%mrange(1,1):this%mrange(2,1) ) )
ENDIF
IF(this%mdim==2) THEN
ALLOCATE(this%mreal8d2p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2) ) )
ENDIF
IF(this%mdim==3) THEN
ALLOCATE(this%mreal8d3p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3) ) )
ENDIF
IF(this%mdim==4) THEN
ALLOCATE(this%mreal8d4p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3), &
this%mrange(1,4):this%mrange(2,4) ) )
ENDIF
IF(this%mdim==5) THEN
ALLOCATE(this%mreal8d5p(this%mrange(1,1):this%mrange(2,1), &
this%mrange(1,2):this%mrange(2,2), &
this%mrange(1,3):this%mrange(2,3), &
this%mrange(1,4):this%mrange(2,4), &
this%mrange(1,5):this%mrange(2,5) ) )
ENDIF
END SELECT
END SUBROUTINE
SUBROUTINE dyd_insert(this,head)
TYPE(dynadata),POINTER,INTENT(INOUT) :: this
TYPE(dynadata),POINTER,INTENT(INOUT) :: head
!
TYPE(dynadata),POINTER :: tail
IF( ASSOCIATED(head) ) THEN
CALL dyd_gettail(head,tail)
!WRITE(*,*) 'Tail name is ',TRIM(tail%mname)
this%mprev => tail
this%mhead => head
tail%mnext => this
NULLIFY(this%mnext)
ELSE
head => this
this%mprev => NULL()
this%mhead => this
this%mnext => NULL()
ENDIF
END SUBROUTINE
SUBROUTINE dyd_gettail(this,tail)
TYPE(dynadata),POINTER,INTENT(IN) :: this
TYPE(dynadata),POINTER,INTENT(OUT) :: tail
IF(ASSOCIATED(this)) THEN
tail => this
DO WHILE( ASSOCIATED(tail%mnext) )
tail => tail%mnext
END DO
ENDIF
END SUBROUTINE
FUNCTION dyd_isvalidname(head,aname)
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
LOGICAL dyd_isvalidname
!
TYPE(dynadata),POINTER :: iter
dyd_isvalidname = .TRUE.
IF( ASSOCIATED(head) ) THEN
iter=>head
DO WHILE( ASSOCIATED(iter%mnext) )
!WRITE(*,*) TRIM(aname), ' ',TRIM(iter%mname)
!READ(*,*)
IF(TRIM(aname)==TRIM(iter%mname)) THEN
dyd_isvalidname = .FALSE.
EXIT
ENDIF
iter => iter%mnext
END DO
ENDIF
END FUNCTIOn
FUNCTION dyd_findnode(head,aname)
TYPE(dynadata),POINTER,INTENT(IN) :: head
CHARACTER(LEN=*),INTENT(IN) :: aname
TYPE(dynadata),POINTER :: dyd_findnode
!
TYPE(dynadata),POINTER :: iter=>NULL()
NULLIFY(dyd_findnode)
IF(.NOT.ASSOCIATED(head)) RETURN
iter=>head
loop: DO
IF(TRIM(aname)==TRIM(iter%mname)) THEN
dyd_findnode => iter ! 找到
EXIT loop
ENDIF
IF(.NOT.ASSOCIATED(iter%mnext)) THEN
EXIT loop ! 队列结束，没找到
ELSE
iter=>iter%mnext
ENDIF
END DO loop
END FUNCTION
! other
elemental subroutine upper_case(word)
! convert a word to lower case
character (len=*) , intent(in out) :: word
integer :: i,ic,nlen
nlen = len(word)
do i=1,nlen
ic = ichar(word(i:i))
if (ic >= 97 .and. ic < 122) word(i:i) = char(ic-32)
end do
end subroutine upper_case
END MODULE
测试代码如下
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE TEST_int4d0
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),POINTER :: X
INTEGER ::ierr
!CALL dyd_new(head=list, atype=int4, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
WRITE(*,*) X
X = 5
WRITE(*,*)iter%mint4d0p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_int4d1
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),DIMENSION(:),POINTER :: X
INTEGER ::i, irange(2,2),ierr
irange(1,1) =0; irange(2,1) =9
!CALL dyd_new(head=list, atype=int4, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
range11=0,range12=9)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
WRITE(*,*) X
WRITE(*,*) SHAPE(X),LBOUND(X),UBOUND(X)
X = 5
DO i=0,9; X(i) = i; ENDDO
WRITE(*,*)iter%mint4d1p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_int4d2
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),DIMENSION(:,:),POINTER :: X
INTEGER ::i,j,k, irange(2,2),ierr
irange(1,1) =0; irange(2,1) =9
irange(1,2) =1; irange(2,2) =5
!CALL dyd_new(head=list, atype=int4, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
range11=0,range12=9,range21=1,range22=5)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO j=1,5;DO i=0,9; X(i,j) = k; k=k+1; ENDDO;ENDDO
WRITE(*,*)iter%mint4d2p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_int4d3
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),DIMENSION(:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,3),ierr
irange(1,1) =0; irange(2,1) =9
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
!CALL dyd_new(head=list, atype=int4, adim=3,arange=irange, aname='X')
CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9,1:5,1:3)', &
dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO l=1,3; DO j=1,5;DO i=0,9
X(i,j,l) = k; k=k+1
ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mint4d3p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_int4d4
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),DIMENSION(:,:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,4),ierr
irange(1,1) =0; irange(2,1) =2
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
irange(1,4) =-1; irange(2,4) =1
CALL dyd_new(head=list, atype=int4, adim=4,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:9,1:5,1:3)', &
! dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO m=irange(1,4), irange(2,4)
DO l=irange(1,3), irange(2,3)
DO j=irange(1,2), irange(2,2)
DO i=irange(1,1), irange(2,1)
X(i,j,l,m) = k; k=k+1
ENDDO;ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mint4d4p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_int4d5
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
INTEGER(4),DIMENSION(:,:,:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,5),ierr
irange(1,1) =0; irange(2,1) =2
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
irange(1,4) =-1; irange(2,4) =1
irange(1,5) =5; irange(2,5) =8
!CALL dyd_new(head=list, atype=int4, adim=5,arange=irange, aname='X')
CALL dyd_new(head=list, elementtype='integer(4)', arrayrange='(0:2,1:5,1:3,-1:1,5:8)', &
dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='integer(4)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mint4d1p
!WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO n=irange(1,5), irange(2,5)
DO m=irange(1,4), irange(2,4)
DO l=irange(1,3), irange(2,3)
DO j=irange(1,2), irange(2,2)
DO i=irange(1,1), irange(2,1)
X(i,j,l,m,n) = k; k=k+1
ENDDO;ENDDO;ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mint4d5p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d0
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),POINTER :: X
INTEGER ::ierr
!CALL dyd_new(head=list, atype=real8, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
X = 5
WRITE(*,*)iter%mreal8d0p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d1
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),DIMENSION(:),POINTER :: X
INTEGER ::i, irange(2,2),ierr
irange(1,1) =0; irange(2,1) =9
!CALL dyd_new(head=list, atype=real8, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr, &
range11=0,range12=9)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
WRITE(*,*) SHAPE(X),LBOUND(X),UBOUND(X)
X = 5
DO i=0,9; X(i) = i; ENDDO
WRITE(*,*)iter%mreal8d1p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d2
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),DIMENSION(:,:),POINTER :: X
INTEGER ::i,j,k, irange(2,2),ierr
irange(1,1) =0; irange(2,1) =9
irange(1,2) =1; irange(2,2) =5
!CALL dyd_new(head=list, atype=real8, adim=1,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:9)', &
! dataname='X', ierr=ierr)
CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr, &
range11=0,range12=9,range21=1,range22=5)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO j=1,5;DO i=0,9; X(i,j) = k; k=k+1; ENDDO;ENDDO
WRITE(*,*)iter%mreal8d2p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d3
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),DIMENSION(:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,3),ierr
irange(1,1) =0; irange(2,1) =9
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
!CALL dyd_new(head=list, atype=real8, adim=3,arange=irange, aname='X')
CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:9,1:5,1:3)', &
dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO l=1,3; DO j=1,5;DO i=0,9
X(i,j,l) = k; k=k+1
ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mreal8d3p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d4
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),DIMENSION(:,:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,4),ierr
irange(1,1) =0; irange(2,1) =2
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
irange(1,4) =-1; irange(2,4) =1
CALL dyd_new(head=list, atype=real8, adim=4,arange=irange, aname='X')
!CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:9,1:5,1:3)', &
! dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO m=irange(1,4), irange(2,4)
DO l=irange(1,3), irange(2,3)
DO j=irange(1,2), irange(2,2)
DO i=irange(1,1), irange(2,1)
X(i,j,l,m) = k; k=k+1
ENDDO;ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mreal8d4p
CALL dyd_empty(iter)
END SUBROUTINE
SUBROUTINE TEST_real8d5
USE xx
IMPLICIT NONE
TYPE(dynadata),POINTER :: list=>NULL(),iter=>NULL()
REAL(8),DIMENSION(:,:,:,:,:),POINTER :: X
INTEGER ::i,j,k,l,m,n, irange(2,5),ierr
irange(1,1) =0; irange(2,1) =2
irange(1,2) =1; irange(2,2) =5
irange(1,3) =1; irange(2,3) =3
irange(1,4) =-1; irange(2,4) =1
irange(1,5) =5; irange(2,5) =8
!CALL dyd_new(head=list, atype=real8, adim=5,arange=irange, aname='X')
CALL dyd_new(head=list, elementtype='REAL(8)', arrayrange='(0:2,1:5,1:3,-1:1,5:8)', &
dataname='X', ierr=ierr)
!CALL dyd_new(head=list, elementtype='REAL(8)', dataname='X', ierr=ierr, &
! range11=0,range12=9,range21=1,range22=5,range31=1,range32=3)
CALL getdata(list,'X', X)
iter => dyd_findnode(head=list,aname='X'); !X => iter%mreal8d1p
!WRITE(*,*) X
WRITE(*,*) 'SHAPE(X)=',SHAPE(X),'LBOUND(X)=',LBOUND(X),'UBOUND(X)=',UBOUND(X)
k=1
DO n=irange(1,5), irange(2,5)
DO m=irange(1,4), irange(2,4)
DO l=irange(1,3), irange(2,3)
DO j=irange(1,2), irange(2,2)
DO i=irange(1,1), irange(2,1)
X(i,j,l,m,n) = k; k=k+1
ENDDO;ENDDO;ENDDO;ENDDO;ENDDO
WRITE(*,*)iter%mreal8d5p
CALL dyd_empty(iter)
END SUBROUTINE
PROGRAM TEST
IMPLICIT NONE
! CALL TEST_int4d0
! CALL TEST_int4d1
! CALL TEST_int4d2
! CALL TEST_int4d3
! CALL TEST_int4d4
! CALL TEST_int4d5
CALL TEST_real8d0
CALL TEST_real8d1
CALL TEST_real8d2
CALL TEST_real8d3
CALL TEST_real8d4
CALL TEST_real8d5
END
还好,fortran的指针还不算太差
指针可传给c语言做void*
program test implicit none integer,pointer,dimension(:,:) :: dat, ptr integer, target :: buff(10,10) integer(8) :: n,l,leng,buffsize leng = 100 buffsize = leng*sizeof(buff(1,1)) allocate(dat(1:10,1:10)) leng = 1 do n=1,10 do l=1,10 dat(l,n) = leng leng = leng + 1 enddo enddo leng = 100 ptr => dat call cwrite(ptr, buffsize) ptr => buff call cread(ptr,buffsize) do n=1,10 write(*,'(10I5)') buff(:,n) enddo deallocate(dat) end
#include <stdio.h>
void cwrite_(void * buff, long int* len)
{
long int i;
FILE * fp;
fp = fopen("test.dat","wb");
fwrite(buff,1,*len,fp);
fclose(fp);
}
void cread_(void * buff, long int* len)
{
long int i;
FILE * fp;
fp = fopen("test.dat","rb");
fread(buff,1,*len,fp);
fclose(fp);
}