Context Navigation

radlwsw.F @ 3094

Last change on this file since 3094 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

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======================================================================
[2534]	30	use dimphy, only: klon,klev
[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
[2560]	35	use radinc_h, only: ini_radinc_h
[2534]	36
	37	#ifdef CPP_XIOS
	38	use xios_output_mod, only: send_xios_field
	39	#endif
	40
[101]	41	IMPLICIT none
[2534]	42	include "YOMCST.h"
	43	include "clesphys.h"
	44	include "comcstVE.h"
	45	include "nlteparams.h"
[1301]	46
	47	!===========
	48	! Arguments
	49	!===========
[2534]	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)
[1301]	58
	59	!===========
	60	! Local
	61	!===========
[953]	62	INTEGER k, kk, i, j, band
[2560]	63	integer,save :: i_sw
	64	integer,parameter :: subloop=100
[1301]	65
[892]	66	REAL PPB(klev+1)
[1723]	67	REAL PPA(klev)
[1301]	68
[1591]	69	REAL zfract, zrmu0,latdeg
[1301]	70
[892]	71	REAL zheat(klev), zcool(klev)
[1301]	72	real temp(klev),znivs(klev+1)
[892]	73	REAL ZFSNET(klev+1),ZFLNET(klev+1)
[3]	74	REAL ztopsw, ztoplw
	75	REAL zsolsw, zsollw
	76	cIM BEG
	77	REAL zsollwdown
	78	cIM END
[101]	79	real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
[953]	80
[892]	81	real psi(0:klev+1,0:klev+1)
	82	real deltapsi(0:klev+1,0:klev+1)
[953]	83	real pt0(0:klev+1)
	84	real bplck(0:klev+1,nnuve) ! Planck luminances in table layers
[1301]	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
[953]	88	real factp,factz,ksi
[1639]	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 --------------------------------------
[2534]	95	logical,save :: firstcall=.true.
[3]	96
[1639]	97	cERROR ! For checking if the file it's being read
[3]	98	c-------------------------------------------
	99	c Initialisations
	100	c-----------------
	101
	102	if (firstcall) then
[101]	103
	104	c ---------- ksive --------------
[892]	105	allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
[3]	106	call load_ksi(ksive)
	107
[2560]	108	i_sw = subloop
	109	c for sw_venus_corrk
	110	if (solarchoice.eq.2) call ini_radinc_h(klev)
	111
[3]	112	endif ! firstcall
[953]	113	c-------------------------------------------
[3]	114
	115	DO k = 1, klev
[953]	116	DO i = 1, klon
[2560]	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.
[3]	120	cool(i,k)=0.
[953]	121	ENDDO
[3]	122	ENDDO
[953]	123
[1639]	124
[3]	125	c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[973]	126	DO j = 1, klon
[953]	127
	128	c======================================================================
	129	c Initialisations
	130	c ---------------
	131
[3]	132	DO k = 1, klev
	133	zheat(k) = 0.0
	134	zcool(k) = 0.0
	135	ENDDO
[1639]	136	c zheat(1:klev)=0.0 !Explicit loop (no change in performance)
	137	c zcool(1:klev)=0.0
	138
[3]	139	DO k = 1, klev+1
	140	ZFLNET(k) = 0.0
	141	ZFSNET(k) = 0.0
	142	ENDDO
[1639]	143	c ZFLNET(1:klev+1)=0.0
	144	c ZFSNET(1:klev+1)=0.0
	145
[3]	146	ztopsw = 0.0
	147	ztoplw = 0.0
	148	zsolsw = 0.0
	149	zsollw = 0.0
	150	zsollwdown = 0.0
	151
[1639]	152	zfract = fract(j)
	153	zrmu0 = rmu0(j)
	154
[953]	155	DO k = 1, klev+1
[3]	156	PPB(k) = paprs(j,k)/1.e5
[953]	157	ENDDO
[1723]	158	DO k = 1,klev
	159	PPA(k) = pplay(j,k)/1.e5
	160	ENDDO
[953]	161
	162	pt0(0) = tsol(j)
	163	DO k = 1, klev
	164	pt0(k) = t(j,k)
	165	ENDDO
	166	pt0(klev+1) = 0.
[1639]	167
[953]	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
[1639]	174
[3]	175	c======================================================================
[953]	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 --------------------------
[1639]	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
[1301]	278	ips0=0
[1639]	279	ips0A=0
	280	ips0B=0
	281	if (nbpsve(lat).gt.1) then ! Interpolation on ps
