Context Navigation

radlwsw.F @ 1945

Last change on this file since 1945 was 1945, checked in by slebonnois, 7 years ago
Correction of a bug in radlwsw Venus
File size: 19.3 KB

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

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

Download in other formats:

Original Format