Context Navigation

calchim.F @ 1243

Last change on this file since 1243 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 18.3 KB

Rev	Line
[1056]	1	SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
[1126]	2	. ctemp,cplay,cplev,czlay,czlev,
[3]	3	. dqyc)
	4
	5	c-------------------------------------------------
	6	c
[102]	7	c Introduction d une routine chimique
[3]	8	c
	9	c Auteur: S. Lebonnois, 01/2000 \| 09/2003
	10	c adaptation pour Titan 3D: 02/2009
[1056]	11	c adaptation pour // : 04/2013
[1126]	12	c extension chimie jusqu a 1300km : 10/2013
	13	c
[3]	14	c-------------------------------------------------
	15	c
[102]	16	use dimphy
[1056]	17	use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
[1126]	18	. NLEV,NLRT,NC,ND,NR
[1056]	19	USE comgeomphy, only: rlatd
[1126]	20	USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
[102]	21	implicit none
[3]	22	#include "dimensions.h"
	23	#include "clesphys.h"
	24	#include "paramet.h"
	25	#include "YOMCST.h"
	26
	27	c Arguments
	28	c ---------
	29
[1056]	30	INTEGER nlon ! nb of horiz points
[3]	31	INTEGER ny ! nb de composes (nqmax-nmicro)
[1056]	32	REAL qy_c(nlon,klev,NC) ! Especes chimiques apres adv.+diss.
[3]	33	character*10 nomqy_c(NC+1) ! Noms des especes chimiques
	34	REAL declin_rad,ls_rad ! declinaison et long solaire en radians
	35	REAL dtchim ! pas de temps chimie
[1056]	36	REAL ctemp(nlon,klev) ! Temperature
	37	REAL cplay(nlon,klev) ! pression (Pa)
[1126]	38	REAL cplev(nlon,klev+1) ! pression intercouches (Pa)
	39	REAL czlay(nlon,klev) ! altitude (m)
	40	REAL czlev(nlon,klev+1) ! altitude intercouches (m)
[3]	41
[1056]	42	REAL dqyc(nlon,klev,NC) ! Tendances especes chimiques
[3]	43
	44	c Local variables :
	45	c -----------------
[1056]	46
	47	integer i,j,l,ic,jm1
	48
[3]	49	c variables envoyees dans la chimie: double precision
	50
[1126]	51	REAL temp_c(NLEV)
	52	REAL press_c(NLEV) ! T,p(mbar) a 1 lat donnee
	53	REAL cqy(NLEV,NC) ! frac mol qui seront modifiees
	54	REAL cqy0(NLEV,NC) ! frac mol avant modif
	55	REAL surfhaze(NLEV)
	56	REAL cprodaer(NLEV,4),cmaer(NLEV,4)
	57	REAL ccsn(NLEV,4),ccsh(NLEV,4)
	58	! rmil: milieu de couche, grille pour K,D,p,ct,T,y
	59	! rinter: intercouche (grille RA dans la chimie)
	60	REAL rmil(NLEV),rinter(NLEV),nb(NLEV)
	61	REAL,save :: kedd(NLEV)
	62
[3]	63	REAL a,b
	64	character str11,str22
[1126]	65	REAL latit
[3]	66	character*20 formt1,formt2,formt3
	67
	68	c variables locales initialisees au premier appel
	69
	70	LOGICAL firstcal
	71	DATA firstcal/.true./
	72	SAVE firstcal
	73
[1126]	74	integer dec,declinint,ialt
	75	REAL declin_c ! declinaison en degres
	76	real factalt,factdec,krpddec,krpddecp1,krpddecm1
	77	real duree
	78	REAL,save :: mass(NC)
	79	REAL,save,allocatable :: md(:,:)
	80	REAL,save :: botCH4
[3]	81	DATA botCH4/0.05/
[1126]	82	real,save :: r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
[104]	83	REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
[1126]	84	integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
	85	integer,save :: prod(200,NC),perte(2,200,NC)
	86	character dummy30,fich7,ficdec15,curdec15,name*10
	87	real ficalt,oldq,newq,zalt
	88	logical okfic
	89
