Context Navigation

callcorrk_pluto_mod.F90 @ 3698

Last change on this file since 3698 was 3698, checked in by emillour, 9 months ago

Pluto PCM:
Some fixes to enable runnnig with dynamico:

add "strictboundcorrk" flag in callcorrk_pluto to enable running even if outside of kmatrix temperatures (when strictboundcorrk=.true.)
add premature exiting of writediagsoil if not with lon-lat grid
while at it, turned surfprop.F90 into a module
in physiq, enforce the possibility to output subsurface-related field in most cases, not only when "fast=.true."
adapt reference xml files: subsurface quantities need to be defined on a dedicated grid, otherwise XIOS will generate misleading garbage values. Updated files are put in "deftank/dynamico" for now.

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

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

Download in other formats:

Original Format