Context Navigation

calchim.F90 @ 3567

Last change on this file since 3567 was 3318, checked in by slebonnois, 8 months ago
Titan PCM update : optics + microphysics
File size: 23.0 KB

Rev	Line
[1908]	1	SUBROUTINE calchim(ngrid,qy_c,declin,dtchim, &
[2097]	2	ctemp,cpphi,cpphis,cplay,cplev,czlay,czlev,dqyc)
[1672]	3
[1908]	4	!---------------------------------------------------------------------------------------------------------
	5	!
	6	! Purpose : Interface subroutine to photochemical C model for Titan GCM.
	7	! -------
	8	! The subroutine computes the chemical processes for a single vertical column.
	9	!
	10	! - Only tendencies are returned.
	11	! - With moyzon_ch=.true. and input vectors zonally averaged
	12	! the calculation is done only once per lat. band
	13	!
	14	! Authors: + S. Lebonnois : 01/2000 \| 09/2003
	15	! ------- adaptation for Titan 3D : 02/2009
	16	! adaptation for // : 04/2013
	17	! extension chemistry up to 1300km : 10/2013
	18	!
	19	! + J. Vatant d'Ollone
[1946]	20	! + 02/17 - adaptation for the new generic-forked physics
[1908]	21	! + 01/18 - 03/18 - Major transformations :
	22	! - Upper chemistry fields are now stored in startfi
	23	! and defined on a pressure grid from Vervack profile
	24	! - These modifs enables to run chemistry with others resolution than 32x48x55 !
	25	! - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
	26	! - Chemistry can still be done in 2D
	27	! -> Calcul. once per band lat and put same tendency in all longi.
	28	! Check for negs in physiq_mod.
	29	! -> If procs sharing a lat band, no problem, the calcul will just be done twice.
	30	! -> Will not work with Dynamico, where the chemistry will have to be done in 3D.
	31	! ( and there'll be work to do to get rid of averaged fields )
	32	!
[2097]	33	! + 02/19 : To always have correct photodissociations rates, altitudes sent here by physiq are always
	34	! calculated with effective g - and with reference to the body not the local surface -
	35	! even if in physiq we keep altitudes coherent with dynamics !
	36	!
[2099]	37	! + STILL TO DO : + Replug the interaction with haze (cf titan.old) -> to see with JB.
	38	! + Use iso_c_binding for the fortran-C exchanges.
[1908]	39	!---------------------------------------------------------------------------------------------------------
	40
	41	! --------------------------------------------------------------------
	42	! Structure :
	43	! -----------
	44	! 0. Declarations
	45	! I. Init and firstcall
	46	! 1. Read and store Vervack profile
	47	! 2. Compute planetar averaged atm. properties
	48	! 3. Init compound caracteristics
	49	! 4. Init photodissociations rates from actinic fluxes
	50	! 5. Init chemical reactions
	51	! 6. Init eddy diffusion coeff
	52	! II. Loop on latitudes/grid-points
	53	! 0. Check on 2D chemistry
	54	! 1. Compute atm. properties at grid points
	55	! 2. Interpolate photodissociation rates at lat,alt,dec
	56	! 3. Read composition
	57	! 4. Call main solver gptitan C routine
	58	! 5. Calculate output tendencies on advected tracers
	59	! 6. Update upper chemistry fields
	60	! 0bis. If 2D chemsitry, don't recalculate if needed
	61	! -----------------------------------------------------------------
	62
[1956]	63	USE, INTRINSIC :: iso_c_binding
[3318]	64	use tracer_h
[1908]	65	USE comchem_h
	66	USE dimphy
	67	USE datafile_mod, ONLY: datadir
	68	USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
	69	USE geometry_mod, ONLY: latitude
[2368]	70	#ifndef MESOSCALE
[1908]	71	USE logic_mod, ONLY: moyzon_ch
[1947]	72	USE moyzon_mod, ONLY: tmoy, playmoy
[2368]	73	#endif
[1908]	74
	75	IMPLICIT NONE
	76
[1903]	77	! ------------------------------------------
	78	! ********* 0. Declarations ***********
	79	! ------------------------------------------
[1672]	80
[1908]	81	! Arguments
	82	! ---------
[1672]	83
[1908]	84	INTEGER, INTENT(IN) :: ngrid ! Number of atmospheric columns.
	85	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(IN) :: qy_c ! Chemical species on GCM layers after adv.+diss. (mol/mol).
	86	REAL*8, INTENT(IN) :: declin ! Solar declination (rad).
	87	REAL*8, INTENT(IN) :: dtchim ! Chemistry timsetep (s).
	88	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: ctemp ! Mid-layer temperature (K).
	89	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cpphi ! Mid-layer geopotential (m2.s-2).
[2097]	90	REAL*8, DIMENSION(ngrid), INTENT(IN) :: cpphis ! Surface geopotential (m2.s-2).
[1908]	91	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cplay ! Mid-layer pressure (Pa).
	92	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: cplev ! Inter-layer pressure (Pa).
[2097]	93	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: czlay ! Mid-layer effective altitude (m) : ref = geoid.
	94	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: czlev ! Inter-layer effective altitude (m) ref = geoid.
[1672]	95
[1908]	96	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT) :: dqyc ! Chemical species tendencies on GCM layers (mol/mol/s).
[1672]	97
[1908]	98	! Local variables :
	99	! -----------------
