CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCC CCCCCCCCCCC PROGRAMA PARA EXTRAIR VARIAS SECOES BAROCLINICAS 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 calc_secmerid_UTOT_parte_archv_ATLa008_expt020 USE IFPORT implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... DEFINE INDICE DO REGISTRO A SER LIDO NO ARQUIVO DO HYCOM c ... Modelo 1/24 metarea5 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 integer(4), parameter:: KK=185 !No de profs integer(4), parameter:: j1=1 !inicio secao integer(4), parameter:: j2=JDM ! fim secao c ==================================================================== integer(4), parameter:: jjj=j2-j1+1 integer(4), parameter:: n2drec2=((jjj*kk+4095)/4096)*4096 c ----------------------------------------------------------- c integer(4), parameter:: nrolonmax=11 integer(4), parameter:: nrolonmax=1 integer(4):: nn,isec,isec0(nrolonmax),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(jjj,kdm),lt1(jjj,kdm) real(4):: h0(idm,jdm),lon(nrolonmax) real(4):: h2(kdm),lt2(kdm),h_interp1(kk),h_interp2(jjj,kk) real(4):: prof(kk) real(4):: spval=1.e20 real(4):: bb=1/2.0001 c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! prof(1)=0 do k=2,11 prof(k)=prof(k-1)+1 enddo do k=12,17 prof(k)=prof(k-1)+2.5 enddo do k=18,22 prof(k)=prof(k-1)+5.0 enddo do k=23,26 prof(k)=prof(k-1)+10.0 enddo do k=27,44 prof(k)=prof(k-1)+20.0 enddo do k=45,108 prof(k)=prof(k-1)+25.0 enddo do k=109,kk prof(k)=prof(k-1)+50.0 enddo !expt0='ATLa0.08' expt1='GLB025t08' anoini=2004 anoend=2007 diaini=1 c ano=anoini inter=1 c --- ###################################### c --- Indices das latitudes acima descritas para o expt 1/16 lon(1)=333 !35 W isec0(1)=1080 !plon c lon(1)=44.0 !44 W c isec0(1)=902 c --- ##################################### c --- define seção inicial e final n1=1 n2=nrolonmax !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!! inicia loop latitudes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do 15 nn = n1,n2 isec=isec0(nn) do 20 ano = anoini,anoend IF (ano .EQ. 2004 .OR. ano .EQ. 2016) THEN diaend=366 ELSE diaend=365 ENDIF recnum=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... abre log de saida c write(*,*)lon(nn) write(arq_log,'("sec_i2_",i4.4,"_",a9,"_U_",i4.4,"_",i3.3, &"a",i3.3,".log")')isec,expt1,ano,diaini,diaend open(archv_log,file=arq_log,status='unknown') !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ... abre arquivo de saida ascii .b write(arq_b,'("sec_i_",i4.4,"_",a9,"_U_",i4.4,"_",i3.3,"a", &i3.3,".b")')isec,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)lon(nn),isec 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 meridional de velocidade longitude: ',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,'("sec_i2_",i4.4,"_",a9,"_U_",i4.4,"_",i3.3,"a", &i3.3,".a")')isec,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 = UBARO_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+UVEL_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. c --- para extrair seção meridional de U na grade de p devo interpolar c --- isec e isec+1 de U. do j=j1,j2 if(h(isec,j).gt.1.e10) then h1(j-j1+1,k) = spval else h1(j-j1+1,k) = ( h(isec,j) + h0(isec,j) + h(isec+1,j) + & h0(isec+1,j))*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 j=j1,j2 if(lt(isec,j).gt.1.e10) then lt1(j-j1+1,k) = spval else lt1(j-j1+1,k) = lt(isec,j) / 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 j=1,jjj do k=1,kdm h2(k) = h1(j,k) lt2(k) = lt1(j,k) enddo call interpz(h_interp1,lt2,h2,prof,kdm,KK,bb,j) do k=1,kk h_interp2(j,k) = h_interp1(k); enddo enddo ! do de j !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c salva arquivo binario com campo 2D extraido/manipulado call write2d(h_interp2,archv_idout,jjj,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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!