داده های میانگین ماهانه

سلام
آیا هیچ فایل توضیحی برای داده‌ها وجود داره؟ منظورم اینه که میدونید با چه ساختاری نوشته شدن؟ سوال بعد اینکه فرمت اونا ascii (اسکی) هست یا binari (دودویی)؟

به طور مثال، در برنامه‌ای یکی از روش‌های استخراج ماههای  موردنظر را نوشته ام.

http://paste.starnix.ir/phprmwdtn

http://paste.starnix.ir/p9gkxa4lb

در برنامه اول (test-w.f95) برای 36 ماه (3سال) اعداد 1 تا 36 در نظر گرفته و در فایل باینری به اسم mean.dat به صورت خطی و یک ترازه ذخیره می شود.

در برنامه دوم (test.f95) همین فایل، خوانده شده و فقط سه ماه تابستان (tir و mordad و shahrivar) استخراج می‌گردد.

program testo
INTEGER,PARAMETER :: y=3,n=y*12

REAL*4,DIMENSION(n) :: field_in

OPEN(unit=300, file='mean.dat', form='unformatted', status='unknown', access='direct', recl=4*n)

do i=1,n
field_in(i)=i
end do

WRITE(300,rec=1) (field_in(i),i=1,n)

end program

program testo
INTEGER,PARAMETER :: y=3, n=y*12
REAL*4,DIMENSION(n) :: field_in
REAL*4,DIMENSION(y) :: tir,mordad, shahrivar

OPEN(unit=300, file='mean.dat', form='unformatted', status='unknown', access='direct', recl=4*n)

 READ(300,rec=1) field_in
j=0		!( j+1) reffers to number of year 

do i=1,n

if (i == (4+j*12) ) tir(j+1)=field_in(i)	! for example: for  ( j=0)  =>  (i==4)  =>  tir(1)=field_in(4)  _____ 4th month of first year will be choosed and so on...

if (i == (5+j*12)) mordad(j+1)=field_in(i)

if (i == (6+j*12)) then 
shahrivar(j+1)=field_in(i)
j=j+1
end if 

end do

print*,'tir=', tir
print*,'mordad=', mordad
print*,'shahrivar=', shahrivar
end program

اگر ذخیره داده‌های شما دو بعدی باشه با ساختاری مثل لینک زیر (مقیاس قایم بر حسب سال و مقیاس افقی بر حسب ماه)، کار جداسازی داده‌ها بدون استفاده از شرط و ساده‌تر هم می‌شود

http://cri.ac.ir/files/Data/Synoptic%20Data%20up%20to%20l2005/ESFAHAN

دو برنامه زیر داده‌ها با ساختار دو بعدی را پردازش می‌کنند:

program testo
INTEGER,PARAMETER :: y=3

REAL*4,DIMENSION(12,y) :: field_in

OPEN(unit=300, file='mean.dat', form='unformatted', status='unknown', access='direct', recl=4*12*y)
n=0
do j=1,y
  do i=1,12
  n=n+1
field_in(i,j)=n
  end do
end do

WRITE(300,rec=1) ((field_in(i,j),i=1,12),j=1,y)

end program

program testo
INTEGER,PARAMETER :: y=3 
REAL*4,DIMENSION(12,y) :: field_in
REAL*4,DIMENSION(y) ::tir , mordad,shahrivar

OPEN(unit=300, file='mean.dat', form='unformatted', status='unknown', access='direct', recl=4*12*y)

 READ(300,rec=1) field_in

 
 do j=1,y

 write(*,*) (field_in(i,j),i=1,12)

 end do

 do i=1,y
 tir(i)=field_in(4,i)		! data for TIR will be choosed
 mordad(i)=field_in(5,i)
 shahrivar(i)=field_in(6,i)
 end do

print*,'tir=', tir
print*,'mordad=', mordad
print*,'shahrivar=', shahrivar
end program

مشکل در واقع این هست که ما اندازه درست ابعاد آرایه داده‌ها رو نمیدونیم.

شاید بشه با  allocate کردن بفهمیم ابعاد واقعیش چقدره. قدر مسلم اینه که محور x و y  نماینده طول و عرض جغرافیایی هستند.

من با این فایل کنترلی توی گردس می‌تونم داده‌ها رو بخونم، اما هنوز نمیفهمم ابعاد آرایه رو برای فرترن چطور تنظیم کنم.