[1672]	100
[1908]	101	INTEGER :: i , l, ic, ig, igm1
[1672]	102
[1908]	103	INTEGER :: dec, idec, ipres, ialt, klat
[1672]	104
[1908]	105	REAL*8 :: declin_c ! Declination (deg).
	106	REAL*8 :: factp, factalt, factdec, factlat, krpddec, krpddecp1, krpddecm1
	107	REAL*8 :: temp1, logp
[1672]	108
[1908]	109	! Variables sent into chemistry module (must be in double precision)
	110	! ------------------------------------------------------------------
[1903]	111
[1908]	112	REAL*8, DIMENSION(nlaykim_tot) :: temp_c ! Temperature (K).
	113	REAL*8, DIMENSION(nlaykim_tot) :: press_c ! Pressure (Pa).
	114	REAL*8, DIMENSION(nlaykim_tot) :: phi_c ! Geopotential (m2.s-2) - actually not sent in chem. module but used to compute alts.
	115	REAL*8, DIMENSION(nlaykim_tot) :: nb ! Density (cm-3).
[1672]	116
[1908]	117	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy ! Chemical species in whole column (mol/mol) sent to chem. module.
	118	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy0 ! " " " " " " before modifs.
[1672]	119
[1908]	120	REAL*8 :: surfhaze(nlaykim_tot)
	121	REAL*8 :: cprodaer(nlaykim_tot,4), cmaer(nlaykim_tot,4)
	122	REAL*8 :: ccsn(nlaykim_tot,4), ccsh(nlaykim_tot,4)
[1672]	123
[1908]	124	REAL*8, DIMENSION(nlaykim_tot) :: rmil ! Mid-layer distance (km) to planetographic center.
	125	REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
	126	! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
[1672]	127
[1908]	128	! Saved variables initialized at firstcall
	129	! ----------------------------------------
[1672]	130
[1908]	131	LOGICAL, SAVE :: firstcall = .TRUE.
	132	!$OMP THREADPRIVATE(firstcall)
[1672]	133
[1908]	134	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: kedd ! Eddy mixing coefficient for upper chemistry (cm^2.s-1)
	135	!$OMP THREADPRIVATE(kedd)
[1672]	136
[2099]	137	REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: md ! Mean molecular diffusion coefficients (cm^2.s-1)
	138	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: mass ! Molar mass of the compounds (g.mol-1)
