Context Navigation

callcorrk_pluto_mod.F90 @ 3572

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

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

Download in other formats:

Original Format