DSET ^hgt_500_1948_2005.dat
TITLE ghashang
undef 99999
XDEF 144 linear 0 2.5
YDEF 73 linear -90 2.5 
ZDEF 1 linear 1 1
TDEF 684 LINEAR  1jan1948 1mo 
VARS 1
z 1 99 hgt
ENDVARS

خودم با یه برنامه مثل این، داده‌هایی با ساختار مشابه رو خونده‌ام(کل سطح زمین بود در 23 تراز قائم و 25 زمان). اما واسه داده‌های شما هنوز جواب نداده .

program testo
        
        INTEGER,PARAMETER :: nx=144, ny=73, nz=1, y=57, nt=y*12
        INTEGER :: i,j,k,irec,tt,irec1,t
        REAL*4,DIMENSION(nx,ny) :: field_in
        REAL*8,DIMENSION(nx,ny,nz,nt) :: h
        REAL*4,DIMENSION(y):: Hjun500,Hjul500,Haug500
        
OPEN(unit=100, file='hgt_500_1948_2005.dat', form='unformatted', status='unknown', access='direct', recl=4*nx*ny*nz*nt)

  
        irec=1
        irec1=1

        DO tt=1,nt
              write(*,*)'irec=',irec
              print*,'nt=',tt
             DO k=1,nz
                 READ(100,rec=irec) field_in

                 DO j=1,ny
                    DO i=1,nx
                       h(i,j,k,tt)=field_in(i,j)
                    ENDDO                    
                 ENDDO
                 irec=irec+1
             ENDDO
             write(*,*) i,j,k,tt,h(i,j,k,tt)
	ENDDO
      

      WRITE(*,*) "Enter specific LONGITUDE (0.0E to 357.5E : 2.5 degree)"
      READ (*,*) long
      WRITE(*,*) "Enter specific LATITUDE (90.0S {0} to 90.0N {180} : 2.5 degree)"
      READ (*,*) lat
!       WRITE(*,*) "Enter specific LEVEL ....)"					if there would be some pressure levels with respet to your data 
!       READ (*,*) lev    
 lev=1


 j=0										!hereafter ( j+1) reffers to number of year 

DO  t=1,nt

IF (t == (6+j*12) ) Hjun500(j+1)=h(long,lat,lev,t)		

IF (t == (7+j*12)) Hjul500(j+1)=h(long,lat,lev,t)

IF (t == (8+j*12)) THEN 
Haug500(j+1)=h(long,lat,lev,t)
j=j+1
END IF

ENDDO

end program

خداروشکر بلاخره درست شد. مشکل این بود که ما فکر میکردیم باید مثل یک فرمت استاندارد داده‌های چند ترازه خونده بشه. درحالی که این خودش خروجی یک برنامه دیگه بوده و باید در نظر میگرفتیم که نویسنده اون برنامه احتمالا بُعد زمان رو به عنوان بعد سوم در راستای z در نظر گرفته و همه داده‌ها رو در یک رکورد ذخیره کرده. برای همین بود که ما نمیتونستیم بیشتر از یک رکورد رو بخونیم. 1 نکته برنامه نویسی: خوبه زمانی که در یک حلقه قراره داده‌های زیادی خونده بشه، چیزی مثل ضربان قلب برای چک کردن این که حلقه درست کار میکنه و اگر اشتباهه تا کجا پیش رفته، قرار بدیم (مثل خط 17 برنامه ای که اینجا گذاشتم).

مورد بعدی، مسأله جدا کردن یک طول و عرض خاص هست. در جغرافیا طول و عرض از 0 شروع میشه درحالی که در برنامه ما با شماره یک مولفه آرایه‌ها سروکار داریم که از 1 شروع میشن. پس با تقسیم بر 2.5 که در واقع همون تفکیک فضایی داده‌هاست و اضافه کردن عدد 1، میتونیم شماره آرایه متناسب با طول و عرض مورد نظر رو داشته باشیم. اما از اونجایی که ممکنه طول و عرض مورد نظر مضرب صحیحی از 2.5 نباشه این محاسبه منجر به یک عدد اعشاری میشه که برای شماره یک آرایه درست نیست. پس میتونیم با اضافه کردن یک شرط عدد بدست اومده رو به سمت بالا یا به سمت پایین به صورت یک عدد صحیح گرد کنیم.

