Context Navigation

sw_venus_corrk.F90 @ 3803

Last change on this file since 3803 was 3794, checked in by slebonnois, 6 months ago
SL: VenusPCM - update of CIA in solar module (adapted from Generic PCM) + potential output of band-resolved net flux
File size: 24.2 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, wavev, wnov
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 nfluxv_nu(L_NLAYRAD,L_NSPECTV) ! vertical profil, band-resolved VI net flux (W/m2)
102	REAL*8 heatrate_nu(L_NLAYRAD,L_NSPECTV) ! vertical profil, band-resolved VI heating rate (K/s/micron)
103	REAL*8 fmnetv(L_NLAYRAD)
104	REAL*8 fluxupv(L_NLAYRAD)
105	REAL*8 fluxdnv(L_NLAYRAD)
106	REAL*8 acosz
107	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
108
109	INTEGER ig,l,k,nw,iaer
110
111	real*8,allocatable,save :: taugsurf(:,:)
112	!$OMP THREADPRIVATE(taugsurf)
113	real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol).
114
115	! Local aerosol optical properties for each column on RADIATIVE grid.
116	real8,save,allocatable :: QXVAER(:,:,:) ! Extinction coeff (QVISsQREFQREFvis)
117	real*8,save,allocatable :: QSVAER(:,:,:)
118	real*8,save,allocatable :: GVAER(:,:,:)
119	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER)
120	real, dimension(:,:,:), save, allocatable :: QREFvis3d
121	!$OMP THREADPRIVATE(QREFvis3d)
122
123
124	! Miscellaneous :
125	real*8 temp,temp1,temp2,pweight
126	character(len=10) :: tmp1
127	character(len=10) :: tmp2
128	character(len=100) :: message
129	character(len=10),parameter :: subname="sw_corrk"
130
131	! For fixed water vapour profiles.
132	integer i_var
133	real RH
134	real psat,qsat
135
136	logical OLRz
137	real*8 NFLUXGNDV_nu(L_NSPECTV)
138
139	! Included by RW for runaway greenhouse 1D study.
140	real vtmp(nlayer)
141	REAL*8 muvarrad(L_LEVELS)
142
143
144	!===============================================================
145	! I.a Initialization on first call
146	!===============================================================
147
148
149	if(firstcall) then
150
151	call iniaerosol()
152
153	allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
154	allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
155	allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
156
157	ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
158	! Effective radius and variance of the aerosols
159	allocate(reffrad(ngrid,nlayer,naerkind))
160	allocate(nueffrad(ngrid,nlayer,naerkind))
161
162	!--------------------------------------------------
163	! Set up correlated k
164	!--------------------------------------------------
165
166	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
167	print*, "corrkdir = ",corrkdir
168	write( tmp1, '(i3)' ) NBinfrared
169	write( tmp2, '(i3)' ) NBvisible
170	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
171	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
172
173	call su_gases ! reading of gases.def bypassed in this, hardcoded. cf su_gases.F90
174	call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
175	call setspv ! Basic visible properties.
176	call sugas_corrk ! Set up gaseous absorption properties.
177	call suaer_corrk ! Set up aerosol optical properties.
178
179	! now that L_NGAUSS has been initialized (by sugas_corrk)
180	! allocate related arrays
181	ALLOCATE(dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
182	ALLOCATE(cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
183	ALLOCATE(wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS))
184	ALLOCATE(tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS))
185	ALLOCATE(taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS))
186	ALLOCATE(taugsurf(L_NSPECTV,L_NGAUSS-1))
187
188	OSR_nu(:,:) = 0.
189
190	! Surface albedo in the solar range
191	! example of data: Cutler et al 2020
192	! reflectance of basalts in UV and visible varies from a few to 10%
193	! it also depends on mineralogy
194	! in the absence of further data, I take 5% as a mean value... Sensitivity could be tested...
195	albedo(:,:) = 0.05
196
197	end if ! of if (firstcall)
198
199	!=======================================================================
200	! I.b Initialization on every call
201	!=======================================================================
202
203	qxvaer(:,:,:)=0.0
204	qsvaer(:,:,:)=0.0
205	gvaer(:,:,:) =0.0
206
207	! How much light do we get ?
208	do nw=1,L_NSPECTV
209	stel(nw)=stellarf(nw)/(dist_star**2)
210	end do
211
212	! Get 3D aerosol optical properties.
213	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
214	QVISsQREF3d,omegaVIS3d,gVIS3d,QREFvis3d)
215
216	! Get aerosol optical depths.
217	call aeropacity(ngrid,nlayer,pplay,pplev,pt,aerosol, &
218	reffrad,nueffrad,QREFvis3d,tau_col)
219
220	!-----------------------------------------------------------------------
221	!-----------------------------------------------------------------------
222	do ig=1,ngrid ! Starting Big Loop over every GCM column
223	!-----------------------------------------------------------------------
224	!-----------------------------------------------------------------------
225
226	!=======================================================================
227	! II. Transformation of the GCM variables
228	!=======================================================================
229
230	!-----------------------------------------------------------------------
231	! Aerosol optical properties Qext, Qscat and g.
232	! The transformation in the vertical is the same as for temperature.
233	!-----------------------------------------------------------------------
234
235	do iaer=1,naerkind
236	! Shortwave.
237	do nw=1,L_NSPECTV
238
239	do l=1,nlayer
240
241	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
242	*QREFvis3d(ig,nlayer+1-l,iaer)
243
244	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
245	*QREFvis3d(ig,max(nlayer-l,1),iaer)
246
247	qxvaer(2*l,nw,iaer) = temp1
248	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
249
250	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
251	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
252
253	qsvaer(2*l,nw,iaer) = temp1
254	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
255
256	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
257	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
258
259	gvaer(2*l,nw,iaer) = temp1
260	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
261
262	end do ! nlayer
263
264	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
265	qxvaer(2*nlayer+1,nw,iaer)=0.
266
267	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
268	qsvaer(2*nlayer+1,nw,iaer)=0.
269
270	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
271	gvaer(2*nlayer+1,nw,iaer)=0.
272
273	end do ! L_NSPECTV
274	end do ! naerkind
275
276	! Test / Correct for freaky s. s. albedo values.
277	do iaer=1,naerkind
278	do k=1,L_LEVELS
279
280	do nw=1,L_NSPECTV
281	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
282	message='Serious problems with qsvaer values'
283	call abort_physic(subname,message,1)
284	endif
285	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
286	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
287	endif
288	end do
289
290	end do ! L_LEVELS
291	end do ! naerkind
292
293	!-----------------------------------------------------------------------
294	! Aerosol optical depths
295	!-----------------------------------------------------------------------
296
297	do iaer=1,naerkind
298	do k=0,nlayer-1
299
300	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
301	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
302	! As 'aerosol' is at reference (visible) wavelenght we scale it as
303	! it will be multplied by qxi/v in optci/v
304	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
305	tauaero(2k+2,iaer)=max(temppweight,0.d0)
306	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
307
308	end do
309	! boundary conditions
310	tauaero(1,iaer) = tauaero(2,iaer)
311	!tauaero(1,iaer) = 0.
312	!JL18 at time of testing, the two above conditions gave the same results bit for bit.
313
314	end do ! naerkind
315
316	! Albedo
317	DO nw=1,L_NSPECTV ! Short Wave loop.
318	albv(nw)=albedo(ig,nw)
319	ENDDO
320
321	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
322
323	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324	!!! Note by JL13 : In the following, some indices were changed in the interpolations,
325	!!! so that the model results are less dependent on the number of layers !
326	!!!
327	!!! --- The older versions are commented with the comment !JL13index ---
328	!!!
329	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
330
331
332	!-----------------------------------------------------------------------
333	! Variable species... Not used for Venus
334	!-----------------------------------------------------------------------
335
336	do k=1,L_LEVELS
337	qvar(k) = 1.0D-7
338	end do
339
340	DO l=1,nlayer
341	muvarrad(2*l) = RMD
342	muvarrad(2*l+1) = RMD
343	END DO
344
345	muvarrad(1) = muvarrad(2)
346	muvarrad(2*nlayer+1)=RMD
347
348	! Keep values inside limits for which we have radiative transfer coefficients !!!
349	if(L_REFVAR.gt.1)then ! (there was a bug here)
350	message='no variable species for Venus yet'
351	call abort_physic(subname,message,1)
352	endif
353
354	!-----------------------------------------------------------------------
355	! Pressure and temperature
356	!-----------------------------------------------------------------------
357
358	DO l=1,nlayer
359	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
360	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
361	tlevrad(2*l) = pt(ig,nlayer+1-l)
362	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
363	END DO
364
365	plevrad(1) = 0.
366	! 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.
367
368	tlevrad(1) = tlevrad(2)
369	tlevrad(2*nlayer+1)=tsurf(ig)
370
371	pmid(1) = pplay(ig,nlayer)/scalep
372	pmid(2) = pmid(1)
373
374	tmid(1) = tlevrad(2)
375	tmid(2) = tmid(1)
376
377	DO l=1,L_NLAYRAD-1
378	tmid(2l+1) = tlevrad(2l+1)
379	tmid(2l+2) = tlevrad(2l+1)
380	pmid(2l+1) = plevrad(2l+1)
381	pmid(2l+2) = plevrad(2l+1)
382	END DO
383	pmid(L_LEVELS) = plevrad(L_LEVELS)
384	tmid(L_LEVELS) = tlevrad(L_LEVELS)
385
386	!!Alternative interpolation:
387	! pmid(3) = pmid(1)
388	! pmid(4) = pmid(1)
389	! tmid(3) = tmid(1)
390	! tmid(4) = tmid(1)
391	! DO l=2,L_NLAYRAD-1
392	! tmid(2l+1) = tlevrad(2l)
393	! tmid(2l+2) = tlevrad(2l)
394	! pmid(2l+1) = plevrad(2l)
395	! pmid(2l+2) = plevrad(2l)
396	! END DO
397	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
398	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
399
400	! Test for out-of-bounds pressure.
401	if(plevrad(3).lt.pgasmin)then
402	! print*,'Minimum pressure is outside the radiative'
403	! print*,'transfer kmatrix bounds, exiting.'
404	! message="Minimum pressure outside of kmatrix bounds"
405	! call abort_physic(subname,message,1)
406	elseif(plevrad(L_LEVELS).gt.pgasmax)then
407	! print*,'Maximum pressure is outside the radiative'
408	! print*,'transfer kmatrix bounds, exiting.'
409	! message="Maximum pressure outside of kmatrix bounds"
410	! call abort_physic(subname,message,1)
411	endif
412
413	! Test for out-of-bounds temperature.
414	do k=1,L_LEVELS
415
416	if(tlevrad(k).lt.tgasmin)then
417	print*,'Minimum temperature is outside the radiative'
418	print*,'transfer kmatrix bounds'
419	print*,"k=",k," tlevrad(k)=",tlevrad(k)
420	print*,"tgasmin=",tgasmin
421	! if (strictboundcorrk) then
422	! message="Minimum temperature outside of kmatrix bounds"
423	! call abort_physic(subname,message,1)
424	! else
425	print,'**********************************************'
426	print*,'we allow model to continue with tlevrad<tgasmin'
427	print*,' ... we assume we know what you are doing ... '
428	print*,' ... but do not let this happen too often ... '
429	print,'**********************************************'
430	!tlevrad(k)=tgasmin ! Used in the source function !
431	! endif
432	elseif(tlevrad(k).gt.tgasmax)then
433	print*,'Maximum temperature is outside the radiative'
434	print*,'transfer kmatrix bounds, exiting.'
435	print*,"k=",k," tlevrad(k)=",tlevrad(k)
436	print*,"tgasmax=",tgasmax
437	! if (strictboundcorrk) then
438	! message="Maximum temperature outside of kmatrix bounds"
439	! call abort_physic(subname,message,1)
440	! else
441	print,'**********************************************'
442	print*,'we allow model to continue with tlevrad>tgasmax'
443	print*,' ... we assume we know what you are doing ... '
444	print*,' ... but do not let this happen too often ... '
445	print,'**********************************************'
446	!tlevrad(k)=tgasmax ! Used in the source function !
447	! endif
448	endif
449
450	enddo
451
452	do k=1,L_NLAYRAD+1
453	if(tmid(k).lt.tgasmin)then
454	print*,'Minimum temperature is outside the radiative'
455	print*,'transfer kmatrix bounds, exiting.'
456	print*,"k=",k," tmid(k)=",tmid(k)
457	print*,"tgasmin=",tgasmin
458	! if (strictboundcorrk) then
459	! message="Minimum temperature outside of kmatrix bounds"
460	! call abort_physic(subname,message,1)
461	! else
462	print,'**********************************************'
463	print*,'we allow model to continue but with tmid=tgasmin'
464	print*,' ... we assume we know what you are doing ... '
465	print*,' ... but do not let this happen too often ... '
466	print,'**********************************************'
467	tmid(k)=tgasmin
468	! endif
469	elseif(tmid(k).gt.tgasmax)then
470	print*,'Maximum temperature is outside the radiative'
471	print*,'transfer kmatrix bounds, exiting.'
472	print*,"k=",k," tmid(k)=",tmid(k)
473	print*,"tgasmax=",tgasmax
474	! if (strictboundcorrk) then
475	! message="Maximum temperature outside of kmatrix bounds"
476	! call abort_physic(subname,message,1)
477	! else
478	print,'**********************************************'
479	print*,'we allow model to continue but with tmid=tgasmax'
480	print*,' ... we assume we know what you are doing ... '
481	print*,' ... but do not let this happen too often ... '
482	print,'**********************************************'
483	tmid(k)=tgasmax
484	! endif
485	endif
486	enddo
487
488	!=======================================================================
489	! III. Calling the main radiative transfer subroutines
490	!=======================================================================
491
492	!-----------------------------------------------------------------------
493	! Short Wave Part
494	!-----------------------------------------------------------------------
495
496	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
497	do nw=1,L_NSPECTV
498	stel_fract(nw)= stel(nw) * fract(ig)
499	end do
500
501	call optcv(dtauv,tauv,taucumv,plevrad, &
502	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
503	tmid,pmid,taugsurf,qvar,muvarrad)
504
505	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
506	acosz,stel_fract, &
507	nfluxtopv,fluxtopvdn,nfluxv_nu,nfluxgndv_nu, &
508	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
509
510	else ! During the night, fluxes = 0.
511	nfluxtopv = 0.0d0
512	fluxtopvdn = 0.0d0
513	nfluxv_nu(:,:) = 0.0d0
514	nfluxgndv_nu(:) = 0.0d0
515	do l=1,L_NLAYRAD
516	fmnetv(l)=0.0d0
517	fluxupv(l)=0.0d0
518	fluxdnv(l)=0.0d0
519	end do
520	end if
521
522	!-----------------------------------------------------------------------
523	! Transformation of the correlated-k code outputs
524	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
525
526	! Net incident flux tops, positive downward
527	nfluxtop(ig) = -1.*nfluxtopv
528
529	! Net incident flux at surface, positive downward
530	nfluxsurf(ig) = -1.*fmnetv(L_NLAYRAD)
531
532	! Flux incident at the top of the atmosphere
533	fluxtop_dn(ig)= fluxtopvdn
534
535	fluxabs_sw(ig) = real(-nfluxtopv)
536	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
537
538	if(fluxtop_dn(ig).lt.0.0)then
539	print*,'Achtung! fluxtop_dn has lost the plot!'
540	print*,'fluxtop_dn=',fluxtop_dn(ig)
541	print*,'acosz=',acosz
542	print*,'aerosol=',aerosol(ig,:,:)
543	print*,'temp= ',pt(ig,:)
544	print*,'pplay= ',pplay(ig,:)
545	message="Achtung! fluxtop_dn has lost the plot!"
546	call abort_physic(subname,message,1)
547	endif
548
549	! Net solar flux, positive downward
550
551	do l=1,L_NLAYRAD
552	netflux(ig,L_NLAYRAD+1-l)= -1.*fmnetv(l)
553	enddo
554	netflux(ig,L_NLAYRAD+1) = -1.*nfluxtopv
555
556	! Finally, the heating rates
557
558	DO l=2,L_NLAYRAD
559	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
560	RG/(cpdet(tmid(l))scalep(plevrad(2l+1)-plevrad(2*l-1)))
561	END DO
562
563	! These are values at top of atmosphere
564	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
565	RG/(cpdet(tmid(1))scalep*(plevrad(3)-plevrad(2)))
566
567
568	!-----------------------------------------------------------------------
569	!-----------------------------------------------------------------------
570	end do ! End of big loop over every GCM column.
571	!-----------------------------------------------------------------------
572	!-----------------------------------------------------------------------
573
574	! Diagnostic in 1D:
575	! output the vertical profil of band-resolved heating rates, in K/s/microns
576	if ((1.eq.1).and.is_master.and.(ngrid.eq.1).and.firstcall) then
577	open(unit=11,file="nfluxnu.dat")
578	write(11,*) "lambda(microns)",wavev(:)
579	write(11,) "dlbda(microns)",1.e4dwnv(:)/wnov(:)**2
580	write(11,*) "pressure(Pa) dT_nu(K/s/micron)"
581	DO l=2,L_NLAYRAD
582	heatrate_nu(L_NLAYRAD+1-l,:)=(nfluxv_nu(l,:)-nfluxv_nu(l-1,:)) &
583	RG/(cpdet(tmid(l))scalep(plevrad(2l+1)-plevrad(2*l-1))) &
584	! dlbda(microns)=-dnu(cm-1)*1E4/nu(cm-1)^2
585	/ (1.e4dwnv(:)/wnov(:)*2)
586	write(11,*) pplay(1,L_NLAYRAD+1-l),heatrate_nu(L_NLAYRAD+1-l,:)
587	END DO
588	close(11)
589	endif
590
591	end subroutine sw_venus_corrk
592

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

Download in other formats:

Original Format