[3]	90	c-----------------------------------------------------------------------
	91	c***********************************************************************
	92	c
	93	c Initialisations :
	94	c ----------------
	95
	96	duree = dtchim ! passage en real*4 pour appel a routines C
	97
	98	IF (firstcal) THEN
	99
	100	print*,'CHIMIE, premier appel'
	101
	102	c ************************************
	103	c Au premier appel, initialisation de toutes les variables
	104	c pour les routines de la chimie.
	105	c ************************************
	106
[1126]	107	allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
[104]	108
[1126]	109	c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
[3]	110	c ----------------------------------
	111
	112	if (ny.ne.NC) then
	113	print*,'PROBLEME de coherence dans le nombre de composes:'
	114	. ,ny,NC
	115	STOP
	116	endif
	117
[1126]	118	c calcul de temp_c, densites et press_c en moyenne planetaire :
	119	c --------------------------------------------------------------
[3]	120
[1126]	121	print*,'pression, densites et temp (init chimie):'
[3]	122	print*,'level, press_c, nb, temp_c'
	123	DO l=1,klev
[1126]	124	rinter(l) = (zlevmoy(l)+RA)/1000.
	125	rmil(l) = (zlaymoy(l)+RA)/1000.
[3]	126	c temp_c (K):
[1126]	127	temp_c(l) = tmoy(l)
[3]	128	c press_c (mbar):
[1126]	129	press_c(l) = playmoy(l)/100.
[3]	130	c nb (cm-3):
	131	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
[1126]	132	print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
[3]	133	ENDDO
[1126]	134
	135	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
	136	do l=klev+1,NLEV
	137	rinter(l) = rinter(klev)
	138	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
	139	rmil(l) = rmil(klev)
	140	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
	141	enddo
	142
	143	c lecture de tcp.ver, une seule fois
	144	c remplissage r1d,t1d,ct1d,p1d
	145	open(11,file='../../INPUT/tcp.ver',status='old')
	146	read(11,*) dummy
	147	do i=1,131
	148	read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
	149	c print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
	150	enddo
	151	close(11)
	152
	153	c extension pour klev+1 a NLEV avec tcp.ver
	154	c la jonction klev/klev+1 est brutale... Tant pis ?
	155	ialt=1
	156	do l=klev+1,NLEV
	157	zalt=rmil(l)-RA/1000.
	158	do i=ialt,130
	159	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
	160	ialt=i
	161	endif
	162	enddo
	163	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
	164	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
	165	& + log(p1d(ialt+1)) * factalt )
	166	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
	167	& + log(ct1d(ialt+1)) * factalt )
	168	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
	169	print*,l,zalt,press_c(l),nb(l),temp_c(l)
	170	enddo
[3]	171
	172	c caracteristiques des composes:
	173	c -----------------------------
[1126]	174	mass(:) = 0.0
	175	call comp(nomqy_c,nb,temp_c,mass,md)
[3]	176	print*,' Mass'
	177	do ic=1,NC
	178	print*,nomqy_c(ic),mass(ic)
[1126]	179	c if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
[3]	180	enddo
	181
	182
