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

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)
[1301]	57
[1591]	58	REAL zfract, zrmu0,latdeg
[1301]	59
[892]	60	REAL zheat(klev), zcool(klev)
[1301]	61	real temp(klev),znivs(klev+1)
[892]	62	REAL ZFSNET(klev+1),ZFLNET(klev+1)
[3]	63	REAL ztopsw, ztoplw
	64	REAL zsolsw, zsollw
	65	cIM BEG
	66	REAL zsollwdown
	67	cIM END
[101]	68	real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
[953]	69
[892]	70	real psi(0:klev+1,0:klev+1)
	71	real deltapsi(0:klev+1,0:klev+1)
[953]	72	real pt0(0:klev+1)
	73	real bplck(0:klev+1,nnuve) ! Planck luminances in table layers
[1301]	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
[953]	77	real factp,factz,ksi
[1639]	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 --------------------------------------
[3]	84	logical firstcall
	85	data firstcall/.true./
	86	save firstcall
	87
[1639]	88	cERROR ! For checking if the file it's being read
[3]	89	c-------------------------------------------
	90	c Initialisations
	91	c-----------------
	92
	93	if (firstcall) then
[101]	94
	95	c ---------- ksive --------------
[892]	96	allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
[3]	97	call load_ksi(ksive)
	98
	99	endif ! firstcall
[953]	100	c-------------------------------------------
[3]	101
	102	DO k = 1, klev
[953]	103	DO i = 1, klon
[3]	104	heat(i,k)=0.
	105	cool(i,k)=0.
[953]	106	ENDDO
[3]	107	ENDDO
[953]	108
[1639]	109
[3]	110	c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[973]	111	DO j = 1, klon
[953]	112
	113	c======================================================================
	114	c Initialisations
	115	c ---------------
	116
[3]	117	DO k = 1, klev
	118	zheat(k) = 0.0
	119	zcool(k) = 0.0
	120	ENDDO
[1639]	121	c zheat(1:klev)=0.0 !Explicit loop (no change in performance)
	122	c zcool(1:klev)=0.0
	123
[3]	124	DO k = 1, klev+1
	125	ZFLNET(k) = 0.0
	126	ZFSNET(k) = 0.0
	127	ENDDO
[1639]	128	c ZFLNET(1:klev+1)=0.0
	129	c ZFSNET(1:klev+1)=0.0
	130
[3]	131	ztopsw = 0.0
	132	ztoplw = 0.0
	133	zsolsw = 0.0
	134	zsollw = 0.0
	135	zsollwdown = 0.0
	136
[1639]	137	zfract = fract(j)
	138	zrmu0 = rmu0(j)
	139
[953]	140	DO k = 1, klev+1
[3]	141	PPB(k) = paprs(j,k)/1.e5
[953]	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.
[1639]	149
[953]	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
[1639]	156
[3]	157	c======================================================================
[953]	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 --------------------------
[1639]	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
[1301]	260	ips0=0
[1639]	261	ips0A=0
	262	ips0B=0
	263	if (nbpsve(lat).gt.1) then ! Interpolation on ps
[1301]	264	do ips=1,nbpsve(lat)-1
	265	if ( (psurfve(ips,lat).ge.paprs(j,1))
[1639]	266	. .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
[1301]	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))
[953]	271	exit
	272	endif
	273	enddo
[1442]	274	else ! Only one ps, no interpolation
	275	ips0=1
	276	endif
[1639]	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
[1301]	318	isza0=0
[1639]	319	isza0A=0
	320	isza0B=0
[1301]	321	if (nbszave(lat).gt.1) then
	322	do isza=1,nbszave(lat)-1
	323	if ( (szave(isza,lat).ge.zrmu0)
[1639]	324	. .and.(szave(isza+1,lat).lt.zrmu0) ) then
[1301]	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
[1639]	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
[1687]	360	! init to avoid outside values (near midnight so similar compo...)
	361	isza0B = nbszave(latB)
[1639]	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
[1301]	377
[1442]	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
[1639]	385	if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1)))
	386	. then
[1442]	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
[1301]	394	if ((ips0.eq.0).or.(isza0.eq.0)) then
[953]	395	write(,) 'Finding the right matrix in radlwsw'
[1301]	396	print*,'Interpolation problem, grid point ig=',j
	397	print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
[953]	398	stop
	399	endif
[1639]	400
[1301]	401	if (isza0.eq.-99) then
[1639]	402	mat0 = indexve(lat) +ips0
[1642]	403	if (latA.gt.0) then
	404	mat0A = indexve(latA)+ips0A
	405	mat0B = indexve(latB)+ips0B
	406	endif
[1301]	407	else
[1639]	408	mat0 = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
[1642]	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
[1301]	413	endif
[1639]	414
	415	c write(,) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
	416
[953]	417	c interpolation of ksi and computation of psi,deltapsi
	418	c ----------------------------------------------------
	419
[1639]	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
[953]	535	c======================================================================
[3]	536	c LW call
	537	c---------
[973]	538	temp(1:klev)=t(j,1:klev)
[3]	539	CALL LW_venus_ve(
[973]	540	. PPB,temp,psi,deltapsi,
[3]	541	. zcool,
	542	. ztoplw,zsollw,
	543	. zsollwdown,ZFLNET)
	544
	545	c---------
	546	c SW call
	547	c---------
[1301]	548	znivs=zzlev(j,:)
[1591]	549	latdeg=abs(latitude_deg(j))
	550
[1442]	551	c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale
[1301]	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
[1591]	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
[973]	563	S PPB,temp,
[3]	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
[892]	575	DO k = 1, klev+1
[3]	576	lwnet (j,k) = ZFLNET(k)
	577	swnet (j,k) = ZFSNET(k)
	578	ENDDO
	579
	580	c
[1310]	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
[973]	602	ENDDO ! of DO j = 1, klon
[3]	603	c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[1310]	604
[973]	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)
[3]	616
	617	c tests
	618
	619	c j = klon/2
	620	c j = 1
	621	c print*,'mu0=',rmu0(j)
[1301]	622	c print*,' net flux vis HEAT(K/Eday)'
[892]	623	c do k=1,klev
[1301]	624	c print,k,ZFSNET(k),heat(j,k)86400.
[3]	625	c enddo
[1301]	626	c print*,' net flux IR COOL(K/Eday)'
[892]	627	c do k=1,klev
[1301]	628	c print,k,ZFLNET(k),cool(j,k)86400.
[3]	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