Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sw_venus_corrk.F90 @ 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: 23.3 KB

Line
1	subroutine sw_venus_corrk(ngrid,nlayer, &
2	mu0,pplev,pplay,pt,tsurf,fract,dist_star, &
3	dtsw,nfluxsurf,nfluxtop,netflux,firstcall)
4
5	use mod_phys_lmdz_para, only : is_master
6	use radinc_h, only: NBinfrared, NBvisible, L_NSPECTV, naerkind,&
7	L_LEVELS, L_NGAUSS, L_NLEVRAD, L_NLAYRAD, L_REFVAR
8	use radcommon_h, only: fzerov, gasv, &
9	gweight, pfgasref, pgasmax, pgasmin, &
10	pgasref, tgasmax, tgasmin, tgasref, scalep, &
11	stellarf, dwnv, tauray
12	use datafile_mod, only: datadir,banddir,corrkdir
13	use ioipsl_getin_p_mod, only: getin_p
14	use gases_h, only: ngasmx
15	use optcv_mod, only: optcv
16	use cpdet_phy_mod, only: cpdet
17
18	implicit none
19	#include "YOMCST.h"
20	#include "clesphys.h"
21
22	!==================================================================
23	!
24	! Purpose
25	! -------
26	! Solve the radiative transfer using the correlated-k method for
27	! the gaseous absorption and the Toon et al. (1989) method for
28	! scatttering due to aerosols.
29	!
30	! Based on callcorrk (Generic GCM)
31	! adapted for the SW in the Venus GCM (S. Lebonnois, summer 2020)
32	!==================================================================
33
34	!-----------------------------------------------------------------------
35	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
36	! Layer #1 is the layer near the ground.
37	! Layer #nlayer is the layer at the top.
38	!-----------------------------------------------------------------------
39
40
41	! INPUT
42	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
43	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
44	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
45	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
46	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
47	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
48	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
49	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
50	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
51	logical,intent(in) :: firstcall ! Signals first call to physics.
52
53	! OUTPUT
54	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
55	REAL,INTENT(OUT) :: nfluxsurf(ngrid) ! Net SW flux absorbed at surface (W/m2).
56	REAL,INTENT(OUT) :: nfluxtop(ngrid) ! Net incident top of atmosphere SW flux (W/m2).
57	REAL,INTENT(OUT) :: netflux(ngrid,nlayer) ! net SW flux (W/m2).
58
59	! interesting variables
60	! albedo could be an input map... For future adaptation
61	REAL :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
62	! potential outputs
63	REAL :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
64	REAL :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
65	REAL :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
66	REAL :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
67
68	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
69	REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
70	REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
71	REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
72
73	REAL :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
74	REAL :: tau_col(ngrid) ! Diagnostic from aeropacity.
75
76	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
77	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
78	!$OMP THREADPRIVATE(reffrad,nueffrad)
79
80	!-----------------------------------------------------------------------
81	! Declaration of the variables required by correlated-k subroutines
82	! Numbered from top to bottom (unlike in the GCM)
83	!-----------------------------------------------------------------------
84
85	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
86	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
87
88	! Optical values for the optci/cv subroutines
89	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
90	! NB: Arrays below are "save" to avoid reallocating them at every call
91	! not because their content needs be reused from call to the next
92	REAL*8,allocatable,save :: dtauv(:,:,:)
93	REAL*8,allocatable,save :: cosbv(:,:,:)
94	REAL*8,allocatable,save :: wbarv(:,:,:)
95	!$OMP THREADPRIVATE(dtauv,cosbv,wbarv)
96	REAL*8,allocatable,save :: tauv(:,:,:)
97	REAL*8,allocatable,save :: taucumv(:,:,:)
98	!$OMP THREADPRIVATE(tauv,taucumv)
99	REAL*8 tauaero(L_LEVELS,naerkind)
100	REAL*8 nfluxtopv,fluxtopvdn
101	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
102	REAL*8 fmnetv(L_NLAYRAD)
103	REAL*8 fluxupv(L_NLAYRAD)
104	REAL*8 fluxdnv(L_NLAYRAD)
105	REAL*8 acosz
106	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
107
108	INTEGER ig,l,k,nw,iaer
109
110	real*8,allocatable,save :: taugsurf(:,:)
111	!$OMP THREADPRIVATE(taugsurf)
112	real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol).
113
114	! Local aerosol optical properties for each column on RADIATIVE grid.
115	real8,save,allocatable :: QXVAER(:,:,:) ! Extinction coeff (QVISsQREFQREFvis)
116	real*8,save,allocatable :: QSVAER(:,:,:)
117	real*8,save,allocatable :: GVAER(:,:,:)
118	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER)
119	real, dimension(:,:,:), save, allocatable :: QREFvis3d
120	!$OMP THREADPRIVATE(QREFvis3d)
121
122
123	! Miscellaneous :
124	real*8 temp,temp1,temp2,pweight
125	character(len=10) :: tmp1
126	character(len=10) :: tmp2
127	character(len=100) :: message
128	character(len=10),parameter :: subname="sw_corrk"
129
130	! For fixed water vapour profiles.
131	integer i_var
132	real RH
133	real psat,qsat
134
135	logical OLRz
136	real*8 NFLUXGNDV_nu(L_NSPECTV)
137
138	! Included by RW for runaway greenhouse 1D study.
139	real vtmp(nlayer)
140	REAL*8 muvarrad(L_LEVELS)
141
142
143	!===============================================================
144	! I.a Initialization on first call
145	!===============================================================
146
147
148	if(firstcall) then
149
150	call iniaerosol()
151
152	allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
153	allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
154	allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
155
156	ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
157	! Effective radius and variance of the aerosols
158	allocate(reffrad(ngrid,nlayer,naerkind))
159	allocate(nueffrad(ngrid,nlayer,naerkind))
160
161	!--------------------------------------------------
162	! Set up correlated k
163	!--------------------------------------------------
164
165	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
166	print*, "corrkdir = ",corrkdir
167	write( tmp1, '(i3)' ) NBinfrared
168	write( tmp2, '(i3)' ) NBvisible
169	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
170	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
171
172	call su_gases ! reading of gases.def bypassed in this, hardcoded. cf su_gases.F90
173	call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
174	call setspv ! Basic visible properties.
175	call sugas_corrk ! Set up gaseous absorption properties.
176	call suaer_corrk ! Set up aerosol optical properties.
177
178	! now that L_NGAUSS has been initialized (by sugas_corrk)
179	! allocate related arrays
180	ALLOCATE(dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
181	ALLOCATE(cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
182	ALLOCATE(wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
183	ALLOCATE(tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS))
184	ALLOCATE(taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS))
185	ALLOCATE(taugsurf(L_NSPECTV,L_NGAUSS-1))
186
187	OSR_nu(:,:) = 0.
188
189	! Surface albedo in the solar range
190	! example of data: Cutler et al 2020
191	! reflectance of basalts in UV and visible varies from a few to 10%
192	! it also depends on mineralogy
193	! in the absence of further data, I take 5% as a mean value... Sensitivity could be tested...
194	albedo(:,:) = 0.05
195
196	end if ! of if (firstcall)
197
198	!=======================================================================
199	! I.b Initialization on every call
200	!=======================================================================
201
202	qxvaer(:,:,:)=0.0
203	qsvaer(:,:,:)=0.0
204	gvaer(:,:,:) =0.0
205
206	! How much light do we get ?
207	do nw=1,L_NSPECTV
208	stel(nw)=stellarf(nw)/(dist_star**2)
209	end do
210
211	! Get 3D aerosol optical properties.
212	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
213	QVISsQREF3d,omegaVIS3d,gVIS3d,QREFvis3d)
214
215	! Get aerosol optical depths.
216	call aeropacity(ngrid,nlayer,pplay,pplev,pt,aerosol, &
217	reffrad,nueffrad,QREFvis3d,tau_col)
218
219	!-----------------------------------------------------------------------
220	!-----------------------------------------------------------------------
221	do ig=1,ngrid ! Starting Big Loop over every GCM column
222	!-----------------------------------------------------------------------
223	!-----------------------------------------------------------------------
224
225	!=======================================================================
226	! II. Transformation of the GCM variables
227	!=======================================================================
228
229	!-----------------------------------------------------------------------
230	! Aerosol optical properties Qext, Qscat and g.
231	! The transformation in the vertical is the same as for temperature.
232	!-----------------------------------------------------------------------
233
234	do iaer=1,naerkind
235	! Shortwave.
236	do nw=1,L_NSPECTV
237
238	do l=1,nlayer
239
240	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
241	*QREFvis3d(ig,nlayer+1-l,iaer)
242
243	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
244	*QREFvis3d(ig,max(nlayer-l,1),iaer)
245
246	qxvaer(2*l,nw,iaer) = temp1
247	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
248
249	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
250	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
251
252	qsvaer(2*l,nw,iaer) = temp1
253	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
254
255	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
256	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
257
258	gvaer(2*l,nw,iaer) = temp1
259	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
260
261	end do ! nlayer
262
263	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
264	qxvaer(2*nlayer+1,nw,iaer)=0.
265
266	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
267	qsvaer(2*nlayer+1,nw,iaer)=0.
268
269	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
270	gvaer(2*nlayer+1,nw,iaer)=0.
271
272	end do ! L_NSPECTV
273	end do ! naerkind
274
275	! Test / Correct for freaky s. s. albedo values.
276	do iaer=1,naerkind
277	do k=1,L_LEVELS
278
279	do nw=1,L_NSPECTV
280	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
281	message='Serious problems with qsvaer values'
282	call abort_physic(subname,message,1)
283	endif
284	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
285	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
286	endif
287	end do
288
289	end do ! L_LEVELS
290	end do ! naerkind
291
292	!-----------------------------------------------------------------------
293	! Aerosol optical depths
294	!-----------------------------------------------------------------------
295
296	do iaer=1,naerkind
297	do k=0,nlayer-1
298
299	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
300	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
301	! As 'aerosol' is at reference (visible) wavelenght we scale it as
302	! it will be multplied by qxi/v in optci/v
303	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
304	tauaero(2k+2,iaer)=max(temppweight,0.d0)
305	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
306
307	end do
308	! boundary conditions
309	tauaero(1,iaer) = tauaero(2,iaer)
310	!tauaero(1,iaer) = 0.
311	!JL18 at time of testing, the two above conditions gave the same results bit for bit.
312
313	end do ! naerkind
314
315	! Albedo
316	DO nw=1,L_NSPECTV ! Short Wave loop.
317	albv(nw)=albedo(ig,nw)
318	ENDDO
319
320	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
321
322	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
323	!!! Note by JL13 : In the following, some indices were changed in the interpolations,
324	!!! so that the model results are less dependent on the number of layers !
325	!!!
326	!!! --- The older versions are commented with the comment !JL13index ---
327	!!!
328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329
330
331	!-----------------------------------------------------------------------
332	! Variable species... Not used for Venus
333	!-----------------------------------------------------------------------
334
335	do k=1,L_LEVELS
336	qvar(k) = 1.0D-7
337	end do
338
339	DO l=1,nlayer
340	muvarrad(2*l) = RMD
341	muvarrad(2*l+1) = RMD
342	END DO
343
344	muvarrad(1) = muvarrad(2)
345	muvarrad(2*nlayer+1)=RMD
346
347	! Keep values inside limits for which we have radiative transfer coefficients !!!
348	if(L_REFVAR.gt.1)then ! (there was a bug here)
349	message='no variable species for Venus yet'
350	call abort_physic(subname,message,1)
351	endif
352
353	!-----------------------------------------------------------------------
354	! Pressure and temperature
355	!-----------------------------------------------------------------------
356
357	DO l=1,nlayer
358	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
359	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
360	tlevrad(2*l) = pt(ig,nlayer+1-l)
361	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
362	END DO
363
364	plevrad(1) = 0.
365	! plevrad(2) = 0. !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all.
366
367	tlevrad(1) = tlevrad(2)
368	tlevrad(2*nlayer+1)=tsurf(ig)
369
370	pmid(1) = pplay(ig,nlayer)/scalep
371	pmid(2) = pmid(1)
372
373	tmid(1) = tlevrad(2)
374	tmid(2) = tmid(1)
375
376	DO l=1,L_NLAYRAD-1
377	tmid(2l+1) = tlevrad(2l+1)
378	tmid(2l+2) = tlevrad(2l+1)
379	pmid(2l+1) = plevrad(2l+1)
380	pmid(2l+2) = plevrad(2l+1)
381	END DO
382	pmid(L_LEVELS) = plevrad(L_LEVELS)
383	tmid(L_LEVELS) = tlevrad(L_LEVELS)
384
385	!!Alternative interpolation:
386	! pmid(3) = pmid(1)
387	! pmid(4) = pmid(1)
388	! tmid(3) = tmid(1)
389	! tmid(4) = tmid(1)
390	! DO l=2,L_NLAYRAD-1
391	! tmid(2l+1) = tlevrad(2l)
392	! tmid(2l+2) = tlevrad(2l)
393	! pmid(2l+1) = plevrad(2l)
394	! pmid(2l+2) = plevrad(2l)
395	! END DO
396	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
397	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
398
399	! Test for out-of-bounds pressure.
400	if(plevrad(3).lt.pgasmin)then
401	! print*,'Minimum pressure is outside the radiative'
402	! print*,'transfer kmatrix bounds, exiting.'
403	! message="Minimum pressure outside of kmatrix bounds"
404	! call abort_physic(subname,message,1)
405	elseif(plevrad(L_LEVELS).gt.pgasmax)then
406	! print*,'Maximum pressure is outside the radiative'
407	! print*,'transfer kmatrix bounds, exiting.'
408	! message="Maximum pressure outside of kmatrix bounds"
409	! call abort_physic(subname,message,1)
410	endif
411
412	! Test for out-of-bounds temperature.
413	do k=1,L_LEVELS
414
415	if(tlevrad(k).lt.tgasmin)then
416	print*,'Minimum temperature is outside the radiative'
417	print*,'transfer kmatrix bounds'
418	print*,"k=",k," tlevrad(k)=",tlevrad(k)
419	print*,"tgasmin=",tgasmin
420	! if (strictboundcorrk) then
421	! message="Minimum temperature outside of kmatrix bounds"
422	! call abort_physic(subname,message,1)
423	! else
424	print,'**********************************************'
425	print*,'we allow model to continue with tlevrad<tgasmin'
426	print*,' ... we assume we know what you are doing ... '
427	print*,' ... but do not let this happen too often ... '
428	print,'**********************************************'
429	!tlevrad(k)=tgasmin ! Used in the source function !
430	! endif
431	elseif(tlevrad(k).gt.tgasmax)then
432	print*,'Maximum temperature is outside the radiative'
433	print*,'transfer kmatrix bounds, exiting.'
434	print*,"k=",k," tlevrad(k)=",tlevrad(k)
435	print*,"tgasmax=",tgasmax
436	! if (strictboundcorrk) then
437	! message="Maximum temperature outside of kmatrix bounds"
438	! call abort_physic(subname,message,1)
439	! else
440	print,'**********************************************'
441	print*,'we allow model to continue with tlevrad>tgasmax'
442	print*,' ... we assume we know what you are doing ... '
443	print*,' ... but do not let this happen too often ... '
444	print,'**********************************************'
445	!tlevrad(k)=tgasmax ! Used in the source function !
446	! endif
447	endif
448
449	enddo
450
451	do k=1,L_NLAYRAD+1
452	if(tmid(k).lt.tgasmin)then
453	print*,'Minimum temperature is outside the radiative'
454	print*,'transfer kmatrix bounds, exiting.'
455	print*,"k=",k," tmid(k)=",tmid(k)
456	print*,"tgasmin=",tgasmin
457	! if (strictboundcorrk) then
458	! message="Minimum temperature outside of kmatrix bounds"
459	! call abort_physic(subname,message,1)
460	! else
461	print,'**********************************************'
462	print*,'we allow model to continue but with tmid=tgasmin'
463	print*,' ... we assume we know what you are doing ... '
464	print*,' ... but do not let this happen too often ... '
465	print,'**********************************************'
466	tmid(k)=tgasmin
467	! endif
468	elseif(tmid(k).gt.tgasmax)then
469	print*,'Maximum temperature is outside the radiative'
470	print*,'transfer kmatrix bounds, exiting.'
471	print*,"k=",k," tmid(k)=",tmid(k)
472	print*,"tgasmax=",tgasmax
473	! if (strictboundcorrk) then
474	! message="Maximum temperature outside of kmatrix bounds"
475	! call abort_physic(subname,message,1)
476	! else
477	print,'**********************************************'
478	print*,'we allow model to continue but with tmid=tgasmax'
479	print*,' ... we assume we know what you are doing ... '
480	print*,' ... but do not let this happen too often ... '
481	print,'**********************************************'
482	tmid(k)=tgasmax
483	! endif
484	endif
485	enddo
486
487	!=======================================================================
488	! III. Calling the main radiative transfer subroutines
489	!=======================================================================
490
491	!-----------------------------------------------------------------------
492	! Short Wave Part
493	!-----------------------------------------------------------------------
494
495	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
496	do nw=1,L_NSPECTV
497	stel_fract(nw)= stel(nw) * fract(ig)
498	end do
499
500	call optcv(dtauv,tauv,taucumv,plevrad, &
501	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
502	tmid,pmid,taugsurf,qvar,muvarrad)
503
504	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
505	acosz,stel_fract, &
506	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
507	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
508
509	else ! During the night, fluxes = 0.
510	nfluxtopv = 0.0d0
511	fluxtopvdn = 0.0d0
512	nfluxoutv_nu(:) = 0.0d0
513	nfluxgndv_nu(:) = 0.0d0
514	do l=1,L_NLAYRAD
515	fmnetv(l)=0.0d0
516	fluxupv(l)=0.0d0
517	fluxdnv(l)=0.0d0
518	end do
519	end if
520
521	!-----------------------------------------------------------------------
522	! Transformation of the correlated-k code outputs
523	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
524
525	! Net incident flux tops, positive downward
526	nfluxtop(ig) = -1.*nfluxtopv
527
528	! Net incident flux at surface, positive downward
529	nfluxsurf(ig) = -1.*fmnetv(L_NLAYRAD)
530
531	! Flux incident at the top of the atmosphere
532	fluxtop_dn(ig)= fluxtopvdn
533
534	fluxabs_sw(ig) = real(-nfluxtopv)
535	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
536
537	if(fluxtop_dn(ig).lt.0.0)then
538	print*,'Achtung! fluxtop_dn has lost the plot!'
539	print*,'fluxtop_dn=',fluxtop_dn(ig)
540	print*,'acosz=',acosz
541	print*,'aerosol=',aerosol(ig,:,:)
542	print*,'temp= ',pt(ig,:)
543	print*,'pplay= ',pplay(ig,:)
544	message="Achtung! fluxtop_dn has lost the plot!"
545	call abort_physic(subname,message,1)
546	endif
547
548	! Net solar flux, positive downward
549
550	do l=1,L_NLAYRAD
551	netflux(ig,L_NLAYRAD+1-l)= -1.*fmnetv(l)
552	enddo
553	netflux(ig,L_NLAYRAD+1) = -1.*nfluxtopv
554
555	! Finally, the heating rates
556
557	DO l=2,L_NLAYRAD
558	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
559	RG/(cpdet(tmid(l))scalep(plevrad(2l+1)-plevrad(2*l-1)))
560	END DO
561
562	! These are values at top of atmosphere
563	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
564	RG/(cpdet(tmid(1))scalep*(plevrad(3)-plevrad(2)))
565
566	!-----------------------------------------------------------------------
567	!-----------------------------------------------------------------------
568	end do ! End of big loop over every GCM column.
569	!-----------------------------------------------------------------------
570	!-----------------------------------------------------------------------
571
572	end subroutine sw_venus_corrk
573

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

Download in other formats:

Original Format