[1126]	183	c Stockage des composes utilises dans la prod d aerosols
[3]	184	c (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
	185	c ! decalage de 1 car utilise dans le c !
	186
	187	do ic=1,NC
	188
	189	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	190	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
	191	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	192	c if (nomqy_c(ic).eq."CH4") then
	193	c do l=1,llm
	194	c do j=1,ip1jmp1
	195	c if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015
	196	c enddo
	197	c enddo
	198	c endif
	199	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	200
	201	if (nomqy_c(ic).eq."C4H2") then
	202	utilaer(10) = ic-1
	203	endif
	204	if (nomqy_c(ic).eq."HCN") then
	205	utilaer(6) = ic-1
	206	endif
	207	if (nomqy_c(ic).eq."HC3N") then
	208	utilaer(7) = ic-1
	209	endif
	210	if (nomqy_c(ic).eq."NCCN") then
	211	utilaer(14) = ic-1
	212	endif
	213	if (nomqy_c(ic).eq."CH3CN") then
	214	utilaer(15) = ic-1
	215	utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle
	216	endif
	217	if (nomqy_c(ic).eq."H") then
	218	utilaer(1) = ic-1
	219	endif
	220	if (nomqy_c(ic).eq."H2") then
	221	utilaer(2) = ic-1
	222	endif
	223	if (nomqy_c(ic).eq."C2H2") then
	224	utilaer(3) = ic-1
	225	endif
	226	if (nomqy_c(ic).eq."AC6H6") then
	227	utilaer(13) = ic-1
	228	endif
	229	if (nomqy_c(ic).eq."C2H3CN") then
	230	utilaer(16) = ic-1
	231	endif
	232	if (nomqy_c(ic).eq."C2") then
	233	utilaer(4) = ic-1
	234	endif
	235	if (nomqy_c(ic).eq."C2H") then
	236	utilaer(5) = ic-1
	237	endif
	238	if (nomqy_c(ic).eq."C3N") then
	239	utilaer(8) = ic-1
	240	endif
	241	if (nomqy_c(ic).eq."H2CN") then
	242	utilaer(9) = ic-1
	243	endif
	244	if (nomqy_c(ic).eq."C4H3") then
	245	utilaer(11) = ic-1
	246	endif
	247	if (nomqy_c(ic).eq."AC6H5") then
	248	utilaer(12) = ic-1
	249	endif
	250
	251	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	252	c if ((nomqy_c(ic).eq."HC3N").or.
	253	c $ (nomqy_c(ic).eq."C3N")) then
	254	c DO j=1,ip1jmp1
	255	c do l=1,34 ! p>~ 1 mbar
	256	c qy_c(j,l,ic) = 1.e-30
	257	c enddo
	258	c ENDDO
	259	c endif
	260	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	261
	262	enddo
	263
	264	c taux de photodissociations:
	265	c --------------------------
	266	call disso(krpd,jjp1)
	267
	268	c reactions chimiques:
	269	c -------------------
	270	call chimie(nomqy_c,nb,temp_c,krate,reactif,
	271	. nom_perte,nom_prod,perte,prod)
	272	c print*,'nom_prod, nom_perte:'
	273	c do ic=1,NC
	274	c print*,nom_prod(ic),nom_perte(ic)
	275	c enddo
	276	c print*,'premieres prod, perte(1:reaction,2:compagnon):'
	277	c do ic=1,NC
	278	c print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
	279	c enddo
	280
	281	c l = klev-3
	282	c print*,'krate a p=',press_c(l),' reactifs et produits:'
	283	c do ic=1,ND+1
	284	c print,ic,krpd(7,ic,l,4)nb(l)," ",
	285	c . nomqy_c(reactif(1,ic)+1),
	286	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
	287	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
	288	c enddo
	289	c do ic=ND+2,NR
	290	c print*,ic,krate(l,ic)," ",
	291	c . nomqy_c(reactif(1,ic)+1),
	292	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
	293	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
	294	c enddo
	295
[1126]	296
	297	c init kedd
	298	c ---------
	299	c kedd stays constant with time and space
	300	c => linked to pressure rather than altitude...
	301
	302	kedd(:) = 5e2 ! valeur mise par defaut
	303	! UNITE ? doit etre ok pour gptitan
	304	do l=1,NLEV
	305	zalt=rmil(l)-RA/1000. ! z en km
	306	if ((zalt.ge.250.).and.(zalt.le.400.)) then
	307	kedd(l) = 10.**(3.+(zalt-250.)/50.)
	308	! 1E3 at 250 km
	309	! 1E6 at 400 km
	310	elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
	311	kedd(l) = 10.*(6.+1.3(zalt-400.)/200.)
	312	! 2E7 at 600 km
	313	elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
	314	kedd(l) = 10.*(7.3+0.7(zalt-600.)/300.)
	315	! 1E8 above 900 km
	316	elseif ( zalt.gt.900. ) then
	317	kedd(l) = 1.e8
	318	endif
	319	enddo
	320
[3]	321	ENDIF ! premier appel
	322
	323	c***********************************************************************
	324	c-----------------------------------------------------------------------
	325
	326	c calcul declin_c (en degres)
	327	c ---------------------------
	328
	329	declin_c = declin_rad*180./RPI
	330	c print*,'declinaison en degre=',declin_c
	331
	332	c***********************************************************************
	333	c***********************************************************************
	334	c
	335	c BOUCLE SUR LES LATITUDES
	336	c
[1056]	337	DO j=1,nlon
[3]	338
[1056]	339	if (j.eq.1) then
	340	jm1=1
	341	else
	342	jm1=j-1
	343	endif
	344
	345	if((j.eq.1).or.(klat(j).ne.klat(jm1))) then
	346
[3]	347	c***********************************************************************
	348	c***********************************************************************
	349
	350	c-----------------------------------------------------------------------
	351	c
[1126]	352	c Distance radiale (intercouches, en km)
	353	c ----------------------------------------
	354
	355	do l=1,klev
	356	rinter(l) = (RA+czlev(j,l))/1000.
	357	rmil(l) = (RA+czlay(j,l))/1000.
	358	c print*,rinter(l)
	359	enddo
	360
	361	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
	362	do l=klev+1,NLEV
	363	rinter(l) = rinter(klev)
	364	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
	365	rmil(l) = rmil(klev)
	366	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
	367	enddo
	368
	369	c-----------------------------------------------------------------------
	370	c
[3]	371	c Temperature, pression (mbar), densite (cm-3)
	372	c -------------------------------------------
	373
	374	DO l=1,klev
	375	c temp_c (K):
	376	temp_c(l) = ctemp(j,l)
	377	c press_c (mbar):
	378	press_c(l) = cplay(j,l)/100.
	379	c nb (cm-3):
	380	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
	381	ENDDO
[1126]	382	c extension pour klev+1 a NLEV avec tcp.ver
	383	ialt=1
	384	do l=klev+1,NLEV
	385	zalt=rmil(l)-RA/1000.
	386	do i=ialt,130
	387	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
	388	ialt=i
	389	endif
	390	enddo
	391	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
	392	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
	393	& + log(p1d(ialt+1)) * factalt )
	394	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
	395	& + log(ct1d(ialt+1)) * factalt )
	396	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
	397	enddo
	398
