Context Navigation

lw_venus_ve.1mat @ 777

Last change on this file since 777 was 3, checked in by slebonnois, 14 years ago

Creation de repertoires:

chantiers : pour communiquer sur nos projets de modifs
documentation : pour stocker les docs

Ajout de:

libf/phytitan : physique de Titan
libf/chimtitan: chimie de Titan
libf/phyvenus : physique de Venus

File size: 7.8 KB

Line
1	SUBROUTINE LW_venus_ve(
2	S PPB, PT, PTSURF,
3	S PCOOL,
4	S PTOPLW,PSOLLW,PSOLLWDN,
5	S ZFLNET)
6
7	IMPLICIT none
8
9	#include "dimensions.h"
10	#include "dimphy.h"
11	#include "raddim.h"
12	#include "YOMCST.h"
13	C
14	C ------------------------------------------------------------------
15	C
16	C PURPOSE.
17	C --------
18	C
19	c This routine loads the longwave matrix of factors Ksi,
20	c used to build the Net Exchange Rates matrix Psi.
21	c Psi(i,j,nu) = Ksi(i,j,nu) * ( B(i,nu)-B(j,nu) )
22	c
23	c This Ksi matrix has been computed by Vincent Eymet
24	C
25	c The NER matrix is then integrated in frequency, and the output
26	c are calculated.
27	c
28	C AUTHOR.
29	C -------
30	C Sebastien Lebonnois
31	C
32	C MODIFICATIONS.
33	C --------------
34	C ORIGINAL : 27/07/2005
35	C ------------------------------------------------------------------
36	C
37	C* ARGUMENTS:
38	C
39	c inputs
40
41	REAL PPB(KFLEV+1) ! inter-couches PRESSURE (bar)
42	REAL PT(KFLEV) ! Temperature in layer (K)
43	REAL PTSURF ! Surface temperature
44	C
45	c output
46
47	REAL PCOOL(KFLEV) ! LONGWAVE COOLING (K/VENUSDAY) within each layer
48	REAL PTOPLW ! LONGWAVE FLUX AT T.O.A. (net, + vers le haut)
49	REAL PSOLLW ! LONGWAVE FLUX AT SURFACE (net, + vers le haut)
50	REAL PSOLLWDN ! LONGWAVE FLUX AT SURFACE (down, + vers le bas)
51	REAL ZFLNET(KFLEV+1) ! net thermal flux at ppb levels (+ vers le haut)
52
53	C
54	C* LOCAL VARIABLES:
55	C
56	integer nlve,nnuve
57	parameter (nlve=kflev) ! fichiers Vincent: same grid than GCM
58	parameter (nnuve=68) ! fichiers Vincent et Bullock
59	real dureejour
60	parameter (dureejour=10.087e6)
61
62	integer i,j,p,nl0,nnu0,band
63	real lambda(nnuve) ! wavelenght in table (mu->m, middle of interval)
64	real ksive(0:nlve+1,0:nlve+1,nnuve) ! ksi factors
65	real bplck(0:nlve+1,nnuve) ! Planck luminances in table layers
66	real al(nnuve),bl(nnuve) ! for Planck luminances calculations
67	real psive(0:nlve+1,0:nlve+1,nnuve) ! NER in W/m**2 per wavelength band
68	real psi_1(0:nlve+1,0:nlve+1) ! NER in W/m**2 (sum on lambda)
69
70	real ztemp(0:nlve) ! GCM temperature in table layers
71	real zlnet(nlve+1) ! net thermal flux (W/m**2)
72	real dzlnet(0:nlve) ! Radiative budget (W/m**2)
73	character*22 nullchar
74	real lambdamin,lambdamax ! in microns
75	real dlambda ! cm-1
76
77	real y(0:nlve,nnuve) ! intermediaire Planck
78	real pdp(nlve) ! epaisseur de la couche en pression (Pa)
79	real zdblay(nlve,nnuve) ! gradient en temperature de planck
80
81	logical firstcall
82	data firstcall/.true./
83
84	save lambda,ksive,al,bl,firstcall
85
86	c ------------------------
87	c Loading the files
88	c ------------------------
89
90	if (firstcall) then
91
92	c Matrice Ksi
93	c------------
94	open(13,file='ksi_gccr.txt')
95	read(13,*) nl0,nnu0
96	if (nl0.ne.nlve) then
97	print*,'Probleme de dimension entre ksi.txt et lw'
98	print*,'N levels = ',nl0,nlve
99	stop
100	endif
101	if (nnu0.ne.nnuve) then
102	print*,'Probleme de dimension entre ksi.txt et lw'
103	print*,'N freq = ',nnu0,nnuve
104	stop
105	endif
106	do band=1,nnuve
107	read(13,*) lambdamin,lambdamax ! en microns
108	lambda(band)=(lambdamin+lambdamax)/2.*1.e-6 ! en m
109	dlambda =(1./lambdamin-1./lambdamax)*1.e4 ! en cm-1
110	c print*,band,lambdamin,dlambda,lambdamax
111	do i=0,nlve+1
112	read(13,'(52e17.9)') (ksive(i,j,band),j=0,nlve+1)
113	c ecart-type MC sur les ksi: pas utilise
114	c read(13,'(52e17.9)') (psive(i,j,band),j=0,nlve+1)
115	c changement de convention (signe) pour ksi,
116	c et prise en compte de la largeur de bande (en cm-1):
117	do j=0,nlve+1
118	ksive(i,j,band) = -ksive(i,j,band)*dlambda
119	enddo
120	enddo
121	c calcul des coeff al et bl pour luminance Planck
122	al(band) = 2.RHPLARCLUMRCLUM/(lambda(band))*5.
123	c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul
124	c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6
125	. * 1.e-6
126	bl(band) = RHPLARCLUM/(RKBOLlambda(band))
127	enddo
128	close(13)
129
130	endif ! firstcall
131
132	c --------------------------------------
133	c Calculation of the Psi matrix
134	c --------------------------------------
135
136	c temperature in the table layers
137	c -------------------------------
138
139	do j=1,nlve
140	ztemp(j) = PT(j)
141	enddo
142
143	ztemp(0) = PTSURF
144
145	c Planck function
146	c ---------------
147
148	do band=1,nnuve
149	y(0,band) = exp(bl(band)/ztemp(0))-1.
150	bplck(0,band) = al(band)/(y(0,band))
151	do j=1,nlve
152	c Developpement en polynomes ?
153	c bplck(j,band) = xp(1,band)
154	c . +ztemp(j)*(xp(2,band)
155	c . +ztemp(j)*(xp(3,band)
156	c . +ztemp(j)*(xp(4,band)
157	c . +ztemp(j)*(xp(5,band)
158	c . +ztemp(j)*(xp(6,band) )))))
159
160	c B(T,l) = al/(exp(bl/T)-1)
161	y(j,band) = exp(bl(band)/ztemp(j))-1.
162	bplck(j,band) = al(band)/(y(j,band))
163	zdblay(j,band) = al(band)bl(band)exp(bl(band)/ztemp(j))/
164	. ((ztemp(j)*2)(y(j,band)**2))
165	enddo
166	bplck(nlve+1,band) = 0.0
167	enddo
168
169	c Calculation of Psi
170	c ------------------
171
172	do band=1,nnuve
173	do j=0,nlve+1
174	do i=0,nlve+1
175	psive(i,j,band)=ksive(i,j,band)*(bplck(i,band)-bplck(j,band))
176	enddo
177	enddo
178	enddo
179
180	do j=0,nlve+1
181	do i=0,nlve+1
182	psi_1(i,j) = 0.0 ! positif quand nrj de i->j
183	enddo
184	enddo
185
186	do band=1,nnuve
187	do j=0,nlve+1
188	do i=0,nlve+1
189	psi_1(i,j) = psi_1(i,j)+psive(i,j,band)
190	enddo
191	enddo
192	enddo
193
194	c Verif...-----------------------
195	c open(11,file="psi.dat")
196	c do i=0,nlve+1
197	c write(11,'(I3,83E17.9)') i,(psi_1(j,i),j=0,nlve+1)
198	c enddo
199	c close(11)
200	c stop
201	c -------------------------------
202
203	c --------------------------
204	c Calculation of the fluxes
205	c --------------------------
206
207	c flux aux intercouches:
208	c zlnet(i+1) est le flux net traversant le plafond de la couche i (+ vers le haut)
209	do p=0,nlve ! numero de la couche
210	zlnet(p+1) = 0.0
211	do j=p+1,nlve+1
212	do i=0,p
213	zlnet(p+1) = zlnet(p+1)+psi_1(i,j)
214	enddo
215	enddo
216	enddo
217
218	c flux net au sol, + vers le haut:
219	PSOLLW = zlnet(1)
220	c flux vers le bas au sol, + vers le bas:
221	PSOLLWDN = 0.0
222	do i=1,nlve+1
223	PSOLLWDN = PSOLLWDN+max(psi_1(i,0),0.0)
224	enddo
225
226	c dfluxnet = radiative budget (W m-2)
227	do p=0,nlve ! numero de la couche
228	dzlnet(p) = 0.0
229	do j=0,nlve+1
230	dzlnet(p) = dzlnet(p)+psi_1(p,j)
231	enddo
232	enddo
233
234
235	c --------------------------------------
236	c Interpolation in the GCM vertical grid
237	c --------------------------------------
238
239	c Flux net
240	c --------
241
242	do j=1,kflev+1
243	ZFLNET(j) = zlnet(j)
244	enddo
245
246	PTOPLW = ZFLNET(kflev+1)
247
248	c Heating rates
249	c -------------
250
251	c cool (K/s) = dfluxnet (W/m2) ! positif quand nrj sort de la couche
252	c *g (m/s2)
253	c /(-dp) (epaisseur couche, en Pa=kg/m/s2)
254	c /cp (J/kg/K)
255
256
257	do j=1,kflev
258	pdp(j)=(PPB(j)-PPB(j+1))*1.e5
259	enddo
260
261	c calcul direct OU calcul par schema implicit
262	if (1.eq.0) then
263	do j=1,kflev
264	PCOOL(j) = dzlnet(j) *RG/RCPD / pdp(j)
265	enddo
266	else
267	call lwi(nlve,nnuve,dzlnet,zdblay,pdp,ksive,PCOOL)
268	endif
269	c print*,dzlnet
270	c print*,pdp
271	c print*,PCOOL
272
273	do j=1,kflev
274	PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday
275	enddo
276
277	firstcall = .false.
278	return
279	end
280

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format