Context Navigation

calchim_asis.F90 @ 2725

Last change on this file since 2725 was 2542, checked in by aslmd, 3 years ago

Generic GCM:

Large update of the chemical modules

Read chemical network from input files
Init chemistry with ModernTrac?
Photolysis online calculation

Property svn:executable set to *

File size: 20.1 KB

Rev	Line
[1796]	1	subroutine calchim_asis(ngrid,nlayer,nq, &
	2	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,fract, &
	3	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim, &
	4	tauref,co2ice, &
[2542]	5	pu,pdu,pv,pdv,surfdust,surfice,icount,zdtchim)
[1796]	6
[2542]	7	use tracer_h, only: mmol, noms, nesp, is_chim
[1796]	8
	9	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
	10	USE comcstfi_mod
	11	use callkeys_mod
[2542]	12	use time_phylmdz_mod, only: ecritphy, iphysiq ! output rate set by ecritphy
	13	use types_asis, only: v_phot_3d, jlabel, nb_phot_hv_max
	14	use chimiedata_h
	15	use radcommon_h, only: glat
[1796]	16
	17	implicit none
	18
	19	!=======================================================================
	20	!
	21	! subject:
	22	! --------
	23	!
	24	! Prepare the call for the photochemical module, and send back the
	25	! tendencies from photochemistry in the chemical species mass mixing ratios
	26	!
	27	! Author: Sebastien Lebonnois (08/11/2002)
	28	! -------
	29	! update 12/06/2003 for water ice clouds and compatibility with dust
	30	! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
	31	! update 03/05/2005 cosmetic changes (Franck Lefevre)
	32	! update sept. 2008 identify tracers by their names (Ehouarn Millour)
	33	! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
	34	! update 16/03/2012 optimization (Franck Lefevre)
	35	! update 11/12/2013 optimization (Franck Lefevre)
[1811]	36	! update 20/10/2017 adapted to LMDZ GENERIC+cosmetic changes (Benjamin Charnay)
[2542]	37	! update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
[1796]	38	!
	39	! Arguments:
	40	! ----------
	41	!
	42	! Input:
	43	!
	44	! ptimestep timestep (s)
	45	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	46	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
	47	! pt(ngrid,nlayer) Temperature (K)
	48	! pdt(ngrid,nlayer) Temperature tendency (K)
	49	! pu(ngrid,nlayer) u component of the wind (ms-1)
	50	! pdu(ngrid,nlayer) u component tendency (K)
	51	! pv(ngrid,nlayer) v component of the wind (ms-1)
	52	! pdv(ngrid,nlayer) v component tendency (K)
	53	! dist_sol distance of the sun (AU)
	54	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
	55	! fract(ngrid) day fraction (for diurnal=false)
	56	! pq(ngrid,nlayer,nq) Advected fields, ie chemical species here
	57	! pdq(ngrid,nlayer,nq) Previous tendencies on pq
	58	! tauref(ngrid) Optical depth at 7 hPa
	59	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
	60	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
	61	! surfice(ngrid,nlayer) ice surface area (m2/m3)
	62	!
	63	! Output:
	64	!
[2542]	65	! dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry
	66	! dqschim(ngrid,nesp) ! tendencies on qsurf
[1796]	67	!
	68	!=======================================================================
	69
	70	! input:
	71
	72	integer,intent(in) :: ngrid ! number of atmospheric columns
	73	integer,intent(in) :: nlayer ! number of atmospheric layers
	74	integer,intent(in) :: nq ! number of tracers
	75	real :: ptimestep
	76	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
[2542]	77	real :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers
[1796]	78	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
	79	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
	80	real :: pt(ngrid,nlayer) ! temperature
	81	real :: pdt(ngrid,nlayer) ! temperature tendency
	82	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
	83	real :: pdu(ngrid,nlayer) ! u component tendency
	84	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
	85	real :: pdv(ngrid,nlayer) ! v component tendency
	86	real :: dist_sol ! distance of the sun (AU)
	87	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
	88	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
	89	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
	90	real :: zday ! date (time since Ls=0, in martian days)
	91	real :: tauref(ngrid) ! optical depth at 7 hPa
	92	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
	93	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
	94	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
	95
	96	! output:
	97
[2542]	98	real :: dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry
	99	real :: dqschim(ngrid,nesp) ! tendencies on qsurf
[1796]	100
	101	! local variables:
	102
[2542]	103	integer :: iloc(1) ! index of major species
	104	integer :: ig,l,i,iq,iqmax,iesp
	105	integer :: foundswitch, lswitch ! to switch between photochem and thermochem ? (YJ)
