Context Navigation

radlwsw.F @ 1704

Last change on this file since 1704 was 1687, checked in by slebonnois, 8 years ago
SL: bugs corrections outputs/clouds_AStol/upper_atmosphere/start2archive
File size: 18.3 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 Ksi matrices latitudinaly interpolated (I. Garate-Lopez, 2016)
28
29	c======================================================================
30	use dimphy
31	USE geometry_mod, ONLY: latitude_deg
32	USE phys_state_var_mod, only: heat,cool,radsol,
33	. topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet
34	use write_field_phy
35	IMPLICIT none
36	#include "YOMCST.h"
37	#include "clesphys.h"
38	#include "comcstVE.h"
39	#include "nlteparams.h"
40
41	!===========
42	! Arguments
43	!===========
44	real rmu0(klon), fract(klon), dist
45
46	REAL zzlev(klon,klev+1)
47	real paprs(klon,klev+1), pplay(klon,klev)
48	real tsol(klon)
49	real t(klon,klev)
50
51	!===========
52	! Local
53	!===========
54	INTEGER k, kk, i, j, band
55
56	REAL PPB(klev+1)
57
58	REAL zfract, zrmu0,latdeg
59
60	REAL zheat(klev), zcool(klev)
61	real temp(klev),znivs(klev+1)
62	REAL ZFSNET(klev+1),ZFLNET(klev+1)
63	REAL ztopsw, ztoplw
64	REAL zsolsw, zsollw
65	cIM BEG
66	REAL zsollwdown
67	cIM END
68	real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
69
70	real psi(0:klev+1,0:klev+1)
71	real deltapsi(0:klev+1,0:klev+1)
72	real pt0(0:klev+1)
73	real bplck(0:klev+1,nnuve) ! Planck luminances in table layers
74	real y(0:klev,nnuve) ! temporary variable for Planck
75	real zdblay(0:klev+1,nnuve) ! temperature gradient of planck function
76	integer mat0,lat,ips,isza,ips0,isza0
77	real factp,factz,ksi
78	c ------- for lat-interp ----------------
79	integer mat0A, mat0B, latA, latB, kasua
80	integer ipsA, ipsB, iszaA, iszaB, ips0A, ips0B, isza0A, isza0B
81	real lat_deg, latA_deg, latB_deg
82	real factlat, k1, k2, k3, k4
83	c --------------------------------------
84	logical firstcall
85	data firstcall/.true./
86	save firstcall
87
88	cERROR ! For checking if the file it's being read
89	c-------------------------------------------
90	c Initialisations
91	c-----------------
92
93	if (firstcall) then
94
95	c ---------- ksive --------------
96	allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
97	call load_ksi(ksive)
98
99	endif ! firstcall
100	c-------------------------------------------
101
102	DO k = 1, klev
103	DO i = 1, klon
104	heat(i,k)=0.
105	cool(i,k)=0.
106	ENDDO
107	ENDDO
108
109
110	c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
111	DO j = 1, klon
112
113	c======================================================================
114	c Initialisations
115	c ---------------
116
117	DO k = 1, klev
118	zheat(k) = 0.0
119	zcool(k) = 0.0
120	ENDDO
121	c zheat(1:klev)=0.0 !Explicit loop (no change in performance)
122	c zcool(1:klev)=0.0
123
124	DO k = 1, klev+1
125	ZFLNET(k) = 0.0
126	ZFSNET(k) = 0.0
127	ENDDO
128	c ZFLNET(1:klev+1)=0.0
129	c ZFSNET(1:klev+1)=0.0
130
131	ztopsw = 0.0
132	ztoplw = 0.0
133	zsolsw = 0.0
134	zsollw = 0.0
135	zsollwdown = 0.0
136
137	zfract = fract(j)
138	zrmu0 = rmu0(j)
139
140	DO k = 1, klev+1
141	PPB(k) = paprs(j,k)/1.e5
142	ENDDO
143
144	pt0(0) = tsol(j)
145	DO k = 1, klev
146	pt0(k) = t(j,k)
147	ENDDO
148	pt0(klev+1) = 0.
149
150	DO k = 0,klev+1
151	DO i = 0,klev+1
152	psi(i,k) = 0. ! positif quand nrj de i->k
153	deltapsi(i,k) = 0.
154	ENDDO
155	ENDDO
156
157	c======================================================================
158	c Getting psi and deltapsi
159	c ------------------------
160
161	c Planck function
162	c ---------------
163	do band=1,nnuve
164	do k=0,klev
165	c B(T,l) = al/(exp(bl/T)-1)
166	y(k,band) = exp(bl(band)/pt0(k))-1.
167	bplck(k,band) = al(band)/(y(k,band))
168	zdblay(k,band)= al(band)bl(band)exp(bl(band)/pt0(k))/
169	. ((pt0(k)pt0(k))(y(k,band)*y(k,band)))
170	enddo
171	bplck(klev+1,band) = 0.0
172	zdblay(klev+1,band)= 0.0
173	enddo
174
175	c finding the right matrixes
176	c --------------------------
177
178	mat0 = 0
179	mat0A = 0
180	mat0B = 0
181
182	c Latitude
183	c --------
184	lat = 0
185	latA = 0
186	latB = 0
187
188	c write(,) 'nlatve:', nlatve
189
190	lat_deg = abs(latitude_deg(j))
191
192	c if (nlatve.eq.1) then ! clouds are taken as uniform
193	if ((nlatve.eq.1).or.(lat_deg.le.25.)) then
194	lat = 1
195	elseif (lat_deg.le.50.) then
196	lat = 1
197	latA = 1
198	latB = 2
199	latA_deg = 25.0
200	latB_deg = 55.0
201	elseif (lat_deg.le.55.) then
202	lat = 2
203	latA = 1
204	latB = 2
205	latA_deg = 25.0
206	latB_deg = 55.0
207	elseif (lat_deg.le.60.) then
208	lat = 2
209	latA = 2
210	latB = 3
211	latA_deg = 55.0
212	latB_deg = 65.0
213	elseif (lat_deg.le.65.) then
214	lat = 3
215	latA = 2
216	latB = 3
217	latA_deg = 55.0
218	latB_deg = 65.0
219	elseif (lat_deg.le.70.) then
220	lat = 3
221	latA = 3
222	latB = 4
223	latA_deg = 65.0
224	latB_deg = 75.0
225	elseif (lat_deg.le.75.) then
226	lat = 4
227	latA = 3
228	latB = 4
229	latA_deg = 65.0
230	latB_deg = 75.0
231	elseif (lat_deg.le.80.) then
232	lat = 4
233	latA = 4
234	latB = 5
235	latA_deg = 75.0
236	latB_deg = 85.0
237	elseif (lat_deg.le.85.) then
238	lat = 5
239	latA = 4
240	latB = 5
241	latA_deg = 75.0
242	latB_deg = 85.0
243	else
244	lat = 5
245	endif
246
247	c write(,) 'Lat',lat,'LatA',latA,'LatB',latB
248
249	factlat = 0
250	if (latA.gt.0) then
251	factlat = (lat_deg - latA_deg) / (latB_deg - latA_deg)
252	endif
253
254	c write (,) 'Factor de correccion:', factlat
255
256
257	c Pressure at Surface
258	c -------------------
259
260	ips0=0
261	ips0A=0
262	ips0B=0
263	if (nbpsve(lat).gt.1) then ! Interpolation on ps
264	do ips=1,nbpsve(lat)-1
265	if ( (psurfve(ips,lat).ge.paprs(j,1))
266	. .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
267	ips0 = ips
268	c print*,'ig=',j,' ips0=',ips
269	factp = (paprs(j,1) -psurfve(ips0,lat))
270	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
271	exit
272	endif
273	enddo
274	else ! Only one ps, no interpolation
275	ips0=1
276	endif
277
278	if (latA.eq.lat) then
279	ips0A=ips0
280	else
281	if (latA.gt.0) then
282	if (nbpsve(latA).gt.1) then
283	do ipsA=1,nbpsve(latA)-1
284	if ( (psurfve(ipsA,latA).ge.paprs(j,1))
285	. .and.(psurfve(ipsA+1,latA).lt.paprs(j,1)) ) then
286	ips0A = ipsA
287	exit
288	endif
289	enddo
290	else ! Only one ps, no interpolation
291	ips0A=1
292	endif ! nbpsve(latA).gt.1
293	endif ! latA.gt.0 (if latA=0 ips0A is not used, so it doesn't matter)
294	endif ! latA.eq.lat
295
296	if (latB.eq.lat) then
297	ips0B=ips0
298	else
299	if (latB.gt.0) then
300	if (nbpsve(latB).gt.1) then
301	do ipsB=1,nbpsve(latB)-1
302	if ( (psurfve(ipsB,latB).ge.paprs(j,1))
303	. .and.(psurfve(ipsB+1,latB).lt.paprs(j,1)) ) then
304	ips0B = ipsB
305	exit
306	endif
307	enddo
308	else
309	ips0B=1
310	endif ! nbpsve(latB).gt.1
311	endif ! latB.gt.0 (if latB=0 ips0B is not used, so it doesn't matter)
312	endif ! latB.eq.lat
313
314
315	c Solar Zenith Angle
316	c ------------------
317
318	isza0=0
319	isza0A=0
320	isza0B=0
321	if (nbszave(lat).gt.1) then
322	do isza=1,nbszave(lat)-1
323	if ( (szave(isza,lat).ge.zrmu0)
324	. .and.(szave(isza+1,lat).lt.zrmu0) ) then
325	isza0 = isza
326	c print*,'ig=',j,' isza0=',isza
327	factz = (zrmu0 -szave(isza0,lat))
328	. /(szave(isza0+1,lat)-szave(isza0,lat))
329	exit
330	endif
331	enddo
332	else ! Only one sza, no interpolation
333	isza0=-99
334	endif
335
336
337	if (latA.eq.lat) then
338	isza0A=isza0
339	else
340	if (latA.gt.0) then
341	if (nbszave(latA).gt.1) then
342	do iszaA=1,nbszave(latA)-1
343	if ( (szave(iszaA,latA).ge.zrmu0)
344	. .and.(szave(iszaA+1,latA).lt.zrmu0) ) then
345	isza0A = iszaA
346	exit
347	endif
348	enddo
349	else ! Only one sza, no interpolation
350	isza0A=-99
351	endif ! nbszave(latA).gt.1
352	endif ! latA.gt.0 (if latA=0 isza0A is not used, so it doesn't matter)
353	endif ! latA.eq.lat
354
355	if (latB.eq.lat) then
356	isza0B=isza0
357	else
358	if (latB.gt.0) then
359	if (nbszave(latB).gt.1) then
360	! init to avoid outside values (near midnight so similar compo...)
361	isza0B = nbszave(latB)
362	do iszaB=1,nbszave(latB)-1
363	if ( (szave(iszaB,latB).ge.zrmu0)
364	. .and.(szave(iszaB+1,latB).lt.zrmu0) ) then
365	isza0B = iszaB
366	exit
367	endif
368	enddo
369	else ! Only one sza, no interpolation
370	isza0B=-99
371	endif ! nbszave(latB).gt.1
372	endif ! latB.gt.0 (if latB=0 isza0B is not used, so it doesn't matter)
373	endif ! latB.eq.lat
374
375	c write(,) 'nbszave', nbszave(lat),'nbpsve(lat)',nbpsve(lat)
376
377
378	c -------- Probleme aux bords
379	if ((ips0.eq.0).and.(psurfve(1,lat).gt.paprs(j,1))) then
380	ips0 = 1
381	print*,'Extrapolation! ig=',j,' ips0=',ips0
382	factp = (paprs(j,1) -psurfve(ips0,lat))
383	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
384	endif
385	if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1)))
386	. then
387	ips0 = nbpsve(lat)-1
388	print*,'Extrapolation! ig=',j,' ips0=',ips0
389	factp = (paprs(j,1) -psurfve(ips0,lat))
390	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
391	endif
392	c ---------
393
394	if ((ips0.eq.0).or.(isza0.eq.0)) then
395	write(,) 'Finding the right matrix in radlwsw'
396	print*,'Interpolation problem, grid point ig=',j
397	print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
398	stop
399	endif
400
401	if (isza0.eq.-99) then
402	mat0 = indexve(lat) +ips0
403	if (latA.gt.0) then
404	mat0A = indexve(latA)+ips0A
405	mat0B = indexve(latB)+ips0B
406	endif
407	else
408	mat0 = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
409	if (latA.gt.0) then
410	mat0A = indexve(latA)+(isza0A-1)*nbpsve(latA)+ips0A
411	mat0B = indexve(latB)+(isza0B-1)*nbpsve(latB)+ips0B
412	endif
413	endif
414
415	c write(,) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
416
417	c interpolation of ksi and computation of psi,deltapsi
418	c ----------------------------------------------------
419
420	if (isza0.eq.-99) then
421	if (latA.gt.0) then ! Not being in the two extremal bins
422
423	do band=1,nnuve
424	do k=0,klev+1
425	do i=k+1,klev+1
426	k1 = ksive(i,k,band,mat0A)*(1-factlat)
427	. + ksive(i,k,band,mat0B)*factlat
428	k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
429	. + ksive(i,k,band,mat0B+1)*factlat
430	ksi = k1(1-factp) + k2factp
431	psi(i,k) = psi(i,k) +
432	. RPIksi(bplck(i,band)-bplck(k,band))
433	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
434	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
435	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
436
437	kasua=1
438	enddo
439	enddo
440	enddo
441	do k=0,klev+1
442	do i=k+1,klev+1
443	psi(k,i) = -psi(i,k)
444	enddo
445	enddo
446
447	else ! latA=0 --> extremal bins
448
449	do band=1,nnuve
450	do k=0,klev+1
451	do i=k+1,klev+1
452	ksi = ksive(i,k,band,mat0)*(1-factp)
453	. + ksive(i,k,band,mat0+1)*factp
454	psi(i,k) = psi(i,k) +
455	. RPIksi(bplck(i,band)-bplck(k,band))
456	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
457	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
458	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
459
460	kasua=2
461	enddo
462	enddo
463	enddo
464	do k=0,klev+1
465	do i=k+1,klev+1
466	psi(k,i) = -psi(i,k)
467	enddo
468	enddo
469
470	endif ! latA.gt.0
471
472	else ! isza0=!-99
473
474	if (latA.gt.0) then ! Not being in the two extremal bins
475
476	do band=1,nnuve
477	do k=0,klev+1
478	do i=k+1,klev+1
479	k1 = ksive(i,k,band,mat0A)*(1-factlat)
480	. + ksive(i,k,band,mat0B)*factlat
481	k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
482	. + ksive(i,k,band,mat0B+1)*factlat
483	k3 = ksive(i,k,band,mat0A+nbpsve(latA))*(1-factlat)
484	. + ksive(i,k,band,mat0B+nbpsve(latB))*factlat
485	k4 = ksive(i,k,band,mat0A+nbpsve(latA)+1)*(1-factlat)
486	. + ksive(i,k,band,mat0B+nbpsve(latB)+1)*factlat
487	ksi = ( k1(1-factp) + k2factp )*(1-factz)
488	. + ( k3(1-factp) + k4factp )*factz
489	psi(i,k) = psi(i,k) +
490	. RPIksi(bplck(i,band)-bplck(k,band))
491	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
492	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
493	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
494
495	kasua=3
496	enddo
497	enddo
498	enddo
499	do k=0,klev+1
500	do i=k+1,klev+1
501	psi(k,i) = -psi(i,k)
502	enddo
503	enddo
504
505	else ! latA=0 --> extremal bins
506
507	do band=1,nnuve
508	do k=0,klev+1
509	do i=k+1,klev+1
510	ksi = ksive(i,k,band,mat0)(1-factp)(1-factz)
511	. + ksive(i,k,band,mat0+1)factp (1-factz)
512	. + ksive(i,k,band,mat0+nbpsve(lat))(1-factp)factz
513	. + ksive(i,k,band,mat0+nbpsve(lat)+1)factp factz
514	psi(i,k) = psi(i,k) +
515	. RPIksi(bplck(i,band)-bplck(k,band))
516	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
517	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
518	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
519
520	kasua=4
521	enddo
522	enddo
523	enddo
524	do k=0,klev+1
525	do i=k+1,klev+1
526	psi(k,i) = -psi(i,k)
527	enddo
528	enddo
529
530	endif ! latA.gt.0
531	endif ! isza0.eq.-99
532
533	c write(,) 'Kasua:', kasua
534
535	c======================================================================
536	c LW call
537	c---------
538	temp(1:klev)=t(j,1:klev)
539	CALL LW_venus_ve(
540	. PPB,temp,psi,deltapsi,
541	. zcool,
542	. ztoplw,zsollw,
543	. zsollwdown,ZFLNET)
544
545	c---------
546	c SW call
547	c---------
548	znivs=zzlev(j,:)
549	latdeg=abs(latitude_deg(j))
550
551	c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale
552	c CALL SW_venus_ve(zrmu0, zfract,
553	c S PPB,temp,znivs,
554	c S zheat,
555	c S ztopsw,zsolsw,ZFSNET)
556
557	c CALL SW_venus_cl_1Dglobave(zrmu0,zfract, ! pour moy globale
558	c CALL SW_venus_cl(zrmu0,zfract,
559	c CALL SW_venus_dc_1Dglobave(zrmu0,zfract, ! pour moy globale
560	c CALL SW_venus_dc(zrmu0,zfract,
561	CALL SW_venus_rh(zrmu0,zfract,latdeg,
562	c CALL SW_venus_rh_1Dglobave(zrmu0,zfract, ! pour moy globale
563	S PPB,temp,
564	S zheat,
565	S ztopsw,zsolsw,ZFSNET)
566
567	c======================================================================
568	radsol(j) = zsolsw - zsollw ! + vers bas
569	topsw(j) = ztopsw ! + vers bas
570	toplw(j) = ztoplw ! + vers haut
571	solsw(j) = zsolsw ! + vers bas
572	sollw(j) = -zsollw ! + vers bas
573	sollwdown(j) = zsollwdown ! + vers bas
574
575	DO k = 1, klev+1
576	lwnet (j,k) = ZFLNET(k)
577	swnet (j,k) = ZFSNET(k)
578	ENDDO
579
580	c
581	C heat/cool with upper atmosphere
582	C
583	IF(callnlte) THEN
584	DO k = 1,nlaylte
585	heat(j,k) = zheat(k)
586	cool(j,k) = zcool(k)
587	ENDDO
588	c Zero tendencies for any remaining layers between nlaylte and klev
589	if (klev.gt.nlaylte) then
590	do k = nlaylte+1, klev
591	heat(j,k) = 0.
592	cool(j,k) = 0.
593	enddo
594	endif
595	ELSE
596	DO k = 1, klev
597	heat(j,k) = zheat(k)
598	cool(j,k) = zcool(k)
599	ENDDO
600	ENDIF ! callnlte
601
602	ENDDO ! of DO j = 1, klon
603	c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
604
605	! for tests: write output fields...
606	! call writefield_phy('radlwsw_heat',heat,klev)
607	! call writefield_phy('radlwsw_cool',cool,klev)
608	! call writefield_phy('radlwsw_radsol',radsol,1)
609	! call writefield_phy('radlwsw_topsw',topsw,1)
610	! call writefield_phy('radlwsw_toplw',toplw,1)
611	! call writefield_phy('radlwsw_solsw',solsw,1)
612	! call writefield_phy('radlwsw_sollw',sollw,1)
613	! call writefield_phy('radlwsw_sollwdown',sollwdown,1)
614	! call writefield_phy('radlwsw_swnet',swnet,klev+1)
615	! call writefield_phy('radlwsw_lwnet',lwnet,klev+1)
616
617	c tests
618
619	c j = klon/2
620	c j = 1
621	c print*,'mu0=',rmu0(j)
622	c print*,' net flux vis HEAT(K/Eday)'
623	c do k=1,klev
624	c print,k,ZFSNET(k),heat(j,k)86400.
625	c enddo
626	c print*,' net flux IR COOL(K/Eday)'
627	c do k=1,klev
628	c print,k,ZFLNET(k),cool(j,k)86400.
629	c enddo
630
631	firstcall = .false.
632	RETURN
633	END

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

Download in other formats:

Original Format