program degrib2model_driver c implicit none c integer nx,ny,nz !Degrib grid dimensions c character(LEN=255):: degrib_dir,degrib_file,outdir,gribfile c real esat,es integer datsav,n,l,kpds(200) c common /estab/esat(15000:45000),es(15000:45000) common /gdsinfo/datsav(11) c_______________________________________________________________________________ c c *** Initialize tables. c call es_ini c c *** Read file names. c read(5,'(a)') gribfile read(5,'(a)') outdir n=index(gribfile,' ')-1 write(6,*) 'GRIBFILE is ', gribfile(1:n) CALL GET_GDS(gribfile(1:n),datsav,kpds) c C X and Y dimensions now obtained from the degribbed GRIB file c nx=datsav(1) ny=datsav(2) nz=19 c write(6,*) 'calling degrib2model ', nx,ny,nz call degrib2model(gribfile,outdir,nx,ny,nz) c end c c=============================================================================== c subroutine degrib2model(gribfile,outdir,nx,ny,nz) c implicit none c !Niv logical contAss !Se for falso, nao ha assimilacao... !Niv real cp,kappa parameter(cp=1004.,kappa=287./1004.) c integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld . ,iyear,imonth,iday,ifcsthr . ,datsav,ii,jj integer kgds(200) real crot(nx*ny),srot(nx*ny),tbar,old,fact,biassum real z1000t,z975t c real ht(nx,ny,nz) !Isobaric heights (m) . ,tp(nx,ny,nz) !Isobaric temps (K) . ,th(nx,ny,nz) !Isobaric theta (K) . ,uw(nx,ny,nz) !Isobaric u-wind (m/s) . ,vw(nx,ny,nz) !Isobaric v-wind (m/s) . ,rh(nx,ny,nz) !Isobaric rh,mr (%,kg/kg) . ,pr(nz),pri(nz) !Isobaric pressures (mb) . ,lat1,lat2,lon0,sw(2),ne(2) . ,etalevs(19) . ,xe,mrsat,esat,es . ,slp(nx,ny,12) !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC c !Niv Declaracao das variaveis utilizadas na assimilacao ! integer dia,hora,minuto integer istatAss1,istatAss2,istatAss3,istatAss4 integer kmObs,nObs parameter(kmObs=100) character horAss*2 !identifica arquivo que sera modificado character dadosAssim*53,nomeAer(kmObs)*4 character aero*4 integer ind(nx,ny) real latGr(ny),lonGr(nx) real deltay ,deltax real altiObs(kmObs),latiObs(kmObs),longObs(kmObs) real alti,lati,longi real pnmmObs(kmObs),tempObs(kmObs) real uObs(kmObs),vObs(kmObs) real dir,dirRad,NPI,veloc real dif ! !Niv character(LEN=255):: degrib_dir,degrib_file character(LEN=255):: outdir,outfile,gribfile character(LEN=255):: gdsfile character(LEN=2) :: gproj character(LEN=12) :: atime character(LEN=6) :: model real urlat,urlon c data ETALEVS/1000., 925., 900., 850., 800., + 750., 700., 650., 600., 550., + 500., 450., 400., 350., 300, + 250., 200., 150., 100./ common /estab/esat(15000:45000),es(15000:45000) common /gdsinfo/datsav(11) c_______________________________________________________________________________ c c *** Fill pressure levels. c Cmp nsfcfld=12 Cmp pressure levels are hardwired here do k=1,nz pr(k)=ETALEVS(k) write(6,*) 'PRESS ',pr(k) enddo c c *** Read in degrib data. c n=index(gribfile,' ')-1 horAss(1:1)=gribfile(n-5:n-5) horAss(2:2)=gribfile(n-4:n-4) write(6,*) 'calling read_degrib' call read_degrib(gribfile(1:n) . ,nx*ny,nz,pr,ht,tp,rh,uw,vw,slp,atime) write(6,*) 'BACK FROM read_degrib' write(*,*) 'horAss ',horAss(1:2) contAss=.true. contAss=.false. if(horAss(1:2).eq.'00'.and.(contAss)) then write(99,*) 'gribfile ',gribfile(1:n) write(99,*)'horAss ',horAss(1:2) !Niv ! 21 out 2004 ! Esquema de Assimilacao - Vou inserir a partir deste trecho. ! Assimilando: Metar ! kmObs - quantidade de observacoes ! localidades: Os aeroportos observados estao relacionados ! na pasta: /home2/nivaldo/worketa_all/data/metar/saida ! arquivo: locais.txt ! nomAer(kmObs) ! longObs(kmObs) ! latiObs(kmObs) ! altiObs(kmObs) ! Os dados estao na pasta: /home2/nivaldo/worketa_all/data/metar/2004/mm/dd/hh.txt ! mm - mes ! dd - dia ! hh.txt - hora da observacao ! 00.txt - corresponde a zero hora ! Variaveis: ! PNMM - pnmmObs(kmObs) ! TEMP - tempObs(kmObs) ! VENTO - dirObs uObs(kmObs) ! velObs vObs(kmObs) ! ! Para assimilar, preciso criar as variaveis contendo as coordenadas dos ! pontos de grade. Utilizarei, tambem, as valores da variavel poderosa ! slp(i,j,2) - pressao ao nivel medio do mar - Esquema de Messinger ! slp(i,j,3) - pressao a superficie ! slp(i,j,4) - altitude da superficie ! Estes campos sao utilizados para o ajuste do geopotencial a temperatura. ! A segunda variavel poderosa, utilizada, e datsav(1) - pontos na direcao x ! (2) - pontos na direcao y ! (3) - latitude inicial da grade ! (4) - longitude inicial da grade ! (5) - latitude final ! (6) - longitude final ! Com as informacoes da variavel poderosa 2 vou criar as seguintes variaveis: ! latGr(jGr),lonGr(iGr) latGr( 1)=float(datsav(3))/1000. latGr(ny)=float(datsav(5))/1000. deltay=(latGr(ny)-latGr(1))/float(ny-1) do j=2,ny latGr(j)=latGr(1)+float(j-1)*deltay enddo lonGr( 1)=float(datsav(4))/1000. lonGr(nx)=float(datsav(6))/1000. deltax=(lonGr(nx)-lonGr(1))/float(nx-1) do i=2,nx lonGr(i)=lonGr(1)+float(i-1)*deltax enddo ! print* ! print*,'Longitude' ! print*,lonGr ! print* ! print*,'Latitude' ! print*,latGr ! quero identificar o nivel p para cada ponto ! de grade do i=1,nx do j=1,ny k=1 do dif=slp(i,j,3)/100.0-pr(k) if(dif.gt.0.) then ind(i,j)=k !if(ind(i,j).gt.nz) then ! print*,ind(i,j),slp(i,j,3),i,j,k !endif exit else k=k+1 endif enddo enddo enddo dadosAssim(1:44)='/dados/mestre/etampi/worketa_all/data/metar/' dadosAssim(45:53)='assim.txt' open(10,file=dadosAssim,status='old',iostat=istatAss1) print*,dadosAssim,istatAss1 k=1 if(istatAss1.ne.0) then print*,'Problemas no arquivo a ser assimilado...' print*,'Nao havera assimilacao...' else NPI=4.*atan(1.0) do read(10,'(1x,3(i2,1x),a4,1x,f5.2,1x,f7.2,1x,f7.2,1x,f6.1)' & ,iostat=istatAss2) & dia,hora,minuto,nomeAer(k), & tempObs(k),pnmmObs(k),veloc,dir if(istatAss2.ne.0) exit ! write(*,'(1x,3(i2,1x),a4,4(1x,f8.2))',iostat=istatAss2) ! & dia,hora,minuto,nomeAer(k), ! & tempObs(k),pnmmObs(k),veloc,dir dirRad=dir*NPI/180. uObs(k)=-veloc*sin(dirRad) vObs(k)=-veloc*cos(dirRad) tempObs(k)=tempObs(k)+273.16 pnmmObs(k)=pnmmObs(k)*100.00 k=k+1 if(k.gt.kmObs) then print*,'quantidade de pontos assimilados alcancou valor' print*,'acima do valor permitido kmObs ',kmObs exit endif enddo nObs=k-1 endif close(10) open(10,file='/dados/mestre/etampi/worketa_all/data/metar/' & //'locais.txt',status='old',iostat=istatAss3) if(istatAss3.ne.0) then print*,'Nao encontrei o arquivo com as coordenadas' print*,'das observacoes...' print*,'Teremos problemas...' else k=1 do if(k.gt.nObs)exit read(10,'(a4,1x,3(f7.2,1x))',iostat=istatAss4) & aero,longi,lati,alti if(istatAss4.ne.0) exit k=1 write(*,'(a4,1x,3(f7.2,1x))',iostat=istatAss4) & aero,longi,lati,alti do if(aero.eq.nomeAer(k)) then longObs(k)=longi if(longObs(k) < 0.0) longObs(k)=360.0+longObs(k) latiObs(k)=lati altiObs(k)=alti !print*,nomeAer(k) !print*,longObs(k),latiObs(k),altiObs(k),k !print*,longi,lati,alti exit endif k=k+1 if(k.gt.nObs)exit enddo enddo endif close(10) print*,'chequei aqui... Niv 1' print*,'entrando na ajuste dados...',nObs call AjusteDados(nObs,altiObs,pnmmObs,tempObs,uObs,vObs, & slp, tp, uw, vw, & nx,ny,nz,lonGr,latGr,longObs,latiObs,pr,ind) !Niv endif 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) if (datsav(11) .ne. 104) then print *,'Calling RH patch.' call rh_fix(nx,ny,nz,rh,pr) else Cmp moisture every 50 for grid 104...handle in special way if (mod(k,2) .eq. 0 .or. k.eq.nz) then if (k.lt.(nz-1)) then write(6,*) 'ave value from ', pr(K+1),pr(K-1),k+1,k-1 DO J=1,NY DO I=1,NX RH(I,J,K)=(RH(I,J,K+1)+RH(I,J,K-1))*0.5 ENDDO ENDDO else write(6,*) 'copying value from ', pr(K-1),k-1 DO J=1,NY DO I=1,NX RH(I,J,K)=RH(I,J,K-1) ENDDO ENDDO endif endif endif elseif (uw(1,1,k) .eq. -99999.) then print *,'U-wind data missing at level: ',pr(k) stop elseif (vw(1,1,k) .eq. -99999.) then print *,'V-wind data missing at level: ',pr(k) stop endif enddo C*** C*** Change 1000 hPa heights to match lowest temps C*** biassum=0. do J=1,NY do I=1,NX C old=ht(i,j,1) fact=287.*(pr(2)-pr(3))/(9.81*pr(2)) z975t=(ht(i,j,3)-ht(i,j,2))/fact fact=287.*(pr(1)-pr(2))/(9.81*pr(1)) z1000t=(ht(i,j,2)-ht(i,j,1))/fact C C this criteria empirically derived. Was found that where 950-975 thick > C about 7 m greater than 975-1000 thick, lowest-lev temperatures became hosed. C if (z975t-z1000t .gt. 2.5) then C tbar=(1./4.)*(tp(i,j,1)+tp(i,j,2)+tp(i,j,3)+tp(i,j,4)) tbar=amax1(tp(i,j,1),tp(i,j,2),tp(i,j,3),tp(i,j,4)) ht(i,j,1)=ht(i,j,2)-fact*tbar biassum=biassum+1 endif 633 format(7(f6.2,1x)) enddo enddo C call smooth(ht,nx,ny,nz,2) write(6,*) 'modified Z1000 at ', biassum, ' points' C*** C*** C*** Cmp Cmp Handle missing surface data.... Cmp do k=1,12 if (slp(1,1,k) .eq. -99999.) then write(6,*) 'SURFACE DATA MISSING... ', K if (datsav(11) .ne. 104) then 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 elseif (datsav(11) .eq. 104) then if (k.eq.7 .or. 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,5) enddo enddo elseif(k.eq.8 .or. 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,6) enddo enddo endif endif endif 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 c print *,'ua :',ht(nx/2,ny/2,nz),th(nx/2,ny/2,nz),rh(nx/2,ny/2,nz) . ,uw(nx/2,ny/2,nz),vw(nx/2,ny/2,nz) print *,'ua :',ht(nx/2,ny/2,1),th(nx/2,ny/2,1),rh(nx/2,ny/2,1) . ,uw(nx/2,ny/2,1),vw(nx/2,ny/2,1) Cmp Cmp Write out every 10th grid point in both directions? Cnot goto 1076 write(6,*) 'level 10', nx,ny write(6,*) 'mixr * 1000.' do J=ny,1,-(ny/15) write(6,744) (rh(I,J,10)*1000.,I=1,nx,nx/9) enddo write(6,*) 'HT' do J=ny,1,-(ny/15) write(6,744) (ht(I,J,10),I=1,nx,nx/9) enddo write(6,*) 'UW' do J=ny,1,-(ny/15) write(6,744) (uw(I,J,10),I=1,nx,nx/9) enddo write(6,*) 'VW' do J=ny,1,-(ny/15) write(6,744) (vw(I,J,10),I=1,nx,nx/9) enddo write(6,*) 'TH',' nx ',nx,' ny ',ny,nx/9 do J=ny,1,-(ny/15) write(6,744) (th(I,J,10),I=1,nx,nx/9) enddo 1076 continue 743 format(30(e9.3,1x)) 744 format(30(f7.1,1x)) c c *** Create output file name. c n=index(outdir,' ')-1 if (datsav(11) .eq. 221) model='ETA221' if (datsav(11) .eq. 212) model='ETA212' if (datsav(11) .eq. 104) model='ETA104' if (datsav(11) .eq. 255) model='ETAwrk' outfile=outdir(1:n)//'/'//atime//'.'//model n=index(outfile,' ')-1 print *,model,' data ---> ',outfile(1:n) open(1,file=outfile(1:n),status='unknown',form='unformatted') c c *** Write header stuff. c C nsfcfld=12 if (datsav(10) .eq. 3) gproj='LC' if (datsav(10) .eq. 5) gproj='PS' if (datsav(10) .eq. 0) gproj='CE' write(1) nz n=index(gribfile,' ')-1 write(6,*) n write(6,*) gribfile write(6,*) 'calling get_fullgds with ', + gribfile(1:n) call get_fullgds(gribfile(1:n),datsav(1),datsav(2),kgds) write(6,*) 'back from get_fullgds' C new stuff for GDS file C n=index(outdir,' ')-1 gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model n=index(gdsfile,' ')-1 write(6,*) 'GDS written to ', gdsfile(1:n) open(unit=14,file=gdsfile(1:n),form='unformatted', + access='sequential') write(6,*) 'writing kgds ', (kgds(I),i=1,14) write(14) kgds close(14) C C end new stuff write(6,*) '********************************************' write(6,*) 'writing out the grid with the following info' write(6,*) 'dimensions nx,ny,nz ', nx,ny,nz write(6,*) '********************************************' c c *** Write isobaric upper air data. c write(1) uw write(1) vw write(1) ht write(1) rh write(1) pr write(1) slp write(6,*) 'writing these sfc values to the ETA file...' write(6,*) 'at center of input grid' do k=1,nsfcfld write(6,*) 'field, value ', k, slp(nx/2,ny/2,k) enddo c close(1) c return end c c=============================================================================== c subroutine read_degrib(etafile,nxny,nz,pr . ,ht,tp,rh,uw,vw,slp,atime) c implicit none c integer nxny,nz,i,datsav c real pr(nz),ht(nxny,nz),tp(nxny,nz) . ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz) . ,dummy,slp(nxny,12),tmp(nxny) real crot(nxny),srot(nxny),rmagb,rmagaft,ubef,vbef,urlat, . urlon c integer kgds(200),kpds(50),len,kerr,jpds(200),jgds(200) . ,lenpds,lenkgds,nwords,kpdsl . ,j,k,nx,IRETO,KNUM,IRET1 c character*(*) etafile character(LEN=1):: pds(50) character(LEN=12):: atime character(LEN=2):: gproj logical bitmap(nxny) common /gdsinfo/datsav(11) c_______________________________________________________________________________ c len=index(etafile//' ',' ')-1 c c *** Check that the input dimensions match grib file dimensions. c write(6,*) 'inside read_degrib ' CCC CCC Should be able to derive all header info in the DEGRIB file CCC from the GDS/PDS of the GRIB file. CCC CCC c c *** Fill time stamp (assume all records are for same time). call get_gds(etafile(1:len),datsav,kpds) write(6,*) 'KPDS vals ', kpds(10),kpds(9),kpds(8) if (kpds(8) .ge. 100) kpds(8)=kpds(8)-100 write(atime,'(i2.2,i2.2,i2.2,i2.2,i4.4)') . kpds(8),kpds(9),kpds(10),kpds(11),kpds(14) c c c *** Fill a missing value flag into first space of each variable. c do k=1,nz ht(1,k)=-99999. tp(1,k)=-99999. rh(1,k)=-99999. uw(1,k)=-99999. vw(1,k)=-99999. enddo Cmp initialize surface fields to -99999. so the interp code can handle Cmp appropriately do k=1,12 do j=1,nxny slp(j,k)=-99999. enddo enddo Cmp c c *** Now put the data into the corresponding arrays. c Cmp C add something to read in surface fields in here C call baopen(11,etafile,IRETO) if (IRETO .ne. 0) write(6,*) 'BAOPEN TROUBLE!!!! ', IRETO jpds=-1 jgds=-1 jpds(5)=81 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,1),IRET1) write(6,*) 'first GETGB, IRET1= ', IRET1 write(6,*) 'LAND/SEA READ!!!!! ' write(6,*) (slp(j,1),j=nwords/2,nwords/2+5) jpds(5)=136 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,2),IRET1) slp(:,2)=slp(:,2)*100 write(6,*) 'PMSL READ!!!!! ' write(6,*) (slp(j,2),j=nwords/2,nwords/2+8) jpds(5)=135 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,3),IRET1) slp(:,3)=slp(:,3)*100 write(6,*) 'PSFC READ!!!!! ' write(6,*) (slp(j,3),j=nwords/2,nwords/2+8) jpds(5)=132 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,4),IRET1) write(6,*) 'ZSFC READ!!!!! ' write(6,*) (slp(j,4),j=nwords/2,nwords/2+8) C Cmp SOIL FIELDS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C if (mod(kpds(7)-i,256).eq.0) then Cmp write(6,*) 'kpds(5)= ', kpds(5) Cmp write(6,*) 'lower is ', i Cmp write(6,*) 'upper is ', (kpds(7)-i)/256. C endif jpds(5)=191 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,5),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil temp over ',jpds(7), 'layer!!' write(6,*) (slp(j,5),j=nwords/2,nwords/2+8) endif jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,7),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil temp over ',jpds(7), 'layer!!', IRET1 write(6,*) (slp(j,7),j=nwords/2,nwords/2+8) endif jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,9),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil temp over ',jpds(7), 'layer!!',IRET1 write(6,*) (slp(j,9),j=nwords/2,nwords/2+5) endif jpds(5)=175 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,11),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil temp over ',jpds(7), 'layer!!',IRET1 write(6,*) (slp(j,11),j=nwords/2,nwords/2+5) endif C MOISTURE jpds(5)=182 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,6),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil wet over ',jpds(7), 'layer!!', IRET1 write(6,*) (slp(j,6),j=nwords/2,nwords/2+8) endif jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,8),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil wet over ',jpds(7), 'layer!!', IRET1 write(6,*) (slp(j,8),j=nwords/2,nwords/2+8) endif jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,10),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil wet over ',jpds(7), 'layer!!',IRET1 write(6,*) (slp(j,10),j=nwords/2,nwords/2+5) endif jpds(5)=183 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,12),IRET1) if (IRET1. eq. 0) then write(6,*) 'found soil wet over ',jpds(7), 'layer!!',IRET1 write(6,*) (slp(j,12),j=nwords/2,nwords/2+5) endif C jpds(5)=7 jpds(6)=100 Cmp doing this in proper order?????? Corig do k=1,nz do k=nz,1,-1 jpds(7)=nint(pr(k)) write(6,*) 'LEVEL ', jpds(7) jpds(5)=52 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,rh(1,k),IRET1) if (rh(20,k).gt.1.e20) write(6,*) 'getgbd(1) rh of ',k,rh(20,k) jpds(5)=7 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,ht(1,k),IRET1) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=11 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tp(1,k),IRET1) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) ! jpds(5)=52 ! call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, ! & BITMAP,rh(1,k),IRET1) ! if (rh(20,k).gt.1.e20) write(6,*) 'getgbd(2) rh of ',k,rh(20,k) ! if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=33 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,uw(1,k),IRET1) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=34 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,vw(1,k),IRET1) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) write(6,*) 'Z,T,Q,U,V ', ht(nxny/2+55,k),tp(nxny/2+55,k), + rh(nxny/2+55,k),uw(nxny/2+55,k),vw(nxny/2+55,k) enddo c c *** Normal finish. c 1000 continue C rotate the winds at this point nx=KGDS(2) write(6,*) 'using datsav(10)= ', datsav(10) if (datsav(10) .eq. 3) gproj='LC' if (datsav(10) .eq. 5) gproj='PS' if (gproj .eq. 'LC') then write(6,*) 'rotating a lambert projection' call rotate_lcc(kgds,crot,srot) elseif (gproj .eq. 'PS') then call rotate_str(kgds,crot,srot,urlat,urlon) else write(6,*) 'not doing a wind rotation...' goto 1075 endif do K=1,nz do I=1,nxny rmagb=(uw(I,K)**2. + vw(I,K)**2.)**(0.5) ubef=uw(I,K) vbef=vw(I,K) uw(I,K)=crot(I)*ubef+srot(I)*vbef vw(I,K)=crot(I)*vbef-srot(I)*ubef rmagaft=(uw(I,K)**2. + vw(I,K)**2.)**(0.5) if (abs(rmagaft-rmagb).gt.3.) then write(6,*) 'MAG, I,K,old,new==> ',I,K,rmagb,rmagaft write(6,*) 'original components..', ubef,vbef write(6,*) 'new components..', uw(I,k),vw(I,K) write(6,*) 'rotation cosines ', crot(I),srot(I) write(6,*) '.................................' endif ENDDO ENDDO 1075 continue return c c *** Premature end of file. c 1100 continue print *,'Premature end of file.' print *,'Abort...' stop c end c c=============================================================================== c subroutine rh_fix(nx,ny,nz,rh,pr) c implicit none c integer nx,ny,nz,i,j,k,kk c real 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 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 write(6,*) 'checking RH at lev= ', pr(k) 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=============================================================================== subroutine rotate_lcc(kgds,crot,srot) Cmp stolen/adapted from iplib code gdswiz03 C C SUBPROGRAM: GDSWIZ03 GDS WIZARD FOR LAMBERT CONFORMAL CONICAL C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63) C AND RETURNS ONE OF THE FOLLOWING: C (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES C (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES C FOR LAMBERT CONFORMAL CONICAL PROJECTIONS. C IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT C BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT C OUTPUT ELEMENTS ARE SET TO FILL VALUES. C THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO. C LAMBERT CONFORMAL GRIDS C (2) - NX NR POINTS ALONG X-AXIS C (3) - NY NR POINTS ALONG Y-AXIS C (4) - LA1 LAT OF ORIGIN (LOWER LEFT) C (5) - LO1 LON OF ORIGIN (LOWER LEFT) C (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) C (7) - LOV - ORIENTATION OF GRID C (8) - DX - X-DIR INCREMENT C (9) - DY - Y-DIR INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) C (12) - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER C (13) - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER parameter (NPTS=100000) INTEGER KGDS(200) REAL RLON(NPTS),RLAT(NPTS) REAL CROT(NPTS),SROT(NPTS) real DLON,AN real DE,DR REAL PI,DPR PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - FILL=-999 LROT=1 IROT=1 IM=KGDS(2) JM=KGDS(3) RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) ORIENT=KGDS(7)*1.E-3 DX=KGDS(8) DY=KGDS(9) IPROJ=MOD(KGDS(10)/128,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) RLATI1=KGDS(12)*1.E-3 RLATI2=KGDS(13)*1.E-3 H=(-1.)**IPROJ HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DXS=DX*HI DYS=DY*HJ IF(RLATI1.EQ.RLATI2) THEN AN=SIN(H*RLATI1/DPR) ELSE AN=LOG(COS(RLATI1/DPR)/COS(RLATI2/DPR))/ & LOG(TAN((H*RLATI1+90)/2/DPR)/TAN((H*RLATI2+90)/2/DPR)) ENDIF DE=RERTH*COS(RLATI1/DPR)*TAN((H*RLATI1+90)/2/DPR)**AN/AN IF(H*RLAT1.EQ.90) THEN XP=1 YP=1 ELSE DR=DE/TAN((H*RLAT1+90)/2/DPR)**AN DLON1=MOD(RLON1-ORIENT+180+3600,360.)-180 C atmp=RLON1-ORIENT+180.+3600. C DLON1= atmp - INT (atmp/360.)*360. - 180. XP=1-H*SIN(AN*DLON1/DPR)*DR/DXS YP=1+COS(AN*DLON1/DPR)*DR/DYS ENDIF ANTR=1/(2*AN) DE2=DE**2 XMIN=0 XMAX=IM+1 YMIN=0 YMAX=JM+1 NRET=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE GRID COORDINATES TO EARTH COORDINATES C XP=1 C YP=1 DO N=1,IM*JM J=INT((N-1)/IM)+1 I=N-(J-1)*IM IF(I.GE.XMIN.AND.I.LE.XMAX.AND. & J.GE.YMIN.AND.J.LE.YMAX) THEN DI=(I-XP)*DXS DJ=(J-YP)*DYS DR2=DI**2+DJ**2 IF(DR2.LT.DE2*1.E-6) THEN RLON(N)=0. RLAT(N)=H*90. ELSE RLON(N)=MOD(ORIENT+H/AN*DPR*ATAN2(DI,-DJ)+3600,360.) C atmp=ORIENT+H/AN*DPR*ATAN2(DI,-DJ)+3600 C RLON(N)= atmp - INT(atmp/360.) * 360. RLAT(N)=H*(2*DPR*ATAN((DE2/DR2)**ANTR)-90) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN C atmp=RLON(N)-ORIENT+180.+3600. C DLON=atmp - INT(atmp/360.)*360.-180. DLON=MOD(RLON(N)-ORIENT+180+3600,360.)-180 CROT(N)=H*COS(AN*DLON/DPR) SROT(N)=SIN(AN*DLON/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE RLON(N)=FILL RLAT(N)=FILL ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine get_gds(etafile,gdsinfo,kpds) character*(*) etafile character(LEN=1):: pds(50) integer kgds(200),kpds(200),len,kerr . ,lenpds,lenkgds,nwords,kpdsl . ,j,k,gdsinfo(11) . ,gdsav,IRETO,JGDS(200),JPDS(200) real tmp(100000) logical bitmap(100000) nxny=100000 JPDS=-1 JGDS=-1 len=index(etafile//' ',' ')-1 call baopen(11,etafile(1:len),IRETO) write(6,*) 'BAOPEN in get_gds: ', IRETO if (IRETO .ne. 0) then print *,'Error opening unit=11, file name = ',etafile(1:len) . ,' iostat = ',kerr stop endif jpds(5)=11 jpds(6)=100 jpds(7)=500 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tmp,IRET1) write(6,*) 'back from getgb in get_gds: ', IRET1 gdsinfo(1)=KGDS(2) gdsinfo(2)=KGDS(3) gdsinfo(3)=KGDS(4) gdsinfo(4)=KGDS(5) gdsinfo(5)=KGDS(7) gdsinfo(6)=KGDS(8) gdsinfo(7)=KGDS(9) gdsinfo(8)=KGDS(12) gdsinfo(9)=KGDS(13) gdsinfo(10)=KGDS(1) gdsinfo(11)=KPDS(3) write(6,*) gdsinfo write(6,*) 'KPDSinfo ', (kpds(I),i=1,10) return print *,'GETGDS PROBLEM' stop end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine lambert(gdslatur,gdslonur) C C Subroutine written 9 March 1999 by M. Pyle to support tiled 221 input. C Code adapted from GEMPAK routine gblamb.c. Whole purpose is to get the C UR corner lat and lon for use in workstation Eta. integer latin1,latin2,nx,ny,la1,lo1,lov integer dx,dy,datsav real(kind=8):: earth_rad, const_cone, xll,yll,xur,yur,lat1, + lon1,loncnt,angle1,angle2,x1,x2,y1,y2,alpha,rtemp real gdslatur,gdslonur + ,gdslatll,gdslonll parameter(rpi=3.141592654) parameter(d2r=rpi/180.) parameter(r2d=180./rpi) parameter(radius=6370000.) common /gdsinfo/ datsav(11) latin1=datsav(8) latin2=datsav(9) la1=datsav(3) lo1=datsav(4) dx=datsav(6) dy=datsav(7) nx=datsav(1) ny=datsav(2) lov=datsav(5) write(6,*) 'values in lambert ' write(6,*) latin1,latin2,la1,lo1,dx,dy,nx,ny,lov lat1= (la1/1000.0)*d2r if (lo1 .eq. 180000) lo1=lo1-360000 if (lo1 .lt. -180000) lo1=lo1+360000 lon1= (lo1/1000.0)*d2r Cmp now have LL corner in radians, W is negative if (lov .eq. 180000) lov=lov-360000 if (lov .lt. -180000) lov=lov+360000 loncnt= (lov/1000.0)*d2r angle1= (rpi/2.) - ( abs(latin1/1000.0) * d2r ) angle2= (rpi/2.) - ( abs(latin2/1000.0) * d2r ) if (latin1 .eq. latin2) then const_cone=cos(angle1) else const_cone= ( log ( sin(angle2) ) - log ( sin ( angle1 ) ) )/ + ( log ( tan ( angle2/2 ) ) - log ( tan ( angle1/2 ) ) ) endif write(6,*) 'const_cone= ', const_cone earth_rad=radius/const_cone cmp assuming NH x1 = earth_rad * tan( (rpi/2.-lat1) / 2 )**(const_cone)* + sin (const_cone * ( lon1 - loncnt ) ) y1 = -earth_rad * tan( (rpi/2.-lat1) / 2 )**(const_cone)* + cos (const_cone * ( lon1 - loncnt ) ) alpha= (tan(angle1 / 2 )**const_cone)/sin (angle1) x2=x1 + ( nx - 1 ) * alpha * dx y2=y1 + ( ny - 1 ) * alpha * dy xll=min(x1,x2) xur=max(x1,x2) yll=min(y1,y2) yur=max(y1,y2) gdslatll= ( rpi/2. - 2 * + atan ( ( (abs(xll)**2.+abs(yll)**2.)**(0.5)/earth_rad)** + (1/const_cone) ) ) * r2d write(6,*) 'xll, yll, xll**2., yll**2. ', xll, yll, xll**2., + yll**2. write(6,*) 'lat pieces: ', rpi/2., (xll**2.+yll**2.), earth_rad, + r2d rtemp= atan2 ( xll, -yll ) * ( 1 / const_cone ) + loncnt if ( rtemp .gt. rpi ) then gdslonll = ( rtemp - 2.*rpi ) * r2d else if ( rtemp .lt. -rpi ) then gdslonll = ( rtemp + 2.*rpi ) * r2d else gdslonll = rtemp * r2d endif gdslatur= ( rpi/2. - 2 * + atan ( ( (abs(xur)**2.+abs(yur)**2.)**(0.5)/earth_rad)** + (1/const_cone) ) ) * r2d rtemp= atan2 ( xur, -yur ) * ( 1 / const_cone ) + loncnt if ( rtemp .gt. rpi ) then gdslonur = ( rtemp - 2.*rpi ) * r2d else if ( rtemp .lt. -rpi ) then gdslonur = ( rtemp + 2.*rpi ) * r2d else gdslonur = rtemp * r2d endif write(6,*) 'output==> ' write(6,*) 'LL points ', gdslatll,gdslonll write(6,*) 'UR points ', gdslatur,gdslonur return END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine rotate_str(kgds,crot,srot,urlat,urlon) C C Program adapted 26 March 99 for workstation Eta. Original code from C iplib program GDSWZD05 C C C SUBPROGRAM: GDSWZD05 GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63) C AND RETURNS ONE OF THE FOLLOWING: C (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES C (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES C FOR POLAR STEREOGRAPHIC AZIMUTHAL PROJECTIONS. C IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT C BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT C OUTPUT ELEMENTS ARE SET TO FILL VALUES. C THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO. C OPTIONALLY, THE VECTOR ROTATIONS AND THE MAP JACOBIANS C FOR THIS GRID MAY BE RETURNED AS WELL. C POLAR STEREOGRAPHIC GRIDS C (2) - N(I) NR POINTS ALONG LAT CIRCLE C (3) - N(J) NR POINTS ALONG LON CIRCLE C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) C (7) - LOV GRID ORIENTATION C (8) - DX - X DIRECTION INCREMENT C (9) - DY - Y DIRECTION INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28) parameter (NPTS=80000) INTEGER KGDS(200) REAL RLON(NPTS),RLAT(NPTS) REAL CROT(NPTS),SROT(NPTS) real DLON,AN real DE,DR,urlat,urlon REAL PI,DPR PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LROT=1 IROT=1 IM=KGDS(2) JM=KGDS(3) RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) ORIENT=KGDS(7)*1.E-3 DX=KGDS(8) DY=KGDS(9) IPROJ=MOD(KGDS(10)/128,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) H=(-1.)**IPROJ HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DXS=DX*HI DYS=DY*HJ DE=(1.+SIN(60./DPR))*RERTH DR=DE*COS(RLAT1/DPR)/(1+H*SIN(RLAT1/DPR)) XP=1-H*SIN((RLON1-ORIENT)/DPR)*DR/DXS YP=1+COS((RLON1-ORIENT)/DPR)*DR/DYS DE2=DE**2 XMIN=0 XMAX=IM+1 YMIN=0 YMAX=JM+1 NRET=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE GRID COORDINATES TO EARTH COORDINATES DO N=1,IM*JM J=INT((N-1)/IM)+1 I=N-(J-1)*IM IF(I.GE.XMIN.AND.I.LE.XMAX.AND. & J.GE.YMIN.AND.J.LE.YMAX) THEN DI=(I-XP)*DXS DJ=(J-YP)*DYS DR2=DI**2+DJ**2 IF(DR2.LT.DE2*1.E-6) THEN RLON(N)=0. RLAT(N)=H*90. ELSE RLON(N)=MOD(ORIENT+H*DPR*ATAN2(DI,-DJ)+3600,360.) RLAT(N)=H*DPR*ASIN((DE2-DR2)/(DE2+DR2)) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN CROT(N)=H*COS((RLON(N)-ORIENT)/DPR) SROT(N)=SIN((RLON(N)-ORIENT)/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - urlat=rlat(im*jm) urlon=rlon(im*jm) if (urlon .gt. 180) urlon=urlon-360. if (urlon .lt. -180) urlon=urlon+360. write(6,*) 'in stereo sub found lat,lon ', urlat,urlon C write(6,*) 'lat/lon in stereo ' do N=1,IM*JM J=INT((N-1)/IM)+1 I=N-(J-1)*IM C write(6,69) I,J,RLAT(N),RLON(N) enddo 69 format(I3,1x,I3,1x,2(f7.3,1x)) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine get_fullgds(etafile,nx,ny,kgds) character*(*) etafile character(LEN=1):: pds(50) integer kgds(200),kpds(200),len,kerr,jpds(200) . ,lenpds,lenkgds,nwords,kpdsl,jgds(200) . ,j,k,KNUM,nx,ny logical bitmap(nx*ny) real tmp(nx*ny) write(6,*) 'inside get_fullgds...', etafile nxny=nx*ny len=index(etafile//' ',' ')-1 jpds=-1 jgds=-1 jpds(5)=11 jpds(6)=100 jpds(7)=500 write(6,*) 'calling getgb ' C write(6,*) 'jpds = ', jpds C write(6,*) 'jgds = ', jgds call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tmp,IRET1) write(6,*) 'return from getgb ', IRET1 if (IRET1 .ne. 0) then print *,'Error getting GDS in get_fullgds ', IRET1 stop endif write(6,*) 'leaving get_fullgds' return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine smooth(a,idim,jdim,ldim,npass) integer idim,jdim,ldim,npass real a(idim,jdim,ldim) DO N=1,NPASS do L=1,LDIM do J=2,JDIM-1 do I=2,IDIM-1 a(i,j,l)=1./8.*(4.*a(i,j,l)+a(i+1,j,l)+a(i-1,j,l) + +a(i,j+1,l)+a(i,j-1,l)) enddo enddo enddo ENDDO return end !Niv ! subroutine AjusteDados(nObs,altO,pnmmObs,tempObs,uObs,vObs, & slp, tp, uw, vw, & im,jm,km, x, y, xo, yo,pr,ind) implicit none logical contAlt integer i,i1,ifator,im,j,j1,jfator,jm,k,karq,km integer kmObs,nObs parameter(kmObs=100) integer irec real Rdef,tan,yVal real x0,y0 integer ind(im,jm) real altP(im,jm),dist(im,jm,kmObs) real altO(kmObs) real W(im,jm,kmObs) real aux(im,jm),aux1(im,jm),aux2(im,jm) real x(im),xo(kmObs),y(jm),yo(kmObs) real pnmmObs(kmObs),tempObs(kmObs),uObs(kmObs),vObs(kmObs) real tp(im,jm,km),uw(im,jm,km),vw(im,jm,km) real pr(km) real slp(im,jm,12) do i=1,im do j=1,jm altP(i,j)=slp(i,j,4) enddo enddo print*,'altitudes das observacoes...' !do k=1,nObs ! write(*,*) altO(k) !enddo W=9999.9 Rdef=2.5 Rdef=5.0 Rdef=7.5 contAlt=.true. if(contAlt) then ! Testando o efeito de altitude... do k=1,nObs do j=1,jm do i=1,im dist(i,j,k)=sqrt((yo(k)-y(j))**2+(xo(k)-x(i))**2) !write(99,*)dist(i,j,k),xo(k),yo(k),i,j,k !write(99,*)x(i),y(j) !write(99,*)altP(i,j) !write(99,*) if(dist(i,j,k) < Rdef) then if(xo(k) /= x(i)) then tan=(yo(k)-y(j))/(xo(k)-x(i)) else tan=9999.99 endif if((tan > 0.).and.(tan < 9999.99))then if(xo(k) > x(i)) then i1=i+1 if(i1 > im) i1=im ifator=1 jfator=1 else i1=i-1 if(i1 < 1) i1=1 ifator=-1 jfator=-1 endif j1=j do !if(x(i1) < xo(k)) then if((x(i1)-xo(k))*ifator < 0) then yval=(x(i1)-x(i))*tan + y(j) if((yval > y(j1)).and.(yval <= y(j1+1)))then if((altP(i,j) < altP(i1,j1 )).and. & (altP(i,j) < altP(i1,j1+1*jfator)).and. & (altO(k) < altP(i1,j1 )).and. & (altO(k) < altP(i1,j1+1*jfator))) then !O ponto i,j nao e influenciado pela obs. W(i,j,k)=0.0 exit else i1=i1+1*ifator endif else j1=j1+1*jfator endif !if(i1 >= im) exit else ! !O ponto e influenciado. !Preciso W(i,j,k)=1.0 exit endif enddo else if(tan < 0.0) then !case(:0.0) if(x(i) < xo(k)) then i1=i+1 if(i1 > im) i1 = im ifator=1 jfator=-1 else i1=i-1 if(i1 < 1) i1 = 1 ifator=-1 jfator=1 endif j1=j do !read* !if(x(i1) > xo(k))then if((x(i1)-xo(k))*ifator < 0) then yval=(x(i1)-x(i))*tan+y(j) if((yval > y(j1)).and.(yval<=y(j1+1))) then if((altP(i,j) < altP(i1,j1 )).and. & (altP(i,j) < altP(i1,j1+1*jfator)).and. & (altO(k) < altP(i1,j1 )).and. & (altO(k) < altP(i1,j1+1*jfator))) then !Este ponto nao e influenciado W(i,j,k)=0.0 exit else i1=i1+1*ifator endif else j1=j1+1*jfator !if(j1 >= jm) then ! W(i,j,k)=1.0 ! exit !endif endif !if(i1 <=1 ) then ! W(i,j,k)=1.0 ! exit !endif else !O ponto e influenciado, vou calcular !o peso. W(i,j,k)=1.0 exit endif enddo else if(tan == 0.) then if(xo(k) < x(i)) then j1=j i1=i-1 if(i1 < 1) i1 =1 ifator=-1 do if((altP(i,j) < altP(i1,j1)).and. & (altO(k) < altP(i1,j-1))) then !Este ponto nao e influenciado W(i,j,k)=0.0 exit else i1=i1+1*ifator if((i1 < 1).or.(xo(k) > x(i1))) then W(i,j,k)=1.0 exit endif endif enddo else j1=j i1=i+1 if(i1 > im) i1=im ifator=1 do if((altP(i,j) < altP(i1,j1)).and. & (altO(k) < altP(11,j1))) then !Este ponto nao e influenciado W(i,j,k)=0.0 exit else i1=i1+1*ifator if((i1 > im).or.(xo(k) < x(i1))) then W(i,j,k)=1.0 exit endif endif enddo endif endif if(tan == 9999.99) then if(yo(k) < y(j)) then i1=i j1=j-1 if(j1 < 1) j1 =1 jfator=-1 do if(altP(i,j) < altP(i1,j1)) then !.and. !Este ponto nao e influenciado W(i,j,k)=0.0 exit else j1=j1+1*jfator if((j1 < 1).or.(yo(k) > y(j1))) then W(i,j,k)=1.0 exit endif endif enddo else if((yo(k)==y(j)).and.(xo(k)==x(i))) then W(i,j,k)=1.0 !exit else i1=i j1=j+1 if(j1 > jm) j1=jm jfator=1 do if(altP(i,j) < altP(i1,j1)) then !.and. !Este ponto nao e influenciado W(i,j,k)=0.0 exit else j1=j1+1*jfator if((j1 > jm).or.(yo(k) < y(j1))) then W(i,j,k)=1.0 exit endif endif enddo endif endif endif endif endif else W(i,j,k)=0.0 endif enddo enddo enddo else W=1.0 endif print*,'vou entrar na Assimila...' do i=1,im do j=1,jm aux(i,j)=uw(i,j,ind(i,j)) aux2(i,j)=aux(i,j) enddo enddo irec=0 call Assimila(uObs,aux,aux1,W,altP,altO,im,jm,nObs, & x,y,xo,yo,irec) do i=1,im do j=1,jm write(88,*) aux2(i,j)-aux1(i,j) uw(i,j,ind(i,j))=aux1(i,j) enddo enddo print*,'assimilei u' do i=1,im do j=1,jm aux(i,j)=vw(i,j,ind(i,j)) enddo enddo call Assimila(vObs,aux,aux1,W,altP,altO,im,jm,nObs, & x,y,xo,yo,irec) do i=1,im do j=1,jm vw(i,j,ind(i,j))=aux1(i,j) enddo enddo print*,'assimilei v' do i=1,im do j=1,jm aux(i,j)=tp(i,j,ind(i,j)) enddo enddo call Assimila(tempObs,aux,aux1,W,altP,altO,im,jm,nObs, & x,y,xo,yo,irec) do i=1,im do j=1,jm tp(i,j,ind(i,j))=aux1(i,j) enddo enddo print*,'assimilei Temp' do i=1,im do j=1,jm aux(i,j)=slp(i,j,2) enddo enddo call Assimila(pnmmObs,aux,aux1,W,altP,altO,im,jm,nObs, & x,y,xo,yo,irec) do i=1,im do j=1,jm slp(i,j,2)=aux1(i,j) enddo enddo print*,'assimilei pnmm' print*,'saindo da Ajuste...' return end !Niv subroutine Assimila(tsmObs,tsm,tsmNova,W1,altP,altO,nx,ny,nObs, & lon,lat,lonObs,latObs,irec) !definicao das variaveis implicit none integer irec integer i,istat1,istat2,istat3,im,j,jm,k,km,kmax integer nx,ny integer kmObs,nObs parameter(kmObs=100) real alt,angZen real delta,difLon,difLat,dist,distV real fatorCor,fatorCorV real latMax,latMin,lonMax,lonMin real rdef,rVert,somaPeso,somaPesoV real tsm1,tsm2,tsm3,tsm4 real lon(nx),lat(ny),latObs(kmObs),lonObs(kmObs) real tsm(nx,ny),tsmNova(nx,ny) real tsmInt(kmObs),tsmObs(kmObs) real altP(nx,ny),altO(kmObs) real w1(nx,ny,kmObs) real w(kmObs),wVert(kmObs) rdef=2.5 rdef=5.0 rdef=7.5 rVert=25. rVert=50. rVert=75. tsmNova=tsm latMax=-999.9 latMin= 999.9 lonMax=-999.9 lonMin= 999.9 do k=1,nObs if(lonObs(k) < 0.0) lonObs(k)=360.0+lonObs(k) if(latMax < latObs(k)) latMax=latObs(k) if(latMin > latObs(k)) latMin=latObs(k) if(lonMax < lonObs(k)) lonMax=lonObs(k) if(lonMin > lonObs(k)) lonMin=lonObs(k) enddo close(10) !Os dados ja foram lidos... !Geracao das latitudes e longitudes dos campos do NMC do k=1,nObs i=1 do !1 difLon=lonObs(k)-lon(i) j=1 if(difLon < delta) then do !2 difLat=latObs(k)-lat(j) if(difLat < delta) then tsm1=tsm(i,j)+(tsm(i+1,j)-tsm(i,j))* & (lonObs(k)-lon(i))/(lon(i+1)-lon(i)) tsm2=tsm(i,j+1)+(tsm(i+1,j+1)-tsm(i,j+1))* & (lonObs(k)-lon(i))/(lon(i+1)-lon(i)) tsm3=tsm(i,j)+(tsm(i,j+1)-tsm(i,j))* & (latObs(k)-lat(j))/(lat(j+1)-lat(j)) tsm4=tsm(i+1,j)+(tsm(i+1,j+1)-tsm(i+1,j))* & (latObs(k)-lat(j))/(lat(j+1)-lat(j)) tsmInt(k)=(tsm1+tsm2+tsm3+tsm4)*0.25 i=nx+1 exit else j=j+1 endif if(j>ny)exit enddo!2 else i=i+1 endif if(i>nx)exit enddo!1 enddo !read* !Calculo do peso das observacoe e ajustes... tsmNova=tsm do i=1,nx if((lon(i) <= lonMax).and.(lon(i)>=lonMin)) then do j=1,ny if((lat(j) <= latMax).and.(lat(j)>=latMin)) then somaPeso=0.0 do k=1,nObs if(W1(i,j,k)==1.0) then dist=(lonObs(k)-lon(i))**2+(latObs(k)-lat(j))**2 dist=sqrt(dist) if(dist < Rdef) then w(k)=(Rdef**2-dist**2)/(Rdef**2+dist**2) w(k)=exp(-1.0*dist*dist/(2.0*Rdef*Rdef)) else w(k)=0.0 endif distV=abs(altP(i,j)-altO(k)) ! if(distV < RVert) then ! wVert(k)=(RVert**2-distV**2)/(RVert**2+distV**2) ! else ! wVert(k)=0.0 ! endif if(distV > RVert) then wVert(k)=0.0 else if(distV <= RVert/4.0) then wVert(k)=1.0 else wVert=1.0-distV/(altP(i,j)+altO(k)) endif endif somaPeso=somaPeso+w(k) somaPesoV=somaPesoV+wVert(k) else wVert(k)=0.0 W(k)=0.0 endif enddo fatorCor=0.0 fatorCorV=0.0 ! wVert=1.0 !SITUACAO ESPECIAL... ESTOU FAZENDO TESTES... ! somaPesoV=0.0 !Estou fazendo testes... do k=1,nObs fatorCor=fatorCor+w(k)**2*wVert(k)**2 & *(tsmObs(k)-tsmInt(k)) enddo if((somaPeso+somaPesoV) > 0.0) then fatorCor=fatorCor/(somaPeso+somaPesoV) endif tsmNova(i,j)=tsm(i,j)+fatorCor write(99,'(4(e18.8,1x),2(i4,1x))') & tsmNova(i,j),tsm(i,j), & fatorCor,somaPeso,i,j endif enddo endif enddo close(10) OPEN(10,file='dadAjus.dat',form='unformatted',status='unknown', & access='direct',recl=nx*ny*4) irec=irec+1 write(10,REC=irec)tsm irec=irec+1 write(10,REC=irec)tsmNova irec=irec+1 write(10,REC=irec)altP close(10) open(10,file='dados1.dat',status='unknown') do k=1,nObs write(10,'(i4,1x,2(e12.4,1x))')k,tsmObs(k),tsmInt(k) enddo close(10) print*,'saindo da Assimila...' return end !Niv