[1796]	106	integer,save :: chemthermod
	107
[2542]	108	real :: zq(ngrid,nlayer,nesp) ! pq+pdq*ptimestep before chemistry
[1796]	109	! new mole fraction after
	110	real :: zt(ngrid,nlayer) ! temperature
	111	real :: zu(ngrid,nlayer) ! u component of the wind
	112	real :: zv(ngrid,nlayer) ! v component of the wind
	113	real :: taucol ! optical depth at 7 hPa
[2542]	114	real :: xmmol(nlayer,nesp) ! mmol/mmean but only for chemical species
	115
[1796]	116	logical,save :: firstcall = .true.
	117
	118	! for each column of atmosphere:
	119
	120	real :: zpress(nlayer) ! Pressure (mbar)
	121	real :: zdens(nlayer) ! Density (cm-3)
	122	real :: ztemp(nlayer) ! Temperature (K)
	123	real :: zlocal(nlayer) ! Altitude (km)
[2542]	124	real :: zycol(nlayer,nesp) ! Composition (mole fractions)
[1796]	125	real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1)
	126	real :: szacol ! Solar zenith angle
	127	real :: fract(ngrid) ! day fraction
	128	real :: fractcol ! day fraction
	129	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
	130	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
	131
	132	integer :: iter(nlayer) ! Number of chemical iterations
	133	! within one physical timestep
[2542]	134	integer :: icount
[1796]	135	! for output:
	136
[2542]	137	logical :: output ! to issue calls to writediagfi and stats
[1796]	138	parameter (output = .true.)
[2542]	139	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
	140	! within one physical timestep
	141	integer :: ierr
	142	real :: zdtchim(ngrid,nlayer) ! temperature modification by chemistry
	143	real :: dEzchim(ngrid,nlayer) ! energy modification by chemistry
	144	real :: zdtchim_output(ngrid) ! flux modification by chemistry in W.m-2
[1796]	145
	146	!=======================================================================
	147	! initialization of the chemistry (first call only)
	148	!=======================================================================
	149
	150	if (firstcall) then
	151
[2542]	152	!! Moved to routine indice in photochemistry_asis
	153	!! because nb_phot_hv_max value needed in order
	154	!! to choose if we call read_phototable or not.
	155	!! A cleaner solution need to be find.
	156	! if (photochem .and. .not. jonline) then
	157	! print*,'calchim: Read photolysis lookup table'
	158	! call read_phototable
	159	! end if
	160
	161	if (.not.allocated(SF_mode)) allocate(SF_mode(nesp))
	162	if (.not.allocated(SF_value)) allocate(SF_value(nesp))
	163	if (.not.allocated(prod_rate)) allocate(prod_rate(nesp))
	164	if (.not.allocated(surface_flux)) allocate(surface_flux(ngrid,nesp))
	165	if (.not.allocated(surface_flux2)) allocate(surface_flux2(ngrid,nesp))
	166	if (.not.allocated(escape)) allocate(escape(ngrid,nesp))
	167	if (.not.allocated(chemnoms)) allocate(chemnoms(nesp))
	168
	169	surface_flux(:,:) = 0.
	170	surface_flux2(:,:) = 0.
	171	escape(:,:) = 0.
	172	SF_mode(:) = 2
	173	SF_value(:) = 0.
	174	prod_rate(:) = 0.
	175	iter_3d(:,:) = 0.
	176	iter(:) = 0.
	177
	178	call ini_tracchim
	179
	180	! Sanity check mmol /= 0. in chemistry
	181	do iq = 1,nq
	182	if (is_chim(iq).eq.1 .and. mmol(iq).eq.0.) then
	183	write(,) 'Error in calchim:'
	184	write(,) 'Mmol cannot be equal to 0 for chemical species'
	185	stop
[1796]	186	end if
[2542]	187	end do
[1796]	188
	189	firstcall = .false.
	190	end if ! if (firstcall)
	191
	192	! Initializations
	193
	194	zycol(:,:) = 0.
	195	dqchim(:,:,:) = 0.
	196	dqschim(:,:) = 0.
	197
	198	!=======================================================================
	199	! loop over grid
	200	!=======================================================================
	201
	202	do ig = 1,ngrid
	203
	204	foundswitch = 0
	205	do l = 1,nlayer
[2542]	206	iesp = 0
	207	do iq = 1,nq
	208	if (is_chim(iq).eq.1) then
	209	iesp = iesp + 1
	210	zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
	211	xmmol(l,iesp) = mmol(iq)/mmean(ig,l)
	212	zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp)
	213	end if
