Context Navigation

nirco2abs.F @ 690

Last change on this file since 690 was 690, checked in by acolaitis, 12 years ago
Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration
File size: 8.5 KB

Rev	Line
[414]	1	SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq,
[38]	2	$ mu0,fract,declin,pdtnirco2)
	3
	4	IMPLICIT NONE
	5	c=======================================================================
	6	c subject:
	7	c --------
	8	c Computing heating rate due to
	9	c absorption by CO2 in the near-infrared
	10	c This version includes NLTE effects
	11	c
	12	c (Scheme to be described in Forget et al., JGR, 2003)
	13	c (old Scheme described in Forget et al., JGR, 1999)
	14	c
	15	c This version updated with a new functional fit,
	16	c see NLTE correction-factor of Lopez-Valverde et al (1998)
	17	c Stephen Lewis 2000
[414]	18	c
	19	c jul 2011 malv+fgg New corrections for NLTE implemented
[38]	20	c 08/2002 : correction for bug when running with diurnal=F
	21	c
	22	c author: Frederic Hourdin 1996
	23	c ------
	24	c Francois Forget 1999
	25	c
	26	c input:
	27	c -----
	28	c ngrid number of gridpoint of horizontal grid
	29	c nlayer Number of layer
	30	c dist_sol sun-Mars distance (AU)
	31	c mu0(ngridmx)
	32	c fract(ngridmx) day fraction of the time interval
	33	c declin latitude of subslar point
	34	c
	35	c output:
	36	c -------
	37	c
	38	c pdtnirco2(ngrid,nlayer) Heating rate (K/s)
	39	c
	40	c
	41	c=======================================================================
	42	c
	43	c 0. Declarations :
	44	c ------------------
	45	c
	46	#include "dimensions.h"
	47	#include "dimphys.h"
	48	#include "comcstfi.h"
	49	#include "callkeys.h"
	50	#include "comdiurn.h"
[498]	51	#include "nirdata.h"
[552]	52	#include "tracer.h"
[38]	53
	54	c-----------------------------------------------------------------------
	55	c Input/Output
	56	c ------------
[552]	57	integer,intent(in) :: ngrid ! number of (horizontal) grid points
	58	integer,intent(in) :: nlayer ! number of atmospheric layers
	59	real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
	60	real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
	61	integer,intent(in) :: nq ! number of tracers
	62	real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
	63	real,intent(in) :: mu0(ngridmx) ! solar angle
	64	real,intent(in) :: fract(ngridmx) ! day fraction of the time interval
	65	real,intent(in) :: declin ! latitude of sub-solar point
	66
	67	real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
[38]	68	c
	69	c Local variables :
	70	c -----------------
[414]	71	INTEGER l,ig, n, nstep,i
[38]	72	REAL co2heat0, zmu(ngridmx)
	73
	74	c special diurnal=F
	75	real mu0_int(ngridmx),fract_int(ngridmx),zday_int
	76	real ztim1,ztim2,ztim3,step
	77
	78	c
	79	c local saved variables
	80	c ---------------------
[552]	81	logical,save :: firstcall=.true.
[575]	82	integer,save :: ico2=0 ! index of "co2" tracer
	83	integer,save :: io=0 ! index of "o" tracer
[38]	84	c p0noonlte is a pressure below which non LTE effects are significant.
	85	c REAL p0nonlte
	86	c DATA p0nonlte/7.5e-3/
	87	c SAVE p0nonlte
	88
	89	c parameters for CO2 heating fit
	90	real n_a, n_p0, n_b
	91	parameter (n_a=1.1956475)
	92	parameter (n_b=1.9628251)
	93	parameter (n_p0=0.0015888279)
	94
[414]	95	c Variables added to implement NLTE correction factor (feb 2011)
	96	real pyy(nlayer)
	97	real cor1(nlayer),oldoco2(nlayer),alfa2(nlayer)
	98	real p2011,cociente1,merge
	99	real cor0,oco2gcm
	100
[38]	101	c----------------------------------------------------------------------
	102
	103	c Initialisation
	104	c --------------
[552]	105	if (firstcall) then
	106	if (nircorr.eq.1) then
	107	! we will need co2 and o tracers
	108	ico2=igcm_co2
	109	if (ico2==0) then
	110	write(,) "nirco2abs error: I need a CO2 tracer"
	111	write(,) " when running with nircorr==1"
	112	stop
	113	endif
	114	io=igcm_o
	115	if (io==0) then
	116	write(,) "nirco2abs error: I need an O tracer"
	117	write(,) " when running with nircorr==1"
	118	stop
	119	endif
	120	endif
	121	firstcall=.false.
	122	endif
	123
	124
