program degrib2model_driver c implicit none c integer nx,ny,nz,kpds(200),datsav(9),n,k !Degrib grid dimensions c character(LEN=255):: gribfile,outfile c real esat,es c common /estab/esat(15000:45000),es(15000:45000) c_______________________________________________________________________________ c c *** Initialize tables. c call es_ini c c *** Read file names. c read(5,'(a)') gribfile read(5,'(a)') outfile c! gribfile='../data_in/init/gfsdata'//c_k//'.grb' c! outfile ='../data_in/init/gfsdata'//c_k//'.dat' c n=index(gribfile,' ')-1 print *,n write(6,*) 'first call of GET_GDS' c CALL GET_GDS(gribfile(1:n),datsav,kpds)!!!!!!!!!!!!!!!!!!! DRAGAN c nx=datsav(1) !!!!!!!!!!!!!!! c ny=datsav(2) !!!!!!!!!!!!!!!!!! cGSM nx=360 !!!!!!!!!!!!!!!!!! cGSM ny=181 !!!!!!!!!!!!!!!!!! nx=1440 !!!!!!!!!!!!!!!!!! ny=721 !!!!!!!!!!!!!!!!!! if (nx .eq. 360 .and. ny .eq. 181) then nz=26 elseif (nx .eq. 144 .and. ny .eq. 73) then nz=17 elseif (nx .eq. 1440 .and. ny .eq. 721) then nz=26 else write(6,*) 'what grid is this??? ', nx,ny write(6,*) ' ' write(6,*) 'I do not know the vertical levels for this ' write(6,*) 'grid...please add some code to dgeta2model.f' write(6,*) ' ' STOP 999 endif print *,"nx,ny,nz,",nx,ny,nz,gribfile,outfile c call degrib2model(gribfile,outfile,nx,ny,nz) c end c c=============================================================================== c subroutine degrib2model(gribfile,outfile,nx,ny,nz) c implicit none c real cp,kappa,dlat,dlon,sum,avg parameter(cp=1004.,kappa=287./1004.) c integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld . ,iyear,imonth,iday,ifcsthr . ,count,kgds(200),ireto,n1,n2,n3,LenStr c cGSM real ht(nx,ny,nz) !Isobaric heights (m) cGSM . ,tp(nx,ny,nz) !Isobaric temps (K) cGSM . ,th(nx,ny,nz) !Isobaric theta (K) cGSM . ,uw(nx,ny,nz) !Isobaric u-wind (m/s) cGSM . ,vw(nx,ny,nz) !Isobaric v-wind (m/s) cGSM . ,rh(nx,ny,nz) !Isobaric rh,mr (%,kg/kg) cGSM . ,pr(nz),pri(nz) !Isobaric pressures (mb) real lat1,lat2,lon0,sw(2),ne(2) . ,xe,mrsat,esat,es . ,avnlevs(26),reanlevs(17) cGSM . ,slp(nx,ny,14) !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL !@@@ . ,temp(nx,ny) !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 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 REAL, ALLOCATABLE:: htout(:,:,:), + uwout(:,:,:),vwout(:,:,:),rhout(:,:,:), + slpout(:,:,:) c character*255 gribfile,outdir,outfile,gdsfile character*2 gproj character*11 atime character*8 atime_r character*7 model character*5 timex data AVNLEVS/1000,975,950,925,900,850,800,750,700,650, + 600,550,500,450,400,350,300,250,200,150, + 100, 70, 50, 30, 20, 10/ data REANLEVS/1000,925,850,700,600,500, + 400,300,250,200,150,100, + 70, 50, 30, 20, 10/ c common /estab/esat(15000:45000),es(15000:45000) 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,14)) allocate(pr(nz)) allocate(pri(nz)) c_______________________________________________________________________________ c c *** Fill pressure levels. c Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... print *,'degrib2model...' nsfcfld=14 c nsfcfld=12 !@@@ Cmp write(6,*) 'using hardwired stuff! ', nx,ny,nz do k=1,nz if (nz .eq. 26) then pr(k)=AVNLEVS(k) elseif (nz .eq. 17) then pr(k)=reanlevs(k) endif enddo write(6,*) 'return read_degrib' cGSM OPEN(11,file='../INIT_DATA/feb1996/1996020900.bin', cGSM & form='unformatted', cGSM & access='direct',status='old',recl=360*181*4) n=index(gribfile,' ')-1 open(11,file=gribfile(1:n),form='unformatted',status='old', & access='direct',recl=1440*721*4) write(0,*) 'OPEN FILE ',gribfile(1:n) c############################################# !GSM NOVA LEITURA PARA OS ARQUIVOS GFS !GSM do l=1,nz !GSM read(11) ht(:,:,l) !GSM print*,' ht l ',l,ht(241,521,l) !GSM enddo !GSM do l=1,nz !GSM read(11) uw(:,:,l) !GSM print*,'uw l ',l,uw(241,521,l) !GSM enddo !GSM do l=1,nz !GSM read(11) vw(:,:,l) !GSM print*,'vw l ',l, vw(241,521,l) !GSM enddo !GSM do l=1,nz !GSM read(11) tp(:,:,l) !GSM print*,'tp l ',l, tp(241,521,l) !GSM enddo !GSM do l=1,nz !GSM read(11) rh(:,:,l) !GSM print*,'rh l ',l, rh(241,521,l) !GSM enddo ! Missing RH at 20hPa ! !GSM rh(:,:,nz)=-99999. !GSM rh(:,:,nz-1)=-99999. !GSM do l=1,12 !GSM if (l.eq.2.or.l.eq.3) then !GSM !GSM read(11) slp(:,:,l) !GSM slp(:,:,l)=slp(:,:,l)/100 !GSM else !GSM read(11) slp(:,:,l) !GSM endif !GSM print*,'slp l ',l, slp(241,521,l) !GSM enddo !GSM close(11) c############################################# do 20 l=1,nz k=l read(11,rec=k) ((uw(i,j,l),i=1,nx),j=ny,1,-1) c if (l .eq. 20) then c print*,'uw=',((uw(1,1,l))) c endif 20 continue c print*,'uw=',uw(1,1,20) do 30 l=1,nz k=l+26 read(11,rec=k) ((vw(i,j,l),i=1,nx),j=ny,1,-1) 30 continue c print*,'vw=',vw(1,1,20) do 40 l=1,nz k=l+52 read(11,rec=k) ((ht(i,j,l),i=1,nx),j=ny,1,-1) 40 continue c print*,'ht=',ht(1,1,20) do 50 l=1,nz k=l+78 read(11,rec=k) ((tp(i,j,l),i=1,nx),j=ny,1,-1) 50 continue c print*,'tp=',tp(1,1,20) do 60 l=1,nz k=l+104 read(11,rec=k) ((rh(i,j,l),i=1,nx),j=ny,1,-1) 60 continue c print*,'rh=',rh(1,1,20) read(11,rec=131) ((slp(i,j,1),i=1,nx),j=ny,1,-1) read(11,rec=132) ((slp(i,j,2),i=1,nx),j=ny,1,-1) read(11,rec=133) ((slp(i,j,3),i=1,nx),j=ny,1,-1) read(11,rec=134) ((slp(i,j,4),i=1,nx),j=ny,1,-1) read(11,rec=135) ((slp(i,j,5),i=1,nx),j=ny,1,-1) read(11,rec=136) ((slp(i,j,6),i=1,nx),j=ny,1,-1) read(11,rec=137) ((slp(i,j,7),i=1,nx),j=ny,1,-1) read(11,rec=138) ((slp(i,j,8),i=1,nx),j=ny,1,-1) read(11,rec=139) ((slp(i,j,9),i=1,nx),j=ny,1,-1) read(11,rec=140) ((slp(i,j,10),i=1,nx),j=ny,1,-1) read(11,rec=141) ((slp(i,j,11),i=1,nx),j=ny,1,-1) read(11,rec=142) ((slp(i,j,12),i=1,nx),j=ny,1,-1) ccc read(11,rec=143) ((slp(i,j,13),i=1,nx),j=ny,1,-1) ccc read(11,rec=144) ((slp(i,j,14),i=1,nx),j=ny,1,-1) OPEN(77,file='../INIT_DATA/feb1998/ICEC_WEASD.bin', & form='unformatted',access='direct',status='old',recl=360*181*4) read(77,rec=1) ((slp(i,j,13),i=1,360),j=181,1,-1) read(77,rec=2) ((slp(i,j,14),i=1,360),j=181,1,-1) do i=1,360 do j=1,181 if (slp(i,j,13).eq.9.9990003E+20) then slp(i,j,13)=0. endif if (slp(i,j,14).eq.9.9990003E+20) then slp(i,j,14)=0. endif enddo enddo do i=1,360 do j=1,181 ! print*,slp(i,j,13) enddo enddo c######################################### print *,'SLP',1,slp(234,70,1) print *,'SLP',2,slp(234,70,2) print *,'SLP',3,slp(234,70,3) print *,'SLP',4,slp(234,70,4) print *,'SLP',5,slp(234,70,5) print *,'SLP',6,slp(234,70,6) print *,'SLP',7,slp(234,70,7) print *,'SLP',8,slp(234,70,8) print *,'SLP',9,slp(234,70,9) print *,'SLP',10,slp(234,70,10) print *,'SLP',11,slp(234,70,11) print *,'SLP',12,slp(234,70,12) print *,'SLP',13,slp(234,70,13) print *,'SLP',14,slp(234,70,14) do l=5,12 c! print *,'SLP',14,slp(:,:,l) enddo c print*,slp(:,:,5) c print*,rh(22,:,23) close(11) close(77) ccc close(12) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! print *,'rh(1,1,1),',rh(1,1,1) c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !GSM !GSM do i=1,nx !GSM do j=1,ny !GSM if (slp(i,j,1)<0.1) then !GSM slp(i,j,1)=0 !GSM else !GSM slp(i,j,1)=1 !GSM endif !GSM enddo !GSM enddo do i=1,nx do j=ny,1,-1 do l=1,nz if (ht(i,j,l).eq.9.9990003E+20) ht(i,j,l)=-99999. if (uw(i,j,l).eq.9.9990003E+20) uw(i,j,l)=-99999. if (vw(i,j,l).eq.9.9990003E+20) vw(i,j,l)=-99999. if (tp(i,j,l).eq.9.9990003E+20) tp(i,j,l)=-99999. if (rh(i,j,l).eq.9.9990003E+20) rh(i,j,l)=-99999. enddo !GSM do l=1,12 if (slp(i,j,l).eq.9.9990003E+20) slp(i,j,l)=0. enddo !GSM c############################################## DRAGAN do l=1,21 if (rh(i,j,l).gt.100) rh(i,j,l)=100 enddo do l=1,nz if (rh(i,j,l).lt.0) rh(i,j,l)=1 enddo enddo enddo c########################################################################### print*,slp(360,:,3) c c *** Check for any missing data. c do k=1,nz if (ht(1,1,k) .eq. -99999.) then print *,'Height data missing at level: ',pr(k) stop elseif (tp(1,1,k) .eq. -99999.) then print *,'Temperature data missing at level: ',pr(k) stop elseif (rh(1,1,k) .eq. -99999.) then print *,'RH data missing at level: ',k,pr(k) print *,'Calling RH patch.' call rh_fix(nx,ny,nz,rh,pr) elseif (uw(1,1,k) .eq. -99999.) then print *,'U-wind data missing at level: ',pr(k),uw(1,1,k) stop elseif (vw(1,1,k) .eq. -99999.) then print *,'V-wind data missing at level: ',pr(k) stop endif enddo do k=1,nsfcfld if (slp(1,1,k) .eq. 9.9990003E+20) then print *,'SLP data missing at level: ',pr(k) endif enddo Cmp Cmp Handle missing surface data.... Cmp C CFor reanalysis runs these data will not be used. Just bogus Cto hold a place. C write(6,*) 'going to use nsfcfld= ', nsfcfld CIf rean put bogus data into slp if (nz .eq. 17) then do K=5,8 do j=1,ny do i=1,nx if (mod(k,2) .eq. 0.and.k.ne.2) then !!!!!!!!!!!????????????????? slp(i,j,k)=0.14 else slp(i,j,k)=275. endif enddo enddo enddo endif !!! do k=5,nsfcfld do k=5,12 if (slp(1,1,k) .eq. -99999.) then write(6,*) 'SURFACE DATA MISSING... ', K if (k.eq.9.or.k.eq.11) then write(6,*) 'filling soil temp data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,7) enddo enddo elseif(k.eq.10.or.k.eq.12.) then write(6,*) 'filling soil moisture data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,8) enddo enddo endif endif Cmp if (mod(k,2) .eq. 0.and.k.ne.2) then !!!!!!!!!!!!!!!! do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=0.14 endif enddo enddo else do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=273.15 endif enddo enddo endif Cmp enddo c c *** Convert 3d temp to theta. c *** Compute Exner function. c *** Convert 3d rh to mr. c do k=1,nz pri(k)=1./pr(k) enddo c do k=1,nz do j=1,ny do i=1,nx th(i,j,k)=tp(i,j,k)*(1000.*pri(k))**kappa C ex(i,j,k)=cp*tp(i,j,k)/th(i,j,k) it=tp(i,j,k)*100 it=min(45000,max(15000,it)) xe=esat(it) mrsat=0.00622*xe/(pr(k)-xe) rh(i,j,k)=rh(i,j,k)*mrsat enddo enddo enddo Cmp Cmp for some reason AVN data appears to be "upside down" in the N-S Cmp sense. Flip the arrays of data. If reanalysis data, skip this part Cmp c! if (nz .eq. 17) goto 969 do k=1,nz do j=1,ny do i=1,nx temp(i,j)=th(i,j,k) enddo enddo do J=1,ny do i=1,nx th(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=uw(i,j,k) enddo enddo do J=1,ny do i=1,nx uw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=vw(i,j,k) enddo enddo do J=1,ny do i=1,nx vw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=ht(i,j,k) enddo enddo if (k.eq.1) write(6,*) 'flipping isobaric data - ' do J=1,ny do i=1,nx ht(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=rh(i,j,k) enddo enddo do J=1,ny do i=1,nx rh(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx C temp(i,j)=ex(i,j,k) enddo enddo do J=1,ny do i=1,nx C ex(i,j,k)=temp(i,ny-J+1) enddo enddo C----------------------------------------- enddo Cmp generalize eventually! ??????????????????????????????????????????????????? do k=1,nsfcfld C write(6,*) 'flipping surface fields... ' do j=1,ny do i=1,nx temp(i,j)=slp(i,j,k) enddo enddo do J=1,ny do i=1,nx slp(i,j,k)=temp(i,ny-J+1) enddo enddo enddo CC CC CC Add bogus column if reanalysis data CC CC if (nx .eq. 144) then ALLOCATE(HTOUT(145,ny,nz),UWOUT(145,ny,nz),VWOUT(145,ny,nz)) ALLOCATE(RHOUT(145,ny,nz),SLPOUT(145,ny,14)) DO J=1,NY DO I=1,NX DO K=1,NSFCFLD SLPOUT(I,J,K)=SLP(I,J,K) ENDDO DO K=1,NZ HTOUT(I,J,K)=HT(I,J,K) UWOUT(I,J,K)=UW(I,J,K) VWOUT(I,J,K)=VW(I,J,K) RHOUT(I,J,K)=RH(I,J,K) ENDDO ENDDO DO K=1,NSFCFLD SLPOUT(145,J,K)=SLP(1,J,K) ENDDO DO K=1,NZ HTOUT(145,J,K)=HT(1,J,K) UWOUT(145,J,K)=UW(1,J,K) VWOUT(145,J,K)=VW(1,J,K) RHOUT(145,J,K)=RH(1,J,K) ENDDO ENDDO endif c c c *** Create output file name. c n=index(outfile,' ')-1 Cmp make ETA_avn to differentiate from files based on eta grib input if (nz .eq. 26) then model='ETA_avn' elseif (nz .eq. 17) then model='ETA_rnl' nx=145 endif n=index(outfile,' ')-1 cccccccccccccc outfile=outdir(1:n)//'/'//atime(1:LenStr)//'.'//model print *,model,' data ---> ',outfile(1:n) !@@@ open(1,file=outfile(1:n),status='unknown',form='unformatted') c c *** Write header stuff. c ip=1 jp=1 Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... c nsfcfld=12 !@@@ nsfcfld=14 gproj='LL' write(1) nz C new stuff for GDS file C C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C new stuff for GDS file C n=index(gribfile,' ')-1 write(6,*) 'calling get_fullgds with ', + gribfile(1:n) c call get_fullgds(gribfile(1:n),nx,ny,kgds) cGSM kgds(2)=360 cGSM kgds(3)=181 kgds(2)=1440 kgds(3)=721 kgds(4)=-90000 kgds(6)=128 kgds(7)=90000 kgds(8)=-1000 kgds(9)=1000 kgds(10)=1000 kgds(20)=255 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! write(6,*) 'back from get_fullgds' c call baclose(11,ireto) n=index(outfile,' ')-1 gdsfile='gdsinfo.ETA_avn' n=index(gdsfile,' ')-1 write(6,*) 'GDS written to ', gdsfile(1:n) open(unit=14,file=gdsfile(1:n),form='unformatted', + access='sequential') C C MODIFY GDS TO REFLECT FLIPPING OF DATA IN N-S SENSE C if (nz .eq. 26) then kgds(4)= -kgds(4) kgds(7)= -kgds(7) endif write(6,*) 'writing kgds ', (kgds(I),i=1,14) write(14) kgds close(14) ! write(6,*) 'past 14 close' C C end new stuff c c *** Write isobaric upper air data. c write(6,*) 'starting writes to unit1' if (nx .ne. 145) then write(1) uw write(1) vw write(1) ht write(1) rh write(1) pr write(1) slp else WRITE(1) UWOUT WRITE(1) VWOUT WRITE(1) HTOUT WRITE(1) RHOUT WRITE(1) PR WRITE(1) SLPOUT endif if (allocated(uwout)) deallocate(uwout) if (allocated(vwout)) deallocate(vwout) if (allocated(htout)) deallocate(htout) if (allocated(rhout)) deallocate(rhout) if (allocated(slpout)) deallocate(slpout) Cmp write(6,*) 'writing these sfc values to the ETA_avn file...' do k=1,nsfcfld write(6,*) 'field, value ', k, slp(1,1,k) enddo cGSM do k=1,nsfcfld cGSM write(6,*) 'field, value ', k, slp(81,261,k) cGSM enddo c close(1) c return end c c=============================================================================== c subroutine rh_fix(nx,ny,nz,rh,pr) c implicit none c integer*4 nx,ny,nz,i,j,k,kk c real*4 rh(nx,ny,nz),pr(nz) c_______________________________________________________________________________ c c *** Fix bottom levels if necessary. c if (rh(1,1,1) .eq. -99999.) then do k=2,nz do j=1,ny do i=1,nx ! print *,i,j,k,rh(i,j,k) enddo enddo if (rh(1,1,k) .ne. -99999.) then do kk=k-1,1,-1 print *,'Copying',nint(pr(kk+1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk+1) enddo enddo enddo goto 10 endif enddo print *,'RH patch did not work.' stop endif c c *** Fix upper levels if necessary. c 10 continue if (rh(1,1,nz) .eq. -99999.) then do k=nz-1,1,-1 if (rh(1,1,k) .ne. -99999.) then do kk=k+1,nz print *,'Copying',nint(pr(kk-1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk-1) enddo enddo enddo goto 20 endif enddo endif c 20 continue do k=1,nz if (rh(1,1,k) .eq. -99999.) then print *,'RH patch did not work.' stop endif enddo c return end c c=============================================================================== c subroutine es_ini c common /estab/esat(15000:45000),es(15000:45000) c c *** Create tables of the saturation vapour pressure with up to c two decimal figures of accuraccy: c do it=15000,45000 t=it*0.01 p1 = 11.344-0.0303998*t p2 = 3.49149-1302.8844/t c1 = 23.832241-5.02808*alog10(t) esat(it) = 10.**(c1-1.3816E-7*10.**p1+ . 8.1328E-3*10.**p2-2949.076/t) es(it) = 610.78*exp(17.269*(t-273.16)/(t-35.86)) enddo c return end c c===============================================================================