Context Navigation

radlwsw.F @ 3803

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

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

Download in other formats:

Original Format