[38]	125	c co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
	126	co2heat0=n_a(1.52/dist_sol)*2/daysec
	127
	128	c Simple calcul for a given sun incident angle (if diurnal=T)
	129	c --------------------------------------------
	130
	131	IF (diurnal) THEN
	132	do ig=1,ngrid
	133	zmu(ig)=sqrt(1224.mu0(ig)mu0(ig)+1.)/35.
[414]	134
	135	if(nircorr.eq.1) then
	136	do l=1,nlayer
	137	pyy(l)=pplay(ig,l)
	138	enddo
	139
	140	call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
	141	call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
	142	call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
	143	endif
	144
	145	do l=1,nlayer
	146	! Calculations for the O/CO2 correction
	147	if(nircorr.eq.1) then
	148	cor0=1./(1.+n_p0/pplay(ig,l))**n_b
[552]	149	if(pq(ig,l,ico2).gt.1.e-6) then
	150	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
[414]	151	else
	152	oco2gcm=1.e6
	153	endif
	154	cociente1=oco2gcm/oldoco2(l)
	155	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
	156	$ (1.-alfa2(l))
	157	merge=10**merge
	158	p2011=sqrt(merge)*cor0
	159	else if (nircorr.eq.0) then
	160	p2011=1.
	161	cor1(l)=1.
	162	endif
	163
	164	if(fract(ig).gt.0.) pdtnirco2(ig,l)=
[38]	165	& co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
	166	& /(1.+n_p0/pplay(ig,l))**n_b
[414]	167	! Corrections from tabulation
	168	$ * cor1(l) * p2011
[38]	169	c OLD SCHEME (forget et al. 1999)
	170	c s co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
	171	c s / (1+p0nonlte/pplay(ig,l))
	172	enddo
	173	enddo
[414]	174
[38]	175
	176	c Averaging over diurnal cycle (if diurnal=F)
	177	c -------------------------------------------
	178	c NIR CO2 abs is slightly non linear. To remove the diurnal
	179	c cycle, it is better to average the heating rate over 1 day rather
	180	c than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)
	181
	182	ELSE ! if (.not.diurnal) then
	183
	184	nstep = 20 ! number of integration step /sol
	185	do n=1,nstep
	186	zday_int = (n-1)/float(nstep)
	187	ztim2=COS(declin)COS(2.pi*(zday_int-.5))
	188	ztim3=-COS(declin)SIN(2.pi*(zday_int-.5))
	189	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	190	s ztim1,ztim2,ztim3,
	191	s mu0_int,fract_int)
	192	do ig=1,ngrid
	193	zmu(ig)=sqrt(1224.mu0_int(ig)mu0_int(ig)+1.)/35.
[414]	194
	195	if(nircorr.eq.1) then
	196	do l=1,nlayer
	197	pyy(l)=pplay(ig,l)
	198	enddo
[690]	199	call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
	200	call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
	201	call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
[414]	202	endif
	203
	204	do l=1,nlayer
	205	if(nircorr.eq.1) then
	206	cor0=1./(1.+n_p0/pplay(ig,l))**n_b
[552]	207	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
[414]	208	cociente1=oco2gcm/oldoco2(l)
	209	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
	210	$ (1.-alfa2(l))
	211	merge=10**merge
	212	p2011=sqrt(merge)*cor0
	213	else if (nircorr.eq.0) then
	214	p2011=1.
	215	cor1(l)=1.
	216	endif
	217
[38]	218	if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
	219	& pdtnirco2(ig,l) + (1/float(nstep))*
	220	& co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
	221	& /(1.+n_p0/pplay(ig,l))**n_b
[414]	222	! Corrections from tabulation
	223	$ * cor1(l) * p2011
[38]	224	enddo
	225	enddo
	226	end do
	227	END IF
	228
	229	return
	230	end
	231
[414]	232
	233
	234	subroutine interpnir(escout,p,nlayer,escin,pin,nl)
	235	C
	236	C subroutine to perform linear interpolation in pressure from 1D profile
	237	C escin(nl) sampled on pressure grid pin(nl) to profile
	238	C escout(nlayer) on pressure grid p(nlayer).
	239	C
	240	real escout(nlayer),p(nlayer)
	241	real escin(nl),pin(nl),wm,wp
	242	integer nl,nlayer,n1,n,nm,np
	243	do n1=1,nlayer
	244	if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
	245	escout(n1) = 0.0
	246	else
	247	do n = 1,nl-1
	248	if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
	249	nm=n
	250	np=n+1
	251	wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
	252	wp=1.0 - wm
	253	endif
	254	enddo
	255	escout(n1) = escin(nm)wm + escin(np)wp
	256	endif
	257	enddo
	258	return
	259	end

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

Download in other formats:

Original Format