Context Navigation

nirco2abs.F @ 3094

Last change on this file since 3094 was 2580, checked in by slebonnois, 3 years ago
SL: Version used for VCD 1.1 (tuneupperatm => key+photochemistry, nirco2abs) add prod and loss from photochem new nonorographic routine (not used for VCD 1.1) update of deftank files
File size: 15.6 KB

Rev	Line
[1310]	1	SUBROUTINE nirco2abs(nlon,nlev,nplay,dist_sol,nq,pq,
[1442]	2	$ mu0,fract,pdtnirco2)
[1310]	3
	4	use dimphy
[1545]	5	use geometry_mod, only: longitude_deg, latitude_deg
[1442]	6	use chemparam_mod, only: i_co2, i_o
[1310]	7	c use compo_hedin83_mod2
	8
	9
	10	IMPLICIT NONE
	11	c=======================================================================
	12	c subject:
	13	c --------
	14	c Computing heating rate due to
	15	c absorption by CO2 in the near-infrared
	16	c This version includes NLTE effects
	17	c
	18	c (Scheme to be described in Forget et al., JGR, 2003)
	19	c (old Scheme described in Forget et al., JGR, 1999)
	20	c
	21	c This version updated with a new functional fit,
	22	c see NLTE correction-factor of Lopez-Valverde et al (1998)
	23	c Stephen Lewis 2000
	24	c
[2198]	25	c apr 2019 d.quirino Improving NLTE params, SOIR/SPICAV Temp comparison
[1442]	26	c oct 2014 g.gilli Coupling with photochemical model
	27	C jan 2014 g.gilli Revision (following martian non-lte param)
[1310]	28	C jun 2013 l.salmi First adaptation to Venus and NIR NLTE param
	29	c jul 2011 malv+fgg New corrections for NLTE implemented
	30	c 08/2002 : correction for bug when running with diurnal=F
	31	c
	32	c author: Frederic Hourdin 1996
	33	c ------
	34	c Francois Forget 1999
	35	c
	36	c input:
	37	c -----
	38	c nlon number of gridpoint of horizontal grid
	39	c nlev Number of layer
	40	c dist_sol sun-Venus distance (AU)
	41	c mu0(nlon)
	42	c fract(nlon) day fraction of the time interval
	43	c declin latitude of subslar point
	44	c
	45	c output:
	46	c -------
	47	c
	48	c pdtnirco2(nlon,nlev) Heating rate (K/sec)
	49	c
	50	c
	51	c=======================================================================
	52	c
	53	c 0. Declarations :
	54	c ------------------
	55	c
	56
	57	#include "YOMCST.h"
	58	#include "clesphys.h"
	59	c#include "comdiurn.h"
	60	#include "nirdata.h"
	61	c#include "tracer.h"
[1442]	62	#include "mmol.h"
[1310]	63	c-----------------------------------------------------------------------
	64	c Input/Output
	65	c ------------
	66	integer,intent(in) :: nlon ! number of (horizontal) grid points
	67	integer,intent(in) :: nlev ! number of atmospheric layers
[1442]	68
[1310]	69	real,intent(in) :: nplay(nlon,nlev) ! Pressure
	70	real,intent(in) :: dist_sol ! Sun-Venus distance (in AU)
	71	integer,intent(in) :: nq ! number of tracers
[1442]	72	real,intent(in) :: pq(nlon,nlev,nq) ! mass mixing ratio tracers
[1310]	73	real,intent(in) :: mu0(nlon) ! solar angle
	74	real,intent(in) :: fract(nlon) ! day fraction of the time interval
	75	c real,intent(in) :: declin ! latitude of sub-solar point
[1442]	76	real :: co2vmr_gcm(nlon,nlev), o3pvmr_gcm(nlon,nlev)
	77
[1310]	78	real,intent(out) :: pdtnirco2(nlon,nlev) ! heating rate (K/sec)
	79
	80	c
	81	c Local variables :
	82	c -----------------
	83	INTEGER l,ig, n, nstep,i
	84	REAL co2heat0, zmu(nlon)
	85
	86	c special diurnal=F
	87	real mu0_int(nlon),fract_int(nlon),zday_int
	88	real ztim1,ztim2,ztim3,step
	89