[1908]	139	!$OMP THREADPRIVATE(mass,md)
	140
	141	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: r1d, ct1d, p1d, t1d ! Vervack profile
	142	! JVO 18 : No threadprivate for those as they'll be read in tcp.ver by master
	143
	144	REAL*8, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: krpd ! Photodissociations rate table
	145	REAL*8, DIMENSION(:,:) , ALLOCATABLE, SAVE :: krate ! Reactions rate ( photo + chem )
	146	!$OMP THREADPRIVATE(krpd,krate)
	147
	148	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nom_prod, nom_perte
	149	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: reactif
	150	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: prod
	151	INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: perte
	152	!$OMP THREADPRIVATE(nom_prod,nom_perte,reactif,prod,perte)
	153
	154	! TEMPORARY : Dummy parameters without microphysics
	155	! Here to keep the whole stuff running without modif chem. module
	156	! ---------------------------------------------------------------
	157
	158	INTEGER :: utilaer(16)
	159	INTEGER :: aerprod = 0
	160	INTEGER :: htoh2 = 0
	161
	162	! -----------------------------------------------------------------------
	163	! *************** I. Initialisations and Firstcall ******************
	164	! -----------------------------------------------------------------------
	165
	166	IF (firstcall) THEN
	167
	168	PRINT*, 'CHIMIE, premier appel'
	169
[2326]	170	if (ngrid .eq. 1) then ! if 1D no dynamic mixing, we set the kedd in all column
	171	call check(nlaykim_tot,0,nlrt_kim,nkim)
	172	else
	173	call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
	174	endif
[1908]	175
	176	ALLOCATE(r1d(131))
	177	ALLOCATE(ct1d(131))
	178	ALLOCATE(p1d(131))
	179	ALLOCATE(t1d(131))
	180
	181	ALLOCATE(md(nlaykim_tot,nkim))
	182	ALLOCATE(mass(nkim))
	183
	184	ALLOCATE(kedd(nlaykim_tot))
	185
	186	ALLOCATE(krate(nlaykim_tot,nr_kim))
[1950]	187	ALLOCATE(krpd(nd_kim+1,nlrt_kim,15,nlat_actfluxes))
[1908]	188
	189	ALLOCATE(nom_prod(nkim))
	190	ALLOCATE(nom_perte(nkim))
	191	ALLOCATE(reactif(5,nr_kim))
	192	ALLOCATE(prod(200,nkim))
	193	ALLOCATE(perte(2,200,nkim))
	194
[2099]	195	! 0. Deal with characters for C-interoperability
	196	! ----------------------------------------------
	197	! NB ( JVO 19 ) : Using iso_c_binding would do things in an even cleaner way !
	198	DO ic=1,nkim
	199	nomqy_c(ic) = trim(cnames(ic))//char(0) ! Add the C null terminator
	200	ENDDO
	201	nomqy_c(nkim+1)="HV"//char(0) ! For photodissociations
	202
