Autor: Marcelo Barbio Rosa. Data: Set/2014 DOS TERMOS DE SUTCLIFF A teoria de Sutcliff descreve a evolução de ciclones extratropicais e foi desenvolvida na década de 1950 por Sutcliff e Pettersen (1956). Esta teoria tenta explicar, usando-se apenas as informações de superfície e do nível de 500 hPa, tradicionalmente conhecido como Nível de Não Divergência (NND), o desenvolvimento destes sistemas com base em cinco termos que são: Advecção de Vorticidade em 500 hPa (Advrt), Laplaciano do Campo de Advecção de Espessura (LAdesp), Gradiente de Vorticidade Planetária devido a Translação do Sistema (Cqel), Laplaciano devido ao Aquecimento Diabatico (Ldia), Laplaciano do Movimento Vertical Adiabático (Lsig). DADOS DE ENTRADA: Estes devem ser em pontos de grade regular com resolução espacial definida pelo parâmetros reslon e relat dados em graus. A matriz deve estar gravada de sul para norte e a latitude inicial, em graus, deve ser dada. Vento Zonal (uvel) em m/s Vento Meridional (vvel) em m/s Velocidade vertical (omeg) em Pa/s Altura geopotencial (geop) em m Temperatura do bulbo seco (tempo) em K Umidade Relativa (urel) adimensional Topografia (topo) em metros. Caso não tenha faça topo=0.0 Pressao (Pa) dos níveis padrões dos dados de entrada Undef deve ser dado, usualmente é definido como -9999.9 DeltaT é o tempo em segundos entre dois campos de dados Sup é o número do primeiro nível considerado como superfície. Normalmente = 1 Nnd é número do nível de não divergência (i.e. de 500 hPa) da matriz de entrada. ! ! ********************************************************************** SUBORUTINE SUTCLIFF (im,jm,km,undef,latin,reslat,reslon,deltaT,sup,nnd& uvel,vvel,omeg,geop,temp,urel,topo,pres,& Advrt,Lsig,LAdesp,Cqel,Ldia) IMPLICIT NONE integer,intent(in) :: im,jm,km,sup,nnd real,intent(in) :: latin,reslat,reslon,undef,deltaT real,dimension(im,jm,km),intent(in) :: uvel,vvel,omeg,geop,temp,urel real,dimension(km), intent(in) :: pres(km),topo(im,jm,1) real,dimension(im,jm,km),intent(out) :: Advrt,Lsig,LAdesp,Cqel,Ldia real, parameter :: pi=3.14159265358 !pi real, parameter :: rad=pi/180.0 !radianos real, parameter :: Lv=2.5200e6 !Calor latente condensacao(j/kg) real, parameter :: Ra=6378000.0 !Raio da terra (m) real, parameter :: Rg=287.05 !Cte do ar seco (j/kg/k) real, parameter :: grav=9.81 !acel da gravidade (m/s2) real, parameter :: Re=287.05 !Cte do ar seco (j/kg/k) real, parameter :: cp=1004.6 !calor especifico (j/kg/k) real, parameter :: gamaD=grav/cp !adiabatica seca real, parameter :: epsilon=0.622 real, parameter :: T0=273.16 real, parameter :: omg=7.2921159e-5 !velocidade angular da terra (rad/s) real, dimension(jm) :: fcori real, dimension(im,jm,1) :: pnm real, dimension(im,jm,km) :: uesp,uespSt real, dimension(im,jm,km) :: ugeo,vgeo,Raio,DeltaS,DeltaD,DeltaN,umod real, dimension(im,jm,km) :: vort1,vort2 real, dimension(im,jm,km) :: DvelS,DvelN,ang1,DdegN real, dimension(im,jm,km) :: DgeoN,DgeoS,DtmpS,tempN,tempS real, dimension(im,jm,km) :: uvelS,uvelN,vvelS,vvelN,geopN,geopS real, dimension(im,jm,km) :: var,sgg,uespN,uespS,Aduesp real, dimension(im,jm,1) :: qel,varN,varS real, dimension(im,jm,1) :: umed,Adesp,dltN,dltS,Ltemp real, dimension(im,jm,1) :: LapN,LapS,heat,sig real :: cte,cte1,cte2,ang2,res real :: lat(jm),pres(km) integer :: erro,i,j,k,h,t1,cont res=(reslat+reslon)*0.5 do j=1,jm lat(j)=latinc+reslat*float(j-1) enddo ! ************************************************************************ ! Forca de Coriolis (s^-1) do j=1,jm cte=lat(j) if (abs(cte).lt.5) cte=5*abs(lat(j))/lat(j) fcori(j)=2.0*omg*sin(cte*rad) enddo ! ************************************************************************ ! Calculo da pressao em superficie (hPa) e da Uesp (g/g) do j=1,jm do i=1,im cte=(temp(i,j,1)-temp(i,j,2))/(geop(i,j,1)-geop(i,j,2)) cte=temp(i,j,1)-0.5*cte*geop(i,j,2) pnm(i,j,1)=1e-2*pres(1)*exp(grav*geop(i,j,1)/(Re*cte)) pnm(i,j,1)=1e-2*pres(1)*exp(geop(i,j,1)/(29.2898*temp(i,j,1))) enddo enddo call Hurrel(im,jm,km,T0,Rg,epsilon,temp,pres,urel,uesp,uespSt) ! ************************************************************************ ! Calculo do modulo do vento real (m/s) do k=1,km do j=1,jm do i=1,im umod(i,j,k)=sqrt(uvel(i,j,k)**2.0+vvel(i,j,k)**2.0) enddo enddo enddo ! ************************************************************************ ! Vento geostrofico em coordenadas cartesianas (m/s) do k=1,km do j=1,jm-1 do i=1,im ugeo(i,j,k)=-grav*(geop(i,j+1,k)-geop(i,j,k))/(reslat*Ra*rad)/fcori(j) ugeo(i,jm,k)=ugeo(i,jm-1,k) enddo enddo do j=1,jm do i=1,im-1 cte=reslon*Ra*rad*cos(lat(j)*rad) vgeo(i,j,k)=grav*(geop(i+1,j,k)-geop(i,j,k))/cte/fcori(j) vgeo(im,j,k)=vgeo(im-1,j,k) enddo enddo enddo ! ************************************************************************ ! Angulo do vento tangencial (graus) ang1=0.0 do k=1,km do j=1,jm do i=1,im if(ugeo(i,j,k).ne.0.0) ang1(i,j,k)=atan2(vgeo(i,j,k),ugeo(i,j,k))/rad if(ang1(i,j,k).lt.0.0) ang1(i,j,k)=360+ang1(i,j,k) enddo enddo enddo ! ************************************************************************ ! Calculo das Coordenadas Naturais call normal (im,jm,km,ang1,undef,uvel,uvelN,uvelS) call normal (im,jm,km,ang1,undef,vvel,vvelN,vvelS) do k=1,km do j=1,jm do i=1,im DvelN(i,j,k)=sqrt(uvelN(i,j,k)**2.0+vvelN(i,j,k)**2.0)-umod(i,j,k) enddo enddo enddo call normal (im,jm,km,ang1,undef,ugeo,uvelN,uvelS) call normal (im,jm,km,ang1,undef,vgeo,vvelN,vvelS) do k=1,km do j=1,jm do i=1,im if(uvelS(i,j,k).ne.0.0) ang2=atan2(vvelS(i,j,k),uvelS(i,j,k))/rad !graus if(ang2.lt.0.0) ang2=360+ang2 DeltaD(i,j,k)=ang2-ang1(i,j,k) enddo enddo enddo call normal (im,jm,km,ang1,undef,geop,geopN,geopS) do k=1,km do j=1,jm do i=1,im DgeoN(i,j,k)=geopN(i,j,k)-geop(i,j,k) DgeoS(i,j,k)=geopS(i,j,k)-geop(i,j,k) enddo enddo enddo DeltaS=0.0 DeltaN=0.0 do k=1,km do j=1,jm do i=1,im if(abs(cos(ang1(i,j,k)*rad)).ge.0.5*sqrt(2.0))then DeltaS(i,j,k)=abs(reslat*Ra*rad*cos(lat(j)*rad)/cos(ang1(i,j,k)*rad)) DeltaN(i,j,k)=abs(reslat*Ra*rad/cos(ang1(i,j,k)*rad)) else DeltaS(i,j,k)=abs(reslon*Ra*rad/sin(ang1(i,j,k)*rad)) DeltaN(i,j,k)=abs(reslon*Ra*rad*cos(lat(j)*rad)/sin(ang1(i,j,k)*rad)) endif enddo enddo enddo ! ************************************************************************ ! Raio de Orbita da particula (m) Raio=undef do k=1,km do j=1,jm do i=1,im Raio(i,j,k)=DeltaS(i,j,k)/(DeltaD(i,j,k)*rad) if(abs(Raio(i,j,k)*1e-3).gt.1.0e4) Raio(i,j,k)=undef enddo enddo enddo ! ************************************************************************ ! Suavizacao e preenchimento de buracos do k=1,2 t1=1 call media (im,jm,km,t1,undef,res,rad,Raio) call media (im,jm,km,t1,undef,res,rad,DeltaD) call media (im,jm,km,t1,undef,res,rad,ang1) call media (im,jm,km,t1,undef,res,rad,DvelN) call media (im,jm,km,t1,undef,res,rad,DeltaS) call media (im,jm,km,t1,undef,res,rad,DeltaN) enddo ! ************************************************************************ ! Vorticidade em Coordenadas Naturais (s^-1) vort1=undef; vort2=undef do k=1,km do j=1,jm do i=1,im if(abs(Raio(i,j,k)).gt.50.0.and.Raio(i,j,k).ne.undef) then vort1(i,j,k)=umod(i,j,k)/Raio(i,j,k) if(abs(vort1(i,j,k)).gt.1e-4) vort1(i,j,k)=undef endif if(deltaN(i,j,k).ne.undef) vort2(i,j,k)=DvelN(i,j,k)/DeltaN(i,j,k) if(abs(vort2(i,j,k)).gt.1e-4) vort2(i,j,k)=undef enddo enddo enddo ! ************************************************************************ ! Gradiente de Temperatura em Coordenadas Naturais (K/m) call normal (im,jm,km,ang1,undef,temp,tempN,tempS) do k=1,km do j=1,jm do i=1,im if(DeltaS(i,j,k).ne.undef) then DtmpS(i,j,k)=(tempS(i,j,k)-temp(i,j,k))/DeltaS(i,j,k) endif enddo enddo enddo ! ************************************************************************ ! Equacao de Desenvolvimento de Sutcliff (s^-1) ! Adveccao de Vorticidade em 500 hPa em Coordenadas Naturais qel=undef do j=1,jm do i=1,im if(vort1(i,j,nnd).ne.undef.and.vort2(i,j,nnd).ne.undef) then qel(i,j,1)=vort1(i,j,nnd)-vort2(i,j,nnd)+fcori(j) endif if(vort1(i,j,nnd).ne.undef.and.vort2(i,j,nnd).eq.undef) then qel(i,j,1)=vort1(i,j,nnd)+fcori(j) endif if(vort1(i,j,nnd).eq.undef.and.vort2(i,j,nnd).ne.undef) then qel(i,j,1)=vort2(i,j,nnd)+fcori(j) endif enddo enddo call normal (im,jm,1,ang1,undef,qel,varN,varS) Advrt=undef do j=1,jm do i=1,im if(qel(i,j,1).ne.undef) then cte=(varS(i,j,1)-qel(i,j,1))/deltaS(i,j,nnd) Advrt(i,j,1)=-umod(i,j,nnd)*cte endif enddo enddo ! ------------------------------------------------------------------------ ! Delta normal e tangencial medios entre 1000 e 500 call mediapond (im,jm,km,undef,sup,nnd,pres,deltaS,dltS) call mediapond (im,jm,km,undef,sup,nnd,pres,deltaN,dltN) ! ------------------------------------------------------------------------ ! Laplaciano do campo de adveccao de espessura call mediapond (im,jm,km,undef,sup,nnd,pres,umod,umed) ! Adveccao de Espessura Adesp=undef do j=1,jm do i=1,im Adesp(i,j,1)=-umed(i,j,1)*(DgeoS(i,j,nnd)/deltaS(i,j,nnd)- & DgeoS(i,j,sup)/deltaS(i,j,sup)) enddo enddo call laplaciano (im,jm,1,undef,ang1,dltN,dltS,Adesp,LapN,LapS) LAdesp=undef do j=1,jm do i=1,im LAdesp(i,j,1)=-grav*(LapN(i,j,1)+LapS(i,j,1))/fcori(j) enddo enddo ! ------------------------------------------------------------------------ ! Gradiente de vorticidade planetaria devido a translacao do sistema qel=undef do j=1,jm ! Adveccao de Vorticidade em 1000 hPa do i=1,im if(vort1(i,j,sup).ne.undef.and.vort2(i,j,sup).ne.undef) then qel(i,j,1)=vort1(i,j,sup)-vort2(i,j,sup)+fcori(j) endif if(vort1(i,j,sup).ne.undef.and.vort2(i,j,sup).eq.undef) then qel(i,j,1)=vort1(i,j,sup)+fcori(j) endif if(vort1(i,j,sup).eq.undef.and.vort2(i,j,sup).ne.undef) then qel(i,j,1)=vort2(i,j,sup)+fcori(j) endif enddo enddo call normal (im,jm,1,ang1,undef,qel,varN,varS) do j=1,jm ! Gradiente de vorticidade em 1000 hPa do i=1,im if(qel(i,j,1).ne.undef) then varN(i,j,1)=(varN(i,j,1)-qel(i,j,1))/deltaN(i,j,sup) varS(i,j,1)=(varS(i,j,1)-qel(i,j,1))/deltaS(i,j,sup) endif enddo enddo ! Na falta do deslocamento do real do sistema, usa-se umed call mediapond (im,jm,km,undef,sup,nnd,pres,uvel,umed) Cqel=undef do j=1,jm ! Variacao devido a translacao do sistema em 1000hPa do i=1,im if(qel(i,j,1).ne.undef) Cqel(i,j,1)=umed(i,j,1)*(varN(i,j,1)+varS(i,j,1)) enddo enddo ! ------------------------------------------------------------------------ ! Laplaciano devido ao Aquecimento Diabatico call normal (im,jm,km,ang1,undef,uesp,uespN,uespS) Aduesp=0.0 do k=1,km ! Adveccao de Uesp na Grande Escala do j=1,jm do i=1,im if(uesp(i,j,k).ne.undef) then cte=(uespS(i,j,k)-uesp(i,j,k))/deltaS(i,j,k) Aduesp(i,j,k)=umod(i,j,k)*cte endif enddo enddo enddo ! Variacao local do Aquecimento Diabatico var=undef do k=1,km-1 do j=1,jm do i=1,im ! if(geop(i,j,k).gt.topo(i,j,1)) then cte1=uesp(i,j,k)+Aduesp(i,j,k) cte2=uesp(i,j,k+1)+Aduesp(i,j,k+1) !Meio representa a eficiencia e -0.000278 o esfriamento radiativo medio (K/s) var(i,j,k)=Lv*(0.5*(cte1-cte2)/deltaT)/cp-0.000278-umod(i,j,k)*DtmpS(i,j,k) ! endif enddo enddo enddo call mediapond (im,jm,km,undef,sup,nnd,pres,var,heat) call laplaciano (im,jm,1,undef,ang1,dltN,dltS,heat,LapN,LapS) Ldia=undef do j=1,jm do i=1,im if(heat(i,j,1).ne.undef) Ldia(i,j,1)=-Re*(LapN(i,j,1)+LapS(i,j,1))/fcori(j) enddo enddo ! ------------------------------------------------------------------------ ! Laplaciano do movimento Vertical Adiabatico ! Sigma da atmosfera media (Holton, Dynamic Meteorology, 2004, pg 51) sgg=undef do k=1,km-1 do j=1,jm do i=1,im ! if(geop(i,j,k).gt.topo(i,j,1)) then sgg(i,j,k)=-omeg(i,j,k)*(Re*temp(i,j,k)/(pres(k)*cp)-& (temp(i,j,k+1)-temp(i,j,k))/(pres(k+1)-pres(k))) ! endif enddo enddo enddo call mediapond (im,jm,km,undef,sup,nnd,pres,sgg,sig) call laplaciano (im,jm,1,undef,ang1,dltN,dltS,sig,LapN,LapS) Lsig=undef do j=1,jm do i=1,im if(sig(i,j,1).ne.undef) Lsig(i,j,1)=-Re*(LapN(i,j,1)+LapS(i,j,1))/fcori(j) enddo enddo ! ------------------------------------------------------------------------ do j=1,jm do i=1,im if(abs(Cqel(i,j,1)) .gt.9e-08) Cqel(i,j,1) =undef if(abs(LAdesp(i,j,1)).gt.9e-08) LAdesp(i,j,1)=undef if(abs(Advrt(i,j,1)) .gt.9e-08) Advrt(i,j,1) =undef if(abs(Ldia(i,j,1)) .gt.9e-08) Ldia(i,j,1) =undef if(abs(Lsig(i,j,1)) .gt.9e-08) Lsig(i,j,1) =undef enddo enddo ! ************************************************************************ ! Caso deseje suavizar os campos de saida desabilite “go to 11” go to 11 do k=1,2 t1=1 call media (im,jm,1,t1,undef,res,rad,Advrt) call media (im,jm,1,t1,undef,res,rad,Lsig) call media (im,jm,1,t1,undef,res,rad,LAdesp) call media (im,jm,1,t1,undef,res,rad,Cqel) call media (im,jm,1,t1,undef,res,rad,Ldia) enddo 11 continue END SUBROUTINE SUTCLIFF ! ************************************************************************ ! ************************************************************************ subroutine media (im,jm,tm,t1,undef,res,rad,var) IMPLICIT NONE integer,intent(in) :: im,jm,tm,t1 real,intent(in) :: undef,res,rad real,intent(inout) :: var(im,jm,tm) real :: cont,soma,var2(im,jm,tm) integer :: i1,i,j,k,t,f,g,j1,j2,f1,f2,f3,f4,ff real :: lat0,lat1,lon0,lon1,dist,cte,raio do k=1,tm do j=1,jm do i=1,im var2(i,j,k)=var(i,j,k) enddo enddo enddo raio=5 do t=1,t1 if (t.gt.1) raio=3 ! em pontos de grade do k=1,tm do j=1,jm lat0=res*float(j-1)-90 ff=raio !+int((2.0*raio-raio)*sin(abs(lat0*rad))) j1=j-ff j2=j+ff if(j1.lt.1) j1=1 if(j2.gt.jm) j2=jm do i=1,im if(t.eq.1.or.(t.gt.1.and.var(i,j,k).eq.undef)) then soma=0 cont=0 lon0=res*float(i-1) f1=i-ff f2=i f3=i+1 f4=i+ff if(f1.lt.1) then f1=im-ff f2=im f3=1 f4=i+ff endif if(f4.gt.im) then f1=i-ff f2=im f3=1 f4=ff endif do g=j1,j2 lat1=res*float(g-1)-90 do f=f1,f2 lon1=res*float(f-1) dist=acos((sin(lat1*rad)*sin(lat0*rad)+cos(lat1*rad)* & cos(lat0*rad)*cos((lon0-lon1)*rad)))/rad if(var(f,g,k).ne.undef.and.dist.le.float(ff)*res) then cte=(1.0-dist/(float(ff)*res))**3.0 soma=var(f,g,k)*cte+soma cont=cont+cte endif enddo do f=f3,f4 lon1=res*float(f-1) dist=acos((sin(lat1*rad)*sin(lat0*rad)+cos(lat1*rad)* & cos(lat0*rad)*cos((lon0-lon1)*rad)))/rad if(var(f,g,k).ne.undef.and.dist.le.float(ff)*res) then cte=(1.0-dist/(float(ff)*res))**3.0 soma=var(f,g,k)*cte+soma cont=cont+cte endif enddo enddo if(cont.gt.0) var2(i,j,k)=soma/cont endif enddo enddo enddo do k=1,tm do j=1,jm do i=1,im if(t.eq.1)then if(var(i,j,k).eq.undef) var(i,j,k)=var2(i,j,k) else var(i,j,k)=var2(i,j,k) endif enddo enddo enddo enddo end subroutine media ! ************************************************************************ subroutine normal (im,jm,km,ang1,undef,var,varN,varS) IMPLICIT NONE integer,intent(in) :: im,jm,km real,intent(in) :: undef real,intent(in ) :: ang1(im,jm,km),var(im,jm,km) real,intent(out ) :: varN(im,jm,km),varS(im,jm,km) integer :: i,j,k real :: angS,angN vars=undef varn=undef do k=1,km do j=2,jm-1 do i=2,im-1 if (var(i,j,k).ne.undef) then angS=ang1(i,j,k) angN=angS+90 if(angN.gt.360) angN=angN-360 ! Buscando para onde o fluxo vai if(angS.lt.45)then varS(i,j,k)=(var(i,j+1,k)*(angS-00)+var(i+1,j+1,k)*( 45-angS))/45.0 varN(i,j,k)=(var(i,j+1,k)*(angN-90)+var(i-1,j+1,k)*(135-angN))/45.0 endif if(angS.ge.45.and.angS.lt.90)then varS(i,j,k)=(var(i+1,j+1,k)*(angS- 45)+var(i,j+1,k)*( 90-angS))/45.0 varN(i,j,k)=(var(i-1,j+1,k)*(angN-135)+var(i-1,j,k)*(180-angN))/45.0 endif if(angS.ge.90.and.angS.lt.135)then varS(i,j,k)=(var(i,j+1,k)*(angS- 90)+var(i-1,j+1,k)*(135-angS))/45.0 varN(i,j,k)=(var(i-1,j,k)*(angN-180)+var(i-1,j-1,k)*(225-angN))/45.0 endif if(angS.ge.135.and.angS.lt.180)then varS(i,j,k)=(var(i-1,j+1,k)*(angS-135)+var(i-1,j,k)*(180-angS))/45.0 varN(i,j,k)=(var(i-1,j-1,k)*(angN-225)+var(i,j-1,k)*(270-angN))/45.0 endif if(angS.ge.180.and.angS.lt.225)then varS(i,j,k)=(var(i-1,j,k)*(angS-180)+var(i-1,j-1,k)*(225-angS))/45.0 varN(i,j,k)=(var(i,j-1,k)*(angN-270)+var(i+1,j-1,k)*(315-angN))/45.0 endif if(angS.ge.225.and.angS.lt.270)then varS(i,j,k)=(var(i-1,j-1,k)*(angS-225)+var(i,j-1,k)*(270-angS))/45.0 varN(i,j,k)=(var(i+1,j-1,k)*(angN-315)+var(i+1,j,k)*(360-angN))/45.0 endif if(angS.ge.270.and.angS.lt.315)then varS(i,j,k)=(var(i,j-1,k)*(angS-270)+var(i+1,j-1,k)*(315-angS))/45.0 varN(i,j,k)=(var(i,j+1,k)*(angN- 00)+var(i+1,j+1,k)*( 45-angN))/45.0 endif if(angS.ge.315.and.angS.le.360)then varS(i,j,k)=(var(i+1,j-1,k)*(angS-315)+var(i+1,j,k)*(360-angS))/45.0 varN(i,j,k)=(var(i+1,j+1,k)*(angN- 45)+var(i,j+1,k)*( 90-angN))/45.0 endif endif enddo enddo enddo do k=1,km do j=1,jm varN(1,j,k) =0.5*(varN(im-1,j,k)+varN(2,j,k)) varN(im,j,k)=varN(1,j,k) varS(1,j,k) =0.5*(varS(im-1,j,k)+varS(2,j,k)) varS(im,j,k)=varS(1,j,k) enddo enddo end subroutine normal ! ************************************************************************ subroutine laplaciano (im,jm,km,undef,ang1,deltaN,deltaS,var,LapN,LapS) IMPLICIT NONE integer,intent(in) :: im,jm,km real,intent(in) :: undef real,dimension(im,jm,km),intent(in ) :: var,ang1,deltaN,deltaS real,dimension(im,jm,km),intent(out) :: LapN,LapS real,dimension(im,jm,km) :: varN,varS,var2N,var2S integer :: i,j,k call normal (im,jm,km,ang1,undef,var,varN,varS) ! Gradiente do k=1,km do j=1,jm do i=1,im if(var(i,j,k).ne.undef) then varN(i,j,k)=(varN(i,j,k)-var(i,j,k))/deltaN(i,j,k) varS(i,j,k)=(varS(i,j,k)-var(i,j,k))/deltaS(i,j,k) endif enddo enddo enddo call normal (im,jm,km,ang1,undef,varN,var2N,var2S) ! Laplaciano da Normal do k=1,km do j=1,jm do i=1,im if(var(i,j,k).ne.undef) then LapN(i,j,k)=(var2N(i,j,k)-varN(i,j,k))/deltaN(i,j,k) endif enddo enddo enddo call normal (im,jm,km,ang1,undef,varS,var2N,var2S) ! Laplaciano do Tangencial do k=1,km do j=1,jm do i=1,im if(var(i,j,k).ne.undef) then LapS(i,j,k)=(var2S(i,j,1)-varS(i,j,1))/deltaS(i,j,k) endif enddo enddo enddo end subroutine Laplaciano ! ************************************************************************ subroutine mediapond (im,jm,km,undef,sup,nnd,pres,var,Mvar) IMPLICIT NONE integer,intent(in) :: im,jm,km,sup,nnd real,intent(in) :: var(im,jm,km),pres(km),undef real,intent(out) :: Mvar(im,jm,1) real :: cte,cte1 integer :: i,j,k Mvar=0.0 do j=1,jm do i=1,im cte =0.0 cte1=0.0 do k=sup,nnd if(pres(k).ge.pres(nnd).and.var(i,j,k).ne.undef) then cte =var(i,j,k)*(pres(k)-pres(k+1))+cte cte1=(pres(k)-pres(k+1))+cte1 endif enddo Mvar(i,j,1)=cte/cte1 enddo enddo end subroutine mediapond ! ***************************************************************************************** subroutine Hurrel(im,jm,km,T0,Rg,epsilon,temp,pres,urel,uesp,uespSt) implicit none integer, intent(in) :: im,jm,km ! Variaveis de Entrada real, intent(in) :: Rg,T0,epsilon real,dimension(im,jm,km),intent(in) :: temp,urel real,dimension(km), intent(in) :: pres real :: es(im,jm,km),cte integer :: i,j,k real,dimension(im,jm,km),intent(out) :: uesp,uespSt ! ************************************************************************ ! Tensao e Umidade especifica do Vapor (Pa)(Iribarne, pg 68) do k=1,km do j=1,jm do i=1,im es(i,j,k)=1e2*10**(-2937.4/temp(i,j,k)-4.9283*log10(temp(i,j,k))+23.547) cte=urel(i,j,k)*epsilon*es(i,j,k)/(pres(k)-es(i,j,k)) uesp(i,j,k)=cte/(cte+1.0) enddo enddo enddo ! Tensao e Umidade especifica do Vapor saturado (Pa)(Iribarne, pg 68) do k=1,km do j=1,jm do i=1,im cte=epsilon*es(i,j,k)/(pres(k)-es(i,j,k)) uespSt(i,j,k)=cte/(cte+1.0) enddo enddo enddo end subroutine Hurrel