[3]	399	c-----------------------------------------------------------------------
	400	c
[1126]	401	c passage krpd => krate
	402	c ---------------------
	403	c initialisation krate pour dissociations
[3]	404
[1126]	405	if ((declin_c*10+267).lt.14.) then
	406	declinint = 0
	407	dec = 0
	408	else
	409	if ((declin_c*10+267).gt.520.) then
	410	declinint = 14
	411	dec = 534
	412	else
	413	declinint = 1
	414	dec = 27
	415	do while( (declin_c*10+267).ge.real(dec+20) )
	416	dec = dec+40
	417	declinint = declinint+1
	418	enddo
	419	endif
	420	endif
	421	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
	422	factdec = ( declin_c - (dec-267)/10. ) / 4.
	423	else
	424	factdec = ( declin_c - (dec-267)/10. ) / 2.7
	425	endif
[3]	426
[1126]	427	do l=1,NLEV
[3]	428
[1126]	429	c INTERPOL EN ALT POUR k (krpd tous les deux km)
	430	ialt = int((rmil(l)-RA/1000.)/2.)+1
	431	factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
[3]	432
[1126]	433	do i=1,ND+1
	434	krpddecm1 = krpd(declinint ,i,ialt ,klat(j)) * (1-factalt)
	435	& + krpd(declinint ,i,ialt+1,klat(j)) * factalt
	436	krpddec = krpd(declinint+1,i,ialt ,klat(j)) * (1-factalt)
	437	& + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
	438	krpddecp1 = krpd(declinint+2,i,ialt ,klat(j)) * (1-factalt)
	439	& + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
[3]	440
[1126]	441	! ND+1 correspond a la dissociation de N2 par les GCR
	442	if ( factdec.lt.0. ) then
	443	krate(l,i) = krpddecm1 * abs(factdec)
	444	& + krpddec * ( 1 + factdec)
	445	endif
	446	if ( factdec.gt.0. ) then
	447	krate(l,i) = krpddecp1 * factdec
	448	& + krpddec * ( 1 - factdec)
	449	endif
	450	if ( factdec.eq.0. ) krate(l,i) = krpddec
	451	enddo
	452	enddo
	453
[3]	454	c-----------------------------------------------------------------------
	455	c
	456	c composition
	457	c ------------
[1126]	458
[3]	459	do ic=1,NC
	460	do l=1,klev
[1126]	461	cqy(l,ic) = qy_c(j,l,ic)
	462	cqy0(l,ic) = cqy(l,ic)
[3]	463	enddo
	464	enddo