[1796]	214	end do
[2542]	215
[1796]	216	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
	217	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
	218	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
	219	zpress(l) = pplay(ig,l)/100.
	220	ztemp(l) = zt(ig,l)
	221	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
	222	zlocal(l) = zzlay(ig,l)/1000.
	223	zmmean(l) = mmean(ig,l)
	224
	225	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
	226
	227	surfdust1d(l) = surfdust(ig,l)*1.e-2
	228	surfice1d(l) = surfice(ig,l)*1.e-2
	229
	230	! search for switch index between regions
	231
	232	! if (photochem .and. thermochem) then
	233	! if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
	234	! lswitch = l
	235	! foundswitch = 1
	236	! end if
	237	! end if
	238	if (.not. photochem) then
	239	lswitch = 22
	240	end if
	241	! if (.not. thermochem) then
	242	lswitch = min(50,nlayer+1)
	243	! end if
	244
	245	end do ! of do l=1,nlayer
	246
[2542]	247	szacol = acos(mu0(ig))*180./pi
	248	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
	249	fractcol = fract(ig)
[1796]	250
	251	!=======================================================================
	252	! call chemical subroutines
	253	!=======================================================================
	254
	255	! chemistry in lower atmosphere
	256
[2542]	257	! if (photochem) then
[1796]	258
[2542]	259	call photochemistry_asis(nlayer,ngrid, &
	260	ig,lswitch,zycol,szacol,fractcol,ptimestep,&
	261	zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
	262	surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:))
[1796]	263
[2542]	264	! diagnostic photochemical heating
	265	zdtchim_output(ig) = 0.
	266	do l = 1,nlayer
	267	dEzchim(ig,l) = zdtchim(ig,l)cpp(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
	268	zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)cpp(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
	269	end do
[1796]	270
	271	do l = 1,nlayer
	272	iter_3d(ig,l) = iter(l)
	273	end do
	274
	275	! chemistry in upper atmosphere
	276
	277	! if (thermochem) then
	278	! call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
	279	! zdens,zpress,zlocal,szacol,ptimestep,zday)
	280	! end if
	281
	282	! dry deposition
	283
[2542]	284	if (depos) then
	285	call deposition_source(ngrid, nlayer, nq, &
	286	ig, zzlay, zzlev, zdens, zycol, ptimestep)
	287	surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:)
	288	if (ngrid==1) then
	289	if(mod(icount,ecritphy).eq.0) then
	290	surface_flux2(ig,:) = surface_flux2(ig,:)/float(ecritphy)
	291	endif
	292	else
	293	if(mod(icount*iphysiq,ecritphy).eq.0) then
	294	surface_flux2(ig,:) = surface_flux2(ig,:)*iphysiq/float(ecritphy)
	295	endif
	296	endif
	297	end if
[1796]	298
	299	!=======================================================================
	300	! tendencies
	301	!=======================================================================
	302
	303	! index of the most abundant species at each level
	304
	305	! major(:) = maxloc(zycol, dim = 2)
	306
	307	! tendency for the most abundant species = - sum of others
	308
	309	do l = 1,nlayer
	310	iloc=maxloc(zycol(l,:))
	311	iqmax=iloc(1)
[2542]	312	do iq = 1,nesp
[1796]	313	if (iq /= iqmax) then
[2542]	314	dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep
[1796]	315	dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep)
	316	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq)
	317	end if
	318	end do
[2542]	319
[1796]	320	end do ! of do l = 1,nlayer
	321
	322
	323
	324
	325
	326
	327	!=======================================================================
	328	! end of loop over grid
	329	!=======================================================================
	330
[2542]	331	end do ! of do ig=1,ngridbidon(ig,:) = 1
[1796]	332
	333	!=======================================================================
	334	! write outputs
	335	!=======================================================================
	336
	337	! value of parameter 'output' to trigger writting of outputs
	338	! is set above at the declaration of the variable.
