CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCC CCCCCCCCCCC PROGRAMA PARA EXTRAIR VARIAS SECOES DE VELOCIDADES CCCCCCCCCCC SIMULTANEAMENTE CCCCCCCCCCC (BINARIO - ACESSO DIRETO), E SALVAR EM UM BINARIO ([.a]) CCCCCCCCCCC CCCCCCCCCCC MARIELA GABIOUX (SETEMBRO DE 2011) CCCCCCCCCCC Baseado na leitura/escrita do hycom (mod_za.f) CCCCCCCCCCC e nos scripts criados pelo Prof. Afonso Paiva e Vladimir CCCCCCCCCCC Santos da Silva CCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCC CCCCCCCCCCC Versao: tira_secao_interp_vsc_mg.f CCCCCCCCCCC le arquivo para varios tempos, extrai uma secao de veloc CCCCCCCCCCC zonal interpolada na vertical e salva em 1 arquivo de saida CCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC program calcula_secao_zonal_VTOT_archv_ATLa008 USE IFPORT implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... DEFINE INDICE DO REGISTRO A SER LIDO NO ARQUIVO DO HYCOM c ... MODELO 1/12 novo atlantico (versao 2.2.18 do HYCOM) integer(4), parameter:: JDM=1100, IDM=1500, KDM=41 integer(4):: MONTG_INDEX=0, SSH_INDEX=1 integer(4):: ONETA_INDEX=2, SURFLX_INDEX=3 integer(4):: WTRFLX_INDEX=4, SALFLX_INDEX=5 integer(4):: BLDEPTH_INDEX=6, MLDEPTH_INDEX=7 c integer(4):: pk14=6 integer(4):: UBARO_INDEX=11, VBARO_INDEX=12 integer(4):: UVEL_INDEX=13, VVEL_INDEX=14, LT_INDEX=15 integer(4):: TEMP_INDEX=16, SALT_INDEX=17 integer(4):: NUM3DF=5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! integer(4), parameter:: n2drec=((idm*jdm+4095)/4096)*4096 integer(4):: k, i, j, INDEX, ios, recnum, ind c ---------------------------------------------------------- c ---- define parametros da secao a extrair CB integer(4), parameter:: KK=33 !185 !276 !33 integer(4), parameter:: i1=1 !-36W integer(4), parameter:: i2=IDM !-20W integer(4), parameter:: iii=i2-i1+1 integer(4), parameter:: n2drec2=((iii*kk+4095)/4096)*4096 c ----------------------------------------------------------- c integer(4), parameter:: nrolatmax=11 integer(4), parameter:: nrolatmax=1 !5S e 10S integer(4):: nn,jsec,jsec0(nrolatmax),n1,n2 character expt0*10,expt1*23 integer(4):: archv_idin=12, archv_idout=13, archv_log=11 integer(4)::archv_b=16 integer(4):: IJDM=IDM*JDM integer(4):: ano,dia,anoini,anoend,diaini,diaend,inter character arq_in*150, arq_out*150, arq_log*150,work*15 character arq_b*150 real(4):: h(idm,jdm),lt(idm,jdm),h1(iii,kdm),lt1(iii,kdm) real(4):: h0(idm,jdm),lat(nrolatmax) real(4):: h2(kdm),lt2(kdm),h_interp1(kk),h_interp2(iii,kk) real(4):: prof(kk) real(4):: spval=1.e20 real(4):: bb=1/2.0001 c ----------------------------------------------------------- c === niveis levitus data (prof(k),k=1,kk) &/0,10,20,30,50,75,100,125,150,200,250,300,400,500,600,700, &800,900,1000,1100,1200,1300,1400,1500,1750,2000,2500,3000,3500, &4000,4500,5000,5500/ !profundidades Levitus c ============================================= !expt0='ATL12' expt1='GLB025t08' anoini=2004 anoend=2007 C diaini=1 inter=1 c --- ###################################### c lat(1)=5 ! 5S lat(1)=-67.95 ! 10S c --- ###################################### c jsec0(1)=757 ! indice de latitude 5 S jsec0(1)=112 !indice de latitude S 10 c --- ##################################### c --- define seção inicial e final n1=1 n2=nrolatmax !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! inicia loop latitudes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c do 15 nn = 1,nrolatmax do 15 nn = n1,n2 jsec=jsec0(nn) do 20 ano = anoini,anoend IF (ano .EQ. 2012 .OR. ano .EQ. 2016) THEN diaini = 1 diaend=366 ELSE diaini = 1 diaend=365 ENDIF recnum=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... abre log de saida write(arq_log,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_", &i3.3,"a",i3.3,".log")')jsec,expt1,ano,diaini,diaend open(archv_log,file=arq_log,status='unknown') !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... abre arquivo de saida ascii .b write(arq_b,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_", &i3.3,"a",i3.3,".b")')jsec,expt1,ano,diaini,diaend open(archv_b,file=arq_b,status='unknown') !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! write(archv_b,*)'*************************** ' write(archv_b,*)' ' write(archv_b,*)' Informacoes do arquivo' write(archv_b,*)' ' write(archv_b,102)lat(nn),jsec write(archv_b,*)' ' write(archv_b,103)ano,diaini,diaend write(archv_b,*)' ' write(archv_b,104)inter write(archv_b,*)' ' write(archv_b,*)'Interpolada na grade de p' write(archv_b,*)' ' write(archv_b,105)kk write(archv_b,*)' ' write(archv_b,*)'*************************** ' 102 format(' Seção zonal de velocidade latitude: ',f7.2, &' indice ',i4.4) 103 format(' Ano: ',i4,' dia inicial: ',i3.3,' dia final: ',i3.3) 104 format(' Informação salva a cada ',i2.2,' dias') 105 format(' total de profundidades variaveis: ',i3.3) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... abre arquivo para salvar resultados (binario) write(arq_out,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_",i3.3, &"a",i3.3,".a")')jsec,expt1,ano,diaini,diaend open(archv_idout,file=arq_out,status='unknown',form='unformatted' &,recl=n2drec2,access='DIRECT',action='WRITE') write(*,'("resultados em: ",a)')arq_out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! inicia loop temporal (ex: todos os dias de 2006) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do 30 dia = diaini,diaend,inter recnum=recnum+1 ! No do registro a ser salvo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!! Experimento 1/16 sem assimilação 2005 !!!!!!!!!!!!!!!!!!!!!!! ! Base Hidrodinâmica REMO 3 write(arq_in,'(" &/pesq/share/monan/ocean_model/hycom/GLBa025/dados/ &data_",i4.4,"/ &archv.",i4.4,"_",i3.3,"_00.a")')ano,ano,dia !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c =========================================================== c =========================================================== c ... abre arquivo open(unit=archv_idin,file=arq_in,status='old',form='unformatted', &action='READ',recl=n2drec,access='direct',iostat=ios) c write(archv_log,*)'ano,dia = ',ano,dia write(archv_log,'(a)')arq_in c write(*,'(a)')arq_in !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... leitura dos campos 2D INDEX = VBARO_INDEX call read2d(h0,archv_idin,idm,jdm,n2drec,index) C ... define o campo a ser lido do 10 k=1,kdm INDEX = ((k-1)*NUM3DF+VVEL_INDEX) c ... le registro - vetor a(n2drec) call read2d(h,archv_idin,idm,jdm,n2drec,index) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... exemplo: modifica spval de terra, mantem TSM como esta c mas se for SSH ou ML_depth, ou dp, entao converte unidades c ... O importante é que neste ponto do programa voce tem uma matriz c h(idm,jdm) contendo o campo 2D lido. Faca o que quiser com ela. do i=i1,i2 if(h0(i,jsec).gt.1.e10) then h1(i-i1+1,k) = spval else h1(i-i1+1,k) = ( h0(i,jsec) + h(i,jsec) + h0(i,jsec+1) + & h(i,jsec+1))*0.5 endif enddo C ... define o campo a ser lido INDEX = ((k-1)*NUM3DF+LT_INDEX) c ... le registro - vetor a(n2drec) call read2d(lt,archv_idin,idm,jdm,n2drec,index) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do i=i1,i2 if(lt(i,jsec).gt.1.e10) then lt1(i-i1+1,k) = spval else lt1(i-i1+1,k) = lt(i,jsec) / 9806. endif enddo 10 continue !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... fazendo a interpolacao vertical para cada ponto c ... entrada: nome do arquivo de saida c ... campo 1D (pontual) de espessura de camada (lt2); c ... variavel do ponto em todas as profundidades (h2); c ... profundidades a ser interpolado (prof); c ... numero total de camadas (KDM) c ... numero de profundidades a ser interpoladas c ... bb= coef para dividir a camada, bb deve ser menor a 1/2 pois nesse caso c ... os dois pontos estariam no centro da camada. c ... Sugestao bb=1/(2.0001) c write(*,*)'interpola' do i=1,iii do k=1,kdm h2(k) = h1(i,k) lt2(k) = lt1(i,k) enddo call interpz(h_interp1,lt2,h2,prof,kdm,KK,bb,i) do k=1,kk h_interp2(i,k) = h_interp1(k); enddo enddo ! do de j !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c salva arquivo binario com campo 2D extraido/manipulado call write2d(h_interp2,archv_idout,iii,kk,n2drec2,recnum) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... fecha arquivo de leitura antes de ir para proximo tempo close(archv_idin) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! fecha loop de tempo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 30 continue c ... fecha arquivo de saida, log e b close(archv_idout) close(archv_log) close(archv_b) 20 continue !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 15 continue !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... encerra programa stop end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!! subrotina read2d !!!!!!!!!! le campo 2D (camada) em vetor a(n2drec), convert endian !!!!!!!!!! e transforma vetor a(n2drec) em matriz aa(idm,jdm), !!!!!!!!!! retornando [aa] para o programa principal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read2d(aa,archv_idin,idm,jdm,n2drec,index) implicit none integer(4):: i,j,idm,jdm,n2drec,index,ios integer(4):: archv_idin real(4):: a(n2drec),aa(idm,jdm) c ... le registro - vetor a(n2drec) read(archv_idin,rec=INDEX+1,iostat=ios)a call zaio_endian(a,n2drec) c ... transforma vetor lido em matriz 2D do j=1,jdm do i=1,idm aa(i,j) = a(i+(j-1)*idm) enddo enddo return end subroutine read2d !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!! subrotina write2d !!!!!!!!!! recebe campo 2D (camada) na matriz aa(idm,jdm), corrige !!!!!!!!!! endian, transforma no vetor a(n2drec) e salva este vetor !!!!!!!!!! no registro recnum, ou seja, salva campo 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine write2d(aa,archv_idout,idm,jdm,n2drec,recnum) implicit none integer(4):: i,j,idm,jdm,n2drec,recnum,ios integer(4):: archv_idout real(4):: a(n2drec),aa(idm,jdm) c ... transforma campo 2D em vetor a(n2drec) para salvar do j=1,jdm do i=1,idm a(i+(j-1)*idm) = aa(i,j) enddo enddo c ... corrige endian e salva arquivo binario call zaio_endian(a,n2drec) write(unit=archv_idout,rec=recnum,iostat=ios)a return end subroutine write2d !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! !!!!!!!!! zaio_endian(a,n) !!!!!!!!! - swap the endian-ness of the array. !!!!!!!!! - assumes integer(kind=1) and integer(kind=4) ocupy !!!!!!!!! one and four bytes respectively. !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine zaio_endian(a,n) implicit none integer, intent(in) :: n integer(kind=4), intent(inout) :: a(n) ! 4-bytes integer k integer(kind=4) ii4, io4 ! 4-bytes integer(kind=1) ii1(4),io1(4) ! 1-byte equivalence (ii4,ii1(1)), (io4,io1(1)) ! non-standard f90 do k= 1,n ii4 = a(k) io1(1) = ii1(4) io1(2) = ii1(3) io1(3) = ii1(2) io1(4) = ii1(1) a(k) = io4 enddo return end subroutine zaio_endian !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine interpz(dadoi,hm,dado,zz,Kt,M,bb,i) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c...function [dadoi]=interpz(hm,dado,zz,K,M,bb); c... c... Interpolacao dos dados de camadas para niveis c... c... Variaveis necessarias: c... hm = espessuras das camadas c... dado = variavel a ser interpolada c... M = n. de profundidades na saida c... Kt = total de camadas c... zz = novo eixo de profundidades; c... bb= coef para dividir a camada, bb deve ser menor a 1/2 pois nesse caso c... os dois pontos estariam no centro da camada. c... Sugestao bb=1/(2.0001) implicit none integer(4):: mm,M,kk,kt,i,j,n,k real(4):: bb real(4):: hh(kt),hm(kt),p(kt+1),z(2*kt+1),dd(2*kt+1) real(4):: dado(kt),zz(M) real(4):: dadoi(M) c ... retira valores extremos do kk=2,kt if (dado(kk).gt..5e2 .or. dado(kk).lt.-.5e2) & dado(kk)=dado(kk-1) enddo c...Calcula profundidades acumuladas p(1)=0.0 do kk=1,Kt hh(kk)=max( hm(kk),0.1 ) p(kk+1)=p(kk)+hh(kk) enddo c...localizando prof dos resultados e interpolando n=1 z(n)=0.0 dd(1)=dado(1) do kk=1,Kt n=n+1; z(n)=p(kk)+hh(kk)*bb dd(n)=dado(kk) n=n+1 z(n)=p(kk+1)-hh(kk)*bb dd(n)=dado(kk) enddo z(n+1)=p(kt+1) dd(n+1)=dd(n) if (zz(M) .gt. z(2*kt+1)) then z(2*kt+1)=zz(M)+1 endif c do k=1,2*kt+1 c if (i.eq.361) write (101,*)k,dd(k),z(k) c enddo call interp_1D(dadoi,z,dd,zz,kt,M) c do k=1,33 c if (i.eq.361) write (102,*)zz(k),dadoi(k) c enddo return end subroutine interpz cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine interp_1D(aa,z,dd,zz,kt,M) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none integer(4):: i,j,kt,M,diff real(4):: bb,cc real(4):: aa(M),dd(2*kt+1),zz(M),z(2*kt+1) do i=1,M j=1 diff = zz(i) - z(j) do while (diff .ge. 0 .and. j .le. 2*kt) j=j+1 diff = zz(i) - z(j) enddo if (j.eq.1) then aa(i) = dd(j) elseif (j.eq. 2*kt+1 .and. diff.ge.0) then aa(i) = dd(j) else bb = abs(zz(i) - z(j-1)) cc = abs(zz(i) - z(j)) aa(i) = (dd(j-1)*cc + dd(j)*bb)/(bb+cc) endif enddo return end subroutine interp_1D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!