[2580]	90	logical onepeak
	91	parameter (onepeak=.false.)
	92	c parameter (onepeak=.true.)
[1310]	93	c
	94	c local saved variables
	95	c ---------------------
	96	logical,save :: firstcall=.true.
	97	integer,save :: ico2=0 ! index of "co2" tracer
	98	integer,save :: io=0 ! index of "o" tracer
	99
[2580]	100	ccc=================================================
[1442]	101	cccc parameters for CO2 heating fit
[2580]	102	ccc=================================================
	103
	104	c--------------------------------------------------
	105	c One-peak martian-type fit => Gabriella (2014+)
	106	c--------------------------------------------------
[1310]	107	c n_a = heating rate for Venusian day at p0, r0, mu =0 [K day-1]
	108	c Here p0 = p_cloud top [Pa]
	109	c n_p0 = is a pressure below which non LTE effects are significant [Pa]
	110	c n_a Solar heating [K/Eday] at the cloud top, taken from Crisps table
	111
	112	real n_a, n_p0, n_b, p_ctop
	113
[2580]	114	cc "Nominal" values used in Gilli+2017
[2198]	115	c parameter (n_a = 18.13/86400.0) !c K/Eday ---> K/sec
	116	c parameter (p_ctop=13.2e2)
	117	c parameter (n_p0=0.008)
	118
	119	cc "New" values used to improve SPICAV/SOIR Temperature comparision (D.Quirino)
[2580]	120	cc Gilli+2021
[2198]	121	parameter (n_a = 15.92/86400.0) !c K/Eday ---> K/sec
	122	parameter (p_ctop=19.85e2)
	123	parameter (n_p0=0.1)
[1442]	124	parameter (n_b=1.362)
[1310]	125
[1442]	126	c -- NLTE Param v2 --
	127	C parameter (n_p0=0.01)
	128	c parameter (n_b = 1.3)
[1310]	129
[2580]	130	c--------------------------------------------------
	131	c Multi-peaks Roldan-type fit => Laura (2013)
	132	c New paramaters (Param9*0.5) => Enora (2021)
	133	c--------------------------------------------------
	134	c ENORA FINE TUNING used for VCD 1.1
	135	c (fit to fig 12 Roldan-2000)
	136	real n_coFB, n_aFB, n_bFB, n_p0FB, n_eFB
	137	real n_coISO, n_aISO, n_bISO, n_p0ISO, n_eISO
	138	real n_coFH, n_aFH, n_bFH, n_p0FH, n_eFH
	139	real n_co43, n_a43, n_b43, n_p043, n_e43
	140	real n_co43b, n_a43b, n_b43b, n_p043b, n_e43b
	141	real n_conir, n_anir, n_bnir, n_p0nir, n_enir
	142
	143	parameter (n_coFB=119./86400.0) !c K/Eday ---> K/sec
	144	parameter (n_aFB=0.185)
	145	parameter (n_bFB=3.7)
	146	parameter (n_p0FB=2.9e-4)
	147	parameter (n_eFB=0.76)
	148
	149	parameter (n_coISO=265./86400.0) !c K/Eday ---> K/sec
	150	parameter (n_aISO=0.313)
	151	parameter (n_bISO=1.65)
	152	parameter (n_p0ISO=0.076)
	153	parameter (n_eISO=0.99)
	154
	155	parameter (n_coFH=2.5/86400.0) !c K/Eday ---> K/sec
	156	parameter (n_aFH=3.98)
	157	parameter (n_bFH=2.9)
	158	parameter (n_p0FH=0.17)
	159	parameter (n_eFH=2.16)
	160
	161	parameter (n_co43=55./86400.0) !c K/Eday ---> K/sec
	162	parameter (n_a43=0.625)
	163	parameter (n_b43=2.6)
	164	parameter (n_p043=0.043)
	165	parameter (n_e43=1.654)
	166
	167	! parameter (n_co43b=100./86400.0) !c K/Eday ---> K/sec
	168	! => fine tuning: not affected by the *0.5 below (see ENORA FINE TUNING)
	169	parameter (n_co43b=200./86400.0) !c K/Eday ---> K/sec
	170	parameter (n_a43b=5.5)
	171	parameter (n_b43b=2.3)
	172	parameter (n_p043b=1.)
	173	parameter (n_e43b=0.4)
	174
	175	parameter (n_conir=6.5/86400.0) !c K/Eday ---> K/sec
	176	parameter (n_anir=35.65)
	177	parameter (n_bnir=2.1)
	178	parameter (n_p0nir=0.046)
	179	parameter (n_enir=0.9)
	180
	181	real :: picFB(nlon,nlev), picISO(nlon,nlev), picFH(nlon,nlev)
	182	real :: pic43(nlon,nlev), pic43b(nlon,nlev), picnir(nlon,nlev)
	183
	184	ccc=================================================
	185
