Context Navigation

radlwsw.F @ 3556

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

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

Download in other formats:

Original Format