効率の良いfortranコードを書きたい

プログラマーではないのだが研究の大部分を占めているのがプログラミング。しかも海洋モデルはfotranで書かれてあるのでいやでもfortranを習得しなければならない。最近並列コンピューティングでシミュレーションをするが大きなデータを扱うのでかなり時間がかかる。そこで少しでも改善の余地があればと思いfortranを効率よく書く秘訣みたいなものを色々ググってみた。以下、その中で理解でき且つ簡単に使えそうなのを書き留めた。

１。do loopは後ろから

do k = 1, kt
  do j = 1, jt
    do i = 1, it
      A(i, j, k) = ...
end do
end do
end do

とする方が

do i = 1, it
  do j = 1, jt
    do k = 1, kt
      A(i, j, k) = ...
end do
end do
end do

より効率いいらしい。


２。同じサイズのarrayはまとめて

real, dimension(3, n) :: a

とする方が

real, dimension(n) :: a, b, c

と三つに分けるよりいいらしい。


1＆２のソース：https://jblevins.org/log/efficient-code

３。コロン（：）は左側に

max(A(:, :, i))

とする方が

max(A(i, :, :))

より計算が早いらしい。

ソース：http://www.fortran90.org/src/best-practices.html


４。割り算より掛け算

(A / 2) / 4

より

(A * 0.5) * 0.25

の方が早いらしい。



ソース：https://www.allinea.com/blog/201407/four-simple-tips-optimizing-your-code


最後に。

コンパイラーに最適化のフラッグがあることに気づいた。
-Oというフラッグで-O3がベスト（コンパイルされる時間は増えるが）。
上記にあげたこともしてくれたりするのかな？

interesting segfault