[1310]	186	c Variables added to implement NLTE correction factor (feb 2011)
	187	real pyy(nlev)
	188	real cor1(nlev),oldoco2(nlev),alfa2(nlev)
	189	real p2011,cociente1,merge
	190	real cor0,oco2gcm
	191
[2580]	192	c----------------------------------------------------------------------
[1310]	193	c Initialisation
	194	c --------------
[1442]	195	if (firstcall) then
	196	if (nircorr.eq.1) then
	197	c ! we will need co2 and o tracers
	198	ico2= i_co2
	199	if (ico2==0) then
	200	write(,) "nirco2abs error: I need a CO2 tracer"
	201	write(,) " when running with nircorr==1"
	202	stop
	203	endif
	204	io=i_o
	205	if (io==0) then
	206	write(,) "nirco2abs error: I need an O tracer"
	207	write(,) " when running with nircorr==1"
	208	stop
	209	endif
	210	endif
	211	firstcall=.false.
	212	endif
[2580]	213	c --------------
	214	c co2heat0 is correction for dist_sol (is 1 for Venus)
	215	co2heat0=(0.7233/dist_sol)**2
[1310]	216
[2580]	217	pdtnirco2(:,:)=0.
	218	c----------------------------------------------------------------------
	219
[1310]	220	c
	221	c Simple calcul for a given sun incident angle (if cycle_diurne=T)
	222	c --------------------------------------------
	223
	224	IF (cycle_diurne) THEN
	225
	226	do ig=1,nlon
	227	zmu(ig)=sqrt(1224.mu0(ig)mu0(ig)+1.)/35.
	228
[2580]	229	c---------------------------
	230	if (onepeak) then
	231	c---------------------------
[1310]	232	if(nircorr.eq.1) then
	233	do l=1,nlev
	234	pyy(l)=nplay(ig,l)
	235	enddo
	236	call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
	237	call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
	238	call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
	239	endif
	240	do l=1,nlev
	241	c Calculations for the O/CO2 correction
	242	if(nircorr.eq.1) then
	243	cor0=1./(1.+n_p0/nplay(ig,l))**n_b
[1442]	244	if(pq(ig,l,ico2) .gt. 1.e-6) then
	245	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
[2464]	246	! handle the rare cases when pq(ig,l,io)<0
	247	if (pq(ig,l,io).lt.0) then
	248	write(,) "nirco2abs: warning ig=",ig," l=",l,
	249	& " pq(ig,l,io)=",pq(ig,l,io)
	250	oco2gcm=1.e6
	251	endif
[1310]	252	else
	253	oco2gcm=1.e6
	254	endif
	255	cociente1=oco2gcm/oldoco2(l)
[1442]	256
	257	c WRITE(,) "nirco2abs line 211", l, cociente1
	258
[1310]	259	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
	260	$ (1.-alfa2(l))
	261	merge=10**merge
	262	p2011=sqrt(merge)*cor0
	263
	264	else if (nircorr.eq.0) then
	265	p2011=1.
	266	cor1(l)=1.
	267	endif
	268
	269	if(fract(ig).gt.0.) pdtnirco2(ig,l)=
