Context Navigation

radlwsw.F @ 1624

Last change on this file since 1624 was 1591, checked in by slebonnois, 9 years ago
SL: Mise a jour de la haute atmosphere, du transfert radiatif (solaire=Rainer Haus; IR: reglages de juin 2016), et implementation des variations de temperature potentielle dans la basse atmosphere (variation de la masse molaire moyenne)
File size: 10.4 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
3	!
4	SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
5	. paprs, pplay,tsol, t)
6	c
7	c======================================================================
8	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
9	c Objet: interface entre le modele et les rayonnements
10	c Arguments:
11	c dist-----input-R- distance astronomique terre-soleil
12	c rmu0-----input-R- cosinus de l'angle zenithal
13	c fract----input-R- duree d'ensoleillement normalisee
14	c zzlev----input-R- altitude a inter-couche (m)
15	c paprs----input-R- pression a inter-couche (Pa)
16	c pplay----input-R- pression au milieu de couche (Pa)
17	c tsol-----input-R- temperature du sol (en K)
18	c t--------input-R- temperature (K)
19
20	c MODIFS pour multimatrices ksi SPECIFIQUE VENUS
21	c S. Lebonnois 20/12/2006
22	c corrections 13/07/2007
23	c New ksi matrix: possibility of different cloud model fct of lat 05/2014
24
25	c With extension NLTE (G. Gilli, 2014)
26
27	c======================================================================
28	use dimphy
29	USE geometry_mod, ONLY: latitude_deg
30	USE phys_state_var_mod, only: heat,cool,radsol,
31	. topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet
32	use write_field_phy
33	IMPLICIT none
34	#include "YOMCST.h"
35	#include "clesphys.h"
36	#include "comcstVE.h"
37	#include "nlteparams.h"
38
39	!===========
40	! Arguments
41	!===========
42	real rmu0(klon), fract(klon), dist
43
44	REAL zzlev(klon,klev+1)
45	real paprs(klon,klev+1), pplay(klon,klev)
46	real tsol(klon)
47	real t(klon,klev)
48
49	!===========
50	! Local
51	!===========
52	INTEGER k, kk, i, j, band
53
54	REAL PPB(klev+1)
55
56	REAL zfract, zrmu0,latdeg
57
58	REAL zheat(klev), zcool(klev)
59	real temp(klev),znivs(klev+1)
60	REAL ZFSNET(klev+1),ZFLNET(klev+1)
61	REAL ztopsw, ztoplw
62	REAL zsolsw, zsollw
63	cIM BEG
64	REAL zsollwdown
65	cIM END
66	real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
67
68	real psi(0:klev+1,0:klev+1)
69	real deltapsi(0:klev+1,0:klev+1)
70	real pt0(0:klev+1)
71	real bplck(0:klev+1,nnuve) ! Planck luminances in table layers
72	real y(0:klev,nnuve) ! temporary variable for Planck
73	real zdblay(0:klev+1,nnuve) ! temperature gradient of planck function
74	integer mat0,lat,ips,isza,ips0,isza0
75	real factp,factz,ksi
76
77	logical firstcall
78	data firstcall/.true./
79	save firstcall
80
81	c-------------------------------------------
82	c Initialisations
83	c-----------------
84
85	if (firstcall) then
86
87	c ---------- ksive --------------
88	allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
89	call load_ksi(ksive)
90
91	endif ! firstcall
92	c-------------------------------------------
93
94	DO k = 1, klev
95	DO i = 1, klon
96	heat(i,k)=0.
97	cool(i,k)=0.
98	ENDDO
99	ENDDO
100
101	c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
102	DO j = 1, klon
103
104	c======================================================================
105	c Initialisations
106	c ---------------
107
108	DO k = 1, klev
109	zheat(k) = 0.0
110	zcool(k) = 0.0
111	ENDDO
112	DO k = 1, klev+1
113	ZFLNET(k) = 0.0
114	ZFSNET(k) = 0.0
115	ENDDO
116	ztopsw = 0.0
117	ztoplw = 0.0
118	zsolsw = 0.0
119	zsollw = 0.0
120	zsollwdown = 0.0
121
122	zfract = fract(j)
123	zrmu0 = rmu0(j)
124
125	DO k = 1, klev+1
126	PPB(k) = paprs(j,k)/1.e5
127	ENDDO
128
129	pt0(0) = tsol(j)
130	DO k = 1, klev
131	pt0(k) = t(j,k)
132	ENDDO
133	pt0(klev+1) = 0.
134
135	DO k = 0,klev+1
136	DO i = 0,klev+1
137	psi(i,k) = 0. ! positif quand nrj de i->k
138	deltapsi(i,k) = 0.
139	ENDDO
140	ENDDO
141
142	c======================================================================
143	c Getting psi and deltapsi
144	c ------------------------
145
146	c Planck function
147	c ---------------
148	do band=1,nnuve
149	do k=0,klev
150	c B(T,l) = al/(exp(bl/T)-1)
151	y(k,band) = exp(bl(band)/pt0(k))-1.
152	bplck(k,band) = al(band)/(y(k,band))
153	zdblay(k,band)= al(band)bl(band)exp(bl(band)/pt0(k))/
154	. ((pt0(k)pt0(k))(y(k,band)*y(k,band)))
155	enddo
156	bplck(klev+1,band) = 0.0
157	zdblay(klev+1,band)= 0.0
158	enddo
159
160	c finding the right matrixes
161	c --------------------------
162	mat0 = 0
163	if (nlatve.eq.1) then ! clouds are taken as uniform
164	lat=1
165	else
166	if (abs(latitude_deg(j)).le.50.) then
167	lat=1
168	elseif (abs(latitude_deg(j)).le.60.) then
169	lat=2
170	elseif (abs(latitude_deg(j)).le.70.) then
171	lat=3
172	elseif (abs(latitude_deg(j)).le.80.) then
173	lat=4
174	else
175	lat=5
176	endif
177	endif
178
179	ips0=0
180	if (nbpsve(lat).gt.1) then
181	do ips=1,nbpsve(lat)-1
182	if ( (psurfve(ips,lat).ge.paprs(j,1))
183	. .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
184	ips0 = ips
185	c print*,'ig=',j,' ips0=',ips
186	factp = (paprs(j,1) -psurfve(ips0,lat))
187	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
188	exit
189	endif
190	enddo
191	else ! Only one ps, no interpolation
192	ips0=1
193	endif
194	isza0=0
195	if (nbszave(lat).gt.1) then
196	do isza=1,nbszave(lat)-1
197	if ( (szave(isza,lat).ge.zrmu0)
198	. .and.(szave(isza+1,lat).lt.zrmu0) ) then
199	isza0 = isza
200	c print*,'ig=',j,' isza0=',isza
201	factz = (zrmu0 -szave(isza0,lat))
202	. /(szave(isza0+1,lat)-szave(isza0,lat))
203	exit
204	endif
205	enddo
206	else ! Only one sza, no interpolation
207	isza0=-99
208	endif
209
210	c -------- Probleme aux bords
211	if ((ips0.eq.0).and.(psurfve(1,lat).gt.paprs(j,1))) then
212	ips0 = 1
213	print*,'Extrapolation! ig=',j,' ips0=',ips0
214	factp = (paprs(j,1) -psurfve(ips0,lat))
215	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
216	endif
217	if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1)))
218	. then
219	ips0 = nbpsve(lat)-1
220	print*,'Extrapolation! ig=',j,' ips0=',ips0
221	factp = (paprs(j,1) -psurfve(ips0,lat))
222	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
223	endif
224	c ---------
225
226	if ((ips0.eq.0).or.(isza0.eq.0)) then
227	write(,) 'Finding the right matrix in radlwsw'
228	print*,'Interpolation problem, grid point ig=',j
229	print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
230	stop
231	endif
232
233	if (isza0.eq.-99) then
234	mat0 = indexve(lat)+ips0
235	else
236	mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0
237	endif
238
239	c interpolation of ksi and computation of psi,deltapsi
240	c ----------------------------------------------------
241	do band=1,nnuve
242	do k=0,klev+1
243	do i=0,klev+1
244	if (isza0.eq.-99) then
245	ksi = ksive(i,k,band,mat0)*(1-factp)
246	. +ksive(i,k,band,mat0+1)*factp
247	else
248	ksi = ksive(i,k,band,mat0)(1-factp)(1-factz)
249	. +ksive(i,k,band,mat0+1)factp (1-factz)
250	. +ksive(i,k,band,mat0+nbpsve(lat))(1-factp)factz
251	. +ksive(i,k,band,mat0+nbpsve(lat)+1)factp factz
252	endif
253
254	psi(i,k) = psi(i,k) +
255	. RPIksi(bplck(i,band)-bplck(k,band))
256	deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
257	enddo
258	enddo
259	enddo
260
261	c======================================================================
262	c LW call
263	c---------
264	temp(1:klev)=t(j,1:klev)
265	CALL LW_venus_ve(
266	. PPB,temp,psi,deltapsi,
267	. zcool,
268	. ztoplw,zsollw,
269	. zsollwdown,ZFLNET)
270
271	c---------
272	c SW call
273	c---------
274	znivs=zzlev(j,:)
275	latdeg=abs(latitude_deg(j))
276
277	c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale
278	c CALL SW_venus_ve(zrmu0, zfract,
279	c S PPB,temp,znivs,
280	c S zheat,
281	c S ztopsw,zsolsw,ZFSNET)
282
283	c CALL SW_venus_cl_1Dglobave(zrmu0,zfract, ! pour moy globale
284	c CALL SW_venus_cl(zrmu0,zfract,
285	c CALL SW_venus_dc_1Dglobave(zrmu0,zfract, ! pour moy globale
286	c CALL SW_venus_dc(zrmu0,zfract,
287	CALL SW_venus_rh(zrmu0,zfract,latdeg,
288	c CALL SW_venus_rh_1Dglobave(zrmu0,zfract, ! pour moy globale
289	S PPB,temp,
290	S zheat,
291	S ztopsw,zsolsw,ZFSNET)
292
293	c======================================================================
294	radsol(j) = zsolsw - zsollw ! + vers bas
295	topsw(j) = ztopsw ! + vers bas
296	toplw(j) = ztoplw ! + vers haut
297	solsw(j) = zsolsw ! + vers bas
298	sollw(j) = -zsollw ! + vers bas
299	sollwdown(j) = zsollwdown ! + vers bas
300
301	DO k = 1, klev+1
302	lwnet (j,k) = ZFLNET(k)
303	swnet (j,k) = ZFSNET(k)
304	ENDDO
305
306	c
307	C heat/cool with upper atmosphere
308	C
309	IF(callnlte) THEN
310	DO k = 1,nlaylte
311	heat(j,k) = zheat(k)
312	cool(j,k) = zcool(k)
313	ENDDO
314	c Zero tendencies for any remaining layers between nlaylte and klev
315	if (klev.gt.nlaylte) then
316	do k = nlaylte+1, klev
317	heat(j,k) = 0.
318	cool(j,k) = 0.
319	enddo
320	endif
321	ELSE
322	DO k = 1, klev
323	heat(j,k) = zheat(k)
324	cool(j,k) = zcool(k)
325	ENDDO
326	ENDIF ! callnlte
327
328	ENDDO ! of DO j = 1, klon
329	c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
330
331	! for tests: write output fields...
332	! call writefield_phy('radlwsw_heat',heat,klev)
333	! call writefield_phy('radlwsw_cool',cool,klev)
334	! call writefield_phy('radlwsw_radsol',radsol,1)
335	! call writefield_phy('radlwsw_topsw',topsw,1)
336	! call writefield_phy('radlwsw_toplw',toplw,1)
337	! call writefield_phy('radlwsw_solsw',solsw,1)
338	! call writefield_phy('radlwsw_sollw',sollw,1)
339	! call writefield_phy('radlwsw_sollwdown',sollwdown,1)
340	! call writefield_phy('radlwsw_swnet',swnet,klev+1)
341	! call writefield_phy('radlwsw_lwnet',lwnet,klev+1)
342
343	c tests
344
345	c j = klon/2
346	c j = 1
347	c print*,'mu0=',rmu0(j)
348	c print*,' net flux vis HEAT(K/Eday)'
349	c do k=1,klev
350	c print,k,ZFSNET(k),heat(j,k)86400.
351	c enddo
352	c print*,' net flux IR COOL(K/Eday)'
353	c do k=1,klev
354	c print,k,ZFLNET(k),cool(j,k)86400.
355	c enddo
356
357	firstcall = .false.
358	RETURN
359	END
360

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

Download in other formats:

Original Format