I had a segmentation fault which is due to NaN. However, this was interesting because although the IF statement was always false, it gave me segfault because NaN would result if the IF statement was true.. 

        write(numout,*) 'hhpass: nancheck ',(1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon30)
        IF(frld(ji,jj)-pfrld(ji,jj) > 0.) THEN ! In case of accretion, newly formed ice will have a tracer concentration of an underlying seawater value.
          write(numout,*) 'hhpass: inside if statement'
          icechl(ji,jj) = icechl(ji,jj)/(trn(ji,jj,1,jpdch)+epsilon30)!((1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon30)*icechl(ji,jj) + ((1.-frld(ji,jj))-(1.-pfrld(ji,jj)))/(1.-frld(ji,jj)+epsilon30)*trn(ji,jj,1,jpdch)
        END IF


 Update: This seems to have to do with floating points.

Writing out the output for the variables in dinomenator and numerator show a particular floating point (i.e. 14):

 hhpass: pfrld    1.00000000000000
 hhpass: frld    1.00000000000000





In the example above, I tried avoid the division by zero by introducing epsilon number 30 (i.e. 1.e-30), but the division produced NaN:

hhpass: nancheck    NaN

After this, I instead introduced a epsilon number that is smaller than the floating point of the two variables (i.e. 4). Then the division produced zero:

 hhpass: nancheck   0.000000000000000E+000


Epsilon30 worked for biogeochemical tracers. But not for ice physical tracers (probably a specic floating point number is defined in LIM_SRC_2 for these tracers?)

macportsでgfortranとnetcdf-fortranのバージョンをマッチさせる手順

前回の日記では具体的な方法を明記していなかったので今回また同じ機会に遭遇したのでメモしておくことに。

まずはgfortran

以下はバージョン4.6をインストールしたい場合のコマンド

sudo port install gcc46

で、これだけではgcc46しかインストールされていないので以下のコマンドでgfortranを使用可能にする

sudo port select gcc mp-gcc46

これでgfortran 4.6がインストール完了。

次にnetcdf-fortranも同じ4.6に対応したものをインストールする。以下のコマンドで。

sudo port install netcdf-fortran +gcc46

以上！

macportsのnetcdf-fortranとgfortranのバージョンはマッチさせんといかんばい

同僚二人のマックブックにGOTMとFABMモデルをインストールしてコンパイルする作業の手伝いをしたがどうもコンパイルがうまくいかない。

use netcdfだったりuse fabm_typesにエラーがどうのこうのというメッセージがでてしまう。

自分のマックで一年ほど前にやったときは問題なかったのになんで最新のOSでできないのか色々みんなで苦しんだ結果、問題はmacportsでインストールしたnetcdf-fortranのバージョンがgfortranのサイトからインストールしたgfortranのバージョンとマッチしていなかったから。gfortranのバージョンは

gfortran --version

で確認できる。以下が出力例：

GNU Fortran (GCC) 4.8.2

一方macportsでインストールしたnetcdf-fortranのバージョンは

port installed netcdf-fortran

で確認できる。以下がその出力例:
 
The following ports are currently installed:
  netcdf-fortran @4.2_12+gcc48 (active)

この場合、gfortranは4.8.2、netcdf-fortranは4.8なのでオーケー。

これが例えば、4.9と4.8のコンビネーションだとミスマッチでコンパイルでエラーになる。

failed to cmake, but succeeded with make

cmake -L ~/GOTM/gotm-git/src 
.
.
.
NetCDF_INCLUDE_DIRS:PATH=/usr/include
NetCDF_LIBRARIES:STRING=-lnetcdf

"make" produces an error (not locating the appropriate netcdf directory).

head /usr/include/netcdf.inc 

!     NetCDF-3.

!

! netcdf version 3 fortran interface:

.

.

.

Thus, by the current (default) setting, cmake is choosing netcdf3, which cannot be used in FABM-GOTM (we need netcdf4). We need to manually set the location of netcdf libraries (and include as well?).

cmake -L -DNetCDF_LIBRARIES=/usr/local/netcdf-4.3.0_GF/lib ~/GOTM/gotm-git/src
-- Found NetCDF: /usr/local/netcdf-4.3.0_GF/lib
-- Configuring done
WARNING: Target "gotm_exe" requests linking to directory "/usr/local/netcdf-4.3.0_GF/lib".  Targets may link only to libraries.  CMake is dropping the item.
.
.
.
NetCDF_INCLUDE_DIRS:PATH=/usr/include

NetCDF_LIBRARIES:STRING=/usr/local/netcdf-4.3.0_GF/lib

"make" still produces an error. probably /usr/include is not right and I need to specify the path for it as well.

cmake -L -DNetCDF_LIBRARIES=/usr/local/netcdf-4.3.0_GF/lib -DNetCDF_INCLUDE_DIRS=/usr/local/netcdf-4.3.0_GF/include ~/GOTM/gotm-git/src
.
.
.
NetCDF_INCLUDE_DIRS:PATH=/usr/local/netcdf-4.3.0_GF/include

NetCDF_LIBRARIES:STRING=/usr/local/netcdf-4.3.0_GF/lib

make
.
.
.
gotm_fabm_output.F90:(.text+0x217a): undefined reference to `__netcdf_MOD_nf90_put_att_one_fourbyteint'
collect2: error: ld returned 1 exit status
make[2]: *** [gotm] Error 1
make[1]: *** [CMakeFiles/gotm_exe.dir/all] Error 2
make: *** [all] Error 2

Still the same error...

I gave up on "cmake" for now. Instead I tried the good old "make" compilation in ~/GOTM/gotm-git/src.

and I was able to do compile with the following script:

doit.sh
#! /bin/sh
export PATH=/usr/local/netcdf-4.3.0_GF/bin:$PATH
export FORTRAN_COMPILER=GFORTRAN
export NETCDF_VERSION=NETCDF4

make

Providing the path to netcdf bin folder helped.

This can be coupled to FABM which is built using the new cmake system:

doit_couple.sh
#! /bin/sh
export PATH=/usr/local/netcdf-4.3.0_GF/bin:$PATH
export FABMDIR=$HOME/FABM/fabm-git
export FORTRAN_COMPILER=GFORTRAN
export FABM=true
export FABM_PREFIX=$HOME/local/fabm/gotm
export NETCDF_VERSION=NETCDF4

make

For now, this is sufficient to couple FABM and GOTM.

segfault issues solved.

Segmentation fault error

When execute "gotm" and get the following error below, "make distclean" followed by "make" in "GOTM/gotm-git/src" will fix the problem. Not sure why though...





    init_gotm_fabm
 Initializing biogeochemical model "uvic_eco"...
*** glibc detected *** gotm: free(): invalid next size (normal): 0x0000000002139fb0 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x7eb96)[0x7f682983db96]
gotm[0x73ea5b]
gotm[0x414a81]
gotm[0x405866]
gotm[0x403b61]
gotm[0x4035ad]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed)[0x7f68297e076d]
gotm[0x4035dd]
======= Memory map: ========
00400000-008bd000 r-xp 00000000 08:02 13501193                           /HOME/hakase/GOTM/gotm-git/bin/gotm_prod_GFORTRAN
00abc000-00abd000 r--p 004bc000 08:02 13501193                           /HOME/hakase/GOTM/gotm-git/bin/gotm_prod_GFORTRAN
00abd000-00acf000 rw-p 004bd000 08:02 13501193                           /HOME/hakase/GOTM/gotm-git/bin/gotm_prod_GFORTRAN
00acf000-00b5c000 rw-p 00000000 00:00 0
02107000-02143000 rw-p 00000000 00:00 0                                  [heap]
7f682936d000-7f6829382000 r-xp 00000000 08:02 10092768                   /lib/x86_64-linux-gnu/libgcc_s.so.1
7f6829382000-7f6829581000 ---p 00015000 08:02 10092768                   /lib/x86_64-linux-gnu/libgcc_s.so.1
7f6829581000-7f6829582000 r--p 00014000 08:02 10092768                   /lib/x86_64-linux-gnu/libgcc_s.so.1
7f6829582000-7f6829583000 rw-p 00015000 08:02 10092768                   /lib/x86_64-linux-gnu/libgcc_s.so.1
7f6829583000-7f68295be000 r-xp 00000000 08:02 4590874                    /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f68295be000-7f68297bd000 ---p 0003b000 08:02 4590874                    /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f68297bd000-7f68297be000 r--p 0003a000 08:02 4590874                    /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f68297be000-7f68297bf000 rw-p 0003b000 08:02 4590874                    /usr/lib/x86_64-linux-gnu/libquadmath.so.0.0.0
7f68297bf000-7f6829974000 r-xp 00000000 08:02 10092621                   /lib/x86_64-linux-gnu/libc-2.15.so
7f6829974000-7f6829b74000 ---p 001b5000 08:02 10092621                   /lib/x86_64-linux-gnu/libc-2.15.so
7f6829b74000-7f6829b78000 r--p 001b5000 08:02 10092621                   /lib/x86_64-linux-gnu/libc-2.15.so
7f6829b78000-7f6829b7a000 rw-p 001b9000 08:02 10092621                   /lib/x86_64-linux-gnu/libc-2.15.so
7f6829b7a000-7f6829b7f000 rw-p 00000000 00:00 0
7f6829b7f000-7f6829c7a000 r-xp 00000000 08:02 10092632                   /lib/x86_64-linux-gnu/libm-2.15.so
7f6829c7a000-7f6829e79000 ---p 000fb000 08:02 10092632                   /lib/x86_64-linux-gnu/libm-2.15.so
7f6829e79000-7f6829e7a000 r--p 000fa000 08:02 10092632                   /lib/x86_64-linux-gnu/libm-2.15.so
7f6829e7a000-7f6829e7b000 rw-p 000fb000 08:02 10092632                   /lib/x86_64-linux-gnu/libm-2.15.so
7f6829e7b000-7f6829f90000 r-xp 00000000 08:02 4591034                    /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f6829f90000-7f682a190000 ---p 00115000 08:02 4591034                    /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f682a190000-7f682a191000 r--p 00115000 08:02 4591034                    /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f682a191000-7f682a193000 rw-p 00116000 08:02 4591034                    /usr/lib/x86_64-linux-gnu/libgfortran.so.3.0.0
7f682a193000-7f682a1b5000 r-xp 00000000 08:02 10092633                   /lib/x86_64-linux-gnu/ld-2.15.so
7f682a305000-7f682a38a000 rw-p 00000000 00:00 0
7f682a3b2000-7f682a3b5000 rw-p 00000000 00:00 0
7f682a3b5000-7f682a3b6000 r--p 00022000 08:02 10092633                   /lib/x86_64-linux-gnu/ld-2.15.so
7f682a3b6000-7f682a3b8000 rw-p 00023000 08:02 10092633                   /lib/x86_64-linux-gnu/ld-2.15.so
7fff0954d000-7fff0956e000 rw-p 00000000 00:00 0                          [stack]
7fff095d9000-7fff095db000 r-xp 00000000 00:00 0                          [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:
#0  0x7F6829E94547
#1  0x7F6829E94B4E
#2  0x7F68297F549F
#3  0x7F68297F5425
#4  0x7F68297F8B8A
#5  0x7F682983339D
#6  0x7F682983DB95
#7  0x73EA5A in __fabm_MOD_fabm_create_model_from_file
#8  0x414A80 in __gotm_fabm_MOD_init_gotm_fabm
#9  0x405865 in __gotm_MOD_init_gotm
#10  0x403B60 in MAIN__ at main.F90:0
Aborted (core dumped)

fortranの仕組みがよくわからん…

DMS fluxをいろんなmethodを使って計算できるようにfabm.nmlに新しいパラメータ「flux」を取り入れた。

他のパラメータ同様それをモデルの中で読み込みfluxの数値によってmethodを変えられるようにする仕組みのはずが…

下記の図によるとflux=0,1,2のときDMS fluxのシミュレーション結果がそれぞれ異なることがわかる。

０と１の場合はそれぞれモデルに

if (flux.eq.0) then
...
else if (flux.eq.1) then
...
end

と記したので思っていたとおりの結果。

しかし、２の場合はデフォルト値（０に設定したはず）と同じ結果を期待したんだが同じではない。２の場合の条件は何も書いていないのでなぜ０でも１でもない結果（しかもちゃんとDMS fluxを計算した結果）が出力されたのかがわからん…。

とりあえず解決策として０、１以外数値を選んだ場合はDMS fluxは０になるようにした。

「netcdfとgfortranはマッチアップさせんといかん」の巻

netcdfが/opt/local内（/includeと/lib）にうまくインストールされたのでいざコンパイル！っとGOTMでやってみたが以下のエラーが…

use netcdf
             1
Fatal Error: Wrong module version '10' (expected '6') for file 'netcdf.mod' opened at (1)

調べた結果、どうやら自分がインストールしたnetcdfとこのパソコンにインストールされているgfortranのバージョンがマッチしないのが原因のよう。

GFortranのサイトのMacOSXインストール手順に従って、（１）現在インストールされているgfortranを削除→（２）Mavericks用（最新）のgfortranをインストールした。

するとうまくいった！！！

gen1d.f90

↓最近作ったプログラム（２０１２年夏）。

! This script extracts the 1-D data from 3-D offline dynamics.

program generate1d
implicit none
INCLUDE 'netcdf.inc'
integer,parameter::jpi=499,jpj=604,jpk=46,jpt=31
integer::ncid,ncfile,tempid,saliid,windid,wateid,swrnid,turbid,icecid,dim1,dim2,dim3,dim4,it,dtemp,dsali,dwind,dicec,dswrn,dturb,dwate
real,dimension(jpi,jpj,jpk,jpt)::temp,sali
real,dimension(jpi,jpj,jpt)::wind,wate,swrn,turb,icec
integer,dimension(3)::stnx,stny,dim1d
integer,dimension(4)::dim2d
real,dimension(3,1,jpk,jpt)::temp1,sali1
real,dimension(3,1,jpt)::wind1,wate1,swrn1,turb1,icec1

! Open the data file
ncfile=nf_open('Lab_T_94_07.nc',nf_nowrite,ncid)
! Inquire the variables
ncfile=nf_inq_varid(ncid,'votemper',tempid)
ncfile=nf_inq_varid(ncid,'vosaline',saliid)
ncfile=nf_inq_varid(ncid,'sowindsp',windid)
ncfile=nf_inq_varid(ncid,'sowaflcd',wateid)
ncfile=nf_inq_varid(ncid,'soshfldo',swrnid)
ncfile=nf_inq_varid(ncid,'somixhgt',turbid)
ncfile=nf_inq_varid(ncid,'soicecov',icecid)
! Extract 3D data
ncfile=nf_get_var_real(ncid,tempid,temp)
ncfile=nf_get_var_real(ncid,saliid,sali)
ncfile=nf_get_var_real(ncid,windid,wind)
ncfile=nf_get_var_real(ncid,wateid,wate)
ncfile=nf_get_var_real(ncid,swrnid,swrn)
ncfile=nf_get_var_real(ncid,turbid,turb)
ncfile=nf_get_var_real(ncid,icecid,icec)
! Close the data file
ncfile=nf_close(ncid)

stnx=(/ 282, 301, 350 /)
stny=(/ 382, 416, 500 /)
do it=1,3
 temp1(it,1,:,:)=temp(stnx(it),stny(it),:,:)
 sali1(it,1,:,:)=sali(stnx(it),stny(it),:,:)
 wind1(it,1,:)=wind(stnx(it),stny(it),:)
 wate1(it,1,:)=wate(stnx(it),stny(it),:)
 swrn1(it,1,:)=swrn(stnx(it),stny(it),:)
 turb1(it,1,:)=turb(stnx(it),stny(it),:)
 icec1(it,1,:)=icec(stnx(it),stny(it),:)
end do

ncfile=nf_create('Lab_T_1d_94_07.nc',nf_share,ncid)
ncfile=nf_def_dim(ncid,'nlon',3,dim1)
ncfile=nf_def_dim(ncid,'nlat',1,dim2)
ncfile=nf_def_dim(ncid,'ndepth',jpk,dim3)
ncfile=nf_def_dim(ncid,'ntime',jpt,dim4)
dim2d(1)=dim1
dim2d(2)=dim2
dim2d(3)=dim3
dim2d(4)=dim4
dim1d(1)=dim1
dim1d(2)=dim2
dim1d(3)=dim4
ncfile=nf_def_var(ncid,'nvotemper',nf_real,4,dim2d,dtemp)
ncfile=nf_def_var(ncid,'nvosaline',nf_real,4,dim2d,dsali)
ncfile=nf_def_var(ncid,'nsowindsp',nf_real,3,dim1d,dwind)
ncfile=nf_def_var(ncid,'nsoshfldo',nf_real,3,dim1d,dswrn)
ncfile=nf_def_var(ncid,'nsoicecov',nf_real,3,dim1d,dicec)
ncfile=nf_def_var(ncid,'nsowaflcd',nf_real,3,dim1d,dwate)
ncfile=nf_def_var(ncid,'nsomixhgt',nf_real,3,dim1d,dturb)
ncfile=nf_enddef(ncid)
ncfile=nf_put_var_real(ncid,dtemp,temp1)
ncfile=nf_put_var_real(ncid,dsali,sali1)
ncfile=nf_put_var_real(ncid,dwind,wind1)
ncfile=nf_put_var_real(ncid,dswrn,swrn1)
ncfile=nf_put_var_real(ncid,dicec,icec1)
ncfile=nf_put_var_real(ncid,dwate,wate1)
ncfile=nf_put_var_real(ncid,dturb,turb1)
ncfile=nf_close(ncid)

end program generate1d

subroutine handle_err(errcode)
        implicit none
        include 'netcdf.inc'
        integer errcode
       
        print *, 'Error: ', nf_strerror(errcode)
        stop 2
        end

divide_5day.f90

昔作ったプログラム、その参。

! This program divides 5day-mean data
! so that the file size does not exceed the memory size on stokes
! created by Hakase Hayashida on Jan.24.2012

! Replace the ones with "!" symbols to the appropriate ones
program dividenetcdfdata
implicit none
INCLUDE 'netcdf.inc'
integer,parameter::jpi=304,jpj=334,jpk=46,jpt=72
integer::ncid,ncfile,dim1,dim2,dim3,dim4,divider,dimtx(3)!,dimtx(4)
integer::tempid,deftemp1,deftemp2,deftemp3,deftemp4
integer::vartemp
real,dimension(jpi,jpj,jpt)::temp
real,dimension(jpi,jpj,18)::temp1,temp2,temp3,temp4
!real,dimension(jpi,jpj,jpk,jpt)::temp
!real,dimension(jpi,jpj,jpk,18)::temp1,temp2,temp3,temp4

!Opening netcdf data file from grid_t,v,u,w
ncfile=nf_open('grid_t29.nc',nf_nowrite,ncid) !
ncfile=nf_inq_varid(ncid,'sohefldo',tempid) !
ncfile=nf_get_var_real(ncid,tempid,temp)
ncfile=nf_close(ncid)

!Dividing the data
divider=18
temp1=temp(:,:,1:divider)
temp2=temp(:,:,divider+1:divider*2)
temp3=temp(:,:,1+divider*2:divider*3)
temp4=temp(:,:,1+divider*3:divider*4)

!Creating netcdf data file for temperature
ncfile=nf_create('lab29f.nc',nf_share,ncid) !
ncfile=nf_def_dim(ncid,'x',jpi,dim1)
ncfile=nf_def_dim(ncid,'y',jpj,dim2)
!ncfile=nf_def_dim(ncid,'deptht',jpk,dim3)
ncfile=nf_def_dim(ncid,'time_counter',divider,dim4)
dimtx(1)=dim1
dimtx(2)=dim2
!dimtx(3)=dim3
dimtx(3)=dim4
!dimtx(4)=dim4
ncfile=nf_def_var(ncid,'sohefldo1',nf_real,3,dimtx,deftemp1) !
ncfile=nf_def_var(ncid,'sohefldo2',nf_real,3,dimtx,deftemp2) !
ncfile=nf_def_var(ncid,'sohefldo3',nf_real,3,dimtx,deftemp3) !
ncfile=nf_def_var(ncid,'sohefldo4',nf_real,3,dimtx,deftemp4) !
ncfile=nf_enddef(ncid)
ncfile=nf_put_var_real(ncid,deftemp1,temp1)
ncfile=nf_put_var_real(ncid,deftemp2,temp2)
ncfile=nf_put_var_real(ncid,deftemp3,temp3)
ncfile=nf_put_var_real(ncid,deftemp4,temp4)
ncfile=nf_close(ncid)

end program dividenetcdfdata


subroutine handle_err(errcode)
      implicit none
      include 'netcdf.inc'
      integer errcode

      print *, 'Error: ', nf_strerror(errcode)
      stop 2
      end

convert5day2mm.f90

昔作ったプログラム、その弐。

program mmtest
! converts 5-day average data into monthly average data
implicit none
INCLUDE 'netcdf.inc'

integer,parameter::jpi=544,jpj=336,jptcon=1468,jpt=73
integer::i,jan(6),feb(6),mar(6),apr(6),may(6),jun(6)
integer::jul(6),aug(7),sep(6),oct(6),nov(6),dec(6)
character(2)::yrpar(10),monpar(12)
character(100)::filnam
character(5)::dirnam
character(8)::param
integer::ncfile,varid,value,ncid,ncid2,newvarid,r,rrr
! rgyr=regularyear(365days),lpyr=leapyear(366days)
integer::idvar,dim1,dim2,dim3,dim4,dimty(3)
real,dimension(jpi,jpj,46,jpt)::varhu
real,dimension(544,336,46)::janhu,febhu,marhu,aprhu,mayhu,junhu,julhu,sephu,octhu,novhu,dechu
real,dimension(544,336,46)::mmjan,mmfeb,mmmar,mmapr,mmmay,mmjun,mmjul,mmaug,mmsep,mmoct,mmnov,mmdec
real,dimension(544,336,46)::aughu
real,dimension(544,336,46,12)::mm
real,dimension(12,46,336,544)::mmyear_fixed

yrpar=(/'01','02','03','04','05','06','07','08','09','10'/)

jan=(/(i,i=1,6)/)
feb=(/(i,i=7,12)/)
mar=(/(i,i=13,18)/)
apr=(/(i,i=19,24)/)
may=(/(i,i=25,30)/)
jun=(/(i,i=31,36)/)
jul=(/(i,i=37,42)/)
aug=(/(i,i=43,49)/)
sep=(/(i,i=50,55)/)
oct=(/(i,i=56,61)/)
nov=(/(i,i=62,67)/)
dec=(/(i,i=68,73)/)

do rrr=1,10
!!!!!!!!!!!!!!!! change the three lines !!!!!!!!!!!!!!!!
dirnam='P/POS'
filnam='NAO_'//dirnam//'NAO_5d_00'//yrpar(rrr)//'0101_00'//yrpar(rrr)//'1231_grid_T.nc'
param='vosaline'
ncfile=nf_open(filnam,nf_nowrite,ncid)
ncfile=(nf_inq_varid(ncid,param,varid))
ncfile=(nf_get_var_real(ncid,varid,varhu))
ncfile=nf_close(ncid)

!dat1=varhu(:,:,:,1:6)
!janhu=sum(dat1,dim=4)

janhu=sum(varhu(:,:,:,jan),dim=4)
febhu=sum(varhu(:,:,:,feb),4)
marhu=sum(varhu(:,:,:,mar),4)
aprhu=sum(varhu(:,:,:,apr),4)
mayhu=sum(varhu(:,:,:,may),4)
junhu=sum(varhu(:,:,:,jun),4)
julhu=sum(varhu(:,:,:,jul),4)
aughu=sum(varhu(:,:,:,aug),4)
sephu=sum(varhu(:,:,:,sep),4)
octhu=sum(varhu(:,:,:,oct),4)
novhu=sum(varhu(:,:,:,nov),4)
dechu=sum(varhu(:,:,:,dec),4)

!end do

mm(:,:,:,1)=janhu/float(6)
mm(:,:,:,2)=febhu/float(6)
mm(:,:,:,3)=marhu/float(6)
mm(:,:,:,4)=aprhu/float(6)
mm(:,:,:,5)=mayhu/float(6)
mm(:,:,:,6)=junhu/float(6)
mm(:,:,:,7)=julhu/float(6)
mm(:,:,:,8)=aughu/float(7)
mm(:,:,:,9)=sephu/float(6)
mm(:,:,:,10)=octhu/float(6)
mm(:,:,:,11)=novhu/float(6)
mm(:,:,:,12)=dechu/float(6)

!mmyear=(/mmjan,mmfeb,mmmar,mmapr,mmmay,mmjun,mmjul,mmaug,mmsep,mmoct,mmnov,mmdec/)
!mmyear=reshape((/mmjan,mmfeb,mmmar,mmapr,mmmay,mmjun,mmjul,mmaug,mmsep,mmoct,mmnov,mmdec/),shape(mmyear))
monpar=(/'01','02','03','04','05','06','07','08','09','10','11','12'/)

do r=1,12

ncfile=nf_create('NAO_'//dirnam//'NAO_grid_'//param//'_year'//yrpar(rrr)//'_month' //monpar(r)//'.nc',nf_share,ncid2)
ncfile=nf_def_dim(ncid2,'x',544,dim1)
ncfile=nf_def_dim(ncid2,'y',336,dim2)
ncfile=nf_def_dim(ncid2,'z',46,dim3)
!ncfile=nf_def_dim(ncid2,'month',12,dim4)
dimty(1)=dim1
dimty(2)=dim2
dimty(3)=dim3
!dimty(4)=dim4
ncfile=nf_def_var(ncid2,param,nf_real,3,dimty,newvarid)
ncfile=nf_enddef(ncid2)
ncfile=nf_put_var_real(ncid2,newvarid,mm(:,:,:,r))
ncfile=nf_close(ncid2)

end do
end do
end program mmtest


subroutine h(status)
      integer::status
      INCLUDE 'netcdf.inc'
      if (status .ne. NF_NOERR) then
         print *, NF_STRERROR(status)
         stop 'stopped here'
      end if
      return
      end

convert_5day_to_month.f90

昔作ったフォートランプログラムをここに保存。

! This program extracts temperature data from grid_T file
! so that the file size does not exceed the memory size on stokes
! created by Hakase Hayashida on Jan.18.2012

program convertnetcdfdata
implicit none
INCLUDE 'netcdf.inc'
integer,parameter::jpi=304,jpj=334,jpk=46,jpt=72,nmonth=12
integer::ncid,ncfile,dim1,dim2,dim3,dim4,dimtx(4)
integer::tempid,saliid,meriid,vertid,deftemp,defsali,defmeri,defvert
integer::ind1,tt,vartemp,varsali,varmeri,varvert
real,dimension(jpi,jpj,jpk,jpt)::temp,sali,meri,vert
real,dimension(jpi,jpj,jpk,nmonth)::montemp,monsali,monmeri,monvert

!Opening netcdf data file from grid_t,v,w
ncfile=nf_open('grid_t28.nc',nf_nowrite,ncid)
ncfile=nf_inq_varid(ncid,'votemper',tempid)
ncfile=nf_inq_varid(ncid,'vosaline',saliid)
ncfile=nf_get_var_real(ncid,tempid,temp)
ncfile=nf_get_var_real(ncid,saliid,sali)
ncfile=nf_close(ncid)
ncfile=nf_open('grid_v28.nc',nf_nowrite,ncid)
ncfile=nf_inq_varid(ncid,'vomecrty',meriid)
ncfile=nf_get_var_real(ncid,meriid,meri)
ncfile=nf_close(ncid)
ncfile=nf_open('grid_w28.nc',nf_nowrite,ncid)
ncfile=nf_inq_varid(ncid,'vovecrtz',vertid)
ncfile=nf_get_var_real(ncid,vertid,vert)
ncfile=nf_close(ncid)


!Converting the 5day-mean data into monthly-mean data
ind1=1
do tt=1,12
 montemp(:,:,:,tt)=sum(temp(:,:,:,ind1:ind1+5),dim=4)/float(6)
 monsali(:,:,:,tt)=sum(sali(:,:,:,ind1:ind1+5),dim=4)/float(6)
 monmeri(:,:,:,tt)=sum(meri(:,:,:,ind1:ind1+5),dim=4)/float(6)
 monvert(:,:,:,tt)=sum(vert(:,:,:,ind1:ind1+5),dim=4)/float(6)
 ind1=ind1+6
end do

!Creating netcdf data file for temperature
ncfile=nf_create('lab28_mm_tsvw.nc',nf_share,ncid)
ncfile=nf_def_dim(ncid,'x',jpi,dim1)
ncfile=nf_def_dim(ncid,'y',jpj,dim2)
ncfile=nf_def_dim(ncid,'deptht',jpk,dim3)
ncfile=nf_def_dim(ncid,'time_counter',nmonth,dim4)
dimtx(1)=dim1
dimtx(2)=dim2
dimtx(3)=dim3
dimtx(4)=dim4
ncfile=nf_def_var(ncid,'votemper',nf_real,4,dimtx,deftemp)
ncfile=nf_def_var(ncid,'vosaline',nf_real,4,dimtx,defsali)
ncfile=nf_def_var(ncid,'vomecrty',nf_real,4,dimtx,defmeri)
ncfile=nf_def_var(ncid,'vovecrtz',nf_real,4,dimtx,defvert)
ncfile=nf_enddef(ncid)
ncfile=nf_put_var_real(ncid,deftemp,montemp)
ncfile=nf_put_var_real(ncid,defsali,monsali)
ncfile=nf_put_var_real(ncid,defmeri,monmeri)
ncfile=nf_put_var_real(ncid,defvert,monvert)
ncfile=nf_close(ncid)

end program convertnetcdfdata

subroutine handle_err(errcode)
      implicit none
      include 'netcdf.inc'
      integer errcode

      print *, 'Error: ', nf_strerror(errcode)
      stop 2
      end