[1908]	203	! 1. Read Vervack profile "tcp.ver", once for all
	204	! -----------------------------------------------
	205
	206	!$OMP MASTER
	207	OPEN(11,FILE=TRIM(datadir)//'/tcp.ver',STATUS='old')
	208	READ(11,*)
	209	DO i=1,131
	210	READ(11,*) r1d(i), t1d(i), ct1d(i), p1d(i)
	211	! For debug :
	212	! -----------
	213	! PRINT*, "tcp.ver", r1d(i), t1d(i), ct1d(i), p1d(i)
	214	ENDDO
	215	CLOSE(11)
	216	!$OMP END MASTER
	217	!$OMP BARRIER
	218
	219	! 2. Calculation of temp_c, densities and altitudes in planetary average
	220	! ----------------------------------------------------------------------
	221
[1947]	222	! JVO18 : altitudes are no more calculated in firstcall, as I set kedd in pressure grid
[1908]	223
[2097]	224	! a. For GCM layers we just copy-paste ( assuming that physiq always send correct altitudes ! )
[1908]	225
	226	PRINT*,'Init chemistry : pressure, density, temperature ... :'
[1947]	227	PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'
[2368]	228
	229	#ifndef MESOSCALE
	230	!! JVO 20 : I put this CPP key because mesoscale cannot access tmoy and playmoy,
	231	! you should fix this when you will propeerly use chemistry !
	232
[1908]	233	IF (ngrid.NE.1) THEN
	234	DO l=1,klev
	235	temp_c(l) = tmoy(l) ! K
	236	press_c(l) = playmoy(l)/100. ! mbar
	237	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
[1979]	238	PRINT*, l, press_c(l), nb(l), temp_c(l)
[1908]	239	ENDDO
	240	ELSE
[2368]	241	#endif
[1908]	242	DO l=1,klev
	243	temp_c(l) = ctemp(1,l) ! K
	244	press_c(l) = cplay(1,l)/100. ! mbar
	245	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
[1947]	246	PRINT*, l, press_c(l), nb(l), temp_c(l)
[1908]	247	ENDDO
[2368]	248	#ifndef MESOSCALE
[1947]	249	ENDIF
[2368]	250	#endif
[1908]	251
	252	! b. Extension in upper atmosphere with Vervack profile
[2097]	253	! NB : Maybe the transition klev/klev+1 is harsh if T profile different from Vervack ...
[1908]	254
	255	ipres=1
	256	DO l=klev+1,nlaykim_tot
	257	press_c(l) = preskim(l-klev) / 100.0
	258	DO i=ipres,130
	259	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
	260	ipres=i
	261	ENDIF
[1672]	262	ENDDO
[1908]	263	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
[1672]	264
[1908]	265	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
	266	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
[1947]	267	PRINT*, l , press_c(l), nb(l), temp_c(l)
[1908]	268	ENDDO
[1672]	269
[1908]	270	! 3. Compounds caracteristics
	271	! ---------------------------
	272	mass(:) = 0.0
	273	call comp(nomqy_c,nb,temp_c,mass,md)
	274	PRINT*,' Mass'
	275	DO ic=1,nkim
	276	PRINT*, nomqy_c(ic), mass(ic)
	277	ENDDO
[1672]	278
[1908]	279	! 4. Photodissociation rates
	280	! --------------------------
	281	call disso(krpd,nlat_actfluxes)
[1672]	282
[1908]	283	! 5. Init. chemical reactions with planetary average T prof.
	284	! ----------------------------------------------------------
[1672]	285
[1908]	286	! NB : Chemical reactions rate are assumed to be constant within the T range of Titan's atm
	287	! so we fill their krate once for all but krate for photodiss will be filled at each timestep
	288
	289	call chimie(nomqy_c,nb,temp_c,krate,reactif, &
	290	nom_perte,nom_prod,perte,prod)
[1672]	291
[1908]	292	! 6. Eddy mixing coefficients (constant with time and space)
	293	! ----------------------------------------------------------
[2053]	294
[1908]	295	kedd(:) = 1.e3 ! Default value =/= zero
[1672]	296
[1946]	297	! NB : Eddy coeffs (e.g. Lavvas et al 08, Yelle et al 08) in altitude but they're rather linked to pressure
[1908]	298	! Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
[1672]	299
[2097]	300	!! First calculate kedd for upper chemistry layers
	301	!DO l=klev-4,nlaykim_tot
	302	! logp=-log10(press_c(l))
	303	!! 2E6 at 400 km ~ 10-2 mbar
	304	! IF ( logp.ge.2.0 .and. logp.le.3.0 ) THEN
	305	! kedd(l) = 2.e6 * 5.0**(logp-2.0)
	306	!! 1E7 at 500 km ~ 10-3 mbar
	307	! ELSE IF ( logp.ge.3.0 .and. logp.le.4.0 ) THEN
	308	! kedd(l) = 1.e7 * 3.0**(logp-3.0)
	309	!! 3E7 above 700 km ~ 10-4 mbar
	310	! ELSEIF ( logp.gt.4.0 ) THEN
	311	! kedd(l) = 3.e7
	312	! ENDIF
	313	!ENDDO
	314
	315	! Kedd from (E7) in Vuitton 2019
[2099]	316	if (ngrid .eq. 1) then ! if 1D no dynamic mixing, we set the kedd in all column
	317	DO l=1,nlaykim_tot
	318	kedd(l) = 300.0 * ( 1.0E2 / press_c(l) )*1.5 3.0E7 / &
	319	( 300.0 * ( 1.0E2 / press_c(l) )**1.5 + 3.0E7 )
	320	ENDDO
	321	else
	322	DO l=klev-4,nlaykim_tot
	323	! JVO 18 : We keep the nominal profile in the GCM 5 upper layers
	324	! to have a correct vertical mixing in the sponge layer
	325	kedd(l) = 300.0 * ( 1.0E2 / press_c(l) )*1.5 3.0E7 / &
[2097]	326	( 300.0 * ( 1.0E2 / press_c(l) )**1.5 + 3.0E7 )
[2099]	327	ENDDO
	328	endif
[1908]	329
[2099]	330	if (ngrid .gt. 1) then ! not in 1D, no dynamic mixing
	331	! Then adjust 10 layers profile fading to default value depending on kedd(ptop)
	332	DO l=klev-15,klev-5
	333	temp1 = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev-4)/press_c(klev-15)) )
	334	kedd(l) = 10.*( 3.0 + log10(kedd(klev-4)/1.e3) temp1 )
	335	ENDDO
	336	endif
