! ! PROGRAMA DE LEITURA de arquivos tipo ETAGRDhh.tmhh ! program leit ! use variables integer, parameter :: kmax=50 integer :: dt,kgtype,imdlty,igout,jgout,im,jm,NRSTRT,NPINCR real :: dx,lonw,alatvt,lats,dy character(len=20) :: field character(len=6) :: nufile,outype,proj,readco,readll,itag,dataset,datset character(len=12) :: compl logical :: ltsoil, lqsoil,north logical :: infloop ! ! Dimension 2D variables ! real,allocatable,dimension(:,:) :: lon !LONGITUDE real,allocatable,dimension(:,:) :: lat !LATITUDE real,allocatable,dimension(:,:) :: mask !LAND/SEA MASK real,allocatable,dimension(:,:) :: cssf !AVE SFC SENHEAT FX real,allocatable,dimension(:,:) :: cmsf !AVE SFC MOMENTUM FX real,allocatable,dimension(:,:) :: ps !SURFACE PRESSURE real,allocatable,dimension(:,:) :: ttprec !ACM TOTAL PRECIP real,allocatable,dimension(:,:) :: zorl !ROUGHNESS LENGTH real,allocatable,dimension(:,:) :: GRIB real,allocatable,dimension(:,:) :: tsh !SHELTER TEMPERATURE real,allocatable,dimension(:,:) :: zs !SURFACE HEIGHT real,allocatable,dimension(:,:) :: tdsh !SHELTER DEWPOINT real,allocatable,dimension(:,:) :: uanem !U WIND AT ANEMOM HT real,allocatable,dimension(:,:) :: vanem !V WIND AT ANEMOM HT ! ! ! Dimension 3D variables ! real,allocatable,dimension(:,:,:) :: phi !HEIGHT ON ETA SFCS real,allocatable,dimension(:,:,:) :: peta !PRESS ON ETA SFCS real,allocatable,dimension(:,:,:) :: q !SPEC HUM ON ETA SFCS real,allocatable,dimension(:,:,:) :: tpot !POT TEMP ON ETA SFCS real,allocatable,dimension(:,:,:) :: tke !TRBLNT KE ON ETA SFC real,allocatable,dimension(:,:,:) :: u !U WIND ON ETA SFCS real,allocatable,dimension(:,:,:) :: v !V WIND ON ETA SFCS infloop=.true. read(5,*)ITAG,NRSTRT,NPINCR,dataset !GSM dataset=TRIM(dataset) write(6,*)ITAG,NRSTRT,NPINCR,dataset ! READ OUTPUT GRID SPECIFICATIONS. ! OPEN(18,FILE='cntrl.parm_'//dataset,status='old') rewind(18) read(18,1000) kgtype read(18,1000) imdlty read(18,1030) datset read(18,1030) outype read(18,1030) nufile read(18,1030) proj read(18,1010) north read(18,1000) im read(18,1000) jm read(18,1020) dx read(18,1020) lonw read(18,1020) alatvt read(18,1020) lats read(18,1020) dy read(18,1030) readll read(18,1030) readco 1000 format(T28,I5) 1010 format(T28,L1) 1020 format(T28,F11.6) 1030 format(T28,A6) print*,"im: ",im," jm: ",jm,"kmax: ",kmax ! Allocate variables allocate(GRIB(im,jm)) allocate(lon(im,jm)) allocate(lat(im,jm)) allocate(tsh(im,jm)) allocate(tdsh(im,jm)) allocate(zs(im,jm)) allocate(uanem(im,jm)) allocate(vanem(im,jm)) allocate(mask(im,jm)) allocate(ps(im,jm)) allocate(ttprec(im,jm)) allocate(cssf(im,jm)) allocate(cmsf(im,jm)) allocate(zorl(im,jm)) ! Variaveis 3D allocate(peta(im,jm,kmax)) allocate(phi(im,jm,kmax)) allocate(tpot(im,jm,kmax)) allocate(q(im,jm,kmax)) allocate(u(im,jm,kmax)) allocate(v(im,jm,kmax)) allocate(tke(im,jm,kmax)) OPEN(unit=51,file=datset//itag//'.T00S',status='old',form='unformatted') ! get I and J indexes for the smaller domain ! IW= ABS((-83) - (-83))/ 0.40 + 1 ! IE= ABS((-83) - (-25.8))/ 0.40 + 1 ! JS= ABS((-50.2) - (-50.2))/ 0.40 + 1 ! JN= ABS((-50.2) - (12.2))/ 0.40 + 1 IW= 1 IE= im JS= 1 JN= jm ! ! get the nx and ny number of points in x and y direction nx = ie - iw + 1 ny = jn - js + 1 ! ! get the exactly lat/lon of the smaller domain ! clonll= -lonw + ((iw-1) * dx) clonur= -lonw + ((ie-1) * dx) clatll= lats + ((js-1) * dy) clatur= lats + ((jn-1) * dy) print*, iw, ie, js, jn print*, clonll, clonur, clatll, clatur READ(51) IHRST,IMM,IDD,IYY,IHH print*, ' reading header' READ(51) KGTYP,PROJ,NORTH,IMOUT,JMOUT,POLEI,POLEJ,ALATVT,ALONVT,XMESHL kqs=0 kts=0 kt=0 kq=0 kr=0 kp=0 kw=0 ku=0 kv=0 kc=0 ltsoil=.false. lqsoil=.false. do while (infloop) READ(51,END=999) FIELD,SFC print*, field, sfc READ(51) GRIB if (field.eq. 'LONGITUDE ') then; lon=grib elseif (field.eq. 'LATITUDE ') then; lat=grib elseif (field.eq. 'SURFACE HEIGHT ') then; zs=grib elseif (field.eq. 'SHELTER TEMPERATURE ') then; tsh=grib elseif (field.eq. 'SHELTER DEWPOINT ') then; tdsh=grib elseif (field.eq. 'U WIND AT ANEMOM HT ') then; uanem=grib elseif (field.eq. 'V WIND AT ANEMOM HT ') then; vanem=grib elseif (field.eq. 'LAND/SEA MASK ') then; mask=grib elseif (field.eq. 'SURFACE PRESSURE ') then; ps=grib elseif (field.eq.'ACM TOTAL PRECIP ') then; ttprec=grib/1000 elseif (field.eq.'AVE SFC SENHEAT FX ') then; cssf=grib elseif (field.eq.'AVE SFC MOMENTUM FX ') then; cmsf=grib elseif (field.eq.'ROUGHNESS LENGTH ') then; zorl=grib ! ! retrieving 3d fields ! elseif (field.eq.'PRESS ON ETA SFCS ') then; kw=kw+1; peta(1:im,1:jm,kw)=grib elseif (field.eq.'HEIGHT ON ETA SFCS ') then; kp=kp+1; phi(1:im,1:jm,kp)=grib elseif (field.eq.'POT TEMP ON ETA SFCS') then; kt=kt+1; tpot(1:im,1:jm,kt)=grib elseif (field.eq.'SPEC HUM ON ETA SFCS') then; kq=kq+1; q(1:im,1:jm,kq)=grib elseif (field.eq.'U WIND ON ETA SFCS ') then; ku=ku+1; u(1:im,1:jm,ku)=grib elseif (field.eq.'V WIND ON ETA SFCS ') then; kv=kv+1; v(1:im,1:jm,kv)=grib elseif (field.eq.'TRBLNT KE ON ETA SFC') then; kc=kc+1; tke(1:im,1:jm,kc)=grib endif enddo 999 continue ! ! ! save fields for Grads plotting. First 2D, then 3D fields ! OPEN(unit=17,file='latlon_'//itag,STATUS='UNKNOWN', & & form='unformatted',access='direct' , & & recl=4*nx*ny) irec=1 write(17,rec=irec)((lon(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((lat(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((zs(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((tsh(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((tdsh(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((uanem(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((vanem(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((mask(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((ps(i,j),i=iw,ie),j=js,jn) irec=irec+1 where(ttprec<=0.) ttprec=0. write(17,rec=irec)((ttprec(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((cssf(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((cmsf(i,j),i=iw,ie),j=js,jn) irec=irec+1 write(17,rec=irec)((zorl(i,j),i=iw,ie),j=js,jn) print*, ' 2D irec=',irec-1 irec=irec+1 do k=kmax,1,-1 write(17,rec=irec)((peta(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((phi(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((tpot(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((q(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((u(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((v(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo do k=kmax,1,-1 write(17,rec=irec)((tke(i,j,k),i=iw,ie),j=js,jn) irec=irec+1 enddo stop end program leit