Program IRTest
include 'iruph.fh'
integer Nx, Ny
parameter (Nx=1, Ny=1)
integer k, l, i, j
double precision Fintup_f(Nx, Ny, 0:Nlev),
& Fintdown_f(Nx, Ny, 0:Nlev)
C Heat_f -- радиационный нагрев (K/сут)
double precision Heat_f(Nx, Ny, Nlev)
double precision Fup(NBand, 0:Nlev), Fdown(NBand, 0:Nlev)
integer in, out
CCC *********************
call IR_init(in,out)
do j = 1, Ny
do i = 1, Nx
write(*,*) i, j
CC if ((j .ge. 2) .and. (i.ge.2)) debug = .TRUE.
call IR_read(in,i,j)
call IR_vert
call irflux(Fup, Fdown)
do k = 0, Nlev
Fintup_f(i,j,k) = 0.0
Fintdown_f(i,j,k) = 0.0
do l = 1, NBand
Fintup_f(i,j,k) = Fintup_f(i,j,k) +
& Fup(l,k)*DeltaNu(l)
Fintdown_f(i,j,k) = Fintdown_f(i,j,k) +
& Fdown(l,k)*DeltaNu(l)
end do
end do
do k = 1, Nlev
Heat_f(i,j,k) = (Fintup_f(i,j,k-1) - Fintup_f(i,j,k) -
& Fintdown_f(i,j,k-1) + Fintdown_f(i,j,k)) * g /
& (cp * (P(k-1)-P(k)) * mb2cgs) * day2sec
end do
do k = 0, Nlev
write(out,"(6(F10.3))") H(k), P(k), P(k) / P(0),
& Fintup_f(i,j,k) * 1.D-3,
& Fintdown_f(i,j,k) * 1.D-3,
& (Fintup_f(i,j,k) - Fintdown_f(i,j,k)) * 1.D-3
if (k .lt. Nlev) then
write(out,"('#',9X,2F10.3,35X,F10.3)") 0.5D0*(P(k+1) +
& P(k)), 0.5D0*(P(k+1) + P(k))/P(0), Heat_f(i,j,k+1)
end if
end do
end do
end do
end
block data IR_globals
include 'iruph.fh'
integer k
data NuC /70D0, 210D0, 380D0, 530D0, 610D0, 670D0, 730D0,
& 870D0, 1020D0, 1120D0, 1210D0, 1270D0, 1330D0, 1390D0,
& 1550D0, 1940D0, 2230D0/
data DeltaNu /140D0, 140D0, 200D0, 100D0, 60D0, 60D0, 60D0,
& 220D0, 80D0, 120D0, 60D0, 60D0, 60D0, 60D0,
& 260D0, 520D0, 60D0/
C
data debug /.FALSE./, do_H2O /.FALSE./, do_CO2 /.FALSE./,
& do_O3 /.FALSE./, do_cont /.FALSE./, do_cld /.FALSE./
end
subroutine IR_init(in, out)
include 'iruph.fh'
character*12 modname, outname, key
integer in, out
logical EOF
in = 10
out = 11
open(8, file='ir_pars.dat')
10 read(8,*,err=10) modname
15 read(8,*,err=15) outname
do while (.NOT. EOF)
20 read(8,*,END=30,ERR=20) key
if(key .EQ. 'debug') then
debug = .TRUE.
else if (key .EQ. 'h2o') then
do_H2O = .TRUE.
else if (key .EQ. 'co2') then
do_CO2 = .TRUE.
else if (key .EQ. 'o3') then
do_O3 = .TRUE.
else if (key .EQ. 'cont') then
do_cont = .TRUE.
else if (key .EQ. 'cld') then
do_cld = .TRUE.
end if
goto 40
30 EOF=.TRUE.
40 continue
end do
close(8)
open(in, file=modname)
open(out, file=outname)
write(out,"('#',A)") " H2O CO2 O3 cnt cld"
write(out,"('#',5(L3,X))") do_H2O, do_CO2, do_O3, do_cont, do_cld
write(out,"('#')")
write(out,"('#',2A)") " H(км) P(мб) P/Ps F_u(Вт/м2) ",
& "F_d(Вт/м2) F_e(Вт/м2) Нагрев(К/сут)"
end
subroutine IR_read(in,i,j)
include 'iruph.fh'
integer in,i,j
integer k
double precision H_rec, P_rec, T_rec, RhoAir_rec,
& RhoVap_rec, RhoO3_rec
character*12 modname, key
do k = 0, Nlev
50 read(in,*,err=50,end=60) H_rec, P_rec, T_rec, RhoAir_rec,
& RhoVap_rec, RhoO3_rec
T(k) = T_rec
P(k) = P_rec
H(k) = H_rec
cc q_H2O(k) = RhoVap_rec / RhoAir_rec
q_O3(k) = RhoO3_rec / RhoAir_rec
q_H2O(k) = RhoVap_rec / RhoAir_rec ! Increase H2O content on 20%
c ex q_O3(k) = 0.7 * RhoO3_rec / RhoAir_rec ! Decrease O3 content on 30%
C [CO2] = 300 ppmv
q_CO2(k) = 0.456D-3
end do
60 if (debug) write(*,*) '***', k, 'Records read'
end
subroutine IR_vert
C расчет поглощающих масс и связаных величин для данной вертикали (i,j).
include 'iruph.fh'
integer k
double precision w_H2O(0:Nz_cld), Pw_H2O(0:Nz_cld),
& Tw_H2O(0:Nz_cld),
& w_CO2(0:Nlev), Pw_CO2(0:Nlev), Tw_CO2(0:Nlev),
& w_O3(0:Nlev), Pw_O3(0:Nlev), Tw_O3(0:Nlev),
& tau_cont(0:Nz_cld)
C tau_cont -- оптич. длина пути для континуального поглощения водяного пара
common /ir_H2O/ w_H2O, Pw_H2O, Tw_H2O
common /ir_CO2/ w_CO2, Pw_CO2, Tw_CO2
common /ir_O3/ w_O3, Pw_O3, Tw_O3
common /ir_cont/ tau_cont
save /ir_H2O/, /ir_CO2/, /ir_O3/, /ir_cont/
C Ниже -- временные переменные для вычисления пропускания облаков
C Параметры для расчета температурной зависимости конинуального поглощения
C по Roberts, Selby, Biberman, 1977
C gamma_cont выбрана по рекомендации Б.А. Фомина.
double precision T1_cont, T0_cont, gamma_cont
parameter (T1_cont = 296.D0, T0_cont = 1800.D0,
& gamma_cont = 0.002)
C
C
do k = 1, Nlev
DeltaL(k) = (P(k-1) - P(k))/g * mb2CGS * betadiff
end do
C Заполняем массивы из специфических для каждого в-ва common-блоков
w_H2O(0) = 0.D0
Pw_H2O(0) = 0.D0
Tw_H2O(0) = 0.D0
do k=1, Nz_cld
w_H2O(k) = w_H2O(k-1) + 0.5D0*(q_H2O(k-1) + q_H2O(k)) *
& DeltaL(k)
Pw_H2O(k) = Pw_H2O(k-1) + 0.5D0*(q_H2O(k-1) + q_H2O(k)) *
& DeltaL(k) * 0.5D0*(P(k-1) + P(k))
Tw_H2O(k) = Tw_H2O(k-1) + 0.5D0*(q_H2O(k-1) + q_H2O(k)) *
& DeltaL(k) * 0.5D0*(T(k-1) + T(k))
end do
if (debug) write (*,*) 'H2O'
C
w_CO2(0) = 0.D0
Pw_CO2(0) = 0.D0
Tw_CO2(0) = 0.D0
do k=1, Nlev
w_CO2(k) = w_CO2(k-1) + 0.5D0*(q_CO2(k-1) + q_CO2(k)) *
& DeltaL(k)
Pw_CO2(k) = Pw_CO2(k-1) + 0.5D0*(q_CO2(k-1) + q_CO2(k)) *
& DeltaL(k) * 0.5D0*(P(k-1) + P(k))
Tw_CO2(k) = Tw_CO2(k-1) + 0.5D0*(q_CO2(k-1) + q_CO2(k)) *
& DeltaL(k) * 0.5D0*(T(k-1) + T(k))
write(*,*) k, q_CO2(k), w_CO2(k), Pw_CO2(k), Tw_CO2(k)
end do
if (debug) write (*,*) 'CO2'
C
w_O3(0) = 0.D0
Pw_O3(0) = 0.D0
Tw_O3(0) = 0.D0
do k=1, Nlev
w_O3(k) = w_O3(k-1) + 0.5D0*(q_O3(k-1) + q_O3(k)) *
& DeltaL(k)
Pw_O3(k) = Pw_O3(k-1) + 0.5D0*(q_O3(k-1) + q_O3(k)) *
& DeltaL(k) * 0.5D0*(P(k-1) + P(k))
Tw_O3(k) = Tw_O3(k-1) + 0.5D0*(q_O3(k-1) + q_O3(k)) *
& DeltaL(k) * 0.5D0*(T(k-1) + T(k))
end do
if (debug) write (*,*) 'O3'
C
tau_cont(0) = 0.D0
do k = 1, Nz_cld
tau_cont(k) = tau_cont(k-1) + 0.5D0*(q_H2O(k-1) + q_H2O(k)) *
& DeltaL(k) * 0.5D0*(P(k-1) + P(k)) *
& ( gamma_cont + (1.D0 - gamma_cont) *
& 0.5D0*(q_H2O(k-1) + q_H2O(k))/eps_H2O ) *
& exp(T0_cont * ( 2.D0/(T(k-1) + T(k)) - 1.D0/T1_cont ))
end do
if (debug) write (*,*) 'IR_vert done'
end
subroutine IR_vcld
include 'iruph.fh'
integer k, l, ifract
double precision R1(NBand), kappa(NBand)
C R1 -- доля светового потока в рассеянном луче первого порядка
C R1 = (R1_s + R1_p)/2, см. Шифрин, стр. 117-118
C kappa -- мнимая часть показателя преломления
double precision rho_drop
C rho_drop -- плотность капли [г/см3]
parameter(rho_drop = 1.D0)
double precision tau_cld(NBand, 0:Nz_cld)
C c_cld балл облачности, отнесенный к единце (0<=c_cld<=1)
C r_cld(i)[см] -- радиус капель i-й фракции
common /ir_cld/ tau_cld
save /ir_cld/
data R1 /0.16779, 0.09924, 0.09537, 0.09193, 0.08031, 0.06807,
& 0.06098, 0.03934, 0.04757, 0.05615, 0.05871, 0.06047,
& 0.06237, 0.06425, 0.06416, 0.06271, 0.06799/
data kappa /0.37587, 0.47742, 0.45261, 0.43073, 0.42651,
& 0.42184, 0.34696, 0.17827, 0.05494, 0.04997, 0.04860,
& 0.04767, 0.04670, 0.04573, 0.07361, 0.01500, 0.01625/
do l = 1, NBand
tau_cld(l, 0) = 0.0D0
do k = 1, Nz_cld
tau_cld(l, k) = tau_cld(l, k-1)
do ifract = 1, Nfract
tau_cld(l, k) = tau_cld(l, k) +
& 0.5*(q_cld(ifract, k-1) + q_cld(ifract, k)) * DeltaL(k) *
& 0.75*(1.-R1(l)) / (r_cld(ifract,k) * rho_drop) *
& (1. - exp(-8.*pi*kappa(l)*NuC(l)*r_cld(ifract,k))) *
& c_cld(k)
end do
end do
end do
end
double precision function Fi(T)
C compute the ice particle formation probability according to
C H. Sundkvist, Beitr. phys. Atmosph., May 1993, Vol. 66 #1-2.
double precision T, T1, Tci, Tn, A
parameter(T1 = 273., Tci = 232., Tn = 24.04, A = 1.0576936)
if (T .LE. Tci) then
Fi = 1.
return
else if (T .GE. T1) then
Fi = 0.
return
else
Fi = 1. - A * ( 1. - exp( -((T - Tci)/Tn)**2 ) )
return
end if
end
double precision function FPlank(Nufp, Tfp)
C FPlank - значение функции Планка (спектральной плотности потока черного излучения)
C при температуре Tfp для волнового числа Nufp
implicit none
double precision a,b,pi
C константы, входящие в формулу для спектральной плотности
C потока черного излучения B(T). a[эрг*см^2/с], b[см*K]
parameter (a=1.19107D-5, b=1.43875D0, pi=3.14159D0)
double precision Tfp, Nufp
FPlank = pi*(a * Nufp * Nufp * Nufp) /
& (exp(b * Nufp / Tfp) - 1.0D0)
end
subroutine irflux(Fup, Fdown)
include 'iruph.fh'
double precision Fup(NBand,0:Nlev), Fdown(NBand,0:Nlev)
C Fup(j,k), Fup(j,k) (эрг/(см^3*c)) - спектральные плотности восходящего
C и нисходящего соответственно потоков излучения в
C спектральном интервале j на k-м уровне
C "Рабочие" переменные
double precision TrF1(NBand), TrF2(NBand)
C Эти массивы передаются в transm для вычисления функции пропускания
C в каждом из спектральных интервалов
double precision VInt1(NBand), VInt2(NBand), VIntR(NBand)
integer band, k
C
C
do k = 0, Nlev-1
call transm(TrF1, k, Nlev, NOREFL)
call transm(TrF2, k, Nlev, REFL)
call integ_tr(VInt1, VInt2, VIntR, k)
do band = 1, NBand
Fup(band,k) = FPlank(NuC(band), T(k)) -
& (1 - delta) * FPlank(NuC(band), T(Nlev)) * TrF2(band) -
& VInt1(band) + (1 - delta) * VIntR(band)
Fdown(band,k) = FPlank(NuC(band), T(k)) -
& FPlank(NuC(band), T(Nlev)) * TrF1(band) +
& VInt2(band)
end do
end do
C Вычисляются потоки на верхней границе.
C Нисходящий равен 0 -- таковы граничные условия.
call transm(TrF2, Nlev, Nlev, REFL)
call integ_tr(VInt1, VInt2, VIntR, Nlev)
do band = 1, NBand
Fup(band,Nlev) = FPlank(NuC(band), T(Nlev)) * ( 1 -
& (1 - delta) * TrF2(band) ) -
& VInt1(band) + (1 - delta) * VIntR(band)
Fdown(band,Nlev) = 0.D0
end do
end
subroutine transm(TrF, k1, k2, PathKind)
include 'iruph.fh'
integer k1, k2
double precision TrF(NBand)
C
C - N
C |
C - N-1
C |
C - N-2
C |
C - k2 -
C | | Слой, для которого вычисляется функция пропускания
C . |
C . | при PathKind = NOREFL.
C . |
C | |
C - k1 - Порядок k1 и k2 не важен (можно k2>k1).
C |
C - 2
C |
C - 1
C |
C - 0
C
C
C
C - N
C |
C - N-1
C |
C - N-2
C |
C - k2 _
C | | Слой, для которого вычисляется функция пропускания
C . |
C . | при PathKind = REFL.
C . |
C | |
C - k1 | _
C | | |
C - 2 | | Порядок k1 и k2 не важен (можно k2>k1)
C | | |
C - 1 | |
C | | |
C - 0 |_|
C
C
integer kmin, kmax, PathKind
integer band
double precision TrF_H2O(NBand), TrF_cont(NBand), TrF_CO2(NBand),
& TrF_O3(NBand), TrF_cld
data TrF_H2O /NBand * 1.0/
data TrF_cont /NBand * 1.0/
data TrF_CO2 /NBand * 1.0/
data TrF_O3 /NBand * 1.0/
data TrF_cld /1.0/
C
C
CC write(*,*) 'Transm', k1, k2, PathKind
kmin = min(k1,k2)
kmax = max(k1,k2)
if (do_CO2) call TrCO2(TrF_CO2, kmin, kmax, PathKind)
if (do_O3) call TrO3(TrF_O3, kmin, kmax, PathKind)
CC write(*,'(3(i3x)f10.8)') kmin, kmax, PathKind, TrF_CO2(6)
kmin = min(kmin, Nz_cld)
kmax = min(kmax, Nz_cld)
if ( (kmin .eq. kmax) .and. (PathKind .eq. NOREFL) ) goto 10
if (do_H2O) call TrH2O(TrF_H2O, kmin, kmax, PathKind)
if (do_cont) call Trcont(TrF_cont, kmin, kmax, PathKind)
if (do_cld) call Tr_cld(TrF_cld, kmin, kmax, PathKind)
C Конец
10 continue
do band = 1, NBand
TrF(band) = TrF_H2O(band) * TrF_cont(band) *
& TrF_CO2(band) * TrF_O3(band) * TrF_cld
CC if (TrF(band) .ne. 1.) write(*,*)
CC & kmin, kmax, band, PathKind, TrF(band)
end do
CC write(*,'(3(i3x)17f10.8)') kmin, kmax, PathKind,
CC & (TrF(band), band=1, NBand)
end
subroutine Tr_cld(TrF, k1, k2, PathKind)
include 'iruph.fh'
double precision TrF(NBand)
integer k1, k2
integer PathKind
C PathKind = REFL => считается функция пропускания для слоя от k1 до пов-ти +
C от пов-ти до k2; PathKind = NOREFL => рассматривается слой от k1 до k2.
C REFL = 1, NOREFL = -1 -- константы
double precision tau_cld(NBand,0:Nz_cld)
integer band
common /ir_cld/ tau_cld
save /ir_cld/
C
C
do band = 1, NBand
TrF(band) = exp(-(tau_cld(band,k2)+tau_cld(band,k1)*PathKind))
end do
end
subroutine TrH2O(TrF, k1, k2, PathKind)
include 'iruph.fh'
double precision TrF(NBand)
integer k1, k2
integer PathKind
C PathKind = REFL => считается функция пропускания для слоя от k1 до пов-ти +
C от пов-ти до k2; PathKind = NOREFL => рассматривается слой от k1 до k2.
C REFL = 1, NOREFL = -1 -- константы
integer NB_H2O
parameter (NB_H2O = 16)
integer band(NB_H2O)
double precision alpha(0:NT, NB_H2O), beta(0:NT, NB_H2O)
C alpha [см^2*мб/г], beta [см^2/г] Здесь и далее коэффициенты
C alpha и beta взяты из статьи Розанов, Тимофеев, Фролькис 1993, ФАО,т.29#4
double precision w(0:Nz_cld), Pw(0:Nz_cld), Tw(0:Nz_cld)
common /ir_H2O/ w, Pw, Tw
save /ir_H2O/
integer bnd, i_t
C i_t -- номер ближайшего снизу узлового значения температуры, a_t1 и a_t2 -- коэфф. интерполяции
double precision wt, Pwt, Twt, a_t1, a_t2
C Номера спектральных интервалов, в которых поглощает водяной пар.
data band /1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/
CC data alpha /
CCC 200 225 250 275 300
CC & 4.3282D+7, 3.7193D+7, 3.2434D+7, 2.8609D+7, 2.5470D+7, !1
CC & 8.1582D+7, 7.6555D+7, 7.1793D+7, 6.7339D+7, 6.3217D+7, !2
CC & 7.4699D+6, 9.6882D+6, 1.1884D+7, 1.3959D+7, 1.5861D+7, !3
CC & 1.2118D+5, 1.5645D+5, 1.9211D+5, 2.2717D+5, 2.6345D+5, !4
CC & 4.6103D+4, 7.5386D+4, 1.1359D+5, 1.6009D+5, 2.1319D+5, !5
CC & 1.2764D+4, 2.0265D+4, 3.0803D+4, 4.4684D+4, 6.2011D+4, !6
CC & 4.4548D+3, 7.3967D+3, 1.1956D+4, 1.8542D+4, 2.7763D+4, !7
CC & 9.1193D+2, 1.4186D+3, 2.2292D+3, 3.2424D+3, 4.7623D+3, !8
CC & 4.0530D+1, 8.1060D+1, 1.4186D+2, 2.3305D+2, 3.7490D+2, !9
CC & 8.1060D+2, 1.4186D+3, 2.2292D+3, 3.1411D+3, 4.1543D+3, !10
CC & 1.8239D+3, 2.9384D+3, 4.3570D+3, 5.9782D+3, 7.8020D+3, !11
CC & 8.5113D+3, 1.2260D+4, 1.7023D+4, 2.3305D+4, 3.1309D+4, !12
CC & 4.4279D+4, 7.6399D+4, 1.1855D+5, 1.6901D+5, 2.2514D+5, !13
CC & 1.3390D+6, 1.5987D+6, 1.8051D+6, 1.9613D+6, 2.0733D+6, !14
CC & 1.0830D+7, 8.9660D+6, 7.5624D+6, 6.4704D+6, 5.5999D+6, !15
CC & 6.6102D+6, 5.6010D+6, 4.8141D+6, 4.1847D+6, 3.6717D+6/ !16
C
data alpha /
C 200 225 250 275 300
& 4.3282D+7, 3.7193D+7, 3.2434D+7, 2.8609D+7, 2.5470D+7, !1
& 8.1582D+7, 7.6555D+7, 7.1793D+7, 6.7339D+7, 6.3217D+7, !2
& 7.4699D+6, 9.6882D+6, 1.1884D+7, 1.3959D+7, 1.5861D+7, !3
& 1.2118D+5, 1.5645D+5, 1.9211D+5, 2.2717D+5, 2.6345D+5, !4
& 4.6103D+4, 7.5386D+4, 1.1359D+5, 1.6009D+5, 2.1319D+5, !5
& 1.2764D+4, 2.0265D+4, 3.0803D+4, 4.4684D+4, 6.2011D+4, !6
& 4.4548D+3, 7.3967D+3, 1.1956D+4, 1.8542D+4, 2.7763D+4, !7
& 9.1193D+2, 1.4186D+3, 2.2292D+3, 3.2424D+3, 4.7623D+3, !8
& 4.0530D+1, 8.1060D+1, 1.4186D+2, 2.3305D+2, 3.7490D+2, !9
& 8.1060D+2, 1.4186D+3, 2.2292D+3, 3.1411D+3, 4.1543D+3, !10
& 1.8239D+3, 2.9384D+3, 4.3570D+3, 5.9782D+3, 7.8020D+3, !11
& 8.5113D+3, 1.2260D+4, 1.7023D+4, 2.3305D+4, 3.1309D+4, !12
& 4.4279D+4, 7.6399D+4, 1.1855D+5, 1.6901D+5, 2.2514D+5, !13
& 1.3390D+6, 1.5987D+6, 1.8051D+6, 1.9613D+6, 2.0733D+6, !14
& 1.0830D+7, 8.9660D+6, 7.5624D+6, 6.4704D+6, 5.5999D+6, !15
& 6.6102D+6, 5.6010D+6, 4.8141D+6, 4.1847D+6, 3.6717D+6/ !16
C
data beta /
C 200 225 250 275 300
& 4.3095D+3, 4.2589D+3, 4.2041D+3, 4.1438D+3, 4.0803D+3, !1
& 5.1439D+3, 5.6603D+3, 6.1131D+3, 6.5106D+3, 6.8618D+3, !2
& 2.9110D+2, 4.3200D+2, 5.9900D+2, 7.8770D+2, 9.9390D+2, !3
& 4.5000D+0, 7.1000D+0, 1.0400D+1, 1.4600D+1, 1.9900D+1, !4
& 2.0000D+0, 3.8000D+0, 6.4000D+0, 1.0000D+1, 1.4600D+1, !5
& 5.0000D-1, 9.0000D-1, 1.5000D+0, 2.5000D+0, 3.7000D+0, !6
& 1.6000D-1, 3.0000D-1, 6.0000D-1, 1.0000D+0, 1.6000D+0, !7
& 2.0000D-2, 5.0000D-2, 9.0000D-2, 1.5000D-1, 2.5000D-1, !8
& 2.0000D-3, 5.0000D-3, 1.1000D-2, 2.1000D-2, 4.0000D-2, !9
& 4.0000D-2, 7.0000D-2, 1.3000D-1, 2.0000D-1, 3.0000D-1, !10
& 1.4000D-1, 3.0000D-1, 4.0000D-1, 7.0000D-1, 1.0000D+0, !11
& 1.0000D+0, 1.6000D+0, 2.6000D+0, 4.0000D+0, 6.0000D+0, !12
& 6.8000D+0, 1.2000D+1, 1.9300D+1, 2.8800D+1, 4.0500D+1, !13
& 1.2340D+2, 1.6720D+2, 2.1240D+2, 2.5760D+2, 3.0160D+2, !14
& 1.6813D+3, 1.6372D+3, 1.5974D+3, 1.5612D+3, 1.5280D+3, !15
& 2.2340D+2, 2.3100D+2, 2.3650D+2, 2.4100D+2, 2.4480D+2/ !16
C
C
wt = w(k2) + PathKind*w(k1)
Pwt = ( Pw(k2) + PathKind*Pw(k1) ) / wt
Twt = ( Tw(k2) + PathKind*Tw(k1) ) / wt
Twt = max(TLow, min(Twt, THigh * 0.99999D0))
C Температура не должна достигать правой границы отрезка, иначе
C нарушится единообразие формулы линейной интерполяции
i_t = int( (Twt-TLow) / dTtab )
a_t2 = ( Twt - (TLow + i_t * dTtab) ) / dTtab
a_t1 = 1.D0 - a_t2
if (debug) write (*,*) k1, k2, wt, pwt, Twt
do bnd = 1, NB_H2O
TrF(band(bnd)) = exp(
& -( a_t1*beta(i_t, bnd) + a_t2*beta(i_t+1, bnd) ) * wt /
& (
& 1 + ( a_t1*alpha(i_t, bnd) + a_t2*alpha(i_t+1, bnd) ) *
& wt / Pwt
& ) ** 0.5 )
end do
CCC write (91,*) 'H2O', PathKind, ' [ ', k1, k2, ' ]: ', wt, pwt, twt
CCC write(*,*) (TrF(bnd), bnd=1, NBand)
end
subroutine Trcont(TrF, k1, k2, PathKind)
include 'iruph.fh'
double precision TrF(NBand)
integer k1, k2
integer PathKind
C PathKind = REFL => считается функция пропускания для слоя от k1 до пов-ти +
C от пов-ти до k2; PathKind = NOREFL => рассматривается слой от k1 до k2.
C REFL = 1, NOREFL = -1 -- константы
integer NB_cont
CCC parameter(NB_cont = 4)
parameter(NB_cont = 9)
integer band(NB_cont)
double precision C0(NB_cont)
C C0 [см^2/мб*г]
double precision tau(0:Nz_cld)
common /ir_cont/ tau
save /ir_cont/
integer bnd
C Номера спектральных интервалов, в которых происходит континуальное поглощение
CCC data band /8,9,10,11/
data band /3,4,5,6,7,8,9,10,11/
C RSB, a=4.13e-3 см2/(мб*г); b=7.72 см2/(мб*г); beta = 8.30e-3 см.
CCC data C0 /.0135080D0, .0059201D0, .0049407D0, .0045254D0/
data C0 /0.334D0, 0.099D0, 0.053D0, 0.034D0, 0.220D-1, 0.977D-2,
& 0.76D-2, 0.484D-2, 0.447D-2/
C
C
do bnd = 1, NB_cont
TrF(band(bnd)) = exp( -C0(bnd) * (tau(k2) + PathKind*tau(k1)) )
end do
end
subroutine TrCO2(TrF, k1, k2, PathKind)
include 'iruph.fh'
double precision TrF(NBand)
integer k1, k2
integer PathKind
C PathKind = REFL => считается функция пропускания для слоя от k1 до пов-ти +
C от пов-ти до k2; PathKind = NOREFL => рассматривается слой от k1 до k2.
C REFL = 1, NOREFL = -1 -- константы
integer NB_CO2
parameter(NB_CO2 = 9)
integer band(NB_CO2)
double precision w(0:Nlev), Pw(0:Nlev), Tw(0:Nlev)
common /ir_CO2/ w, Pw, Tw
save /ir_CO2/
integer bnd, i_t
double precision wt, Pwt, Twt, a_t1, a_t2
double precision alpha(0:NT, NB_CO2), beta(0:NT, NB_CO2)
C alpha [см^2*мб/г], beta [см^2/г]
C
C Номера спектральных интервалов, в которых поглощает углекислый газ
data band /4,5,6,7,8,9,10,16,17/
data alpha /
C 200 225 250 275 300
& 2.0265D+1, 6.0795D+1, 1.3172D+2, 2.4318D+2, 3.9517D+2, !4
& 5.6134D+4, 6.7888D+4, 7.7412D+4, 8.5316D+4, 9.1193D+4, !5
& 4.4188D+6, 3.5692D+6, 2.9033D+6, 2.3809D+6, 1.9710D+6, !6
& 3.5261D+4, 4.7116D+4, 5.7654D+4, 6.6267D+4, 7.2853D+4, !7
& 3.0398D+1, 8.1060D+1, 1.6212D+2, 2.9384D+2, 4.7623D+2, !8
& 4.0530D+1, 1.0133D+2, 1.9252D+2, 3.1411D+2, 4.8636D+2, !9
& 4.0530D+1, 9.1193D+1, 1.7225D+2, 2.9384D+2, 4.3570D+2, !10
& 8.8153D+2, 7.0928D+2, 5.8769D+2, 5.0663D+2, 4.3570D+2, !16
& 9.6259D+4, 9.3827D+4, 8.8558D+4, 8.1972D+4, 7.4677D+4/ !17
C
data beta /
C 200 225 250 275 300
& 7.0000D-3, 2.0000D-2, 5.0000D-2, 1.2000D-1, 2.3000D-1, !4
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 1.9700D+1, 2.9400D+1, 4.1500D+1, 5.5300D+1, 7.1100D+1, !5
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 1.9326D+3, 1.9461D+3, 1.9659D+3, 1.9925D+3, 2.0263D+3, !6
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 9.2000D+0, 1.5400D+1, 2.3500D+1, 3.3300D+1, 4.4900D+1, !7
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 6.0000D-3, 2.0000D-2, 4.0000D-2, 1.0000D-1, 1.8000D-1, !8
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 4.0000D-3, 1.1000D-2, 2.6000D-2, 5.4000D-2, 1.0000D-1, !9
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 3.0000D-3, 9.0000D-3, 2.0000D-2, 4.0000D-2, 7.0000D-2, !10
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 9.0000D-2, 1.0000D-1, 1.1000D-1, 1.2000D-1, 1.3000D-1, !16
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,
& 1.6000D+1, 2.0800D+1, 2.6000D+1, 3.1600D+1, 3.7400D+1/ !17
CC & 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0/
C
wt = w(k2) + PathKind*w(k1)
Pwt = ( Pw(k2) + PathKind*Pw(k1) ) / wt
Twt = ( Tw(k2) + PathKind*Tw(k1) ) / wt
Twt = max(TLow, min(Twt, THigh * 0.99999D0))
C Температура не должна достигать правой границы отрезка, иначе
C нарушится единообразие формулы линейной интерполяции
i_t = int( (Twt-TLow) / dTtab )
a_t2 = ( Twt - (TLow + i_t * dTtab) ) / dTtab
a_t1 = 1.D0 - a_t2
CCC write (91,*) 'CO2', PathKind, ' [ ', k1, k2, ' ]: ', wt, pwt, twt
do bnd = 1, NB_CO2
TrF(band(bnd)) = exp(
& -( a_t1*beta(i_t, bnd) + a_t2*beta(i_t+1, bnd) ) * wt /
& (
& 1 + ( a_t1*alpha(i_t, bnd) + a_t2*alpha(i_t+1, bnd) ) *
& wt / Pwt
& ) ** 0.5 )
CC if (TrF(band(bnd)) .ne. 1.) then
CC write(*,*) bnd, band(bnd), k1, k2, '*********'
CC write(*,*) Twt, Pwt, wt, i_t
CC write(*,*) a_t1, a_t2, beta(i_t, bnd), beta(i_t+1, bnd)
CC write(*,*) -( a_t1*beta(i_t, bnd) + a_t2*beta(i_t+1, bnd) )
CC & * wt / (
CC & 1 + ( a_t1*alpha(i_t, bnd) + a_t2*alpha(i_t+1, bnd) ) *
CC & wt / Pwt
CC & ) ** 0.5
CC write(*,*)
CC end if
end do
end
subroutine TrO3(TrF, k1, k2, PathKind)
include 'iruph.fh'
double precision TrF(NBand)
integer k1, k2
integer PathKind
C PathKind = REFL => считается функция пропускания для слоя от k1 до пов-ти +
C от пов-ти до k2; PathKind = NOREFL => рассматривается слой от k1 до k2.
C REFL = 1, NOREFL = -1 -- константы
integer NB_O3
parameter(NB_O3 = 9)
integer band(NB_O3)
double precision alpha(0:NT, NB_O3), beta(0:NT, NB_O3)
double precision w(0:Nlev), Pw(0:Nlev), Tw(0:Nlev)
common /ir_O3/ w, Pw, Tw
save /ir_O3/
integer bnd, i_t
double precision wt, Pwt, Twt, a_t1, a_t2
C alpha [см^2*мб/г], beta [см^2/г]
C Номера спектральных интервалов, в которых поглощает озон
data band /1, 5, 6, 7, 8, 9, 10, 11, 16/
data alpha /
C 200 225 250 275 300
& 1.0943D+4, 1.0335D+4, 9.9299D+3, 9.5246D+3, 9.3219D+3, !1
& 2.9384D+3, 3.3437D+3, 3.6477D+3, 3.8504D+3, 4.0530D+3, !5
& 1.6415D+4, 1.4489D+4, 1.2970D+4, 1.1754D+4, 1.0842D+4, !6
& 1.8745D+4, 1.7225D+4, 1.5199D+4, 1.3679D+4, 1.2362D+4, !7
& 3.4451D+3, 3.1411D+3, 3.0398D+3, 3.0398D+3, 3.2424D+3, !8
& 4.1209D+5, 3.3316D+5, 2.7449D+5, 2.3021D+5, 1.9657D+5, !9
& 1.7327D+4, 1.9150D+4, 2.0569D+4, 2.1785D+4, 2.2697D+4, !10
& 2.8371D+3, 3.1411D+3, 3.3437D+3, 3.4451D+3, 3.5464D+3, !11
& 5.7857D+4, 5.0460D+4, 4.4887D+4, 4.0530D+4, 3.7085D+4/ !16
C
data beta /
C 200 225 250 275 300
& 2.9400D+1, 3.2400D+1, 3.5400D+1, 3.8300D+1, 4.1100D+1, !1
& 2.6000D+0, 3.2000D+0, 3.9000D+0, 4.5000D+0, 5.0000D+0, !5
& 6.2800D+1, 6.2100D+1, 6.1300D+1, 6.0400D+1, 5.9400D+1, !6
& 7.2500D+1, 7.0300D+1, 6.8300D+1, 6.6300D+1, 6.4300D+1, !7
& 2.1000D+0, 2.6000D+0, 3.1000D+0, 3.8000D+0, 4.7000D+0, !8
& 2.0027D+3, 1.9733D+3, 1.9477D+3, 1.9266D+3, 1.9101D+3, !9
& 3.1800D+1, 3.7700D+1, 4.3400D+1, 4.8800D+1, 5.3800D+1, !10
& 2.0000D-2, 2.0000D-2, 2.0000D-2, 2.0000D-2, 3.0000D-2, !11
& 3.2900D+1, 3.2800D+1, 3.2600D+1, 3.2400D+1, 3.2100D+1/ !16
C
C
wt = w(k2) + PathKind*w(k1)
Pwt = ( Pw(k2) + PathKind*Pw(k1) ) / wt
Twt = ( Tw(k2) + PathKind*Tw(k1) ) / wt
Twt = max(TLow, min(Twt, THigh * 0.99999D0))
C Температура не должна достигать правой границы отрезка, иначе
C нарушится единообразие формулы линейной интерполяции
i_t = int( (Twt-TLow) / dTtab )
a_t2 = ( Twt - (TLow + i_t * dTtab) ) / dTtab
a_t1 = 1.D0 - a_t2
do bnd = 1, NB_O3
TrF(band(bnd)) = exp(
& -( a_t1*beta(i_t, bnd) + a_t2*beta(i_t+1, bnd) ) * wt /
& (
& 1 + ( a_t1*alpha(i_t, bnd) + a_t2*alpha(i_t+1, bnd) ) *
& wt / Pwt
& ) ** 0.5 )
C -- Так было до учета -- TrF(band(bnd)) = exp( -beta(bnd) * wt /
C -- температ. завис. -- & (1 + alpha(bnd) * wt / Pwt) ** 0.5 )
end do
end
subroutine integ_tr(val1, val2, valR, k)
include 'iruph.fh'
double precision val1(NBand), val2(NBand), valR(NBand)
integer k
C Проинтегрировать TrF(p_k, p')*dB/dT' по T' от T(0) до T(k) --> val1,
C Проинтегрировать TrF(p_k, p')*dB/dT' по T' от T(k) до T(Nlev) --> val2,
C где TrF -- функция пропускания, вычисляемая при PathKind = NOREFL без захода
C к поверхности.
C Проинтегрировать TrF(p_k, p')*dB/dT' по T' от T(0) до T(Nlev) --> valR,
C где TrF -- функция пропускания, вычисляемая при PathKind = REFL с заходом
C к поверхности.
C Интегрируем методом трапеций.
integer bnd, ktmp
double precision TrF(NBand)
C
C Интегрируем с заходом к поверхности
C Инициализация. Начнем вычислять сумму с последнего члена.
C Последний член суммы, слой между (NLev-1)-м и Nlev-м уровнями.
C Функция пропускания вычисляется с заходом к поверхности, а потому ее надо
C Считать даже при k=Nlev
call transm(TrF, Nlev, k, REFL)
do bnd = 1, NBand
valR(bnd) = 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(Nlev)) -
& FPlank(NuC(bnd), T(Nlev-1)) )
end do
C Основной цикл
do ktmp = 1, Nlev-1
call transm(TrF, ktmp, k, REFL)
do bnd = 1, NBand
valR(bnd) = valR(bnd) + 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(ktmp + 1)) -
& FPlank(NuC(bnd), T(ktmp - 1)) )
end do
end do
C Отдельно рассмотрим случай k=0.
C Тогда функция пропускания слоя от 0 до 0 равна 1, и
C кроме того, val1=0, val2=valR.
if ( k .eq. 0) then
do bnd = 1, NBand
valR(bnd) = 0.5D0 * ( FPlank(NuC(bnd), T(1)) -
& FPlank(NuC(bnd), T(0)) ) + valR(bnd)
end do
do bnd = 1, NBand
val2(bnd) = valR(bnd)
end do
do bnd = 1, NBand
val1(bnd) = 0.D0
end do
return
else
call transm(TrF, 0, k, REFL)
do bnd = 1, NBand
valR(bnd) = 0.5D0 * TrF(bnd) * ( FPlank(NuC(bnd), T(1)) -
& FPlank(NuC(bnd), T(0)) ) + valR(bnd)
end do
end if
C
C Теперь интегрируем без захода к поверхности.
C Отдельно надо рассмотреть случай, когда k=Nlev, чтобы не вычилять
C функцию пропускания для слоя нулевой толщины -- а то zero division!
call transm(TrF, 0, k, NOREFL)
C В последнем члене суммы val1 (слой между (k-1)-м и k-м уровнями) не следует
C вычислять функцию пропускания -- она равна 1 (слой от k до k).
C При инициализации вычислим сразу сумму первого и последнего слагаемых.
do bnd = 1, NBand
val1(bnd) = 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(1)) -
& FPlank(NuC(bnd), T(0)) ) +
& 0.5D0 * ( FPlank(NuC(bnd), T(k)) -
& FPlank(NuC(bnd), T(k-1)) )
end do
C Основной цикл для val1. Дойдем до k-1.
do ktmp = 1, k-1
call transm(TrF, ktmp, k, NOREFL)
do bnd = 1, NBand
val1(bnd) = val1(bnd) + 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(ktmp + 1)) -
& FPlank(NuC(bnd), T(ktmp - 1)) )
end do
end do
C val1 вычислен, valR тоже. Остался val2.
C Отдельно рассмотрим случай k=Nlev. Тогда val2=0, и можно возвращаться.
if (k .eq. Nlev) then
do bnd = 1, NBand
val2(bnd) = 0.D0
end do
return
end if
C Инициализируем val2. Сразу учтем сумму первого и последнего слагаемых.
C В слое от k до k+1 функция пропускания равна 1.
call transm(TrF, Nlev, k, NOREFL)
do bnd = 1, NBand
val2(bnd) = 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(Nlev)) -
& FPlank(NuC(bnd), T(Nlev - 1)) )
& + 0.5D0 * ( FPlank(NuC(bnd), T(k+1)) -
& FPlank(NuC(bnd), T(k)) )
end do
C Основной цикл для val2. Дойдем до Nlev-1
do ktmp = k+1, Nlev-1
call transm(TrF, ktmp, k, NOREFL)
do bnd = 1, NBand
val2(bnd) = val2(bnd) + 0.5D0 * TrF(bnd) *
& ( FPlank(NuC(bnd), T(ktmp+1)) -
& FPlank(NuC(bnd), T(ktmp-1)) )
CC if (val2(bnd) .lt. 0.0) then
CC write(*,*) bnd, k, ktmp, TrF(bnd)
CC end if
end do
end do
C При ktmp=Nlev все уже посчитано при инициализации. Так что все сделано.
end
subroutine IR_ITF
include 'iruph.fh'
double precision TrF(NBand), F0, Ftr
integer l
F0 = 0.D0
do l = 1, NBand
F0 = F0 + FPlank(NuC(l), T(0))*DeltaNu(l)
end do
write(*,*) 'F_surface:', F0
do l = 1, NBand
TrF(l) = 1.D0
end do
call TrH2O(TrF, 0, Nlev, NOREFL)
Ftr= 0.D0
write(*,*) 'H2O'
write (*, "('Nu ',17(F6.0,X))") NuC
write(*,"('B(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l), l=1, NBand)
write(*,"('F(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l)*TrF(l), l=1, NBand)
write (*, "('P(Nu)',17(F6.4,X))") TrF
do l = 1, NBand
Ftr = Ftr + FPlank(NuC(l), T(0))*DeltaNu(l)*TrF(l)
end do
write(*,*) 'ITF_H2O:', Ftr/F0
write(*,*)
do l = 1, NBand
TrF(l) = 1.D0
end do
call TrCO2(TrF, 0, Nlev, NOREFL)
Ftr= 0.D0
write(*,*) 'CO2'
write (*, "('Nu ',17(F6.0,X))") NuC
write(*,"('B(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l), l=1, NBand)
write(*,"('F(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l)*TrF(l), l=1, NBand)
write (*, "('P(Nu)',17(F6.4,X))") TrF
do l = 1, NBand
Ftr = Ftr + FPlank(NuC(l), T(0))*DeltaNu(l)*TrF(l)
end do
write(*,*) 'ITF_CO2:', Ftr/F0
write(*,*)
do l = 1, NBand
TrF(l) = 1.D0
end do
call TrO3(TrF, 0, Nlev, NOREFL)
Ftr= 0.D0
write(*,*) 'O3'
write (*, "('Nu ',17(F6.0,X))") NuC
write(*,"('B(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l), l=1, NBand)
write(*,"('F(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l)*TrF(l), l=1, NBand)
write (*, "('P(Nu)',17(F6.4,X))") TrF
do l = 1, NBand
Ftr = Ftr + FPlank(NuC(l), T(0))*DeltaNu(l)*TrF(l)
end do
write(*,*) 'ITF_O3:', Ftr/F0
write(*,*)
do l = 1, NBand
TrF(l) = 1.D0
end do
call Trcont(TrF, 0, Nlev, NOREFL)
Ftr= 0.D0
write (*, "('Nu ',17(F6.0,X))") NuC
write(*,"('B(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l), l=1, NBand)
write(*,"('F(Nu)',17(F6.0,X))") (Fplank(NuC(l),T(0))*
& DeltaNu(l)*TrF(l), l=1, NBand)
write (*, "('P(Nu)',17(F6.4,X))") TrF
do l = 1, NBand
Ftr = Ftr + FPlank(NuC(l), T(0))*DeltaNu(l)*TrF(l)
end do
write(*,*) 'ITF_cont:', Ftr/F0
write(*,*)
end
subroutine Linterpol(NA, LDA, NG, A, G, NG1, A1, G1)
implicit none
C interpolate a field linearly with respect to one of it's dimensions
integer NA, LDA, NG, NG1
double precision A(NA), G(NG), A1((NG1*NA)/NG), G1(NG1)
C Entry parameters:
C A -- the field to be interpolated
C G -- the grid to be interpolated from
C G1 -- the grid to be interpolated to
C G and G1 are supposed to be monotone
C LDA -- the leading dimension of A (and A1 also), e.g., if A is Nx * Ny * Nz,
C and we interpolate along z axis, then LDA = Nx*Ny;
C if we interpolate along y axis, LDA = Nx, etc.
C NA -- the dimension of A, in the case above NA = Nx*Ny*Nz
C NG -- the dimension of the "source" grid, i.e., dimension of A corresponding to the axis
C along which it is interpolated
C Return parameters:
C A1 -- the interpolated field
integer NGm
parameter(NGm = 30)
integer left(NGm), rigt(NGm)
C si*G(left(i)) <= si*G1(i) < si*G(rigt(i))
double precision aleft(NGm), arigt(NGm)
C aleft and arigt -- coefficients of linear interpolation
integer i, i1, i2, si
if ( NG1 .GT. NGm ) stop 'E0'
if (G1(1) .LT. G1(NG1)) then
si = 1
else
si = -1
end if
if ( (si*G1(1) .LT. si*G(1)) .OR.
& (si*G1(NG1) .GT. si*G(NG)) ) stop 'E1'
i1 = 1
do i = 1, NG1 - 1
do while ( si*G(i1+1) .LE. si*G1(i) )
i1 = i1+1
end do
left(i) = i1
rigt(i) = i1+1
arigt(i) = (G1(i) - G(i1))/(G(i1+1) - G(i1))
aleft(i) = 1.0 - arigt(i)
end do
if ( G1(NG1) .EQ. G(NG) ) then
left(NG1) = NG
rigt(NG1) = NG
aleft(NG1) = 1.0
arigt(NG1) = 0.0
else
do while ( si*G(i1+1) .LE. si*G1(i) )
i1 = i1+1
end do
left(NG1) = i1
rigt(NG1) = i1+1
arigt(NG1) = (G1(NG1) - G(i1))/(G(i1+1) - G(i1))
aleft(NG1) = 1.0 - arigt(NG1)
end if
do i = 1, NA/(LDA*NG)
do i1 = 1, NG1
do i2 = 1, LDA
A1(i2 + ( i1-1 + (i-1)*NG1 ) * LDA) =
& aleft(i1) * A(i2 + ( left(i1)-1 + (i-1)*NG ) * LDA) +
& arigt(i1) * A(i2 + ( rigt(i1)-1 + (i-1)*NG ) * LDA)
end do
end do
end do
end
|
