Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sw_venus_corrk.F90

Last change on this file was 3884, checked in by ikovalenko, 4 months ago

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

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

Download in other formats:

Original Format