Context Navigation

callcorrk_pluto_mod.F90 @ 3684

Last change on this file since 3684 was 3683, checked in by debatzbr, 3 months ago
Solves random variations of the SW heating at the model top BBT
File size: 33.1 KB

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

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

Download in other formats:

Original Format