[2580]	270	& co2heat0n_asqrt((p_ctop*zmu(ig))/nplay(ig,l))
[1310]	271	& /(1.+n_p0/nplay(ig,l))**n_b
	272	c Corrections from tabulation
	273	$ * cor1(l) * p2011
[1442]	274
[2580]	275	enddo !nlev
	276	c---------------------------
	277	else ! multipeak
	278	c---------------------------
	279	do l=1,nlev
	280	if(fract(ig).gt.0.) then
	281	picFB(ig,l)=n_coFB
	282	& ((n_aFB/nplay(ig,l))*n_eFB)
	283	& zmu(ig)*0.82
	284	& /(1.+n_p0FB/nplay(ig,l))**n_bFB
	285
	286	picISO(ig,l)=n_coISO
	287	& ((n_aISO/nplay(ig,l))*n_eISO)
	288	& zmu(ig)*0.55
	289	& /(1.+n_p0ISO/nplay(ig,l))**n_bISO
	290
	291	picFH(ig,l)=n_coFH
	292	& ((n_aFH/nplay(ig,l))*n_eFH)
	293	& zmu(ig)*0.55
	294	& /(1.+n_p0FH/nplay(ig,l))**n_bFH
	295
	296	pic43(ig,l)=n_co43
	297	& ((n_a43/nplay(ig,l))*n_e43)
	298	& zmu(ig)*0.55
	299	& /(1.+n_p043/nplay(ig,l))**n_b43
	300
	301	pic43b(ig,l)=n_co43b
	302	& ((n_a43b/nplay(ig,l))*n_e43b)
	303	& zmu(ig)*0.55
	304	& /(1.+n_p043b/nplay(ig,l))**n_b43b
	305
	306	picnir(ig,l)=n_conir
	307	& ((n_anir/nplay(ig,l))*n_enir)
	308	& zmu(ig)*0.55
	309	& /(1.+n_p0nir/nplay(ig,l))**n_bnir
	310
	311	pdtnirco2(ig,l)=co2heat0*
	312	& (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
	313	& +pic43b(ig,l)+picnir(ig,l))0.5 ! 0.5 = ENORA FINE TUNING
	314
	315	endif
	316	enddo !nlev
	317	c---------------------------
	318	endif
	319	c---------------------------
	320	enddo !nlon
	321
	322
[1310]	323	c Averaging over diurnal cycle (if diurnal=F)
	324	c -------------------------------------------
	325	c NIR CO2 abs is slightly non linear. To remove the diurnal
	326	c cycle, it is better to average the heating rate over 1 day rather
	327	c than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)
	328
	329	ELSE ! if (.not.diurnal) then
	330	nstep = 20 ! number of integration step /sol
	331	do n=1,nstep
	332
	333	zday_int = (n-1)/float(nstep)
	334
[1545]	335	CALL zenang(0.,zday_int,RDAY/nstep,
	336	& latitude_deg,longitude_deg,
	337	& mu0_int,fract_int)
[1310]	338
[2580]	339	do ig=1,nlon
[1310]	340	zmu(ig)=sqrt(1224.mu0_int(ig)mu0_int(ig)+1.)/35.
	341
[2580]	342	c---------------------------
	343	if (onepeak) then
	344	c---------------------------
	345	if(nircorr.eq.1) then
[1310]	346	do l=1,nlev
[2580]	347	pyy(l)=nplay(ig,l)
	348	enddo
	349	call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
	350	call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
	351	call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
	352	endif
	353	do l=1,nlev
[1442]	354	c Calculations for the O/CO2 correction
	355	if(nircorr.eq.1) then
	356	cor0=1./(1.+n_p0/nplay(ig,l))**n_b
[2580]	357	if(pq(ig,l,ico2) .gt. 1.e-6) then
	358	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
	359	! handle the rare cases when pq(ig,l,io)<0
	360	if (pq(ig,l,io).lt.0) then
	361	write(,) "nirco2abs: warning ig=",ig," l=",l,
	362	& " pq(ig,l,io)=",pq(ig,l,io)
	363	oco2gcm=1.e6
	364	endif
	365	else
	366	oco2gcm=1.e6
	367	endif