[3318]	337
	338	! [BdBdT : On force la chimie... kedd nul dans le PCM]
	339	kedd(:klev) = 1.e3 ! Default value =/= zero
[2099]	340
[1908]	341	firstcall = .FALSE.
	342	ENDIF ! firstcall
	343
	344	declin_c = declin*180./pi
[1672]	345
[1908]	346	! -----------------------------------------------------------------------
	347	! ********************* II. Loop on latitudes ***********************
	348	! -----------------------------------------------------------------------
	349
	350	DO ig=1,ngrid
[1672]	351
[1908]	352	IF (ig.eq.1) THEN
	353	igm1=1
	354	ELSE
	355	igm1=ig-1
	356	ENDIF
[1672]	357
[1908]	358	! If 2D chemistry, trick to do the calculation only once per latitude band within the chunk
	359	! NB1 : Will be obsolete with DYNAMICO, the chemistry will necessarly be 3D
[2045]	360	! NB2 : Test of same latitude with dlat=0.1 : I think that if you run sims better than 1/10th degree then
[1908]	361	! either it's with Dynamico and doesn't apply OR it is more than enough in terms of "preco / calc time" !
	362	! -------------------------------------------------------------------------------------------------------
[1672]	363
[2368]	364	#ifndef MESOSCALE
[2099]	365	IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.0)) ) .OR. (.NOT. moyzon_ch) ) THEN
[2368]	366	#endif
[1908]	367	! 1. Compute altitude for the grid point with hydrostat. equilib.
	368	! ---------------------------------------------------------------
