Context Navigation

callcorrk_pluto_mod.F90 @ 3695

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

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

Download in other formats:

Original Format