[1126]	465
	466	c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
	467
	468	WRITE(str2,'(i2.2)') klat(j)
	469	fich = "comp_"//str2//".dat"
	470	inquire (file=fich,exist=okfic)
	471	if (okfic) then
	472	open(11,file=fich,status='old')
	473	c premiere ligne=declin
	474	read(11,'(A15)') ficdec
	475	write(curdec,'(E15.5)') declin_c
	476	c si la declin est la meme: ce fichier a deja ete reecrit
	477	c on lit la colonne t-1 au lieu de la colonne t
	478	c (cas d une bande de latitude partagee par 2 procs)
	479	do ic=1,NC
	480	read(11,'(A10)') name
	481	if (name.ne.nomqy_c(ic)) then
	482	print*,"probleme lecture ",fich
	483	print*,name,nomqy_c(ic)
	484	stop
	485	endif
	486	if (ficdec.eq.curdec) then
	487	do l=klev+1,NLEV
	488	read(11,*) ficalt,cqy(l,ic),newq
	489	enddo
	490	else
	491	do l=klev+1,NLEV
	492	read(11,*) ficalt,oldq,cqy(l,ic)
	493	enddo
	494	endif
	495	enddo
	496	close(11)
	497	else ! le fichier n'est pas la
	498	do ic=1,NC
	499	do l=klev+1,NLEV
	500	cqy(l,ic)=cqy(klev,ic)
	501	enddo
	502	enddo
	503	endif
	504	cqy0 = cqy
	505
[3]	506	c-----------------------------------------------------------------------
	507	c
	508	c total haze area (um2/cm3)
	509	c -------------------------
	510
	511	if (htoh2.eq.1) then
	512	do l=1,klev
	513	! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!!
	514	surfhaze(l) = psurfhaze(j,klev+1-l)
	515	c if (j.eq.25)
	516	c . print*,'psurfhaze en um2/cm3:',surfhaze(l)
	517	enddo
	518	endif
	519
	520	c-----------------------------------------------------------------------
	521	c
	522	c Appel de chimietitan
	523	c --------------------
	524
[1056]	525	call gptitan(rinter,temp_c,nb,
[1126]	526	$ nomqy_c,cqy,
	527	$ duree,(klat(j)-1),mass,md,
	528	$ kedd,botCH4,krate,reactif,
[3]	529	$ nom_prod,nom_perte,prod,perte,
	530	$ aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
	531	$ htoh2,surfhaze)
	532
	533	c Tendances composition
	534	c ---------------------
	535
	536	do ic=1,NC
	537	do l=1,klev
[1126]	538	dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! en /s
[3]	539	enddo
	540	enddo
	541
	542	c-----------------------------------------------------------------------
	543	c
	544	c production aer
	545	c --------------
	546
	547	if (aerprod.eq.1) then
	548
	549	do ic=1,4
	550	do l=1,klev
	551	prodaer(j,l,ic) = cprodaer(l,ic)
	552	maer(j,l,ic) = cmaer(l,ic)
	553	csn(j,l,ic) = ccsn(l,ic)
	554	csh(j,l,ic) = ccsh(l,ic)
	555	enddo
	556	enddo
	557
	558	endif
	559
[1126]	560	c-----------------------------------------------------------------------
	561	c
	562	c sauvegarde compo sur NLEV
	563	c -------------------------
	564
	565	c dans fichier compo_klat(j) (01 à 49)
	566
	567	open(11,file=fich,status='replace')
	568	c premiere ligne=declin
	569	write(11,'(E15.5)') declin_c
	570	do ic=1,NC
	571	write(11,'(A10)') nomqy_c(ic)
	572	do l=klev+1,NLEV
	573	write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
	574	. cqy0(l,ic),cqy(l,ic)
	575	enddo
	576	enddo
	577	close(11)
	578
[3]	579	c***********************************************************************
	580	c***********************************************************************
[1056]	581
[3]	582	c FIN: BOUCLE SUR LES LATITUDES
[1056]	583
	584	else ! same latitude, we don't do calculations again
	585	dqyc(j,:,:) = dqyc(jm1,:,:)
	586	if (aerprod.eq.1) then
	587	prodaer(j,:,:) = prodaer(jm1,:,:)
	588	maer(j,:,:) = maer(jm1,:,:)
	589	csn(j,:,:) = csn(jm1,:,:)
	590	csh(j,:,:) = csh(jm1,:,:)
	591	endif
	592	endif
	593
[3]	594	ENDDO
	595
	596	c***********************************************************************
	597	c***********************************************************************
	598
	599
	600	firstcal = .false.
	601	RETURN
	602	END

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

Download in other formats:

Original Format