!****************************************************************************** !* METLIB * !* * !* FUNCOES MATEMATICAS USUAIS EM METEOROLOGIA * !* * !* Copyright (C) 2007 Sergio Henrique Soares Ferreira * !* * !****************************************************************************** !* Este modulo contem uma colecao de funcoes e sub-rotinas que processam * !* calculos comuns em meteorologias, tais como conversao de direcao e velo- * !* cidade de vento em suas componente zonais e meridionais e calculo de * !* variáveis termodinamicas * !* * !* Estas funcoes e sub-rotinas levam em consideracao os intervalos em que as * !* calculos tem representatividade fisica, retornando valores nulos quando * !* um valor espurio e fornecido. * !****************************************************************************** ! HISTORICO MODULE METLIB real, parameter :: pi=3.141596 real, parameter :: Null=-340282300 private Null public interface Tvirt module procedure Tvirt_Ceusius module procedure Tvirtw_Kelvin end interface CONTAINS !----------------------------------------------------------------------------------------- ! Thickness ! Retorna a espessura em metros entre dois niveis isobariocos (Press1 e Press2) com ! uma temperatura virtural media da camada (TV) ! !--------------------------------------------------------------------------------------- Function thickness(press1, press2, tv_Kelvin);real:: thickness !{ Variaveis de entrada real,intent(in)::press1,press2 ! Pressao da base e topo da camada em (hPa) real,intent(in)::tv_kelvin !Temperatura virtual media da camada em Craus kelvin) !} !{ Variaveis locais real :: p0,p !} p0=press1 p=press2 If (p0 == 0) p0 = 0.00001 If (p == 0 ) p = 0.00001 dz = tv_kelvin * 29.3 * Log(p0 / p) thickness = dz End Function !__________________________________________________________________ ! Tvirt ! Dado a temperatura do ar (TT) ,a temperatura do ponto de orvalho (td), ! e a pressao atmosferica (p), Rerona a Temperatura virtual em graus Ceusius !------------------------------------------------------------------ Function Tvirt_Ceusius(tt, TD, p);real :: Tvirt_Ceusius !{ Variaveis de entrada real,intent(in)::tt,td ! Temperaturas do ar e pto orvalho Graus Ceusius real,intent(in):: p ! Pressao atmosferica em hPa !} If ((TD < 99).And.(p < 1100)) Then UE = 0.622 * PVAPOR(TD) / p ! Fun. umidade especifica Tvirt_Ceusius = (tt + 273.16) * (1 + 0.61 * UE) - 273.16 Else Tvirt_Ceusius = tt End If End Function ! ! Temperatura Virtual em Kelvin ! em funcao da temperatura do ar seco e da Razao de Mistura ! ! Caso seja fornecido valor negativo para umidade especidica ! a umidade sera considerada zero e a temperatura virtual igual a temperatura do ar ! ! Caso seja fornecido temperatura do ar inferior a 0 Kelvin a funcao ! retorna o valor null Function Tvirtw_Kelvin(tt, w);real :: Tvirtw_Kelvin !{ Variaveis de entrada real,intent(in)::tt ! Temperaturas do ar (kelvin) real,intent(in)::w ! Umidade especifica Kg/Kg !} rm=w if (w<0) rm=0 if (tt>0) then Tvirtw_Kelvin = tt * (1 + 0.6 * rm) else Tvirtw_Kelvin=Null end if End Function !____________________________________________________________________ !------------------------------------------------------------------------------ !qff ! `Pressao da Atmosfera Padrao da OACI ! SHSF| !------------------------------------------------------------------------------ ! ----------------------------------------------------------------------------- Function pressure_isa(z);real::pressure_isa ! (metros) !{ Variaveis de interface real,intent(in)::z ! Altitude (metros) !} !{ Variaveis locais real,Dimension(22,3):: np real::p,x1,x2 !} Np(1,1)=1013.25 ; np(1,2)=0 ; np(1,3)= 15.0 Np(2,1)=1000 ; Np(2,2)=111 ; Np(2,3)= 14.3 Np(3,1)=950 ; Np(3,2)=540 ; Np(3,3)= 11.5 Np(4,1)=900 ; Np(4,2)=988 ; Np(4,3)= 8.6 Np(5,1)=850 ; Np(5,2)=1457 ; Np(5,3)= 5.5 Np(6,1)=800 ; Np(6,2)=1949 ; Np(6,3)= 2.3 Np(7,1)=700 ; Np(7,2)=3012 ; Np(7,3)= -4.6 Np(8,1)=600 ; Np(8,2)=4206 ; Np(8,3)= -12.3 Np(9,1)=500 ; Np(9,2)=5574 ; Np(9,3)= -21.2 Np(10,1)=400 ; Np(10,2)=7185 ; Np(10,3)=-31.7 Np(11,1)=300 ; Np(11,2)=9164 ; Np(11,3)=-44.5 Np(12,1)=250 ; Np(12,2)=10363 ; Np(12,3)=-52.3 Np(13,1)=226 ; Np(13,2)=11000 ; Np(13,3)=-56.5 Np(14,1)=200 ; Np(14,2)=11784 ; Np(14,3)=-56.5 Np(15,1)=150 ; Np(15,2)=13608 ; Np(15,3)=-56.5 Np(16,1)=100 ; Np(16,2)=16180 ; Np(16,3)=-56.5 Np(17,1)=70 ; Np(17,2)=18442 ; Np(17,3)=-56.5 Np(18,1)=50 ; Np(18,2)=20576 ; Np(18,3)=-55.9 Np(19,1)=40 ; Np(19,2)=22000 ; Np(19,3)=-54.5 Np(20,1)=30 ; Np(20,2)=23849 ; Np(20,3)=-52.7 Np(21,1)=20 ; Np(21,2)=26481 ; Np(21,3)=-50.0 Np(22,1)=10 ; Np(22,2)=31055 ; Np(22,3)=-45.4 !{ Valor Invalido de PRessao if ((z<100).or.(z>100000)) then pressure_isa=null return end if If ((z >= np(1,2)).And.(p <= np(22,2))) Then do n = 1,21 If ((z >= np(n, 2)).And.(z < np(n + 1, 2))) Then if (z==np(n,2)) then p=np(n,1) exit else x1 = log(np(n, 1)) x2 = log(np(n + 1, 1)) p = exp(ITPL(np(n, 2),x1,np(n + 1, 2),x2, z)) exit ! end if End If end do Elseif (znp(22,2)) then ! Pressao acima de 10 hPa x1 = log(np(n, 1)) x2 = log(np(n + 1, 1)) p = exp(ITPL(np(n, 2),x1,np(n + 1, 2),x2, z)) End If pressure_isa = p End Function pressure_isa Function zpadrao(p);real::zpadrao ! (metros) real,intent(in)::p ! Pressao em hPa zpadrao=altitude_isa(p) end function !------------------------------------------------------------------------------ !altitude_isa ! Altitude da Atmosfera Padrao da OACI ! SHSF| !------------------------------------------------------------------------------ ! ----------------------------------------------------------------------------- Function altitude_isa(p);real::altitude_isa ! (metros) !{ Variaveis de interface real,intent(in)::p ! Pressao atmosferica (hPa) !} !{ Variaveis locais real,Dimension(22,3):: np real::x1,x2,z !} Np(1,1)=1013.25 ; np(1,2)=0 ; np(1,3)= 15.0 Np(2,1)=1000 ; Np(2,2)=111 ; Np(2,3)= 14.3 Np(3,1)=950 ; Np(3,2)=540 ; Np(3,3)= 11.5 Np(4,1)=900 ; Np(4,2)=988 ; Np(4,3)= 8.6 Np(5,1)=850 ; Np(5,2)=1457 ; Np(5,3)= 5.5 Np(6,1)=800 ; Np(6,2)=1949 ; Np(6,3)= 2.3 Np(7,1)=700 ; Np(7,2)=3012 ; Np(7,3)= -4.6 Np(8,1)=600 ; Np(8,2)=4206 ; Np(8,3)= -12.3 Np(9,1)=500 ; Np(9,2)=5574 ; Np(9,3)= -21.2 Np(10,1)=400 ; Np(10,2)=7185 ; Np(10,3)=-31.7 Np(11,1)=300 ; Np(11,2)=9164 ; Np(11,3)=-44.5 Np(12,1)=250 ; Np(12,2)=10363 ; Np(12,3)=-52.3 Np(13,1)=226 ; Np(13,2)=11000 ; Np(13,3)=-56.5 Np(14,1)=200 ; Np(14,2)=11784 ; Np(14,3)=-56.5 Np(15,1)=150 ; Np(15,2)=13608 ; Np(15,3)=-56.5 Np(16,1)=100 ; Np(16,2)=16180 ; Np(16,3)=-56.5 Np(17,1)=70 ; Np(17,2)=18442 ; Np(17,3)=-56.5 Np(18,1)=50 ; Np(18,2)=20576 ; Np(18,3)=-55.9 Np(19,1)=40 ; Np(19,2)=22000 ; Np(19,3)=-54.5 Np(20,1)=30 ; Np(20,2)=23849 ; Np(20,3)=-52.7 Np(21,1)=20 ; Np(21,2)=26481 ; Np(21,3)=-50.0 Np(22,1)=10 ; Np(22,2)=31055 ; Np(22,3)=-45.4 !{ Valor Invalido de PRessao if ((p<0).or.(p>1100)) then altitude_isa=null return end if If ((p <= np(1,1)).And.(p >= np(22,1))) Then do n = 1,21 If ((p <= np(n, 1)).And.(p > np(n + 1, 1))) Then if(p==np(n,1)) then z=np(n,2) exit else x1 = Log(np(n, 1)) x2 = Log(np(n + 1, 1)) z = ITPL(x1, np(n, 2), x2, np(n + 1, 2), Log(p)) exit ! end if End If end do Elseif (p>np(1,1)) then ! Pressao abaixo do nivel do mar x1 = Log(np(1, 1)) x2 = Log(np(2, 1)) z = ITPL(x1, np(n, 2), x2, np(n + 1, 2), Log(p)) elseif (p= MAXY).Or.(Y1 >= MAXY)) Then ITLINEAR = NULY Else ITLINEAR = (Y2 * (IT - x1) + Y1 * (x2 - IT)) / (x2 - x1) End If End Function !___________________________________________________________________ !------------------------------------------------------------------- Function PVAPOR(tt);real::PVAPOR real,intent(in)::tt If (tt > -100) Then PVAPOR = 6.1078 * 10 ** (7.5 * tt / (237.3 + tt))! ' Def. fun tensao do vapor AGUA Else PVAPOR = 0 End If End Function !__________________________________________________________________ !------------------------------------------------------------------ Function RMIST(E, p);real::RMIST ! Razao de Mistura saturada em g/Kg real,intent(in)::E ! Tensao do Valor (hPa) real,intent(in)::P ! Pressao Atmosferica (hPa_ RMIST = 622 * E / (p - E) End Function !__________________________________________________________________ !------------------------------------------------------------------ Function q_humid1(E, P);real::q_humid1 !Umidade Especifica em g/Kg real,intent(in)::E ! Tensao do Vapor (hPa) real,intent(in)::P ! Pressao atmosferica (hPa) q_humid1 = 622 * E / (P - 0.378*E) end Function !___________________________________________________________________ !------------------------------------------------------------------ Function Rseca(tt, PP1, PP2);REAL::Rseca real,intent(in)::tt,pp1,pp2 Rseca = (273.16 + tt) * ((PP2 / PP1) ** (2.0 / 7.0)) - 273.16 !: Razao adiab. SECA End Function !_________________________________________________________________ !----------------------------------------------------------------- Function Rumida( pl, plt, tt);real :: rumida !------------------------------------------------ ! SUB ROTINA P/ CALCULO DA RAZAO ADIABATICA UMIDA !----------------------------------------------- ! ** interface ** REAL,INTENT(IN):: PL,PLT! (PL=PRESS NO NIV. ZR PLT=PRES. NIV. Z) REAL,INTENT(IN):: TT ! TT=TEMP. FORNEC. (C) ! ** local variables ** REAL,DIMENSION(3)::T,TEC,DTE! T=TEMP CALCULADA real:: RCP,CP,ALCP,DT2,DT3,dz,es,r,TE ! ** init variables ** RCP = 287.0 / 1004.0 CP = 0.238 alcp = 2500000.0 / 1004.0 DT2 = 0.0 DT3 = 2.0 dz = (altitude_isa(plt) - altitude_isa(pl)) / 1000.0 !Estimativa da espessura da camada TR = tt + 273.16 ! TEMPERATURA NO NIVEL DE REFERENCIA EM KELVIN !-------------------------------------------------- ! CALCULO DA TEMPERATURA POTENCIAL NO PRIMEIRO NIVEL !---------------------------------------------------- TTO = 273.16 / TR ES = 6.11 * (TTO ** 5.31) * Exp(25.22 * (1 - TTO)) r = (622 * ES / (pl - ES)) / 1000.0 TE = TR * ((1000.0 / (pl - ES)) ** RCP) * Exp(alcp * r / TR) !---------------------------------------------------- ! PRIMEIRA E SEGUNDA ESTIMATIVA (RAZAO DE 6.5 GRAUS /KM) !----------------------------------------------------- T(1) = TR - DZ * 6.5 / 1000.0 T(2) = T(1) + 1 DO i = 1,2 TTO = 273.16 / T(i) ES = 6.11 * TTO ** 5.31 * Exp(25.22 * (1 - TTO)) r = (622 * ES / (plt - ES)) / 1000.0 TEC(i) = T(i) * (1000 / (plt - ES)) ** RCP * Exp(alcp * r / T(i)) DTE(i) = TEC(i) - TE END DO !------------------------------------------------ ! TERCEIRA ESTIMATIVA - POR PROPORCIONALIDADE !------------------------------------------------ T(3) = (DTE(1) * T(2) - DTE(2) * T(1)) / (DTE(1) - DTE(2)) 621 TTO = 273.16 / T(3) ES = 6.11 * TTO ** 5.31 * Exp(25.22 * (1 - TTO)) r = (622 * ES / (plt - ES)) / 1000.0 TEC(3) = T(3) * (1000.0 / (plt - ES)) ** RCP * Exp(alcp * r / T(3)) DTE(3) = TEC(3) - TE !-------------------------------------------- ! FIM DA TERCEIRA ESTIMATIVA T(3) ! VERIFCACAO E APURACAO DA PRECISAO DE CALCULO !--------------------------------------------- If (DTE(3) * DT2 < 0) DT3 = DT3 + 1 DT2 = DTE(3) DIFER = Int(Abs(DTE(3) * 10) + 0.5) If (DIFER > 0) Then T(3) = T(3) - (DTE(3) / DT3) GoTo 621 End If Rumida = (T(3) - 273.16) End Function !---------------------------------------------------- ! TDVUV E TUVDV TRANSFORMACAO DE DIRECAO E VELOCIDADE ! PARA COMPONENTES X E Y E VICE-VERSA !--------------------------------------------------- SUBROUTINE TDVUV(DIR,VEL,CX,CY) real,intent(in)::DIR,VEL real,intent(out)::CX,CY CX = -VEL * SIN(DIR*pi/180) CY = -VEL * COS(DIR*pi/180) END subroutine !______________________________________________________ SUBROUTINE TUVDV(U,V,DIR,VEL) real,intent(in)::U,V ! componentes zonal e meridional real,intent(out)::DIR,VEL ! Direcao e velecidade VEL = SQRT(U ** 2 + V ** 2) if (vel/=0.0) then cddd=-V/VEL IF(CDDD==1.0) DIR=0.0 IF(CDDD==-1.0) DIR=180.0 IF((CDDD/=1.0).AND.(CDDD/=-1.0)) THEN DIR=ACOS(CDDD)*180/pi IF (-U<0) DIR=360.0-DIR IF (DIR<=0) DIR=360+DIR END IF ELSE DIR=0.0 END IF END subroutine !_______________________________________________________ ! ! funcao para interpolacao linear de um ponto x em relacao a ! dois outros pontos (x1,y1)-(x2,y2) Function ITPL(x1, Y1, x2, Y2, x);real :: ITPL real,intent(in)::x1,y1,x2,y2,x real ::ddd ddd = x2 - x1 If (ddd /= 0.0) Then ITPL = (Y2 * (x - x1) + Y1 * (x2 - x)) / (x2 - x1) Else ITPL = Y1 End If End Function ITPL !______________________________________________________ end module