[1672]	369
[1908]	370	! a. For GCM layers we just copy-paste
[2097]	371	! JVO 19 : Now physiq always sent correct altitudes with effective g for chemistry ( even if it's not the case in physiq )
[1672]	372
[1908]	373	DO l=1,klev
	374	rinter(l) = (czlev(ig,l)+rad)/1000.0 ! km
	375	rmil(l) = (czlay(ig,l)+rad)/1000.0 ! km
	376	temp_c(l) = ctemp(ig,l) ! K
	377	phi_c(l) = cpphi(ig,l) ! m2.s-2
	378	press_c(l) = cplay(ig,l)/100. ! mbar
	379	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
	380	ENDDO
[2368]	381	rinter(klev+1)=( (czlay(ig,klev) + ( czlay(ig,klev) - czlev(ig,klev))) +rad )/1000. ! You shall not use czlev(zlev+1) whose value is 1.e7 !
[1672]	382
[1908]	383	! b. Extension in upper atmosphere with Vervack profile
[1672]	384
[1908]	385	ipres=1
	386	DO l=klev+1,nlaykim_tot
	387	press_c(l) = preskim(l-klev) / 100.0
	388	DO i=ipres,130
	389	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
	390	ipres=i
	391	ENDIF
	392	ENDDO
	393	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
[1672]	394
[1908]	395	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
	396	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
	397	ENDDO
	398
	399	! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
	400	! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
[1672]	401
[1908]	402	DO l=klev+1,nlaykim_tot
[1672]	403
[2097]	404	! Compute geopotential on the upper grid with effective g to have correct altitudes
	405	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1908]	406
	407	temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
	408	phi_c(l) = phi_c(l-1) + rtemp1log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
	409
[2097]	410	rmil(l) = ( gradrad / (g*rad - ( phi_c(l) + cpphis(ig) ) ) ) / 1000.0 ! z(phi) with g varying with altitude with reference to the geoid
[1908]	411	ENDDO
	412
	413	DO l=klev+2,nlaykim_tot
[2097]	414	rinter(l) = 0.5*(rmil(l-1) + rmil(l)) ! should be balanced with the intermediate pressure rather than 0.5
[1908]	415	ENDDO
	416
	417	! 2. From krpd, compute krate for dissociations (declination-latitude-altitude interpolation)
	418	! -------------------------------------------------------------------------------------------
	419
	420	! a. Calculate declination dependence
	421
	422	if ((declin_c*10+267).lt.14.) then
	423	idec = 0
	424	dec = 0
	425	else
	426	if ((declin_c*10+267).gt.520.) then
	427	idec = 14
	428	dec = 534
	429	else
	430	idec = 1
	431	dec = 27
	432	do while( (declin_c*10+267).ge.real(dec+20) )
	433	dec = dec+40
	434	idec = idec+1
	435	enddo
	436	endif
[1672]	437	endif
[1908]	438	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
	439	factdec = ( declin_c - (dec-267)/10. ) / 4.
[1672]	440	else
[1908]	441	factdec = ( declin_c - (dec-267)/10. ) / 2.7
[1672]	442	endif
	443
[1908]	444	! b. Calculate klat for interpolation on fixed latitudes of actinic fluxes input
[1672]	445
[1908]	446	klat=1
	447	DO i=1,nlat_actfluxes
	448	IF (latitude(ig).LT.lat_actfluxes(i)) klat=i
	449	ENDDO
	450	IF (klat==nlat_actfluxes) THEN ! avoid rounding problems
	451	klat = nlat_actfluxes-1
	452	factlat = 1.0
	453	ELSE
	454	factlat = (latitude(ig)-lat_actfluxes(klat))/(lat_actfluxes(klat+1)-lat_actfluxes(klat))
	455	ENDIF
[1672]	456
[1908]	457	! c. Altitude loop
[1672]	458
[1908]	459	DO l=1,nlaykim_tot
[1672]	460
[1908]	461	! Calculate ialt for interpolation in altitude (krpd every 2 km)
	462	ialt = int((rmil(l)-rad/1000.)/2.)+1
	463	factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
[1672]	464
[1908]	465	! Altitude can go above top limit of UV levels - in this case we keep the 1310km top fluxes
	466	IF (ialt.GT.nlrt_kim-1) THEN
	467	ialt = nlrt_kim-1 ! avoid out-of-bound array
	468	factalt = 1.0
	469	ENDIF
