Skip to content

Instantly share code, notes, and snippets.

@imm

imm/MW_pkkis_L.f

Created May 29, 2010 16:43
Show Gist options
  • Save imm/418362 to your computer and use it in GitHub Desktop.
Save imm/418362 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile2.f90
real z_, T_, P_, v_, q_, W_
open(2, file = 'model_afgl86_midlat_sum.dat')
read (2, *)
read (2, *)
read (2, *)
read (2, *)
do i = 1, 201
read (2,*) z_, P_, T_, q_, v_, W_
print *, alpha_O3(v_,T_,P_,110.0)
end do
close(2)
Raw
MW_pkkis_L.f
C The subroutine of the calculation of the absorbtion coefficient in oxygen
C
C SUBROUTINE O2(MH,F,H,RO,P,T,G1,TTT1,SAZ)
C
C
C
C ppp - давление, мбар
C qqq - абсолютная влажность, г/м**3
C TTT - температура, К
C f - частота, ГГц
C PKKIS - коэффициент поглощения в км-1
REAL FUNCTION pkkis(ppp,QQQ,TTT,F)
PARAMETER (cv=29.97924)
PARAMETER (MH=1)
C ---------------
real*8 wn, F0V(30)
COMMON /ABSLiebe/ F0(44), AK1(44), AK2(44), AK3(44),
* AK4(44), AK5(44), AK6(44),F0V,BK1(30),
* BK2(30), BK3(30), BK4(30), BK5(30), BK6(30)
C
C
C
REAL H(MH),T(MH),P(MH),RO(MH),G1(MH)
REAL*8 F
REAL TTT1
C
C F0, AK1-6 - ready coefficients for calculation
C Input parameters - MH - the number of the atmosphere arrays elements
C T(MH)-Temperature, K; P(MH)-Pressure, mb(,or gPa,or 10**(-2)Pa)
C RO(MH) - absolute humidity, g/(m**3)
C F - current frequency, GHz
C Output parametrer - G1(MH)- absorbtion coefficientin oxygen along the height, Neper
C
C F=wn*cv
C RO=qqq/6.02205E23*18.*1.E6
C Из ppmv в абс. влажность г/м3
C RO(1)=qqq/1.E4/8.31441/ttt*18.*ppp
C Из г/г в абс. влажность г/м3
C RO(1)=qqq*1.E6*28./18./1.E4/8.31441/ttt*18.*ppp
open(61, file='Tb_ini_Liebe_coeff.dat')
read (61,'(7f10.6)') (F0(j),j=1,44)
read (61,'(9f8.2)') (AK1(j),j=1,44)
read (61,'(9f8.3)') (AK2(j),j=1,44)
read (61,'(9f8.2)') (AK3(j),j=1,44)
read (61,'(9f8.1)') (AK4(j),j=1,44)
read (61,'(9f8.3)') (AK5(j),j=1,44)
read (61,'(9f8.3)') (AK6(j),j=1,44)
read (61,'(7f10.6)') (F0V(j),j=1,30)
read (61,'(9f9.4)') (BK1(j),j=1,30)
read (61,'(9f8.3)') (BK2(j),j=1,30)
read (61,'(9f8.2)') (BK3(j),j=1,30)
read (61,'(9f8.2)') (BK4(j),j=1,30)
read (61,'(9f8.2)') (BK5(j),j=1,30)
read (61,'(9f8.2)') (BK6(j),j=1,30)
close (61)
RO(1)=qqq
T(1)=TTT
P(1)=ppp
saz=0.
H(1)=10.
DO J = 1,MH ! cycle for atmospheric levels
C
TETA = 300./T(J)
E = RO(J)/(7.223*TETA)
PE = 0.1*P(J)- E
C
G = 0.
C Resonant coefficient:
DO 1311 I = 1, 44
SI = AK1(I)*0.000001*PE*TETA**3*EXP(AK2(I)*(1.- TETA))
GAI = AK3(I)*0.001*(PE*TETA**(0.8-AK4(I))+1.1*E*TETA)
DELI = (AK5(I)+AK6(I)*TETA)*0.001*PE*TETA**0.8
AI = GAI*F/F0(I)
XXI = (F0(I)-F)**2+GAI**2
YYI = (F0(I)+F)**2+GAI**2
FI = AI/XXI+AI/YYI-DELI*F/F0(I)*((F0(I)-F)/XXI+(F0(I)+F)/YYI)
1311 G = G+SI*FI
C
C NonResonant coefficient:
SD = 6.14*0.0001*PE*TETA**2
GAD = 5.6*0.001*(PE+1.1*E)*TETA
IF (GAD.EQ.0.) THEN
GNR = 0.
GO TO 1313
END IF
C
AP = 1.40*(1-1.2*0.00001*F**1.5)*10.**(-10)
GNR = SD*F/(GAD*(1+(F/GAD)**2))+AP*F*PE**2*TETA**3.5
C
1313 G = G+GNR
G = 0.1820*F*G
C Converting into Neper:
G1(J) =G/4.34
pkkis=G1(J)
C
C
END DO ! cycle for atmospheric levels
C
C Optical depth:
C
C TTT1=0
C
C DO J=1,(MH-1)
C TTT1=TTT1+G1(J+1)*(H(J+1)-H(J))*SAZ
C IF(TTT1.GE.20.) GO TO 98
C END DO
C
C 98 CONTINUE
C
C
RETURN
END
C The End of the Sub O2(...) - Absorbtion in Oxygen
Raw
MW_pkpar_L.f
C___________________________________________________
C
C
C
C The subroutine of the calculation of the absorbtion coefficient in Water Vapour
C
C SUBROUTINE PAR(MH,F,H,RO,P,T,G2,TTT2,SAZ)
C
C ppp - давление, мбар
C qqq - абсолютная влажность, г/м**3
C TTT - температура, К
C f - частота, ГГц
C PKPAR - коэффициент поглощения в км-1
C
FUNCTION PKPAR(ppp,qqq,TTT,f)
parameter (mh=1)
real*8 wn, f, F0V(30)
C ---------------
COMMON /ABSLiebe/ F0(44), AK1(44), AK2(44), AK3(44),
* AK4(44), AK5(44), AK6(44),F0V,BK1(30),
* BK2(30), BK3(30), BK4(30), BK5(30), BK6(30)
C COMMON /absorption/ F0(44), AK1(44), AK2(44), AK3(44),
C * AK4(44), AK5(44), AK6(44),F0V(30),BK1(30),
C * BK2(30), BK3(30), BK4(30), BK5(30), BK6(30),
C * A0,A1,A2,A3,A4,A5,A6
C REAL F,TTT2
REAL H(MH),RO(MH),T(MH),P(MH),G2(MH)
C
C F0V, BK1-6 - ready coefficients for calculation
C Input parameters - MH - the number of the atmosphere arrays elements
C T(MH)-Temperature, K; P(MH)-Pressure, mb (or hPa,or 10**(-2)Pa)
C RO(MH) - absolute humidity, g/(m**3)
C F - current frequency, GHz
C Output parametrer - G2(MH)- absorbtion coefficient in water vapour, Neper
C TTT2 - optical depth
C
C Из ppmv в абс. влажность г/м3
C RO(1)=qqq/1.E4/8.31441/ttt*18.*ppp
C Из г/г в абс. влажность г/м3
C RO(1)=qqq*1.E6*28./18./1.E4/8.31441/ttt*18.*ppp
open(61, file='Tb_ini_Liebe_coeff.dat')
read (61,'(7f10.6)') (F0(j),j=1,44)
read (61,'(9f8.2)') (AK1(j),j=1,44)
read (61,'(9f8.3)') (AK2(j),j=1,44)
read (61,'(9f8.2)') (AK3(j),j=1,44)
read (61,'(9f8.1)') (AK4(j),j=1,44)
read (61,'(9f8.3)') (AK5(j),j=1,44)
read (61,'(9f8.3)') (AK6(j),j=1,44)
read (61,'(7f10.6)') (F0V(j),j=1,30)
read (61,'(9f9.4)') (BK1(j),j=1,30)
read (61,'(9f8.3)') (BK2(j),j=1,30)
read (61,'(9f8.2)') (BK3(j),j=1,30)
read (61,'(9f8.2)') (BK4(j),j=1,30)
read (61,'(9f8.2)') (BK5(j),j=1,30)
read (61,'(9f8.2)') (BK6(j),j=1,30)
close (61)
RO(1)=qqq
T(1)=TTT
p(1)=ppp
CL = 1.0 ! Strength
CW = 1.0 ! Width
CC = 1.0 ! Continuum
BK1(1) = BK1(1)*CL
BK3(1) = BK3(1)*CW
C
* Case (2) ! Crus Pol for 20 - 24 GHz + Liebe
C
IF (F.LE.32.AND.F.GE.20) THEN ! Crus Pol
CL = 1.0639 ! Strength
CW = 1.0658
CC = 1.2
C
DO 1120 J = 1,MH
C
TETA = 300./T(J)
E = RO(J)/(7.223*TETA)
PE = 0.1*P(J)- E
C
C
G = 0.
C Resonant coefficient:
C
SI = 0.109*E*(TETA**3.5)*EXP(2.143*(1.- TETA))*CL
GAI = 27.84*0.001*(PE*(TETA**(0.6))+4.8*E*(TETA**1.1))
GAI = GAI*CW
AI = GAI*F/22.235
XXI = (22.235-F)**2+GAI**2
YYI = (22.235+F)**2+GAI**2
FI = AI/XXI+AI/YYI
G = SI*FI
C
C NonResonant coefficient:
GNR = F*(3.57*(TETA**7.5)*E+0.113*PE)*0.00001*E*(TETA**3)
GNR = GNR*CC
C
G = G+GNR
G = 0.1820*F*G
C Converting into Neper:
G2(J) = G/4.34
pkpar=G2(J)
C
1120 CONTINUE
C
C
ELSE
! Liebe
C
DO 1220 J = 1,MH
C
TETA = 300./T(J)
E = RO(J)/(7.223*TETA)
PE = 0.1*P(J)- E
C
G = 0.
C Resonant coefficient:
DO 1221 I = 1,30
SI = BK1(I)*E*(TETA**3.5)*EXP(BK2(I)*(1.- TETA))*CL
GAI = BK3(I)*0.001*(PE*(TETA**(BK4(I)))+BK5(I)*E*(TETA**BK6(I)))
GAI = GAI*CW
AI = GAI*F/F0V(I)
XXI = (F0V(I)-F)**2+GAI**2
YYI = (F0V(I)+F)**2+GAI**2
FI = AI/XXI+AI/YYI
1221 G = G+SI*FI
C
C NonResonant coefficient:
GNR = F*(3.57*(TETA**7.5)*E+0.113*PE)*0.00001*E*(TETA**3)
GNR = GNR*CC
C
G = G+GNR
G = 0.1820*F*G
C Converting into Neper:
G2(J) = G/4.34
pkpar=G2(J)
C
1220 CONTINUE
C
END IF
C
C
C TTT2=0
C
C DO J=1,(MH-1)
C TTT2=TTT2+G2(J+1)*(H(J+1)-H(J))*SAZ
C IF(TTT2.GE.20.) GO TO 62
C END DO
C
C
C 62 CONTINUE
C
C
RETURN
END
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment