Context Navigation

callcorrk_pluto_mod.F90 @ 3504

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

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

Download in other formats:

Original Format