program testo
        
        INTEGER,PARAMETER :: nx=144, ny=73, y=57, nt=y*12
        INTEGER :: i,j,k,t
        REAL*4,DIMENSION(nx,ny,nt) :: h
        REAL*8,DIMENSION(y):: Hjun500,Hjul500,Haug500			! u dont know about accuracy of data, so you'd better read and write them with low precission but calculate with double precission
        REAL :: lo,la
OPEN(unit=100, file='hgt_500_1948_2005.dat', form='unformatted', status='unknown', access='direct', recl=4*nx*ny*nt)

!-------------------------------------------READING DATA------------------------------------------

        DO t=1,nt
              print*,'nt=',t

                 READ(100,rec=1) h

       write(*,*) t,h(1,1,t)					! this action works as heart of your program and u are able to check the reading process. u can command this line
	ENDDO
!------------------------------------------CHOOSING A PARTICULAR CORDINATE---------

      WRITE(*,*) "Enter specific LONGITUDE (0.0E to 357.5E : 2.5 degree)"
      READ (*,*) lo
      WRITE(*,*) "Enter specific LATITUDE (90.0S {-} to 90.0N {+} : 2.5 degree)"
      READ (*,*) la
	la=la+90							! in order to transfer (- and +) to {0 to 180} 
	
	IF (AMOD(lo,2.5) > 1.25) THEN
	lon=INT(lo/2.5)+1					!rounded to lower
	ELSE
	lon=INT(lo/2.5)+2					!rounded to upper
	END IF
	
	IF (AMOD(la,2.5) > 1.25) THEN
	lat=INT(la/2.5)+1					!rounded to lower
	ELSE
	lat=INT(la/2.5)+2					!rounded to upper
	END IF
 !--------------------------------------------SELECTING SPECIAL MONTHS------------------------
 
 j=0										!hereafter ( j+1) reffers to number of year 

DO  t=1,nt

IF (t == (6+j*12) ) Hjun500(j+1)=h(lon,lat,t)		

IF (t == (7+j*12)) Hjul500(j+1)=h(lon,lat,t)

IF (t == (8+j*12)) THEN 
Haug500(j+1)=h(long,lat,t)
j=j+1
END IF

ENDDO
!---------------------------------------------------------------------------------------------------------------
end program

امیدوارم موفق باشید

راستی شما اگر محاسبه خاصی (مثلا اعمال فیلتر بر داده‌هاتون) ندارید میتونید داده‌ها رو با گردس هم از طریق اسکریپ‌نویسی پردازش کنید. اینطوری به نظر ملموس‌تر هم هست چون میتونید نتیجه رو بلافاصله رسم و مشاهده کنید. در این زمینه میتونید برای کاربست توابع و برنامه‌های آماده گردس، دو وب‌سایت زیر رو بررسی کنید:

http://www.met.wau.nl/education/atd/Practical/scripts/library_local.html

http://www.met.wau.nl/education/atd/Practical/functions/library.html

در خوندن داده‌ها هیچ اشکالی وجود نداره. می‌تونید این رو با اضافه کردن یک خط دستور مثل

write(*,*) Haug500(j+1)

در شرط مربوط به جداسازی ماه آگوست چک کنید.

قسمت محاسبه میانگین ماه‌ها رو به صورت زیر اصلاح میکنم (شما فراموش کردید اینا رو در یک حلقه محاسبه کنید، یک سری چیزای زاید رو پاک کردم و حلقه‌ای که ایجاد کردم رو با حروف کوچک نوشتم)

Hmjun500=0.0

Hmjul500=0.0

Hmaug500=0.0

do j=1,y

Hmjun500=Hmjun500+Hjun500(j)

Hmjul500=Hmjul500+Hjul500(j)

Hmaug500=Hmaug500+Haug500(j)
 
 end do



Hmjun500=Hmjun500/57

Hmjul500=Hmjul500/57

Hmaug500=Hmaug500/57

WRITE(*,*) 'Hmjun',Hmjun500

WRITE(*,*) 'Hmjul',Hmjul500

WRITE(*,*) 'Hmaug',Hmaug500

ولی از محاسبات آنومالی نفهمیدم قصد دارید با چه ضابطه‌ای اون رو حساب کنید ولی احتمالا اونم به یک حلقه نیاز داره