[1442]	368	cociente1=oco2gcm/oldoco2(l)
[2580]	369
	370	c WRITE(,) "nirco2abs line 211", l, cociente1
	371
[1442]	372	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
	373	$ (1.-alfa2(l))
	374	merge=10**merge
	375	p2011=sqrt(merge)*cor0
[1310]	376
[1442]	377	else if (nircorr.eq.0) then
	378	p2011=1.
	379	cor1(l)=1.
	380	endif
[1310]	381
[2580]	382	if(fract(ig).gt.0.) pdtnirco2(ig,l)=
[1442]	383	& pdtnirco2(ig,l) + (1/float(nstep))*
[2580]	384	& co2heat0n_asqrt((p_ctop*zmu(ig))/nplay(ig,l))
	385	& /(1.+n_p0/nplay(ig,l))**n_b
	386	c Corrections from tabulation
	387	$ * cor1(l) * p2011
	388
	389	enddo !nlev
	390	c---------------------------
	391	else ! multipeak
	392	c---------------------------
	393	do l=1,nlev
	394	if(fract(ig).gt.0.) then
	395	picFB(ig,l)=n_coFB
	396	& ((n_aFB/nplay(ig,l))*n_eFB)
	397	& zmu(ig)*0.82
	398	& /(1.+n_p0FB/nplay(ig,l))**n_bFB
[1310]	399
[2580]	400	picISO(ig,l)=n_coISO
	401	& ((n_aISO/nplay(ig,l))*n_eISO)
	402	& zmu(ig)*0.55
	403	& /(1.+n_p0ISO/nplay(ig,l))**n_bISO
	404
	405	picFH(ig,l)=n_coFH
	406	& ((n_aFH/nplay(ig,l))*n_eFH)
	407	& zmu(ig)*0.55
	408	& /(1.+n_p0FH/nplay(ig,l))**n_bFH
	409
	410	pic43(ig,l)=n_co43
	411	& ((n_a43/nplay(ig,l))*n_e43)
	412	& zmu(ig)*0.55
	413	& /(1.+n_p043/nplay(ig,l))**n_b43
	414
	415	pic43b(ig,l)=n_co43b
	416	& ((n_a43b/nplay(ig,l))*n_e43b)
	417	& zmu(ig)*0.55
	418	& /(1.+n_p043b/nplay(ig,l))**n_b43b
	419
	420	picnir(ig,l)=n_conir
	421	& ((n_anir/nplay(ig,l))*n_enir)
	422	& zmu(ig)*0.55
	423	& /(1.+n_p0nir/nplay(ig,l))**n_bnir
	424
	425	pdtnirco2(ig,l)=
	426	& pdtnirco2(ig,l)+(1/float(nstep))co2heat0
	427	& (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
	428	& +pic43b(ig,l)+picnir(ig,l))0.5 ! 0.5 = ENORA FINE TUNING
	429
	430	endif
	431	enddo !nlev
	432	c---------------------------
	433	endif
	434	c---------------------------
	435	enddo !nlon
	436	enddo !nstep
[1310]	437
[2580]	438	END IF ! diurnal cycle
[1310]	439
	440	return
	441	end
	442
	443
	444	subroutine interpnir(escout,p,nlev,escin,pin,nl)
	445	C
	446	C subroutine to perform linear interpolation in pressure from 1D profile
	447	C escin(nl) sampled on pressure grid pin(nl) to profile
	448	C escout(nlev) on pressure grid p(nlev).
	449	C
	450	real escout(nlev),p(nlev)
	451	real escin(nl),pin(nl),wm,wp
	452	integer nl,nlev,n1,n,nm,np
	453	do n1=1,nlev
	454	if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
[1442]	455	c escout(n1) = 0.0
	456	escout(n1) = 1.e-15
[1310]	457	else
	458	do n = 1,nl-1
	459	if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
	460	nm=n
	461	np=n+1
	462	wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
	463	wp=1.0 - wm
	464	endif
	465	enddo
	466	escout(n1) = escin(nm)wm + escin(np)wp
	467	endif
	468	enddo
	469	return
	470	end

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

Download in other formats:

Original Format