[1301]	282	do ips=1,nbpsve(lat)-1
	283	if ( (psurfve(ips,lat).ge.paprs(j,1))
[1639]	284	. .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
[1301]	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))
[953]	289	exit
	290	endif
	291	enddo
[1442]	292	else ! Only one ps, no interpolation
	293	ips0=1
	294	endif
[1639]	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
[1301]	336	isza0=0
[1639]	337	isza0A=0
	338	isza0B=0
[1301]	339	if (nbszave(lat).gt.1) then
	340	do isza=1,nbszave(lat)-1
	341	if ( (szave(isza,lat).ge.zrmu0)
[1639]	342	. .and.(szave(isza+1,lat).lt.zrmu0) ) then
[1301]	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
[1639]	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
[1687]	378	! init to avoid outside values (near midnight so similar compo...)
	379	isza0B = nbszave(latB)
[1639]	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
[1301]	395
[1442]	396	c -------- Probleme aux bords
[1945]	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
[1442]	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
[1945]	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
[1442]	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
[1945]	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
[1442]	441	c ---------
	442
[1301]	443	if ((ips0.eq.0).or.(isza0.eq.0)) then
[953]	444	write(,) 'Finding the right matrix in radlwsw'
[1301]	445	print*,'Interpolation problem, grid point ig=',j
	446	print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
[953]	447	stop
	448	endif
[1639]	449
[1301]	450	if (isza0.eq.-99) then
[1639]	451	mat0 = indexve(lat) +ips0
[1642]	452	if (latA.gt.0) then
	453	mat0A = indexve(latA)+ips0A
	454	mat0B = indexve(latB)+ips0B
	455	endif
[1301]	456	else
[1639]	457	mat0 = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
[1642]	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
[1301]	462	endif
[1639]	463
	464	c write(,) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
	465
[953]	466	c interpolation of ksi and computation of psi,deltapsi
	467	c ----------------------------------------------------
	468
[1639]	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
[953]	584	c======================================================================
[3]	585	c LW call
	586	c---------
[973]	587	temp(1:klev)=t(j,1:klev)
[3]	588	CALL LW_venus_ve(
[973]	589	. PPB,temp,psi,deltapsi,
[3]	590	. zcool,
	591	. ztoplw,zsollw,
	592	. zsollwdown,ZFLNET)
	593	c---------
	594	c SW call
	595	c---------
[1301]	596	znivs=zzlev(j,:)
[1591]	597	latdeg=abs(latitude_deg(j))
	598
[1442]	599	c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale
[1301]	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
[1591]	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,
[1726]	609	c CALL SW_venus_rh_1Dglobave(zrmu0,zfract, ! pour moy globale
	610	c S PPB,temp,
[2560]	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,
[1723]	617	S PPA,PPB,temp,
[3]	618	S zheat,
	619	S ztopsw,zsolsw,ZFSNET)
[2560]	620	endif ! solarchoice.eq.1
	621
[3]	622	c======================================================================
[2560]	623	c lw into klon grid:
	624	radsol(j) = - zsollw ! + vers bas
[3]	625	toplw(j) = ztoplw ! + vers haut
	626	sollw(j) = -zsollw ! + vers bas
	627	sollwdown(j) = zsollwdown ! + vers bas
	628
[892]	629	DO k = 1, klev+1
[2560]	630	lwnet (j,k) = ZFLNET(k) ! + vers haut
[3]	631	ENDDO
	632
[2560]	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
[3]	644	c
[2560]	645	C cool with upper atmosphere
[1310]	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
[2560]	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
[973]	683	ENDDO ! of DO j = 1, klon
[3]	684	c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[1310]	685
[2560]	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
[973]	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)
[3]	725
	726	c tests
	727
	728	c j = klon/2
	729	c j = 1
	730	c print*,'mu0=',rmu0(j)
[1301]	731	c print*,' net flux vis HEAT(K/Eday)'
[892]	732	c do k=1,klev
[1301]	733	c print,k,ZFSNET(k),heat(j,k)86400.
[3]	734	c enddo
[1301]	735	c print*,' net flux IR COOL(K/Eday)'
[892]	736	c do k=1,klev
[1301]	737	c print,k,ZFLNET(k),cool(j,k)86400.
[3]	738	c enddo
	739
	740	firstcall = .false.
[2534]	741
[3]	742	END

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

Download in other formats:

Original Format