Context Navigation

callcorrk.F @ 3700

Last change on this file since 3700 was 3373, checked in by afalco, 18 months ago
Pluto.old PCM: Corrected some runtime issues with gfortran. AF
File size: 31.6 KB

Line
1	subroutine callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
2	& albedo,emis,mu0,pplev,pplay,pt,
3	& tsurf,fract,dist_star,igout,aerosol,cpp3D,
4	& dtlw,dtsw,fluxsurf_lw,
5	& fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn,
6	& reffrad,tau_col,ptime,pday,firstcall,lastcall,zzlay)
7
8	use radinc_h
9	use radcommon_h
10	use ioipsl_getincom
11	use radii_mod
12	use aerosol_mod
13	use datafile_mod, only: datadir
14
15	implicit none
16
17	!==================================================================
18	!
19	! Purpose
20	! -------
21	! Solve the radiative transfer using the correlated-k method for
22	! the gaseous absorption and the Toon et al. (1989) method for
23	! scatttering due to aerosols.
24	!
25	! Authors
26	! -------
27	! Emmanuel 01/2001, Forget 09/2001
28	! Robin Wordsworth (2009)
29	! Modif Pluton Tanguy Bertrand 2017
30	!==================================================================
31
32	#include "dimphys.h"
33	#include "comcstfi.h"
34	#include "callkeys.h"
35	#include "tracer.h"
36
37	!-----------------------------------------------------------------------
38	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
39	! Layer #1 is the layer near the ground.
40	! Layer #nlayermx is the layer at the top.
41
42	! INPUT
43	INTEGER icount
44	INTEGER ngrid,nlayer
45	INTEGER igout
46	REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau
47	REAL albedo(ngrid) ! SW albedo
48	REAL emis(ngrid) ! LW emissivity
49	REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer
50	REAL zzlay(ngrid,nlayermx) ! altitude at the middle of the layers
51	REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries
52
53	REAL pt(ngrid,nlayermx) ! air temperature (K)
54	REAL tsurf(ngrid) ! surface temperature (K)
55	REAL dist_star,mu0(ngrid) ! distance star-planet (AU)
56	REAL fract(ngrid) ! fraction of day
57
58	! Globally varying aerosol optical properties on GCM grid
59	! Not needed everywhere so not in radcommon_h
60	REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
61	REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
62	REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
63
64	REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
65	REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
66	REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
67
68	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
69	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
70
71	! REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
72	! REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
73
74	REAL reffrad(ngrid,nlayer,naerkind)
75	REAL nueffrad(ngrid,nlayer,naerkind)
76	REAL profhaze(ngrid,nlayer) ! TB17 fixed profile of haze mmr
77
78	! OUTPUT
79	REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
80	REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
81	REAL fluxsurf_lw(ngridmx) ! incident LW flux to surf (W/m2)
82	REAL fluxtop_lw(ngridmx) ! outgoing LW flux to space (W/m2)
83	REAL fluxsurf_sw(ngridmx) ! incident SW flux to surf (W/m2)
84	REAL fluxtop_sw(ngridmx) ! outgoing LW flux to space (W/m2)
85	REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2)
86
87	!-----------------------------------------------------------------------
88	! Declaration of the variables required by correlated-k subroutines
89	! Numbered from top to bottom unlike in the GCM!
90
91	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
92	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
93
94	! Optical values for the optci/cv subroutines
95	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
96	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
97	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
98	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
99	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
100	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
101	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
102	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
103	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
104	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
105
106	REAL*8 tauaero(L_LEVELS+1,naerkind)
107	REAL*8 nfluxtopv,nfluxtopi,nfluxtop
108	real*8 NFLUXTOPV_nu(L_NSPECTV)
109	real*8 NFLUXTOPI_nu(L_NSPECTI)
110	real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
111	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
112	real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
113	real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
114	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
115	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
116	REAL*8 albi,albv,acosz
117
118	INTEGER ig,l,k,nw,iaer,irad
119
120	real fluxtoplanet
121	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
122	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
123
124	real*8 qvar(L_LEVELS) ! mixing ratio of variable component
125	REAL pq(ngridmx,nlayer,nq)
126	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
127	integer nq
128
129	! Local aerosol optical properties for each column on RADIATIVE grid
130	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
131	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
132	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
133	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
134	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
135	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
136
137	save qxvaer, qsvaer, gvaer
138	save qxiaer, qsiaer, giaer
139	save QREFvis3d, QREFir3d
140
141	REAL tau_col(ngrid) ! diagnostic from aeropacity
142
143	! Misc.
144	logical firstcall, lastcall, nantest
145	real*8 tempv(L_NSPECTV)
146	real*8 tempi(L_NSPECTI)
147	real*8 temp,temp1,temp2,pweight
148	character(len=10) :: tmp1
149	character(len=10) :: tmp2
150
151	! for fixed vapour profiles
152	real RH
153	real*8 pq_temp(nlayer)
154	real ptemp, Ttemp, qsat
155
156	! for OLR spec
157	integer OLRcount
158	save OLRcount
159	integer OLRcount2
160	save OLRcount2
161	character(len=2) :: tempOLR
162	character(len=30) :: filenomOLR
163	real ptime, pday
164
165	REAL epsi_ch4
166	SAVE epsi_ch4
167
168	logical diagrad_OLRz,diagrad_OLR,diagrad_surf,diagrad_rates
169	real OLR_nu(ngrid,L_NSPECTI)
170
171	! Allow variations in cp with location
172	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
173
174	! NLTE factor for CH4
175	real eps_nlte_sw23(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
176	real eps_nlte_sw33(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
177	real eps_nlte_lw(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
178	REAL dtlw_nu(nlayermx,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands
179	REAL dtsw_nu(nlayermx,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands
180
181	REAL dpp ! intermediate
182
183	REAL dtlw_co(ngridmx, nlayermx) ! cooling rate (K/s) due to CO (diagnostic)
184	REAL dtlw_hcn_c2h2(ngridmx, nlayermx) ! cooling rate (K/s) due to C2H2/HCN (diagnostic)
185
186	!!read altitudes and vmrch4
187	integer Nfine,ifine
188	parameter(Nfine=701)
189	character(len=100) :: file_path
190	real,save :: levdat(Nfine),vmrdat(Nfine)
191	real :: vmrch4(ngridmx,nlayermx) ! vmr ch4 from vmrch4_proffix
192
193	!=======================================================================
194	! Initialization on first call
195
196	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,qxvaer)
197	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,qsvaer)
198	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,gvaer)
199
200	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,qxiaer)
201	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,qsiaer)
202	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,giaer)
203
204
205	if(firstcall) then
206
207	print*, "callcorrk: Correlated-k data folder:",trim(datadir)
208	call getin("corrkdir",corrkdir)
209	print*, "corrkdir = ",corrkdir
210	write( tmp1, '(i3)' ) L_NSPECTI
211	write( tmp2, '(i3)' ) L_NSPECTV
212	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
213	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
214
215	print*,'starting sugas'
216	call sugas_corrk ! set up gaseous absorption properties
217	print*,'starting setspi'
218	call setspi ! basic infrared properties
219	print*,'starting setspv'
220	call setspv ! basic visible properties
221
222	! Radiative Hazes
223	if (aerohaze) then
224
225	print*,'aerohaze: starting suaer_corrk'
226	call suaer_corrk ! set up aerosol optical properties
227	print*,'ending suaer_corrk'
228
229	!--------------------------------------------------
230	! Effective radius and variance of the aerosols
231	!--------------------------------------------------
232	do iaer=1,naerkind
233	if ((iaer.eq.iaero_haze)) then
234	call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer),
235	& nueffrad(1,1,iaer))
236	endif
237	end do !iaer=1,naerkind.
238	if (haze_radproffix) then
239	print*, 'haze_radproffix=T : fixed profile for haze rad'
240	else
241	print*,'reffrad haze:',reffrad(1,1,iaero_haze)
242	print*,'nueff haze',nueffrad(1,1,iaero_haze)
243	endif
244	endif ! radiative haze
245
246	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
247
248	epsi_ch4=mmol(igcm_ch4_gas)/mmol(igcm_n2)
249
250	! If fixed profile of CH4 gas
251	IF (vmrch4_proffix) then
252	file_path=trim(datadir)//'/gas_prop/vmr_ch4.txt'
253	open(115,file=file_path,form='formatted')
254	do ifine=1,Nfine
255	read(115,*) levdat(ifine), vmrdat(ifine)
256	enddo
257	close(115)
258	ENDIF
259
260	end if ! firstcall
261
262	!=======================================================================
263
264	! L_NSPECTV is the number of Visual(or Solar) spectral intervals
265	! how much light we get
266	fluxtoplanet=0
267	DO nw=1,L_NSPECTV
268	stel(nw)=stellarf(nw)/(dist_star**2) !flux
269	fluxtoplanet=fluxtoplanet + stel(nw)
270	END DO
271
272	!-----------------------------------------------------------------------
273	! Get 3D aerosol optical properties.
274	! ici on selectionne les proprietes opt correspondant a reffrad
275	if (aerohaze) then
276	!--------------------------------------------------
277	! Effective radius and variance of the aerosols if profil non
278	! uniform. Called only if aerohaze=true.
279	!--------------------------------------------------
280	if (haze_radproffix) then
281	call haze_reffrad_fix(ngrid,nlayer,zzlay,
282	& reffrad,nueffrad)
283	endif
284
285	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
286	& QVISsQREF3d,omegaVIS3d,gVIS3d,
287	& QIRsQREF3d,omegaIR3d,gIR3d,
288	& QREFvis3d,QREFir3d)
289
290	! Get aerosol optical depths.
291	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,
292	& reffrad,QREFvis3d,QREFir3d,
293	& tau_col)
294	endif
295
296	!-----------------------------------------------------------------------
297	! Prepare CH4 mixing ratio for radiative transfer
298	IF (methane) then
299	vmrch4(:,:)=0.
300
301	if (ch4fix) then
302	if (vmrch4_proffix) then
303	!! Interpolate on the model vertical grid
304	do ig=1,ngridmx
305	! CALL interp_line(levdat,vmrdat,Nfine,
306	! & zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
307	enddo
308	else
309	vmrch4(:,:)=vmrch4fix
310	endif
311	else
312	vmrch4(:,:)=pq(:,:,igcm_ch4_gas)100.
313	& mmol(igcm_n2)/mmol(igcm_ch4_gas)
314	endif
315	ENDIF
316
317	! Prepare NON LTE correction in Pluto atmosphere
318	IF (nlte) then
319	CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
320	& eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
321	ENDIF
322	c Net atmospheric radiative cooling rate from C2H2 (K.s-1):
323	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
324	dtlw_hcn_c2h2=0.
325	if (cooling) then
326	CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,
327	& pt,dtlw_hcn_c2h2)
328	endif
329
330	!-----------------------------------------------------------------------
331	!=======================================================================
332	!-----------------------------------------------------------------------
333	! starting big loop over every GCM column
334	do ig=1,ngridmx
335
336	!=======================================================================
337	! Transformation of the GCM variables
338
339	!-----------------------------------------------------------------------
340	! Aerosol optical properties Qext, Qscat and g on each band
341	! The transformation in the vertical is the same as for temperature
342	if (aerohaze) then
343	! shortwave
344	do iaer=1,naerkind
345	DO nw=1,L_NSPECTV
346	do l=1,nlayermx
347
348	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)
349	$ *QREFvis3d(ig,nlayermx+1-l,iaer)
350
351	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
352	$ *QREFvis3d(ig,max(nlayermx-l,1),iaer)
353	qxvaer(2*l,nw,iaer) = temp1
354	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
355
356	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
357	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
358
359	qsvaer(2*l,nw,iaer) = temp1
360	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
361
362	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
363	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
364
365	gvaer(2*l,nw,iaer) = temp1
366	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
367
368	end do
369	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
370	qxvaer(2*nlayermx+1,nw,iaer)=0.
371
372	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
373	qsvaer(2*nlayermx+1,nw,iaer)=0.
374
375	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
376	gvaer(2*nlayermx+1,nw,iaer)=0.
377	end do
378
379	! longwave
380	DO nw=1,L_NSPECTI
381	do l=1,nlayermx
382
383	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)
384	$ *QREFir3d(ig,nlayermx+1-l,iaer)
385
386	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
387	$ *QREFir3d(ig,max(nlayermx-l,1),iaer)
388
389	qxiaer(2*l,nw,iaer) = temp1
390	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
391
392	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
393	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
394
395	qsiaer(2*l,nw,iaer) = temp1
396	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
397
398	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
399	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
400
401	giaer(2*l,nw,iaer) = temp1
402	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
403
404	end do
405
406	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
407	qxiaer(2*nlayermx+1,nw,iaer)=0.
408
409	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
410	qsiaer(2*nlayermx+1,nw,iaer)=0.
411
412	giaer(1,nw,iaer)=giaer(2,nw,iaer)
413	giaer(2*nlayermx+1,nw,iaer)=0.
414	end do
415	end do ! naerkind
416
417	! Test / Correct for freaky s. s. albedo values.
418	do iaer=1,naerkind
419	do k=1,L_LEVELS
420
421	do nw=1,L_NSPECTV
422	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
423	print*,'Serious problems with qsvaer values'
424	print*,'in callcorrk'
425	call abort
426	endif
427	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
428	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
429	endif
430	end do
431
432	do nw=1,L_NSPECTI
433	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
434	print*,'Serious problems with qsiaer values'
435	print*,'in callcorrk'
436	call abort
437	endif
438	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
439	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
440	endif
441	end do
442	end do ! levels
443
444	end do ! naerkind
445
446	endif ! aerohaze
447
448	!-----------------------------------------------------------------------
449	! Aerosol optical depths
450	IF (aerohaze) THEN
451	do iaer=1,naerkind ! heritage generic
452	do k=0,nlayer-1
453	pweight=
454	$ (pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/
455	$ (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
456	if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then
457	temp=aerosol(ig,L_NLAYRAD-k,iaer)/
458	$ QREFvis3d(ig,L_NLAYRAD-k,iaer)
459	else
460	print*, 'stop corrk',k,QREFvis3d(ig,L_NLAYRAD-k,iaer)
461	stop
462	end if
463	tauaero(2k+2,iaer)=max(temppweight,0.d0)
464	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0) ! tauaero en L_LEVELS soit deux fois plus que nlayer
465	end do
466
467	! generic New boundary conditions
468	tauaero(1,iaer) = tauaero(2,iaer)
469	!tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
470	!tauaero(1,iaer) = 0.
471	!tauaero(L_LEVELS+1,iaer) = 0.
472
473	end do ! naerkind
474	ELSE
475	tauaero(:,:)=0
476	ENDIF
477	!-----------------------------------------------------------------------
478
479	! Albedo and emissivity
480	albi=1-emis(ig) ! longwave
481	albv=albedo(ig) ! shortwave
482	acosz=mu0(ig) ! cosine of sun incident angle
483
484	!-----------------------------------------------------------------------
485	! Methane vapour
486
487	c qvar = mixing ratio
488	c L_LEVELS (51) different de GCM levels (25) . L_LEVELS = 2(llm-1)+3=2(ngridmx-1)+3
489	c L_REFVAR The number of different mixing ratio values in
490	c datagcm/composition.in for the k-coefficients.
491	qvar(:)=0.
492	IF (methane) then
493
494	do l=1,nlayer
495	qvar(2l) = vmrch4(ig,nlayer+1-l)/100.
496	& mmol(igcm_ch4_gas)/mmol(igcm_n2)
497	qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig,
498	& max(nlayer-l,1)))/2.)/100.*
499	& mmol(igcm_ch4_gas)/mmol(igcm_n2)
500	end do
501	qvar(1)=qvar(2)
502
503
504	! Keep values inside limits for which we have radiative transfer coefficients
505	if(L_REFVAR.gt.1)then ! there was a bug here!
506	do k=1,L_LEVELS
507	if(qvar(k).lt.wrefvar(1))then
508	qvar(k)=wrefvar(1)+1.0e-8
509	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
510	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
511	endif
512	end do
513	endif
514
515	! IMPORTANT: Now convert from kg/kg to mol/mol
516	do k=1,L_LEVELS
517	qvar(k) = qvar(k)/(epsi_ch4+qvar(k)*(1.-epsi_ch4))
518	end do
519	ENDIF ! methane
520
521	!-----------------------------------------------------------------------
522	! Pressure and temperature
523
524	! generic updated:
525	DO l=1,nlayer
526	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
527	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
528	tlevrad(2*l) = pt(ig,nlayer+1-l)
529	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
530	END DO
531
532	plevrad(1) = 0.
533	plevrad(2) = 0. !! Trick to have correct calculations of fluxes
534	! in gflux(i/v).F, but the pmid levels are not impacted by this change.
535
536	tlevrad(1) = tlevrad(2)
537	tlevrad(2*nlayer+1)=tsurf(ig)
538
539	pmid(1) = max(pgasmin,0.0001*plevrad(3))
540	pmid(2) = pmid(1)
541
542	tmid(1) = tlevrad(2)
543	tmid(2) = tmid(1)
544
545	! INI
546	! DO l=1,nlayer
547	! plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
548	! plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
549	! tlevrad(2*l) = pt(ig,nlayer+1-l)
550	! tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,
551	! $ max(nlayer-l,1)))/2
552	! END DO
553
554	!! following lines changed in 03/2015 to solve upper atmosphere bug
555	! plevrad(1) = 0.
556	! plevrad(2) = max(pgasmin,0.0001*plevrad(3))
557	!
558	! tlevrad(1) = tlevrad(2)
559	! tlevrad(2*nlayermx+1)=tsurf(ig)
560	!
561	! tmid(1) = tlevrad(2)
562	! tmid(2) = tlevrad(2)
563	!
564	! pmid(1) = plevrad(2)
565	! pmid(2) = plevrad(2)
566
567	DO l=1,L_NLAYRAD-1
568	tmid(2l+1) = tlevrad(2l+1)
569	tmid(2l+2) = tlevrad(2l+1)
570	pmid(2l+1) = plevrad(2l+1)
571	pmid(2l+2) = plevrad(2l+1)
572	END DO
573	! end of changes
574	pmid(L_LEVELS) = plevrad(L_LEVELS)
575	tmid(L_LEVELS) = tlevrad(L_LEVELS)
576
577	!TB
578	if ((PMID(2).le.1.e-5).and.(ig.eq.1)) then
579	if ((TMID(2).le.30.).and.(ig.eq.1)) then
580	write(,) 'Caution! Pres/temp of upper levels lower than'
581	write(,) 'ref pressure/temp: kcoef fixed for upper levels'
582	!! cf tpindex.F
583	endif
584	endif
585
586	! test for out-of-bounds pressure
587	if(plevrad(3).lt.pgasmin)then
588	print*,'Minimum pressure is outside the radiative'
589	print*,'transfer kmatrix bounds, exiting.'
590	! call abort
591	elseif(plevrad(L_LEVELS).gt.pgasmax)then
592	print*,'Maximum pressure is outside the radiative'
593	print*,'transfer kmatrix bounds, exiting.'
594	! call abort
595	endif
596
597	! test for out-of-bounds temperature
598	do k=1,L_LEVELS
599	if(tlevrad(k).lt.tgasmin)then
600	print*,'Minimum temperature is outside the radiative'
601	print*,'transfer kmatrix bounds, exiting.'
602	print*,tlevrad(k),k,tgasmin
603	! call abort
604	elseif(tlevrad(k).gt.tgasmax)then
605	print*,'Maximum temperature is outside the radiative'
606	print*,'transfer kmatrix bounds, exiting.'
607	! call abort
608	endif
609	enddo
610
611	!=======================================================================
612	! Calling the main radiative transfer subroutines
613
614	!-----------------------------------------------------------------------
615	! Shortwave
616
617	IF(fract(ig) .GE. 1.0e-4) THEN ! only during daylight IPM?! flux UV...
618
619	fluxtoplanet=0.
620	DO nw=1,L_NSPECTV
621	stel_fract(nw)= stel(nw) * fract(ig)
622	fluxtoplanet=fluxtoplanet + stel_fract(nw)
623	END DO
624
625	!print*, 'starting optcv'
626	call optcv(dtauv,tauv,taucumv,plevrad,
627	$ qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,
628	$ tmid,pmid,taugsurf,qvar)
629
630	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,
631	$ acosz,stel_fract,gweight,nfluxtopv,nfluxtopv_nu,
632	$ fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf)
633
634	ELSE ! during the night, fluxes = 0
635	nfluxtopv=0.0
636	DO l=1,L_NLAYRAD
637	fmnetv(l)=0.0
638	fluxupv(l)=0.0
639	fluxdnv(l)=0.0
640	END DO
641	END IF
642
643	!-----------------------------------------------------------------------
644	! Longwave
645
646	call optci(plevrad,tlevrad,dtaui,taucumi,
647	$ qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,
648	$ taugsurfi,qvar)
649	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,
650	$ wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu,
651	$ fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
652
653
654	!-----------------------------------------------------------------------
655	! Transformation of the correlated-k code outputs
656	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
657
658	fluxtop_lw(ig) = fluxupi(1)
659	fluxsurf_lw(ig) = fluxdni(L_NLAYRAD)
660	fluxtop_sw(ig) = fluxupv(1)
661	fluxsurf_sw(ig) = fluxdnv(L_NLAYRAD)
662
663	! Flux incident at the top of the atmosphere
664	fluxtop_dn(ig)=fluxdnv(1)
665
666	! IR spectral output from top of the atmosphere
667	if(specOLR)then
668	do nw=1,L_NSPECTI
669	OLR_nu(ig,nw)=nfluxtopi_nu(nw)
670	end do
671	endif
672
673	! **********************************************************
674	! Finally, the heating rates
675	! g/cp*DF/DP
676	! **********************************************************
677
678	DO l=2,L_NLAYRAD
679	dpp = g/(cppscalep(plevrad(2l+1)-plevrad(2l-1)))
680
681	! DTSW :
682	!dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))*dpp !averaged dtlw on each wavelength
683	do nw=1,L_NSPECTV
684	dtsw_nu(L_NLAYRAD+1-l,nw)=
685	& (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
686	end do
687
688	! DTLW : detailed with wavelength so that we can apply NLTE
689	!dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))*dpp !averaged dtlw on each wavelength
690	do nw=1,L_NSPECTI
691	dtlw_nu(L_NLAYRAD+1-l,nw)=
692	& (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
693	end do
694	END DO
695
696	! values at top of atmosphere
697	dpp = g/(cppscalep(plevrad(3)-plevrad(1)))
698
699	! SW
700	!dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)*dpp
701	do nw=1,L_NSPECTV
702	dtsw_nu(L_NLAYRAD,nw)=
703	& (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp
704	end do
705
706	! LW
707	c dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp
708	do nw=1,L_NSPECTI
709	dtlw_nu(L_NLAYRAD,nw)=
710	& (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp
711	end do
712
713	! **********************************************************
714	! NON NLTE correction in Pluto atmosphere
715	! And conversion of LW spectral heating rates to total rates
716	! **********************************************************
717
718	if (.not.nlte) then
719	eps_nlte_sw23(ig,:) =1. ! IF no NLTE
720	eps_nlte_sw33(ig,:) =1. ! IF no NLTE
721	eps_nlte_lw(ig,:) =1. ! IF no NLTE
722	endif
723
724	do l=1,nlayer
725
726	!LW
727	dtlw(ig,l) =0.
728	! dtlw_co(ig,l) =0. ! only for diagnostic
729	do nw=1,L_NSPECTI
730	! wewei : wavelength in micrometer
731	if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
732	dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
733	else
734	!dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
735	dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
736	! dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
737	end if
738	dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
739	end do
740	! adding c2h2 if cooling active
741	dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l)
742
743	!SW
744	dtsw(ig,l) =0.
745
746	if (strobel) then
747
748	do nw=1,L_NSPECTV
749	if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
750	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
751	elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
752	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
753	else
754	dtsw_nu(l,nw)=dtsw_nu(l,nw)
755	end if
756	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
757	end do
758
759	else ! total heating rates multiplied by nlte coef
760
761	do nw=1,L_NSPECTV
762	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
763	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
764	enddo
765
766	endif
767
768
769	end do
770	! **********************************************************
771
772	! Diagnotics for last call for each grid point
773	!if (lastcall) then
774
775	!print*,'albedi vis=',albv
776	!print*,'albedo ir=',albi
777	!print*,'fluxup ir (:)=',fluxupi
778	!print*,'flux ir net (:)=',fluxdni-fluxupi
779	!print*,'cumulative flux net ir (:)=',fmneti
780	!print*,'cumulative flux net vis (:)=',fmnetv
781	!print*,'fluxdn vis (:)=',fluxdnv
782	!print*,'fluxtop vis=',fluxtop_sw
783	!print*,'fluxsurf vis=',fluxsurf_sw
784	!print*,'fluxtop ir=',fluxtop_lw
785	!print*,'fluxsurf ir=',fluxsurf_lw
786
787	c write(,) 'pressure, eps_nlte_sw, eps_nlte_lw'
788	c do l=1,nlayer
789	c write(,)pplay(1,l),eps_nlte_sw(1,l),eps_nlte_lw(1,l)
790	c end do
791
792	!endif
793
794	! ---------------------------------------------------------------
795	end do ! end of big loop over every GCM column (ig = 1:ngrid)
796
797	!-----------------------------------------------------------------------
798	! Additional diagnostics
799
800	! IR spectral output, for exoplanet observational comparison
801	if(specOLR)then
802	if(ngrid.ne.1)then
803	call writediagspec(ngrid,"OLR3D",
804	& "OLR(lon,lat,band)","W m^-2",3,OLR_nu)
805	endif
806	endif
807
808	if(lastcall)then
809
810	! 1D Output
811	if(ngrid.eq.1)then
812
813	! surface diagnotics
814	diagrad_surf=.true.
815	if(diagrad_surf)then
816	open(116,file='surf_vals.out')
817	write(116,*) tsurf(1),pplev(1,1),
818	& fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1)
819	do nw=1,L_NSPECTV
820	write(116,*) wavev(nw),fmnetv_nu(L_NLAYRAD,nw)
821	enddo
822	do nw=1,L_NSPECTI
823	write(116,*) wavei(nw),fmneti_nu(L_NLAYRAD,nw)
824	enddo
825	close(116)
826	endif
827
828	! OLR by band
829	diagrad_OLR=.true.
830	if(diagrad_OLR)then
831	open(117,file='OLRnu.out')
832	write(117,*) 'IR wavel - band width - OLR'
833	do nw=1,L_NSPECTI
834	write(117,*) wavei(nw),
835	& abs(1.e4/bwnv(nw)-1.e4/bwnv(nw+1)),OLR_nu(1,nw)
836	enddo
837	close(117)
838	endif
839
840	! OLR vs altitude: in a .txt file
841	diagrad_OLRz=.true.
842	if(diagrad_OLRz)then
843	open(118,file='OLRz_plevs.out')
844	open(119,file='OLRz.out')
845	do l=1,L_NLAYRAD
846	write(118,) plevrad(2l)
847	do nw=1,L_NSPECTI
848	write(119,*) fluxupi_nu(l,nw)
849	enddo
850	enddo
851	close(118)
852	close(119)
853	endif
854
855	! Heating rates vs altitude in a .txt file
856	diagrad_rates=.true.
857	if(diagrad_rates)then
858	open(120,file='heating_rates.out')
859	write(120,*) "Pressure - Alt - HR tot - Rates (wavel SW)"
860	do l=1,nlayermx
861	write(120,*) pplay(1,l),zzlay(1,l),dtsw(1,l),dtsw_nu(l,:)
862	enddo
863	close(120)
864
865	open(121,file='cooling_rates.out')
866	write(121,*) "Pressure - Alt - CR tot - Rates (wavel LW)"
867	do l=1,nlayermx
868	write(121,*) pplay(1,l),zzlay(1,l),dtlw(1,l),dtlw_nu(l,:)
869	enddo
870	close(121)
871
872	open(122,file='bands.out')
873	write(122,*) "wavel - bands boundaries (microns)"
874	do nw=1,L_NSPECTV
875	write(122,*) wavev(nw),1.e4/bwnv(nw+1),1.e4/bwnv(nw)
876	enddo
877	do nw=1,L_NSPECTI
878	write(122,*) wavei(nw),1.e4/bwni(nw+1),1.e4/bwni(nw)
879	enddo
880	close(122)
881
882	open(123,file='c2h2_rates.out')
883	write(123,*) "Pressure - Alt - CR c2h2"
884	do l=1,nlayermx
885	write(123,*) pplay(1,l),zzlay(1,l),dtlw_hcn_c2h2(1,l)
886	enddo
887	close(123)
888
889	endif
890
891	endif ! ngrid.eq.1
892
893	endif ! lastcall
894
895	end

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

Download in other formats:

Original Format