PROGRAM covar IMPLICIT NONE INTEGER, PARAMETER :: d_x=144, d_y=157, Nlevels=19 REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),19,19 ) :: rstdv_geo REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),13,13 ) :: rstdv_hum REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),1,1 ) :: rstdv_slp REAL, DIMENSION( 22608,19 ) :: asigma_geo REAL, DIMENSION( 22608,13 ) :: asigma_hum REAL, DIMENSION( 22608,1 ) :: asigma_slp INTEGER :: ilev,ix,iy,irec CALL covargeo(rstdv_geo) CALL covarhum(rstdv_hum) CALL covarslp(rstdv_slp) OPEN(UNIT=10,FILE='sigma_geo_hum_slp.asc',STATUS='unknown') cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE escreve_sigma(rsigma,d_x,d_y,d_z) IMPLICIT NONE INTEGER, INTENT(IN) :: d_x,d_y,d_z REAL, DIMENSION( d_x,d_y,d_z ), INTENT(IN) :: rsigma INTEGER :: ix,iy,iz INTEGER :: ld,lb,ilev,iline,LevBeg,Lev0,LevEnd,irec INTEGER :: Nlon,Nlat,Nlevels,Nlines=2 c INTEGER, DIMENSION(19) :: Plevels=(/ 1000,925,900,850,800,750,700, c & 650,600,550,500,450,400, c & 350,300,250,200,150,100 /) INTEGER, DIMENSION(d_z) :: Plevels REAL :: Xmin=277.,DeltaX=0.399994,Ymin=-50.,DeltaY=0.399998 REAL :: Amiss=1.E20 LOGICAL :: back CHARACTER (LEN=12) :: Ctlfile,Binfile,Ascfile CHARACTER (LEN=9) :: DateEnd='01nov2003' CHARACTER (LEN=8) :: lug CHARACTER (LEN=2) :: Hour='00' CHARACTER (LEN=12) :: luavail Nlon=d_x Nlat=d_y Nlevels=d_z Plevels=(/ 1000,925,900,850,800,750,700, & 650,600,550,500,450,400,350, & 300,250,200,150,100 /) lug=luavail() Ctlfile=lug//'.ctl' open(UNIT=10,FILE=Ctlfile,STATUS='unknown') WRITE(UNIT=10,FMT=11) '^'//lug//'.bin' 11 format('DSET ',a) WRITE(UNIT=10,FMT=12) 12 format('OPTIONS yrev') WRITE(UNIT=10,FMT=13) 13 format('TITLE Forecast merw error "standard deviations" for PSAS') WRITE(UNIT=10,FMT=14)Amiss 14 format('UNDEF',e12.5) WRITE(UNIT=10,FMT=15) 'X',Nlon,Xmin,DeltaX WRITE(UNIT=10,FMT=15) 'Y',Nlat,Ymin,DeltaY 15 FORMAT(a1,'DEF ',i3,' LINEAR ',f7.1,1x,f9.6) WRITE(UNIT=10,FMT=16)Nlevels 16 FORMAT('ZDEF ',i2,' LEVELS') DO 20, iline=1,Nlines Lev0=10*(iline-1) LevBeg=Lev0+1 LevEnd=min(Nlevels,Lev0+10) write(UNIT=10,FMT=21)(PLevels(ilev),ilev=LevBeg, LevEnd) 20 END DO 21 FORMAT(10I4) write(UNIT=10,FMT=23)Hour, DateEnd 23 format('TDEF 1 linear ',a2,'Z',a9,' 6hr') WRITE(UNIT=10,FMT=4) 4 format('vars 1') write(UNIT=10,FMT=25) Nlevels 25 format('mwin ',i3,' 0 Meridional Wind', & ' forecast error standard deviation') write(UNIT=10,FMT=28) 28 format('ENDVARS') close(UNIT=10) Binfile=lug//'.bin' Ascfile=lug//'.asc' open(UNIT=11,FILE=Binfile,FORM="unformatted", & ACCESS="DIRECT",RECL=144*157) irec=1 DO ilev=1,Nlevels WRITE(UNIT=11,REC=irec)((rsigma(1,1,ilev), & ix=1,d_x), iy=1,d_y) irec=irec+1 END DO CLOSE(UNIT=11) open(UNIT=12,FILE=Ascfile,STATUS='unknown') DO ilev=1,Nlevels WRITE(UNIT=11,FMT='(1x,E12.5)') & ( (rsigma(1,1,ilev), ix=1,d_x) & , iy=1,d_y ) END DO CLOSE(UNIT=12) RETURN END SUBROUTINE escreve_sigma cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION luavail() IMPLICIT NONE CHARACTER (LEN=8) :: luavail luavail='mersigma' END FUNCTION luavail cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE covargeo(rstdv) IMPLICIT NONE !Definicao de variaveis CHARACTER (LEN=28), PARAMETER :: s_infla_i= & 'gzip -d /home/alexcp/omf/omf' c & 'gzip -d /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=7), PARAMETER :: s_infla_f='.dat.gz' CHARACTER (LEN=25), PARAMETER :: s_reduz_i= & 'gzip /home/alexcp/omf/omf' c & 'gzip /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=4), PARAMETER :: s_reduz_f='.dat' !INTEGER, PARAMETER :: n_arq=91 INTEGER, PARAMETER :: d_x=144, d_y=157, d_z=19 INTEGER, PARAMETER :: n_arq=91 !Definicao de tipo TYPE lugar INTEGER :: local_dx,local_dy,local_dz END TYPE lugar TYPE lugar_xy INTEGER :: local_dx,local_dy END TYPE lugar_xy !Definicao de variaveis auxiliares INTEGER, DIMENSION(d_z) :: local_z=-1 TYPE(lugar_xy), DIMENSION(d_x,d_y) :: local_xy=lugar_xy(-1,-1) TYPE(lugar), DIMENSION(d_z) :: local=lugar(-1,-1,-1) INTEGER :: k,i,ix,iy,iz,iix,iiy,iiz,icont_max INTEGER :: idx,idy,idz REAL :: status, SYSTEM INTEGER, DIMENSION(n_arq) :: n_ob=0 INTEGER :: iob REAL :: rlat,rlon,rlevel INTEGER :: ijulian,itime,ikx,ikt REAL, DIMENSION(d_z) :: rsigma,rsigma_e,rstdv_z_max=-1. REAL :: romf_a REAL, DIMENSION(n_arq,d_x,d_y,d_z) :: romf=0. REAL, DIMENSION(d_x,d_y,d_z) :: rmedia=0.,rsigma_escrito REAL, DIMENSION( d_x,d_y,d_z,d_z ) :: rcovar_zm=0.,rcovar_zk=0. REAL, DIMENSION( d_z,d_z ) :: rcovar_z=0.,rcovar_zn=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z ) :: rcovar_xym=0., & rcovar_xyn,rcovar_xyk=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xy=0.,rcovar_xynn REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xymn=0. INTEGER :: iqc,iqch CHARACTER (LEN=10) :: arq CHARACTER (LEN=45) :: s_infla CHARACTER (LEN=39) :: s_reduz CHARACTER (LEN=18) :: s3 REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z,d_z ), INTENT(OUT) :: & rstdv DO k=1,2 print *, "k=",k ! arquivo de 00h if( k==1) arq='0020030412' if( k==2) arq='0020030413' if( k==3) arq='0020030414' !'0020030415' !'0020030506' if( k==4) arq='0020030417' !'0020030418' if( k==5) arq='0020030419' if( k==6) arq='0020030420' if( k==7) arq='0020030421' !'0020030422' if( k==8) arq='0020030423' if( k==9) arq='0020030424' if( k==10) arq='0020030425' if( k==11) arq='0020030426' if( k==12) arq='0020030427' if( k==13) arq='0020030428' if( k==14) arq='0020030429' if( k==15) arq='0020030430' if( k==16) arq='0020030501' if( k==17) arq='0020030502' if( k==18) arq='0020030503' if( k==19) arq='0020030504' if( k==20) arq='0020030505' if( k==21) arq='0020030506' !arquivos de 06h if( k==22) arq='0620030412' if( k==23) arq='0620030413' if( k==24) arq='0620030414' if( k==25) arq='0620030415' if( k==26) arq='0620030416' if( k==27) arq='0620030417' if( k==28) arq='0620030418' if( k==29) arq='0620030419' if( k==30) arq='0620030420' if( k==31) arq='0620030421' if( k==32) arq='0620030422' if( k==33) arq='0620030423' if( k==34) arq='0620030424' if( k==35) arq='0620030425' if( k==36) arq='0620030426' if( k==37) arq='0620030427' if( k==38) arq='0620030428' if( k==39) arq='0620030429' if( k==40) arq='0620030430' if( k==41) arq='0620030501' if( k==42) arq='0620030502' if( k==43) arq='0620030503' if( k==44) arq='0620030504' if( k==45) arq='0620030504' if( k==46) arq='0620030505' !arquivos de 12h if( k==47) arq='1220030412' if( k==48) arq='1220030413' if( k==49) arq='1220030414' if( k==50) arq='1220030415' !arq='1220030416' if( k==51) arq='1220030417' if( k==52) arq='1220030418' if( k==53) arq='1220030419' if( k==54) arq='1220030420' if( k==55) arq='1220030421' if( k==56) arq='1220030422' if( k==57) arq='1220030423' if( k==58) arq='1220030424' if( k==59) arq='1220030425' if( k==60) arq='1220030426' if( k==61) arq='1220030427' if( k==62) arq='1220030428' if( k==63) arq='1220030429' if( k==64) arq='1220030430' if( k==65) arq='1220030501' if( k==66) arq='1220030502' if( k==67) arq='1220030503' if( k==68) arq='1220030504' !arquivos de 18h if( k==69) arq='1820030412' if( k==70) arq='1820030413' if( k==71) arq='1820030414' if( k==72) arq='1820030415' if( k==73) arq='1820030416' if( k==74) arq='1820030417' if( k==75) arq='1820030418' if( k==76) arq='1820030419' if( k==77) arq='1820030420' if( k==78) arq='1820030421' if( k==79) arq='1820030422' if( k==80) arq='1820030423' if( k==81) arq='1820030424' if( k==82) arq='1820030425' if( k==83) arq='1820030426' if( k==84) arq='1820030427' if( k==85) arq='1820030428' if( k==86) arq='1820030429' if( k==87) arq='1820030430' if( k==88) arq='1820030501' if( k==89) arq='1820030502' if( k==90) arq='1820030503' if( k==91) arq='1820030504' s_infla=s_infla_i//arq//s_infla_f !print *,s_infla status= SYSTEM( s_infla ) s_reduz=s_reduz_i//arq//s_reduz_f !print *,s_reduz s3='omf'//arq//'.dat' print *, s3 open(UNIT=10,FILE=s3,STATUS="OLD") READ ( UNIT=10, FMT="(48X)") READ ( UNIT=10, FMT="(68X)") READ ( UNIT=10, FMT="(37x,I4)") n_ob(k) READ ( UNIT=10, FMT="(46X)") READ ( UNIT=10, FMT="(44X)") READ ( UNIT=10, FMT="(1X)") READ ( UNIT=10, FMT="(76X)") 1000 FORMAT(I7,F9.2,F9.2,F9.2,I8,I5,I5,I5,F9.2,I5,I5) DO i=1,n_ob(k) !PRINT *, "i=",i READ ( UNIT=10, FMT=1000) & iob,rlat,rlon,rlevel,ijulian,itime,ikx,ikt,romf_a,iqc,iqch IF( (rlon >= (-83.)) .AND. (rlon <= (-25.4)) .AND. & (rlat >= (-50.200001)) .AND. (rlat <= (12.8)) & .AND. ( (iqc/=17) .OR. (iqch/=17) ) .AND. & (ikt==6) ) THEN iy=int( (rlat+50.200001)/0.399994 ) ix=int( (rlon+83.)/0.399998 ) iz=0 IF( rlevel==(1000.) ) iz=1 IF( rlevel==(925.) ) iz=2 IF( rlevel==(900.) ) iz=3 IF( rlevel==(850.) ) iz=4 IF( rlevel==(800.) ) iz=5 IF( rlevel==(750.) ) iz=6 IF( rlevel==(700.) ) iz=7 IF( rlevel==(650.) ) iz=8 IF( rlevel==(600.) ) iz=9 IF( rlevel==(550.) ) iz=10 IF( rlevel==(500.) ) iz=11 IF( rlevel==(450.) ) iz=12 IF( rlevel==(400.) ) iz=13 IF( rlevel==(350.) ) iz=14 IF( rlevel==(300.) ) iz=15 IF( rlevel==(250.) ) iz=16 IF( rlevel==(200.) ) iz=17 IF( rlevel==(150.) ) iz=18 IF( rlevel==(100.) ) iz=19 IF( iz==0 ) goto 10 DO iiz=1,d_z DO iiy=1,d_y DO iix=1,d_x IF( (iiz==iz) .AND. (iiy==iy) .AND. (iix==ix) ) THEN romf(k,ix,iy,iz)=romf_a END IF END DO END DO END DO 10 END IF END DO PRINT *, "calcula media" DO iz=1,d_z DO iy=1,d_x DO ix=1,d_y rmedia(ix,iy,iz)=rmedia(ix,iy,iz)* & REAL(n_arq-1)/REAL(n_arq)+ & romf(k,ix,iy,iz)/REAL(n_arq) END DO END DO END DO status = SYSTEM( s_reduz ) close(UNIT=10) END DO print *,"" print *, "n_arq=", n_arq DO k=1,2 print *, "k=",k print *, "calcula covariancia vertical" DO ix=1,d_x DO iy=1,d_y DO iz=1,d_z DO iiz=1,d_z rcovar_zk(ix,iy,iz,iiz)=rcovar_zk(ix,iy,iz,iiz)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix,iy,iiz)-rmedia(ix,iy,iiz) ) END DO END DO END DO END DO !goto 1 print *, "calcula covariancia horizontal" DO iz=1,d_z DO iy=1,d_y DO idy=0,((d_y-1)-(iy-1)) DO ix=1,d_x DO idx=0,((d_x-1)-(ix-1)) rcovar_xyk(idx,idy,iz)=rcovar_xyk(idx,idy,iz) & *REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix+idx,iy+idy,iz)- & rmedia(ix+idx,iy+idy,iz) ) & /REAL(d_x*d_y) END DO END DO END DO END DO END DO 1 rcovar_zm=rcovar_zm*REAL(n_arq-1)/REAL(n_arq) & +rcovar_zk/REAL(n_arq) where (rcovar_zm<0.) rcovar_zm=0. end where rcovar_xym=rcovar_xym*REAL(n_arq-1)/REAL(n_arq) & +rcovar_xyk/REAL(n_arq) where (rcovar_xym<0.) rcovar_xym=0. end where END DO !stop print *, "final" DO iz=1,d_z DO iiz=1,d_z DO iy=1,d_y DO ix=1,d_x rcovar_z(iz,iiz)=rcovar_z(iz,iiz)* & REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & rcovar_zm(ix,iy,iz,iiz)/REAL(d_x*d_y) END DO END DO END DO END DO DO iz=1,d_z rsigma(iz)=sqrt( rcovar_z(iz,iz) ) END DO call estima_sigma_geo(rsigma,rsigma_e,d_z) DO iz=1,d_z rcovar_z(iz,iz)=(rsigma_e(iz))**2. END DO DO iz=1,d_z DO iiz=1,d_z IF( (ABS( rcovar_z(iz,iz)* & rcovar_z(iiz,iiz) )/=0.) ) THEN rcovar_zn(iz,iiz)=rcovar_z(iz,iiz)* & 1./SQRT( ABS(rcovar_z(iiz,iiz)* & rcovar_z(iz,iz)) ) END IF END DO END DO DO iz=1,d_z call estima_z_geo(iz,rcovar_zn(iz,:),d_z) rcovar_z(iz,:)=( rsigma_e(iz)**2. )*rcovar_zn(iz,:) END DO !stop print *, "v normalizado" rcovar_xyn=rcovar_xym DO iz=1,d_z call estima_xy(rcovar_xyn(:,:,iz),d_x,d_y) END DO DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xy(idx,idy)=rcovar_xy(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xyn(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xynn=rcovar_xy/rcovar_xy(0,0) print *, "h normalizado" call prodext(d_x,d_y,d_z,rcovar_z,rcovar_xynn,rstdv) return DO iz=1,d_z DO iiz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) IF( ABS(rstdv(idx,idy,iz,iiz))>rstdv_z_max(iz) ) THEN rstdv_z_max(iz)=ABS(rstdv(idx,idy,iz,iiz)) local(iz)=lugar(idx,idy,iiz) END IF END DO END DO END DO END DO !PRINT *,"stdv_z_max=",rstdv_z_max print *,"local=" DO iz=11,11 print *, "iz=",iz print *, local(iz) DO iiz=1,5 print *, "iiz=",iiz DO idy=0,4 print *, rstdv(0:4,idy,iz,iiz) END DO END DO END DO DO iz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) rsigma_escrito(idx+1,idy+1,iz)=rstdv(idx,idy,iz,iz) END DO END DO END DO CALL escreve_sigma_geo(rsigma_escrito,d_x,d_y,d_z) !stop RETURN c OPEN(UNIT=20,FILE="geo_e_xy",STATUS="UNKNOWN") DO idy=0,9 c WRITE ( UNIT=20, FMT=3000) ( rcovar_xynn(idx,idy), idx=0,9 ) END DO 3000 FORMAT( 10(1x,E12.3) ) c CLOSE(UNIT=20) DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xymn(idx,idy)=rcovar_xymn(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xym(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xymn=sqrt(rcovar_xymn) rcovar_xymn=rcovar_xymn/rcovar_xymn(0,0) c OPEN(UNIT=21,FILE="geo_xy",STATUS="UNKNOWN") DO idy=0,9 c WRITE ( UNIT=21,FMT=3000) ( rcovar_xymn(idx,idy), idx=0,9 ) END DO c CLOSE(UNIT=21) c OPEN(UNIT=30,FILE="geo_e_z",STATUS="UNKNOWN") c WRITE ( UNIT=30, FMT=3001) ( rstdv(0,0,iz,iz), iz=1,d_z ) 3001 FORMAT( 19(1x,E12.3) ) c CLOSE(UNIT=30) c OPEN(UNIT=31,FILE="geo_z",STATUS="UNKNOWN") c WRITE ( UNIT=31, FMT=3001) ( rsigma(iz), iz=1,d_z ) c CLOSE(UNIT=31) !OPEN(UNIT=32,FILE="Covariancia_z",STATUS="UNKNOWN") !WRITE ( UNIT=32, FMT=3001) ( rcovar_z(19,iz), iz=1,d_z ) !CLOSE(UNIT=32) !OPEN(UNIT=33,FILE="Covariancia_e_z",STATUS="UNKNOWN") !WRITE(UNIT=33,FMT=3001)(rcovar_zm(72,78,19,iz),iz=1,d_z) !CLOSE(UNIT=33) RETURN END SUBROUTINE covargeo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE covarslp(rstdv) IMPLICIT NONE !Definicao de variaveis CHARACTER (LEN=28), PARAMETER :: s_infla_i= & 'gzip -d /home/alexcp/omf/omf' c & 'gzip -d /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=7), PARAMETER :: s_infla_f='.dat.gz' CHARACTER (LEN=25), PARAMETER :: s_reduz_i= & 'gzip /home/alexcp/omf/omf' c & 'gzip /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=4), PARAMETER :: s_reduz_f='.dat' !INTEGER, PARAMETER :: n_arq=91 INTEGER, PARAMETER :: d_x=144, d_y=157, d_z=1 INTEGER, PARAMETER :: n_arq=91 !Definicao de tipo TYPE lugar INTEGER :: local_dx,local_dy,local_dz END TYPE lugar TYPE lugar_xy INTEGER :: local_dx,local_dy END TYPE lugar_xy !Definicao de variaveis auxiliares INTEGER, DIMENSION(d_z) :: local_z=-1 TYPE(lugar_xy), DIMENSION(d_x,d_y) :: local_xy=lugar_xy(-1,-1) TYPE(lugar), DIMENSION(d_z) :: local=lugar(-1,-1,-1) INTEGER :: k,i,ix,iy,iz,iix,iiy,iiz,icont_max INTEGER :: idx,idy,idz REAL :: status, SYSTEM INTEGER, DIMENSION(n_arq) :: n_ob=0 INTEGER :: iob REAL :: rlat,rlon,rlevel INTEGER :: ijulian,itime,ikx,ikt REAL, DIMENSION(d_z) :: rsigma,rsigma_e,rstdv_z_max=-1. REAL :: rsigma_max=-1.,romf_a REAL, DIMENSION(n_arq,d_x,d_y,d_z) :: romf=0. REAL, DIMENSION(d_x,d_y,d_z) :: rmedia=0.,rsigma_escrito REAL, DIMENSION( d_x,d_y,d_z,d_z ) :: rcovar_zm=0.,rcovar_zk=0. REAL, DIMENSION( d_z,d_z ) :: rcovar_z=0.,rcovar_zn=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z ) :: rcovar_xym=0., & rcovar_xyn,rcovar_xyk=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xy=0.,rcovar_xynn REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xymn=0. INTEGER :: iqc,iqch CHARACTER (LEN=10) :: arq CHARACTER (LEN=45) :: s_infla CHARACTER (LEN=39) :: s_reduz CHARACTER (LEN=18) :: s3 REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z,d_z ), INTENT(OUT) :: & rstdv DO k=1,5 print *, "k=",k ! arquivo de 00h if( k==1) arq='0020030412' if( k==2) arq='0020030413' if( k==3) arq='0020030414' !'0020030415' !'0020030506' if( k==4) arq='0020030417' !'0020030418' if( k==5) arq='0020030419' if( k==6) arq='0020030420' if( k==7) arq='0020030421' !'0020030422' if( k==8) arq='0020030423' if( k==9) arq='0020030424' if( k==10) arq='0020030425' if( k==11) arq='0020030426' if( k==12) arq='0020030427' if( k==13) arq='0020030428' if( k==14) arq='0020030429' if( k==15) arq='0020030430' if( k==16) arq='0020030501' if( k==17) arq='0020030502' if( k==18) arq='0020030503' if( k==19) arq='0020030504' if( k==20) arq='0020030505' if( k==21) arq='0020030506' !arquivos de 06h if( k==22) arq='0620030412' if( k==23) arq='0620030413' if( k==24) arq='0620030414' if( k==25) arq='0620030415' if( k==26) arq='0620030416' if( k==27) arq='0620030417' if( k==28) arq='0620030418' if( k==29) arq='0620030419' if( k==30) arq='0620030420' if( k==31) arq='0620030421' if( k==32) arq='0620030422' if( k==33) arq='0620030423' if( k==34) arq='0620030424' if( k==35) arq='0620030425' if( k==36) arq='0620030426' if( k==37) arq='0620030427' if( k==38) arq='0620030428' if( k==39) arq='0620030429' if( k==40) arq='0620030430' if( k==41) arq='0620030501' if( k==42) arq='0620030502' if( k==43) arq='0620030503' if( k==44) arq='0620030504' if( k==45) arq='0620030504' if( k==46) arq='0620030505' !arquivos de 12h if( k==47) arq='1220030412' if( k==48) arq='1220030413' if( k==49) arq='1220030414' if( k==50) arq='1220030415' !arq='1220030416' if( k==51) arq='1220030417' if( k==52) arq='1220030418' if( k==53) arq='1220030419' if( k==54) arq='1220030420' if( k==55) arq='1220030421' if( k==56) arq='1220030422' if( k==57) arq='1220030423' if( k==58) arq='1220030424' if( k==59) arq='1220030425' if( k==60) arq='1220030426' if( k==61) arq='1220030427' if( k==62) arq='1220030428' if( k==63) arq='1220030429' if( k==64) arq='1220030430' if( k==65) arq='1220030501' if( k==66) arq='1220030502' if( k==67) arq='1220030503' if( k==68) arq='1220030504' !arquivos de 18h if( k==69) arq='1820030412' if( k==70) arq='1820030413' if( k==71) arq='1820030414' if( k==72) arq='1820030415' if( k==73) arq='1820030416' if( k==74) arq='1820030417' if( k==75) arq='1820030418' if( k==76) arq='1820030419' if( k==77) arq='1820030420' if( k==78) arq='1820030421' if( k==79) arq='1820030422' if( k==80) arq='1820030423' if( k==81) arq='1820030424' if( k==82) arq='1820030425' if( k==83) arq='1820030426' if( k==84) arq='1820030427' if( k==85) arq='1820030428' if( k==86) arq='1820030429' if( k==87) arq='1820030430' if( k==88) arq='1820030501' if( k==89) arq='1820030502' if( k==90) arq='1820030503' if( k==91) arq='1820030504' s_infla=s_infla_i//arq//s_infla_f !print *,s_infla status= SYSTEM( s_infla ) s_reduz=s_reduz_i//arq//s_reduz_f !print *,s_reduz s3='omf'//arq//'.dat' print *, s3 open(UNIT=10,FILE=s3,STATUS="OLD") READ ( UNIT=10, FMT="(48X)") READ ( UNIT=10, FMT="(68X)") READ ( UNIT=10, FMT="(37x,I4)") n_ob(k) READ ( UNIT=10, FMT="(46X)") READ ( UNIT=10, FMT="(44X)") READ ( UNIT=10, FMT="(1X)") READ ( UNIT=10, FMT="(76X)") 1000 FORMAT(I7,F9.2,F9.2,F9.2,I8,I5,I5,I5,F9.2,I5,I5) DO i=1,n_ob(k) !PRINT *, "i=",i READ ( UNIT=10, FMT=1000) & iob,rlat,rlon,rlevel,ijulian,itime,ikx,ikt,romf_a,iqc,iqch IF( (rlon >= (-83.)) .AND. (rlon <= (-25.4)) .AND. & (rlat >= (-50.200001)) .AND. (rlat <= (12.8)) & .AND. ( (iqc/=17) .OR. (iqch/=17) ) .AND. & (ikt==3) ) THEN iy=int( (rlat+50.200001)/0.399994 ) ix=int( (rlon+83.)/0.399998 ) iz=0 IF( rlevel==(1000.) ) iz=1 IF( rlevel==(925.) ) iz=2 IF( rlevel==(900.) ) iz=3 IF( rlevel==(850.) ) iz=4 IF( rlevel==(800.) ) iz=5 IF( rlevel==(750.) ) iz=6 IF( rlevel==(700.) ) iz=7 IF( rlevel==(650.) ) iz=8 IF( rlevel==(600.) ) iz=9 IF( rlevel==(550.) ) iz=10 IF( rlevel==(500.) ) iz=11 IF( rlevel==(450.) ) iz=12 IF( rlevel==(400.) ) iz=13 IF( rlevel==(350.) ) iz=14 IF( rlevel==(300.) ) iz=15 IF( rlevel==(250.) ) iz=16 IF( rlevel==(200.) ) iz=17 IF( rlevel==(150.) ) iz=18 IF( rlevel==(100.) ) iz=19 IF( iz==0 ) goto 10 DO iiz=1,d_z DO iiy=1,d_y DO iix=1,d_x IF( (iiz==iz) .AND. (iiy==iy) .AND. (iix==ix) ) THEN romf(k,ix,iy,iz)=romf_a END IF END DO END DO END DO 10 END IF END DO PRINT *, "calcula media" DO iz=1,d_z DO iy=1,d_x DO ix=1,d_y rmedia(ix,iy,iz)=rmedia(ix,iy,iz)* & REAL(n_arq-1)/REAL(n_arq)+ & romf(k,ix,iy,iz)/REAL(n_arq) END DO END DO END DO status = SYSTEM( s_reduz ) close(UNIT=10) END DO print *,"" print *, "n_arq=", n_arq DO k=1,5 print *, "k=",k print *, "calcula covariancia vertical" DO ix=1,d_x DO iy=1,d_y DO iz=1,d_z DO iiz=1,d_z rcovar_zk(ix,iy,iz,iiz)=rcovar_zk(ix,iy,iz,iiz)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix,iy,iiz)-rmedia(ix,iy,iiz) ) END DO END DO END DO END DO !goto 1 print *, "calcula covariancia horizontal" DO iz=1,d_z DO iy=1,d_y DO idy=0,((d_y-1)-(iy-1)) DO ix=1,d_x DO idx=0,((d_x-1)-(ix-1)) rcovar_xyk(idx,idy,iz)=rcovar_xyk(idx,idy,iz) & *REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix+idx,iy+idy,iz)- & rmedia(ix+idx,iy+idy,iz) ) & /REAL(d_x*d_y) END DO END DO END DO END DO END DO 1 rcovar_zm=rcovar_zm*REAL(n_arq-1)/REAL(n_arq) & +rcovar_zk/REAL(n_arq) where (rcovar_zm<0.) rcovar_zm=0. end where rcovar_xym=rcovar_xym*REAL(n_arq-1)/REAL(n_arq) & +rcovar_xyk/REAL(n_arq) where (rcovar_xym<0.) rcovar_xym=0. end where END DO !stop print *, "final" DO iz=1,d_z DO iiz=1,d_z DO iy=1,d_y DO ix=1,d_x rcovar_z(iz,iiz)=rcovar_z(iz,iiz)* & REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & rcovar_zm(ix,iy,iz,iiz)/REAL(d_x*d_y) END DO END DO END DO END DO DO iz=1,d_z rsigma(iz)=sqrt( rcovar_z(iz,iz) ) IF( rsigma(iz)>rsigma_max ) rsigma_max=rsigma(iz) END DO IF( rsigma_max>101 ) THEN rsigma_e=101 ELSE call estima_sigma_slp(rsigma,rsigma_e,d_z) END IF DO iz=1,d_z rcovar_z(iz,iz)=(rsigma_e(iz))**2. END DO DO iz=1,d_z DO iiz=1,d_z IF( (ABS( rcovar_z(iz,iz)* & rcovar_z(iiz,iiz) )/=0.) ) THEN rcovar_zn(iz,iiz)=rcovar_z(iz,iiz)* & 1./SQRT( ABS(rcovar_z(iiz,iiz)* & rcovar_z(iz,iz)) ) END IF END DO END DO DO iz=1,d_z call estima_z_slp(iz,rcovar_zn(iz,:),d_z) rcovar_z(iz,:)=( rsigma_e(iz)**2. )*rcovar_zn(iz,:) END DO !stop print *, "v normalizado" rcovar_xyn=rcovar_xym DO iz=1,d_z call estima_xy(rcovar_xyn(:,:,iz),d_x,d_y) END DO DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xy(idx,idy)=rcovar_xy(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xyn(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xynn=rcovar_xy/rcovar_xy(0,0) print *, "h normalizado" call prodext(d_x,d_y,d_z,rcovar_z,rcovar_xynn,rstdv) return DO iz=1,d_z DO iiz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) IF( ABS(rstdv(idx,idy,iz,iiz))>rstdv_z_max(iz) ) THEN rstdv_z_max(iz)=ABS(rstdv(idx,idy,iz,iiz)) local(iz)=lugar(idx,idy,iiz) END IF END DO END DO END DO END DO !PRINT *,"stdv_z_max=",rstdv_z_max print *,"local=" DO iz=1,1 print *, "iz=",iz print *, local(iz) DO iiz=1,5 print *, "iiz=",iiz DO idy=0,4 print *, rstdv(0:4,idy,iz,iiz) END DO END DO END DO DO iz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) rsigma_escrito(idx+1,idy+1,iz)=rstdv(idx,idy,iz,iz) END DO END DO END DO CALL escreve_sigma_hum(rsigma_escrito,d_x,d_y,d_z) !stop RETURN c OPEN(UNIT=20,FILE="slp_e_xy",STATUS="UNKNOWN") DO iz=10,10 DO idy=0,9 c WRITE ( UNIT=20, FMT=3000) ( rcovar_xynn(idx,idy), idx=0,9 ) END DO END DO 3000 FORMAT( 10(1x,E12.3) ) c CLOSE(UNIT=20) DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xymn(idx,idy)=rcovar_xymn(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xym(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xymn=sqrt(rcovar_xymn) IF( rcovar_xymn(0,0) < 1.E-6 ) rcovar_xymn=1.E-6 rcovar_xymn=rcovar_xymn/rcovar_xymn(0,0) c OPEN(UNIT=21,FILE="slp_xy",STATUS="UNKNOWN") DO iz=10,10 DO idy=0,9 c WRITE ( UNIT=21,FMT=3000) ( rcovar_xymn(idx,idy), idx=0,9 ) END DO END DO c CLOSE(UNIT=21) c OPEN(UNIT=30,FILE="slp_e_z",STATUS="UNKNOWN") c WRITE ( UNIT=30, FMT=3001) ( rstdv(0,0,11,iz), iz=1,d_z ) 3001 FORMAT( 19(1x,E12.3) ) c CLOSE(UNIT=30) c OPEN(UNIT=31,FILE="slp_z",STATUS="UNKNOWN") c WRITE ( UNIT=31, FMT=3001) ( rsigma(iz), iz=1,d_z ) c CLOSE(UNIT=31) !OPEN(UNIT=32,FILE="Covariancia_z",STATUS="UNKNOWN") !WRITE ( UNIT=32, FMT=3001) ( rcovar_z(19,iz), iz=1,d_z ) !CLOSE(UNIT=32) !OPEN(UNIT=33,FILE="Covariancia_e_z",STATUS="UNKNOWN") !WRITE(UNIT=33,FMT=3001)(rcovar_zm(72,78,19,iz),iz=1,d_z) !CLOSE(UNIT=33) RETURN END SUBROUTINE covarslp cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE covarhum(rstdv) IMPLICIT NONE !Definicao de variaveis CHARACTER (LEN=28), PARAMETER :: s_infla_i= & 'gzip -d /home/alexcp/omf/omf' c & 'gzip -d /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=7), PARAMETER :: s_infla_f='.dat.gz' CHARACTER (LEN=25), PARAMETER :: s_reduz_i= & 'gzip /home/alexcp/omf/omf' c & 'gzip /home/alex/DMD/RPSAS/omf/omf' CHARACTER (LEN=4), PARAMETER :: s_reduz_f='.dat' !INTEGER, PARAMETER :: n_arq=91 INTEGER, PARAMETER :: d_x=144, d_y=157 , d_z=13 INTEGER, PARAMETER :: n_arq=91 !Definicao de tipo TYPE lugar INTEGER :: local_dx,local_dy,local_dz END TYPE lugar TYPE lugar_xy INTEGER :: local_dx,local_dy END TYPE lugar_xy !Definicao de variaveis auxiliares INTEGER, DIMENSION(d_z) :: local_z=-1 TYPE(lugar_xy), DIMENSION(d_x,d_y) :: local_xy=lugar_xy(-1,-1) TYPE(lugar), DIMENSION(d_z) :: local=lugar(-1,-1,-1) INTEGER :: k,i,ix,iy,iz,iix,iiy,iiz,icont_max INTEGER :: idx,idy,idz REAL :: status, SYSTEM INTEGER, DIMENSION(n_arq) :: n_ob=0 INTEGER :: iob REAL :: rlat,rlon,rlevel INTEGER :: ijulian,itime,ikx,ikt REAL, DIMENSION(d_z) :: rsigma,rsigma_e,rstdv_z_max=-1. REAL :: romf_a REAL, DIMENSION(n_arq,d_x,d_y,d_z) :: romf=0. REAL, DIMENSION(d_x,d_y,d_z) :: rmedia=0.,rsigma_escrito REAL, DIMENSION( d_x,d_y,d_z,d_z ) :: rcovar_zm=0.,rcovar_zk=0. REAL, DIMENSION( d_z,d_z ) :: rcovar_z=0.,rcovar_zn=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z ) :: rcovar_xym=0., & rcovar_xyn,rcovar_xyk=0. REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xy=0.,rcovar_xynn REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ) :: rcovar_xymn=0. INTEGER :: iqc,iqch CHARACTER (LEN=10) :: arq CHARACTER (LEN=45) :: s_infla CHARACTER (LEN=39) :: s_reduz CHARACTER (LEN=18) :: s3 REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z,d_z ), INTENT(OUT) :: & rstdv DO k=1,2 print *, "k=",k ! arquivo de 00h if( k==1) arq='0020030412' if( k==2) arq='0020030413' if( k==3) arq='0020030414' !'0020030415' !'0020030506' if( k==4) arq='0020030417' !'0020030418' if( k==5) arq='0020030419' if( k==6) arq='0020030420' if( k==7) arq='0020030421' !'0020030422' if( k==8) arq='0020030423' if( k==9) arq='0020030424' if( k==10) arq='0020030425' if( k==11) arq='0020030426' if( k==12) arq='0020030427' if( k==13) arq='0020030428' if( k==14) arq='0020030429' if( k==15) arq='0020030430' if( k==16) arq='0020030501' if( k==17) arq='0020030502' if( k==18) arq='0020030503' if( k==19) arq='0020030504' if( k==20) arq='0020030505' if( k==21) arq='0020030506' !arquivos de 06h if( k==22) arq='0620030412' if( k==23) arq='0620030413' if( k==24) arq='0620030414' if( k==25) arq='0620030415' if( k==26) arq='0620030416' if( k==27) arq='0620030417' if( k==28) arq='0620030418' if( k==29) arq='0620030419' if( k==30) arq='0620030420' if( k==31) arq='0620030421' if( k==32) arq='0620030422' if( k==33) arq='0620030423' if( k==34) arq='0620030424' if( k==35) arq='0620030425' if( k==36) arq='0620030426' if( k==37) arq='0620030427' if( k==38) arq='0620030428' if( k==39) arq='0620030429' if( k==40) arq='0620030430' if( k==41) arq='0620030501' if( k==42) arq='0620030502' if( k==43) arq='0620030503' if( k==44) arq='0620030504' if( k==45) arq='0620030504' if( k==46) arq='0620030505' !arquivos de 12h if( k==47) arq='1220030412' if( k==48) arq='1220030413' if( k==49) arq='1220030414' if( k==50) arq='1220030415' !arq='1220030416' if( k==51) arq='1220030417' if( k==52) arq='1220030418' if( k==53) arq='1220030419' if( k==54) arq='1220030420' if( k==55) arq='1220030421' if( k==56) arq='1220030422' if( k==57) arq='1220030423' if( k==58) arq='1220030424' if( k==59) arq='1220030425' if( k==60) arq='1220030426' if( k==61) arq='1220030427' if( k==62) arq='1220030428' if( k==63) arq='1220030429' if( k==64) arq='1220030430' if( k==65) arq='1220030501' if( k==66) arq='1220030502' if( k==67) arq='1220030503' if( k==68) arq='1220030504' !arquivos de 18h if( k==69) arq='1820030412' if( k==70) arq='1820030413' if( k==71) arq='1820030414' if( k==72) arq='1820030415' if( k==73) arq='1820030416' if( k==74) arq='1820030417' if( k==75) arq='1820030418' if( k==76) arq='1820030419' if( k==77) arq='1820030420' if( k==78) arq='1820030421' if( k==79) arq='1820030422' if( k==80) arq='1820030423' if( k==81) arq='1820030424' if( k==82) arq='1820030425' if( k==83) arq='1820030426' if( k==84) arq='1820030427' if( k==85) arq='1820030428' if( k==86) arq='1820030429' if( k==87) arq='1820030430' if( k==88) arq='1820030501' if( k==89) arq='1820030502' if( k==90) arq='1820030503' if( k==91) arq='1820030504' s_infla=s_infla_i//arq//s_infla_f !print *,s_infla status= SYSTEM( s_infla ) s_reduz=s_reduz_i//arq//s_reduz_f !print *,s_reduz s3='omf'//arq//'.dat' print *, s3 open(UNIT=10,FILE=s3,STATUS="OLD") READ ( UNIT=10, FMT="(48X)") READ ( UNIT=10, FMT="(68X)") READ ( UNIT=10, FMT="(37x,I4)") n_ob(k) READ ( UNIT=10, FMT="(46X)") READ ( UNIT=10, FMT="(44X)") READ ( UNIT=10, FMT="(1X)") READ ( UNIT=10, FMT="(76X)") 1000 FORMAT(I7,F9.2,F9.2,F9.2,I8,I5,I5,I5,F9.2,I5,I5) DO i=1,n_ob(k) !PRINT *, "i=",i READ ( UNIT=10, FMT=1000) & iob,rlat,rlon,rlevel,ijulian,itime,ikx,ikt,romf_a,iqc,iqch IF( (rlon >= (-83.)) .AND. (rlon <= (-25.4)) .AND. & (rlat >= (-50.200001)) .AND. (rlat <= (12.8)) & .AND. ( (iqc/=17) .OR. (iqch/=17) ) .AND. & (ikt==5) ) THEN iy=int( (rlat+50.200001)/0.399994 ) ix=int( (rlon+83.)/0.399998 ) iz=0 IF( rlevel==(1000.) ) iz=1 IF( rlevel==(925.) ) iz=2 IF( rlevel==(900.) ) iz=3 IF( rlevel==(850.) ) iz=4 IF( rlevel==(800.) ) iz=5 IF( rlevel==(750.) ) iz=6 IF( rlevel==(700.) ) iz=7 IF( rlevel==(650.) ) iz=8 IF( rlevel==(600.) ) iz=9 IF( rlevel==(550.) ) iz=10 IF( rlevel==(500.) ) iz=11 IF( rlevel==(450.) ) iz=12 IF( rlevel==(400.) ) iz=13 IF( rlevel==(350.) ) iz=14 IF( rlevel==(300.) ) iz=15 IF( rlevel==(250.) ) iz=16 IF( rlevel==(200.) ) iz=17 IF( rlevel==(150.) ) iz=18 IF( rlevel==(100.) ) iz=19 IF( iz==0 ) goto 10 DO iiz=1,d_z DO iiy=1,d_y DO iix=1,d_x IF( (iiz==iz) .AND. (iiy==iy) .AND. (iix==ix) ) THEN romf(k,ix,iy,iz)=romf_a END IF END DO END DO END DO 10 END IF END DO PRINT *, "calcula media" DO iz=1,d_z DO iy=1,d_x DO ix=1,d_y rmedia(ix,iy,iz)=rmedia(ix,iy,iz)* & REAL(n_arq-1)/REAL(n_arq)+ & romf(k,ix,iy,iz)/REAL(n_arq) END DO END DO END DO status = SYSTEM( s_reduz ) close(UNIT=10) END DO print *,"" print *, "n_arq=", n_arq DO k=1,2 print *, "k=",k print *, "calcula covariancia vertical" DO ix=1,d_x DO iy=1,d_y DO iz=1,d_z DO iiz=1,d_z rcovar_zk(ix,iy,iz,iiz)=rcovar_zk(ix,iy,iz,iiz)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix,iy,iiz)-rmedia(ix,iy,iiz) ) END DO END DO END DO END DO !goto 1 print *, "calcula covariancia horizontal" DO iz=1,d_z DO iy=1,d_y DO idy=0,((d_y-1)-(iy-1)) DO ix=1,d_x DO idx=0,((d_x-1)-(ix-1)) rcovar_xyk(idx,idy,iz)=rcovar_xyk(idx,idy,iz) & *REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & ( romf(k,ix,iy,iz)-rmedia(ix,iy,iz) )* & ( romf(k,ix+idx,iy+idy,iz)- & rmedia(ix+idx,iy+idy,iz) ) & /REAL(d_x*d_y) END DO END DO END DO END DO END DO 1 rcovar_zm=rcovar_zm*REAL(n_arq-1)/REAL(n_arq) & +rcovar_zk/REAL(n_arq) where (rcovar_zm<0.) rcovar_zm=0. end where rcovar_xym=rcovar_xym*REAL(n_arq-1)/REAL(n_arq) & +rcovar_xyk/REAL(n_arq) where (rcovar_xym<0.) rcovar_xym=0. end where END DO print *, "final" DO iz=1,d_z DO iiz=1,d_z DO iy=1,d_y DO ix=1,d_x rcovar_z(iz,iiz)=rcovar_z(iz,iiz)* & REAL(d_x*d_y-1)/REAL(d_x*d_y)+ & rcovar_zm(ix,iy,iz,iiz)/REAL(d_x*d_y) END DO END DO END DO END DO DO iz=1,d_z rsigma(iz)=sqrt( rcovar_z(iz,iz) ) END DO call estima_sigma_hum(rsigma,rsigma_e,d_z) DO iz=1,d_z c print *, "no principal" c print *,"iz=",iz," desvp=",rsigma(iz)," desvp_e=",rsigma_e(iz) END DO DO iz=1,d_z rcovar_z(iz,iz)=(rsigma_e(iz))**2. END DO DO iz=1,d_z DO iiz=1,d_z IF( (ABS( rcovar_z(iz,iz)* & rcovar_z(iiz,iiz) )/=0.) ) THEN rcovar_zn(iz,iiz)=rcovar_z(iz,iiz)* & 1./SQRT( ABS(rcovar_z(iiz,iiz)* & rcovar_z(iz,iz)) ) END IF END DO END DO DO iz=1,d_z call estima_z_hum(iz,rcovar_zn(iz,:),d_z) rcovar_z(iz,:)=( rsigma_e(iz)**2. )*rcovar_zn(iz,:) END DO !stop do iz=1,19 c print *,"rcovar_z=" c print *,rcovar_z(iz,iz) end do print *, "v normalizado" rcovar_xyn=rcovar_xym DO iz=1,d_z call estima_xy(rcovar_xyn(:,:,iz),d_x,d_y) END DO DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xy(idx,idy)=rcovar_xy(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xyn(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xynn=rcovar_xy/rcovar_xy(0,0) print *, "h normalizado" call prodext(d_x,d_y,d_z,rcovar_z,rcovar_xynn,rstdv) return DO iz=1,d_z DO iiz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) IF( ABS(rstdv(idx,idy,iz,iiz))>rstdv_z_max(iz) ) THEN rstdv_z_max(iz)=ABS(rstdv(idx,idy,iz,iiz)) local(iz)=lugar(idx,idy,iiz) END IF END DO END DO END DO END DO !PRINT *,"stdv_z_max=",rstdv_z_max print *,"local=" DO iz=11,11 print *, "iz=",iz print *, local(iz) DO iiz=1,5 print *, "iiz=",iiz DO idy=0,4 print *, rstdv(0:4,idy,iz,iiz) END DO END DO END DO DO iz=1,d_z DO idy=0,(d_y-1) DO idx=0,(d_x-1) rsigma_escrito(idx+1,idy+1,iz)=rstdv(idx,idy,iz,iz) END DO END DO END DO CALL escreve_sigma_slp(rsigma_escrito,d_x,d_y,d_z) stop c OPEN(UNIT=20,FILE="hum_e_xy",STATUS="UNKNOWN") DO idy=0,9 c WRITE ( UNIT=20, FMT=3000) ( rcovar_xynn(idx,idy), idx=0,9 ) END DO 3000 FORMAT( 10(1x,E12.3) ) c CLOSE(UNIT=20) DO idx=0,(d_x-1) DO idy=0,(d_y-1) DO iz=1,d_z rcovar_xymn(idx,idy)=rcovar_xymn(idx,idy)* & REAL(d_z-1)/REAL(d_z) & +rcovar_xym(idx,idy,iz)/REAL(d_z) END DO END DO END DO rcovar_xymn=sqrt(rcovar_xymn) rcovar_xymn=rcovar_xymn/rcovar_xymn(0,0) c OPEN(UNIT=21,FILE="hum_xy",STATUS="UNKNOWN") DO idy=0,9 c WRITE ( UNIT=21,FMT=3000) ( rcovar_xymn(idx,idy), idx=0,9 ) END DO c CLOSE(UNIT=21) c OPEN(UNIT=30,FILE="hum_e_z",STATUS="UNKNOWN") c WRITE ( UNIT=30, FMT=3001) ( rstdv(0,0,iz,iz), iz=1,d_z ) 3001 FORMAT( 19(1x,E12.3) ) c CLOSE(UNIT=30) c OPEN(UNIT=31,FILE="hum_z",STATUS="UNKNOWN") c WRITE ( UNIT=31, FMT=3001) ( rsigma(iz), iz=1,19 ) c CLOSE(UNIT=31) !OPEN(UNIT=32,FILE="Covariancia_z",STATUS="UNKNOWN") !WRITE ( UNIT=32, FMT=3001) ( rcovar_z(19,iz), iz=1,d_z ) !CLOSE(UNIT=32) !OPEN(UNIT=33,FILE="Covariancia_e_z",STATUS="UNKNOWN") !WRITE(UNIT=33,FMT=3001)(rcovar_zm(72,78,19,iz),iz=1,d_z) !CLOSE(UNIT=33) RETURN END SUBROUTINE covarhum ccccccccccccccccccccccccccccccccccc SUBROUTINE prodext(d_x,d_y,d_z,rz,rxy,rprod) IMPLICIT NONE INTEGER, INTENT(IN) :: d_x,d_y,d_z REAL, DIMENSION(d_z,d_z), INTENT(IN) :: rz REAL, DIMENSION( 0:(d_x-1),0:(d_y-1) ), INTENT(IN) :: rxy REAL, DIMENSION( 0:(d_x-1),0:(d_y-1),d_z,d_z), INTENT(OUT) :: rprod INTEGER :: i,j,k,m print *,"multiplica" DO m=1,d_z DO k=1,d_z DO j=0,(d_y-1) DO i=0,(d_x-1) rprod(i,j,k,m)=rxy(i,j)*rz(k,m) END DO END DO END DO END DO rprod=sqrt(rprod) END SUBROUTINE prodext cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_xy(rcovar_xy,d_x,d_y) IMPLICIT NONE INTEGER :: ix,iy INTEGER, INTENT(IN) :: d_x,d_y REAL :: rix,riy,rsumc,rw,rsum REAL :: rcov_x,rcov_y REAL, DIMENSION(0:(d_x-1)) :: rcovar_x REAL, DIMENSION(0:(d_y-1)) :: rcovar_y REAL, DIMENSION(0:(d_x-1),0:(d_y-1)) :: rcovar_xy rsum=0. rsumc=0. DO ix=0,(d_x-1) rix=REAL(ix)**(4.) rsum=rix+rsum END DO DO iy=0,(d_y-1) DO ix=0,(d_x-1) rix=REAL(ix)**(2.) if( rcovar_xy(ix,iy)>=1.E-6 ) THEN rcov_x=LOG( rcovar_xy(ix,iy) ) ELSE rcov_x=LOG(1.E-6) ENDIF rsumc=rcov_x*rix+rsumc END DO rw=rsumc/rsum DO ix=0,(d_x-1) rcovar_x(ix)=EXP( rw*REAL(ix)**2. ) END DO rcovar_x(0:(d_x-1))=rcovar_x(0:(d_x-1)) & /rcovar_x(0) END DO rsum=0. rsumc=0. DO iy=0,(d_y-1) riy=REAL(iy)**(4.) rsum=riy+rsum END DO DO ix=0,(d_x-1) DO iy=0,(d_y-1) riy=REAL(iy)**(2.) if( rcovar_xy(ix,iy)>=1.E-6 ) THEN rcov_y=LOG( rcovar_xy(ix,iy) ) ELSE rcov_y=LOG(1.E-6) ENDIF rsumc=rcov_y*riy+rsumc END DO !print *, "rsum=",rsum rw=rsumc/rsum DO iy=0,(d_y-1) rcovar_y(iy)=EXP( rw*REAL(iy)**2. ) END DO rcovar_y(0:(d_y-1))=rcovar_y(0:(d_y-1)) & /rcovar_y(0) END DO DO iy=0,(d_y-1) DO ix=0,(d_x-1) rcovar_xy(ix,iy)=rcovar_x(ix)*rcovar_y(iy) END DO END DO !stop RETURN END SUBROUTINE estima_xy cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_z_geo(iz,rcovar_z,d_z) IMPLICIT NONE INTEGER :: iiz INTEGER, INTENT(IN) :: iz,d_z REAL :: riz,rcov_z,s=0. REAL :: rsumc=0.,rw REAL :: rsum=0. REAL, DIMENSION(d_z) :: rcovar_z IF ( d_z==1 ) RETURN DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(2.*1.6) rsum=riz+rsum END DO DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(1.6) IF( rcovar_z(iiz)>=1.E-6 ) THEN rcov_z=LOG( rcovar_z(iiz) ) ELSE rcov_z=0. ENDIF rsumc=rcov_z*riz+rsumc END DO !print *, "rsum=",rsum rw=rsumc/rsum DO iiz=1,d_z rcovar_z(iiz)=EXP( rw*ABS( LOG(REAL(iiz))- & LOG(REAL(iz)) )**1.6 ) END DO !rcovar_z(1:d_z)=rcovar_z(1:d_z)/rcovar_z(iz) c print *,"iz=",iz c print *, rcovar_z !stop RETURN END SUBROUTINE estima_z_geo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_z_hum(iz,rcovar_z,d_z) IMPLICIT NONE INTEGER :: iiz INTEGER, INTENT(IN) :: iz,d_z REAL :: riz,rcov_z,s=0. REAL :: rsumc=0.,rw REAL :: rsum=0. REAL, DIMENSION(d_z) :: rcovar_z IF ( d_z==1 ) RETURN DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(2.*1.6) rsum=riz+rsum END DO DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(1.6) IF( rcovar_z(iiz)>=1.E-6 ) THEN rcov_z=LOG( rcovar_z(iiz) ) ELSE rcov_z=0. ENDIF rsumc=rcov_z*riz+rsumc END DO !print *, "rsum=",rsum rw=rsumc/rsum DO iiz=1,d_z rcovar_z(iiz)=EXP( rw*ABS( LOG(REAL(iiz))- & LOG(REAL(iz)) )**1.6 ) END DO !rcovar_z(1:d_z)=rcovar_z(1:d_z)/rcovar_z(iz) c print *,"iz=",iz c print *, rcovar_z !stop RETURN END SUBROUTINE estima_z_hum cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_z_slp(iz,rcovar_z,d_z) IMPLICIT NONE INTEGER :: iiz INTEGER, INTENT(IN) :: iz,d_z REAL :: riz,rcov_z,s=0. REAL :: rsumc=0.,rw REAL :: rsum=0. REAL, DIMENSION(d_z) :: rcovar_z IF ( d_z==1 ) RETURN DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(2.*1.6) rsum=riz+rsum END DO DO iiz=1,d_z riz=ABS( LOG(REAL(iiz))-LOG(REAL(iz)) )**(1.6) IF( rcovar_z(iiz)>=1.E-6 ) THEN rcov_z=LOG( rcovar_z(iiz) ) ELSE rcov_z=0. ENDIF rsumc=rcov_z*riz+rsumc END DO !print *, "rsum=",rsum rw=rsumc/rsum DO iiz=1,d_z rcovar_z(iiz)=EXP( rw*ABS( LOG(REAL(iiz))- & LOG(REAL(iz)) )**1.6 ) END DO !rcovar_z(1:d_z)=rcovar_z(1:d_z)/rcovar_z(iz) c print *,"iz=",iz c print *, rcovar_z !stop RETURN END SUBROUTINE estima_z_slp cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE escreve_sigma_geo(rsigma,d_x,d_y,d_z) IMPLICIT NONE INTEGER, INTENT(IN) :: d_x,d_y,d_z REAL, DIMENSION( d_x,d_y,d_z ), INTENT(IN) :: rsigma INTEGER, DIMENSION( 22608,d_z ) :: bsigma REAL, DIMENSION( 22608,d_z ) :: asigma INTEGER :: ix,iy,iz INTEGER :: ld,lb,ilev,iline,LevBeg,Lev0,LevEnd,irec INTEGER :: Nlon=144,Nlat=157,Nlevels=19,Nlines=2 INTEGER :: numero=136 INTEGER, DIMENSION(19) :: Plevels=(/ 1000,925,900,850,800,750,700, & 650,600,550,500,450,400, & 350,300,250,200,150,100 /) INTEGER :: Hour=00 REAL :: Xmin=277.,DeltaX=0.399994,Ymin=-50.,DeltaY=0.399998 REAL :: Amiss=1.E20 LOGICAL :: back c INTEGER :: im,jm,km c REAL, DIMENSION(im,jm) :: slp_stdv !SLP f error stdv REAL, DIMENSION(d_x,d_y,d_z) :: height_stdv !geop height f error stdv c REAL, DIMENSION(im,jm,km) :: mixr_stdv !mixing ratio f error stdv CHARACTER (LEN=12) :: Ctlfile,Binfile CHARACTER (LEN=9) :: Ascfile CHARACTER (LEN=8) :: lug,DateEnd='1nov2003' c CHARACTER (LEN=12) :: luavail bsigma=0 height_stdv=rsigma lug='geosigma' Ctlfile=lug//'.ctl' open(UNIT=10,FILE=Ctlfile,STATUS='unknown') WRITE(UNIT=10,FMT=11) '^'//lug//'.grb' 11 format('DSET ',a) WRITE(UNIT=10,FMT=12) 12 format('OPTIONS sequential big-endian') WRITE(UNIT=10,FMT=13) 13 format('TITLE Forecast error "standard deviations" for PSAS') WRITE(UNIT=10,FMT=14)Amiss 14 format('UNDEF',e12.5) WRITE(UNIT=10,FMT=1) 1 FORMAT('DTYPE grib') WRITE(UNIT=10,FMT=2) '^'//Binfile//'.gmp' 2 FORMAT('index ',a) WRITE(UNIT=10,FMT=15) 'X',Nlon,Xmin,DeltaX WRITE(UNIT=10,FMT=15) 'Y',Nlat,Ymin,DeltaY 15 FORMAT(a1,'DEF ',i3,' LINEAR ',f7.1,1x,f9.6) WRITE(UNIT=10,FMT=16)Nlevels 16 FORMAT('ZDEF ',i2,'LEVELS') DO 20, iline=1,Nlines Lev0=10*(iline-1) LevBeg=Lev0+1 LevEnd=min(Nlevels,Lev0+10) write(UNIT=10,FMT=21)(PLevels(ilev),ilev=LevBeg, LevEnd) 20 END DO 21 FORMAT(10I4) write(UNIT=10,FMT=23)Hour, DateEnd 23 format('TDEF 1 linear ',i2,'Z',a9,' 6hr') WRITE(UNIT=10,FMT=4) 4 format('vars 58') write(UNIT=10,FMT=27) numero 27 format('psml 0, ',i3,', 0, 0 Sea Level Pressure', & ' forecast error standard deviation') write(UNIT=10,FMT=25) Nlevels 25 format('zgeo 0, ',i3,', 0, 0 Geopotential Height', & ' forecast error standard deviation') write(UNIT=10,FMT=26) Nlevels 26 format('umrl 0, ',i3,', 0, 0 Mixing Ratio', & ' forecast error standard deviation') write(UNIT=10,FMT=28) 28 format('ENDVARS') close(UNIT=10) Binfile=lug//'.bin' Ascfile='sigma.asc' c open(UNIT=11,FILE=Binfile,STATUS='unknown') open(UNIT=12,FILE=Ascfile,STATUS='unknown') DO ilev=1,Nlevels DO ix=1,d_x DO iy=1,d_y irec=(ix-1)*d_y+iy !irec=(ix-1)*5+iy bsigma(irec,ilev)=INT(rsigma(ix,iy,ilev)*100) asigma(irec,ilev)=rsigma(ix,iy,ilev) END DO END DO c WRITE(UNIT=11,FMT=29)(bsigma(irec,ilev), irec=1,22608) WRITE(UNIT=12,FMT=30)(asigma(irec,ilev), irec=1,10) END DO 29 FORMAT(22608(1x,B15)) 30 FORMAT(10(1x,E10.4)) c CLOSE(UNIT=11) CLOSE(UNIT=12) RETURN END SUBROUTINE escreve_sigma_geo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE escreve_sigma_hum(rsigma,d_x,d_y,d_z) IMPLICIT NONE INTEGER, INTENT(IN) :: d_x,d_y,d_z REAL, DIMENSION( d_x,d_y,d_z ), INTENT(IN) :: rsigma INTEGER, DIMENSION( 22608,d_z ) :: bsigma REAL, DIMENSION( 22608,d_z ) :: asigma INTEGER :: ix,iy,iz INTEGER :: ld,lb,ilev,iline,LevBeg,Lev0,LevEnd,irec INTEGER :: Nlon=144,Nlat=157,Nlevels=19,Nlines=2 INTEGER :: numero=136 INTEGER, DIMENSION(19) :: Plevels=(/ 1000,925,900,850,800,750,700, & 650,600,550,500,450,400, & 350,300,250,200,150,100 /) INTEGER :: Hour=00 REAL :: Xmin=277.,DeltaX=0.399994,Ymin=-50.,DeltaY=0.399998 REAL :: Amiss=1.E20 LOGICAL :: back c INTEGER :: im,jm,km c REAL, DIMENSION(im,jm) :: slp_stdv !SLP f error stdv REAL, DIMENSION(d_x,d_y,d_z) :: height_stdv !geop height f error stdv c REAL, DIMENSION(im,jm,km) :: mixr_stdv !mixing ratio f error stdv CHARACTER (LEN=12) :: Ctlfile,Binfile CHARACTER (LEN=9) :: Ascfile CHARACTER (LEN=8) :: lug,DateEnd='1nov2003' c CHARACTER (LEN=12) :: luavail bsigma=0 height_stdv=rsigma lug='humsigma' Ctlfile=lug//'.ctl' open(UNIT=10,FILE=Ctlfile,STATUS='unknown') WRITE(UNIT=10,FMT=11) '^'//lug//'.grb' 11 format('DSET ',a) WRITE(UNIT=10,FMT=12) 12 format('OPTIONS sequential big-endian') WRITE(UNIT=10,FMT=13) 13 format('TITLE Forecast error "standard deviations" for PSAS') WRITE(UNIT=10,FMT=14)Amiss 14 format('UNDEF',e12.5) WRITE(UNIT=10,FMT=1) 1 FORMAT('DTYPE grib') WRITE(UNIT=10,FMT=2) '^'//Binfile//'.gmp' 2 FORMAT('index ',a) WRITE(UNIT=10,FMT=15) 'X',Nlon,Xmin,DeltaX WRITE(UNIT=10,FMT=15) 'Y',Nlat,Ymin,DeltaY 15 FORMAT(a1,'DEF ',i3,' LINEAR ',f7.1,1x,f9.6) WRITE(UNIT=10,FMT=16)Nlevels 16 FORMAT('ZDEF ',i2,'LEVELS') DO 20, iline=1,Nlines Lev0=10*(iline-1) LevBeg=Lev0+1 LevEnd=min(Nlevels,Lev0+10) write(UNIT=10,FMT=21)(PLevels(ilev),ilev=LevBeg, LevEnd) 20 END DO 21 FORMAT(10I4) write(UNIT=10,FMT=23)Hour, DateEnd 23 format('TDEF 1 linear ',i2,'Z',a9,' 6hr') WRITE(UNIT=10,FMT=4) 4 format('vars 58') write(UNIT=10,FMT=27) numero 27 format('psml 0, ',i3,', 0, 0 Sea Level Pressure', & ' forecast error standard deviation') write(UNIT=10,FMT=25) Nlevels 25 format('zgeo 0, ',i3,', 0, 0 Geopotential Height', & ' forecast error standard deviation') write(UNIT=10,FMT=26) Nlevels 26 format('umrl 0, ',i3,', 0, 0 Mixing Ratio', & ' forecast error standard deviation') write(UNIT=10,FMT=28) 28 format('ENDVARS') close(UNIT=10) Binfile=lug//'.bin' Ascfile='sigma.asc' c open(UNIT=11,FILE=Binfile,STATUS='unknown') open(UNIT=12,FILE=Ascfile,STATUS='unknown') DO ilev=1,Nlevels DO ix=1,d_x DO iy=1,d_y irec=(ix-1)*d_y+iy !irec=(ix-1)*5+iy bsigma(irec,ilev)=INT(rsigma(ix,iy,ilev)*100) asigma(irec,ilev)=rsigma(ix,iy,ilev) END DO END DO c WRITE(UNIT=11,FMT=29)(bsigma(irec,ilev), irec=1,22608) WRITE(UNIT=12,FMT=30)(asigma(irec,ilev), irec=1,10) END DO 29 FORMAT(22608(1x,B15)) 30 FORMAT(10(1x,E10.4)) c CLOSE(UNIT=11) CLOSE(UNIT=12) RETURN END SUBROUTINE escreve_sigma_hum cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE escreve_sigma_slp(rsigma,d_x,d_y,d_z) IMPLICIT NONE INTEGER, INTENT(IN) :: d_x,d_y,d_z REAL, DIMENSION( d_x,d_y,d_z ), INTENT(IN) :: rsigma INTEGER, DIMENSION( 22608,d_z ) :: bsigma REAL, DIMENSION( 22608,d_z ) :: asigma INTEGER :: ix,iy,iz INTEGER :: ld,lb,ilev,iline,LevBeg,Lev0,LevEnd,irec INTEGER :: Nlon=144,Nlat=157,Nlevels=19,Nlines=2 INTEGER :: numero=136 INTEGER, DIMENSION(19) :: Plevels=(/ 1000,925,900,850,800,750,700, & 650,600,550,500,450,400, & 350,300,250,200,150,100 /) INTEGER :: Hour=00 REAL :: Xmin=277.,DeltaX=0.399994,Ymin=-50.,DeltaY=0.399998 REAL :: Amiss=1.E20 LOGICAL :: back c INTEGER :: im,jm,km c REAL, DIMENSION(im,jm) :: slp_stdv !SLP f error stdv REAL, DIMENSION(d_x,d_y,d_z) :: height_stdv !geop height f error stdv c REAL, DIMENSION(im,jm,km) :: mixr_stdv !mixing ratio f error stdv CHARACTER (LEN=12) :: Ctlfile,Binfile CHARACTER (LEN=9) :: Ascfile CHARACTER (LEN=8) :: lug,DateEnd='1nov2003' c CHARACTER (LEN=12) :: luavail bsigma=0 height_stdv=rsigma lug='slpsigma' Ctlfile=lug//'.ctl' open(UNIT=10,FILE=Ctlfile,STATUS='unknown') WRITE(UNIT=10,FMT=11) '^'//lug//'.grb' 11 format('DSET ',a) WRITE(UNIT=10,FMT=12) 12 format('OPTIONS sequential big-endian') WRITE(UNIT=10,FMT=13) 13 format('TITLE Forecast error "standard deviations" for PSAS') WRITE(UNIT=10,FMT=14)Amiss 14 format('UNDEF',e12.5) WRITE(UNIT=10,FMT=1) 1 FORMAT('DTYPE grib') WRITE(UNIT=10,FMT=2) '^'//Binfile//'.gmp' 2 FORMAT('index ',a) WRITE(UNIT=10,FMT=15) 'X',Nlon,Xmin,DeltaX WRITE(UNIT=10,FMT=15) 'Y',Nlat,Ymin,DeltaY 15 FORMAT(a1,'DEF ',i3,' LINEAR ',f7.1,1x,f9.6) WRITE(UNIT=10,FMT=16)Nlevels 16 FORMAT('ZDEF ',i2,'LEVELS') DO 20, iline=1,Nlines Lev0=10*(iline-1) LevBeg=Lev0+1 LevEnd=min(Nlevels,Lev0+10) write(UNIT=10,FMT=21)(PLevels(ilev),ilev=LevBeg, LevEnd) 20 END DO 21 FORMAT(10I4) write(UNIT=10,FMT=23)Hour, DateEnd 23 format('TDEF 1 linear ',i2,'Z',a9,' 6hr') WRITE(UNIT=10,FMT=4) 4 format('vars 58') write(UNIT=10,FMT=27) numero 27 format('psml 0, ',i3,', 0, 0 Sea Level Pressure', & ' forecast error standard deviation') write(UNIT=10,FMT=25) Nlevels 25 format('zgeo 0, ',i3,', 0, 0 Geopotential Height', & ' forecast error standard deviation') write(UNIT=10,FMT=26) Nlevels 26 format('umrl 0, ',i3,', 0, 0 Mixing Ratio', & ' forecast error standard deviation') write(UNIT=10,FMT=28) 28 format('ENDVARS') close(UNIT=10) Binfile=lug//'.bin' Ascfile='sigma.asc' c open(UNIT=11,FILE=Binfile,STATUS='unknown') open(UNIT=12,FILE=Ascfile,STATUS='unknown') DO ilev=1,Nlevels DO ix=1,d_x DO iy=1,d_y irec=(ix-1)*d_y+iy !irec=(ix-1)*5+iy bsigma(irec,ilev)=INT(rsigma(ix,iy,ilev)*100) asigma(irec,ilev)=rsigma(ix,iy,ilev) END DO END DO c WRITE(UNIT=11,FMT=29)(bsigma(irec,ilev), irec=1,22608) WRITE(UNIT=12,FMT=30)(asigma(irec,ilev), irec=1,10) END DO 29 FORMAT(22608(1x,B15)) 30 FORMAT(22608(1x,E12.5)) c CLOSE(UNIT=11) CLOSE(UNIT=12) RETURN END SUBROUTINE escreve_sigma_slp cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION luavail() IMPLICIT NONE CHARACTER (LEN=8) :: luavail luavail='geosigma' END FUNCTION luavail cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_sigma_geo(rsigma,rsigma_e,d_z) IMPLICIT NONE LOGICAL :: Singular INTEGER :: iz,i,j,d_z_v=0 INTEGER :: n INTEGER, INTENT(IN) :: d_z REAL :: rsigma_max=-1.,riz,rij,ri,rsigma_e_max=-1. REAL, DIMENSION(:), ALLOCATABLE :: rw,rsumc REAL, DIMENSION(:,:), ALLOCATABLE :: rsum REAL, DIMENSION(d_z), INTENT(IN) :: rsigma REAL, DIMENSION(d_z) :: rsigman REAL, DIMENSION(d_z), INTENT(OUT) :: rsigma_e DO iz=1,d_z IF( rsigma(iz)>rsigma_max ) THEN rsigma_max=rsigma(iz) n=iz ENDIF END DO rsigman=rsigma/rsigma_max DO iz=1,d_z IF( rsigman(iz)>=0.1 ) d_z_v=d_z_v+1 END DO IF(d_z_v==d_z) THEN rsigma_e=rsigma RETURN ENDIF print *, "d_z_v=",d_z_v," sera necessario parametrizar" ALLOCATE( rsumc(d_z_v),rw(0:(d_z_v-1)),rsum(d_z_v,d_z_v) ) rsum=0. rsumc=0. DO i=0,(d_z_v-1) DO j=i,(d_z_v-1) DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) rij=REAL(i+j) ri=REAL(i) rsum(i+1,j+1)=rsum(i+1,j+1)+riz**rij END IF END DO END DO DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) ri=REAL(i) rsumc(i+1)=rsumc(i+1)+rsigman(iz)*(riz)**ri END IF END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do DO i=2,(d_z_v-1) DO j=1,i rsum(i,j)=rsum(j,i) END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do c print *,"rsumc=",rsumc CALL eliminacao_gaussiana_geo(rsum,rsumc,d_z_v,rw,Singular) IF( Singular==.TRUE. ) THEN PRINT *, 'Parametrizacao incorreta' STOP END IF print *,"rw=",rw rsigma_e=0. DO iz=1,d_z riz=REAL(iz)/REAL(d_z) DO i=0,(d_z_v-1) ri=REAL(i) rsigma_e(iz)=rsigma_e(iz)+rw(i)*riz**ri END DO END DO DO iz=1,d_z IF( rsigma_e(iz)>rsigma_e_max ) THEN rsigma_e_max=rsigma_e(iz) ENDIF END DO rsigma_e=rsigma_e*rsigma_max !/rsigma_e_max c DO iz=1,d_z c print *,"iz=",iz," desvp=",rsigma(iz)," desvp_e=",rsigma_e(iz) c END DO !stop RETURN END SUBROUTINE estima_sigma_geo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_sigma_hum(rsigma,rsigma_e,d_z) IMPLICIT NONE LOGICAL :: Singular INTEGER :: iz,i,j,d_z_v=0 INTEGER :: n INTEGER, INTENT(IN) :: d_z REAL :: rsigma_max=-1.,riz,rij,ri,rsigma_e_max=-1. REAL, DIMENSION(:), ALLOCATABLE :: rw,rsumc REAL, DIMENSION(:,:), ALLOCATABLE :: rsum REAL, DIMENSION(d_z), INTENT(IN) :: rsigma REAL, DIMENSION(d_z) :: rsigman REAL, DIMENSION(d_z), INTENT(OUT) :: rsigma_e DO iz=1,d_z IF( rsigma(iz)>rsigma_max ) THEN rsigma_max=rsigma(iz) n=iz ENDIF END DO rsigman=rsigma/rsigma_max DO iz=1,d_z IF( rsigman(iz)>=0.1 ) d_z_v=d_z_v+1 END DO IF(d_z_v==d_z) THEN rsigma_e=rsigma RETURN ENDIF print *, "d_z_v=",d_z_v," sera necessario parametrizar" ALLOCATE( rsumc(d_z_v),rw(0:(d_z_v-1)),rsum(d_z_v,d_z_v) ) rsum=0. rsumc=0. DO i=0,(d_z_v-1) DO j=i,(d_z_v-1) DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) rij=REAL(i+j) ri=REAL(i) rsum(i+1,j+1)=rsum(i+1,j+1)+riz**rij END IF END DO END DO DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) ri=REAL(i) rsumc(i+1)=rsumc(i+1)+rsigman(iz)*(riz)**ri END IF END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do DO i=2,(d_z_v-1) DO j=1,i rsum(i,j)=rsum(j,i) END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do c print *,"rsumc=",rsumc CALL eliminacao_gaussiana_hum(rsum,rsumc,d_z_v,rw,Singular) IF( Singular==.TRUE. ) THEN PRINT *, 'Parametrizacao incorreta' STOP END IF print *,"rw=",rw rsigma_e=0. DO iz=1,d_z riz=REAL(iz)/REAL(d_z) DO i=0,(d_z_v-1) ri=REAL(i) rsigma_e(iz)=rsigma_e(iz)+rw(i)*riz**ri END DO END DO DO iz=1,d_z IF( rsigma_e(iz)>rsigma_e_max ) THEN rsigma_e_max=rsigma_e(iz) ENDIF END DO rsigma_e=rsigma_e*rsigma_max !/rsigma_e_max c DO iz=1,d_z c print *,"iz=",iz," desvp=",rsigma(iz)," desvp_e=",rsigma_e(iz) c END DO !stop RETURN END SUBROUTINE estima_sigma_hum cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE estima_sigma_slp(rsigma,rsigma_e,d_z) IMPLICIT NONE LOGICAL :: Singular INTEGER :: iz,i,j,d_z_v=0 INTEGER :: n INTEGER, INTENT(IN) :: d_z REAL :: rsigma_max=-1.,riz,rij,ri,rsigma_e_max=-1. REAL, DIMENSION(:), ALLOCATABLE :: rw,rsumc REAL, DIMENSION(:,:), ALLOCATABLE :: rsum REAL, DIMENSION(d_z), INTENT(IN) :: rsigma REAL, DIMENSION(d_z) :: rsigman REAL, DIMENSION(d_z), INTENT(OUT) :: rsigma_e DO iz=1,d_z IF( rsigma(iz)>rsigma_max ) THEN rsigma_max=rsigma(iz) n=iz ENDIF END DO rsigman=rsigma/rsigma_max DO iz=1,d_z IF( rsigman(iz)>=0.1 ) d_z_v=d_z_v+1 END DO IF(d_z_v==d_z) THEN rsigma_e=rsigma RETURN ENDIF print *, "d_z_v=",d_z_v," sera necessario parametrizar" ALLOCATE( rsumc(d_z_v),rw(0:(d_z_v-1)),rsum(d_z_v,d_z_v) ) rsum=0. rsumc=0. DO i=0,(d_z_v-1) DO j=i,(d_z_v-1) DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) rij=REAL(i+j) ri=REAL(i) rsum(i+1,j+1)=rsum(i+1,j+1)+riz**rij END IF END DO END DO DO iz=1,d_z IF( rsigman(iz)>=0.1 ) THEN riz=REAL(iz)/REAL(d_z) ri=REAL(i) rsumc(i+1)=rsumc(i+1)+rsigman(iz)*(riz)**ri END IF END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do DO i=2,(d_z_v-1) DO j=1,i rsum(i,j)=rsum(j,i) END DO END DO c do i=1,4 c print *,"rsum=",rsum(i,1:4) c end do c print *,"rsumc=",rsumc CALL eliminacao_gaussiana_slp(rsum,rsumc,d_z_v,rw,Singular) IF( Singular==.TRUE. ) THEN PRINT *, 'Parametrizacao incorreta' STOP END IF print *,"rw=",rw rsigma_e=0. DO iz=1,d_z riz=REAL(iz)/REAL(d_z) DO i=0,(d_z_v-1) ri=REAL(i) rsigma_e(iz)=rsigma_e(iz)+rw(i)*riz**ri END DO END DO DO iz=1,d_z IF( rsigma_e(iz)>rsigma_e_max ) THEN rsigma_e_max=rsigma_e(iz) ENDIF END DO rsigma_e=rsigma_e*rsigma_max !/rsigma_e_max c DO iz=1,d_z c print *,"iz=",iz," desvp=",rsigma(iz)," desvp_e=",rsigma_e(iz) c END DO !stop RETURN END SUBROUTINE estima_sigma_slp cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE eliminacao_gaussiana_geo(SisLin,B,N,X,Singular) IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N,N), INTENT(IN) :: SisLin REAL, DIMENSION(N), INTENT(IN) :: B REAL, DIMENSION(N,N+1) :: Matriz REAL, DIMENSION(N), INTENT(OUT) :: X LOGICAL, INTENT(OUT) :: Singular INTEGER :: I, J REAL, DIMENSION(N+1) :: Temp REAL :: PivoAbs, Multiplicador REAL, PARAMETER :: Epsilon=1E-7 INTEGER :: ColPivo Singular=.FALSE. Matriz(1:N,1:N)=SisLin(1:N,1:N) Matriz(1:N,N+1)=B(1:N) DO I=1,N PivoAbs=ABS(Matriz(I,I)) ColPivo=I DO J=I+1,N IF (ABS(Matriz(J,I))>PivoAbs) THEN PivoAbs=ABS(Matriz(J,I)) ColPivo=J END IF END DO IF (PivoAbsPivoAbs) THEN PivoAbs=ABS(Matriz(J,I)) ColPivo=J END IF END DO IF (PivoAbsPivoAbs) THEN PivoAbs=ABS(Matriz(J,I)) ColPivo=J END IF END DO IF (PivoAbs