[2542]	339	if (photochem .and. output) then
[1796]	340
[2542]	341	if (callstats) then
	342	! photoloysis
	343	do i=1,nb_phot_hv_max
	344	call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1), &
	345	's-1',3,v_phot_3d(1,1,i))
	346	end do
	347	! iter
	348	call wstats(ngrid,'iter','iterations', &
[1796]	349	' ',3,iter_3d(1,1))
[2542]	350	! phothcemical heating
	351	call wstats(ngrid,'zdtchim','dT chim', &
	352	'K.s-1',3,zdtchim)
	353	call wstats(ngrid,'dEzchim','dE chim', &
	354	'W m-2',3,dEzchim)
	355	call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
	356	! Mean molecular mass
	357	call wstats(ngrid,'mmean','mean molecular mass', &
[1796]	358	'g.mole-1',3,mmean(1,1))
[2542]	359	! Chemical tendencies
	360	do iesp=1,nesp
	361	call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', &
	362	'kg/kg s^-1',3,dqchim(1,1,iesp) )
	363	end do
	364	if (depos) then
	365	! Surface fluxes
	366	do iesp=1,nesp
	367	call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', &
	368	'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
	369	call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', &
	370	'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
	371	end do
	372	endif ! end depos
	373	endif ! end callstats
	374
	375	! photoloysis
	376	do i=1,nb_phot_hv_max
	377	call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1), &
	378	's-1',3,v_phot_3d(1,1,i))
	379	end do
	380	! iter
	381	call writediagfi(ngrid,'iter','iterations', &
	382	' ',3,iter_3d(1,1))
	383	! phothcemical heating
	384	call writediagfi(ngrid,'zdtchim','dT chim', &
	385	'K.s-1',3,zdtchim)
	386	call writediagfi(ngrid,'dEzchim','dE chim', &
	387	'W m-2',3,dEzchim)
	388	call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
	389	! Mean molecular mass
	390	call writediagfi(ngrid,'mmean','mean molecular mass', &
	391	'g.mole-1',3,mmean(1,1))
	392	! Chemical tendencies
	393	do iesp=1,nesp
	394	call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', &
	395	'kg/kg s^-1',3,dqchim(1,1,iesp) )
	396	end do
	397	if (depos) then
	398	! Surface fluxes
	399	do iesp=1,nesp
	400	call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', &
	401	'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
	402	call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', &
	403	'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
	404	end do
	405	! Restart flux average
	406	if (ngrid == 1) then
	407	if(mod(icount,ecritphy).eq.0) then
	408	surface_flux2(:,:) = 0.0
	409	endif
	410	else
	411	if(mod(icount*iphysiq,ecritphy).eq.0) then
	412	surface_flux2(:,:) = 0.0
	413	endif
[1796]	414	endif
[2542]	415	endif ! end depos
	416
[1796]	417	end if ! of if (output)
	418
	419	return
[2542]	420
	421
	422	contains
	423
	424
	425	!======================================================================
	426
	427	subroutine ini_tracchim
	428
	429	!======================================================================
	430	! YJ Modern tracer.def
	431	! Get in tracer.def chemical parameters
	432	!======================================================================
	433
	434	use chimiedata_h
	435	use tracer_h, only: is_chim
	436	implicit none
	437
	438	! local
	439	character(len=500) :: tracline ! to store a line of text
	440	integer :: nq ! number of tracers
	441	integer :: iesp, iq
	442
	443	! Get nq
	444	open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
	445	if (ierr.eq.0) then
	446	READ(90,'(A)') tracline
	447	IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
	448	READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def
	449	ELSE
	450	DO
	451	READ(90,'(A)',iostat=ierr) tracline
	452	IF (ierr.eq.0) THEN
	453	IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
	454	READ(tracline,*,iostat=ierr) nq
	455	EXIT
	456	ENDIF
	457	ELSE ! If pb, or if reached EOF without having found number of tracer
	458	write(,) "calchim: error reading line of tracers"
	459	write(,) " (first lines of traceur.def) "
	460	stop
	461	ENDIF
	462	ENDDO
	463	ENDIF ! if modern or standard traceur.def
	464	else
	465	write(,) "calchim: error opening traceur.def"
	466	stop
	467	endif
	468
	469	! Get data of tracers
	470	iesp = 0
	471	do iq=1,nq
	472	read(90,'(A)') tracline
	473	if (is_chim(iq).eq.1) then
	474	iesp = iesp + 1
	475	read(tracline,*) chemnoms(iesp)
	476	write(,)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp))
	477	if (index(tracline,'SF_mode=' ) /= 0) then
	478	read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp)
	479	write(,) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp)
	480	else
	481	write(,) ' Parameter value (default) : SF_mode=', SF_mode(iesp)
	482	end if
	483	if (index(tracline,'SF_value=' ) /= 0) then
	484	read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp)
	485	write(,) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp)
	486	else
	487	write(,) ' Parameter value (default) : SF_value=', SF_value(iesp)
	488	end if
	489	if (index(tracline,'prod_rate=' ) /= 0) then
	490	read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp)
	491	write(,) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp)
	492	else
	493	write(,) ' Parameter value (default) : prod_rate=', prod_rate(iesp)
	494	end if
	495	end if
	496	end do
	497
	498	close(90)
	499
	500	end subroutine ini_tracchim
	501
	502	end subroutine calchim_asis
	503

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

Download in other formats:

Original Format