program read_gfs implicit none character(len=255) :: gribfile,outdir !_____________________________________________________________________________________________________________ integer :: nx,ny,nz,np integer :: dx,dy integer :: proj,LatS,LonW,LatN,LonE integer :: i,j,k,l,n,n2 REAL, PARAMETER :: Undef=9.9990003E+20 real, dimension(:,:,:) ,allocatable :: ht !Isobaric heights (m) real, dimension(:,:,:) ,allocatable :: tp !Isobaric temps (K) real, dimension(:,:,:) ,allocatable :: th !Isobaric theta (K) real, dimension(:,:,:) ,allocatable :: uw !Isobaric u-wind (m/s) real, dimension(:,:,:) ,allocatable :: vw !Isobaric v-wind (m/s) real, dimension(:,:,:) ,allocatable :: rh !Isobaric rh,mr (%,kg/kg) real, dimension(:) ,allocatable :: pr,pri !Isobaric pressures (mb) !_____________________________________________________________________________________________________________ real, dimension(:,:,:) ,allocatable :: slp !slp(i,j,1)=LCMSK;slp(i,j,2)=PSLM !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC !slp(i,j,5 & 6) are 0-10 cm STC and SMC !slp(i,j,7 & 8) are 10-200 cm STC and SMC !_____________________________________________________________________________________________________________ character(len=15) :: ModelDriver !_____________________________________________________________________________________________________________ open(1,file='InputModelInf.txt',form='formatted',status='old') rewind 1 read(1,'(a15)')ModelDriver read(1,*)nx,ny,nz read(1,*)proj read(1,*)LatS read(1,*)LonW read(1,*)LatN read(1,*)LonE read(1,*)dx read(1,*)dy allocate(ht(nx,ny,nz)) allocate(tp(nx,ny,nz)) allocate(th(nx,ny,nz)) allocate(uw(nx,ny,nz)) allocate(vw(nx,ny,nz)) allocate(rh(nx,ny,nz)) allocate(slp(nx,ny,12)) allocate(pr(nz)) allocate(pri(nz)) read(5,'(a)') gribfile do k=1,nz read(1,*)pr(k) print*,pr(k) enddo close(1) ! *** Fill pressure levels. write(6,*) 'using hardwired stuff! ', nx,ny,nz write(6,*) (pr(k),k=1,nz) ! *** Read in degrib data. ! n=index(outdir,' ')-1 n2=index(ModelDriver,' ')-1 write(6,*)n2 write(6,*)ModelDriver(1:n2) print*, "MODELO UTILIZADO=",ModelDriver(1:n2) n=index(gribfile,' ')-1 write(*,*) 'OPEN FILE ',gribfile(1:n) open(12,file=gribfile(1:n),form='unformatted',status='old') !!! Dragan 12/07/2018 open(15,form='unformatted', access='direct',recl=1440*721*4) open(15,form='unformatted', access='direct',recl=nx*ny*4) do l=1,nz read(12) ht(:,:,l) print*,' ht l ',l,ht(481,241,l),ht(241,481,l) enddo do l=1,nz read(12) uw(:,:,l) print*,'uw l ',l,uw(481,241,l),ht(241,481,l) enddo do l=1,nz read(12) vw(:,:,l) print*,'vw l ',l, vw(481,241,l),ht(241,481,l) enddo do l=1,nz read(12) tp(:,:,l) print*,'tp l ',l, tp(481,241,l),ht(241,481,l) enddo do l=1,nz read(12) rh(:,:,l) print*,'rh l ',l, rh(481,241,l),ht(241,481,l) enddo ! Missing RH at 20hPa ! ! rh(:,:,nz)=-99999. ! rh(:,:,nz-1)=-99999. do l=1,12 if (l.eq.2.or.l.eq.3) then read(12) slp(:,:,l) slp(:,:,l)=slp(:,:,l)/100 else read(12) slp(:,:,l) endif print*,'slp l ',l, slp(481,241,l),ht(241,481,l) enddo close(12) !############################################################################################################# do i=1,nx do j=1,ny do l=6,12 if(l.eq.6.or.l.eq.8.or.l.eq.10.or.l.eq.12) then if (slp(i,j,l).eq.Undef) slp(i,j,l)=0 if (slp(i,j,l).gt.0.5) slp(i,j,l)=0.3 endif enddo enddo enddo !############################################################################################################# !------Over the Ocean, the value of Temp/Soil data are equal 9.9990003E+20 leaving shadows at the edges, !------so, to solve the problem we fill with neighbors value do l=5,11 if(l.eq.5.or.l.eq.7.or.l.eq.9.or.l.eq.11) then do np=1,200 do I=1,nx do J=1,ny if ( slp(I,J,l).eq.Undef.and.slp(I-1,J,l).ne.Undef) then slp(I,J,l)=slp(I-1,J,l) endif enddo enddo do I=nx-1, 1, -1 do J=1,ny if ( slp(I,J,l).eq.Undef.and.slp(I+1,J,l).ne.Undef) then slp(I,J,l)=slp(I+1,J,l) endif enddo enddo enddo do np=1,200 do I=1,nx do J=1,ny if ( slp(I,J,l).eq.Undef.and.slp(I,J-1,l).ne.Undef) then slp(I,J,l)=slp(I,J-1,l) endif enddo enddo do I=1,nx do J=ny-1, 1, -1 if ( slp(I,J,l).eq.Undef.and.slp(I,J+1,l).ne.Undef) then slp(I,J,l)=slp(I,J+1,l) endif enddo enddo enddo endif enddo !!################################################## WRITE ################################################## do l=1,nz write(15) uw(:,:,l) print*,' uww l ',l,uw(481,241,l),uw(241,481,l) enddo do l=1,nz write(15) vw(:,:,l) print*,' vww l ',l,vw(481,241,l),vw(241,481,l) enddo do l=1,nz write(15) ht(:,:,l) print*,' htw l ',l,ht(481,241,l),ht(241,481,l) enddo do l=1,nz write(15) tp(:,:,l) print*,' tpw l ',l,tp(481,241,l),tp(241,481,l) enddo do l=1,nz write(15) rh(:,:,l) print*,' rhw l ',l,rh(481,241,l),rh(241,481,l) enddo do l=1,12 write(15) slp(:,:,l) print*,' slpw l ',l,slp(481,241,l),slp(241,481,l) enddo end program read_gfs