[1672]	470
[1924]	471	DO i=1,nd_kim+1 ! nd_kim+1 is dissociation of N2 by GCR
[1672]	472
[1950]	473	krpddec = ( krpd(i,ialt ,idec+1,klat) * (1.0-factalt) &
	474	+ krpd(i,ialt+1,idec+1,klat) * factalt ) * (1.0-factlat) &
	475	+ ( krpd(i,ialt ,idec+1,klat+1) * (1.0-factalt) &
	476	+ krpd(i,ialt+1,idec+1,klat+1) * factalt ) * factlat
[1908]	477
[1924]	478	if ( factdec.lt.0. ) then
[1950]	479	krpddecm1 = ( krpd(i,ialt ,idec ,klat) * (1.0-factalt) &
	480	+ krpd(i,ialt+1,idec ,klat) * factalt ) * (1.0-factlat) &
	481	+ ( krpd(i,ialt ,idec ,klat+1) * (1.0-factalt) &
	482	+ krpd(i,ialt+1,idec ,klat+1) * factalt ) * factlat
[1908]	483	krate(l,i) = krpddecm1 * abs(factdec) + krpddec * ( 1.0 + factdec)
[1924]	484	else if ( factdec.gt.0. ) then
[1950]	485	krpddecp1 = ( krpd(i,ialt ,idec+2,klat) * (1.0-factalt) &
	486	+ krpd(i,ialt+1,idec+2,klat) * factalt ) * (1.0-factlat) &
	487	+ ( krpd(i,ialt ,idec+2,klat+1) * (1.0-factalt) &
	488	+ krpd(i,ialt+1,idec+2,klat+1) * factalt ) * factlat
[1908]	489	krate(l,i) = krpddecp1 * factdec + krpddec * ( 1.0 - factdec)
[1924]	490	else if ( factdec.eq.0. ) then
	491	krate(l,i) = krpddec
[1908]	492	endif
	493
	494	ENDDO ! i=1,nd_kim+1
	495	ENDDO ! l=1,nlaykim_tot
	496
	497	! 3. Read composition
	498	! -------------------
	499
	500	DO ic=1,nkim
	501	DO l=1,klev
	502	cqy(l,ic) = qy_c(ig,l,ic) ! advected tracers for the GCM part converted to molar frac.
	503	ENDDO
	504
	505	DO l=1,nlaykim_up
	506	cqy(klev+l,ic) = ykim_up(ic,ig,l) ! ykim_up for the upper atm.
	507	ENDDO
	508	ENDDO
	509
	510	cqy0(:,:) = cqy(:,:) ! Stores compo. before modifs
	511
	512	! 4. Call main Titan chemistry C routine
	513	! --------------------------------------
	514
	515	call gptitan(rinter,temp_c,nb, &
	516	nomqy_c,cqy, &
	517	dtchim,latitude(ig)*180./pi,mass,md, &
[1958]	518	kedd,krate,reactif, &
[1908]	519	nom_prod,nom_perte,prod,perte, &
	520	aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, &
	521	htoh2,surfhaze)
	522
	523	! 5. Calculates tendencies on composition for advected tracers
	524	! ------------------------------------------------------------
	525	DO ic=1,nkim
	526	DO l=1,klev
[3318]	527	dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
[1908]	528	ENDDO
	529	ENDDO
	530
	531	! 6. Update ykim_up
	532	! -----------------
	533	DO ic=1,nkim
	534	DO l=1,nlaykim_up
	535	ykim_up(ic,ig,l) = cqy(klev+l,ic)
	536	ENDDO
	537	ENDDO
	538	! NB: The full vertical composition grid will be created only for the outputs
	539
[2368]	540	#ifndef MESOSCALE
[1956]	541	ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
[1908]	542	dqyc(ig,:,:) = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
	543	ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations
	544	ENDIF
[2368]	545	#endif
[1908]	546
	547	ENDDO
	548
	549	END SUBROUTINE calchim

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

Download in other formats:

Original Format