Context Navigation

phytrac.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: 25.6 KB

Rev	Line
[3]	1	SUBROUTINE phytrac (firstcall,lastcall,
	2	. nqmax,nmicro,ptimestep,appkim,dtkim,
	3	. pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
[175]	4	. lonsol,
	5	. pu,pv,pzlev,pzlay,ftsol,
	6	. tr_seri,qaer,d_tr_mph,d_tr_kim,
	7	. fclat,reservoir)
[3]	8
	9	c======================================================================
	10	c S. Lebonnois, mai 2008
	11	c
	12	c Arguments:
	13	c
	14	c firstcall----input-L-variable logique indiquant le premier passage
	15	c lastcall-----input-L-variable logique indiquant le dernier passage
	16	c nqmax--------input-I-nombre de traceurs (total)
	17	c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
[1056]	18	c ptimestep----input-R-pas d integration pour la physique (seconde)
[3]	19	c appkim-------input-I-appel a la chimie
	20	c dtkim--------input-R-pas de temps chimique (seconde)
	21	c pplev--------input-R-pression pour chaque inter-couche (en Pa)
	22	c pplay--------input-R-pression pour chaque couche (en Pa)
	23	c delp---------input-R-epaisseur d'une couche (en Pa)
	24	c ptemp--------input-R-temperature (K)
	25	c pmu0---------input-R-cos angle zenithal
	26	c pfract-------input-R-fractional day
	27	c pdecli-------input-R-declinaison en radian
	28	c lonsol-------input-R-longitude solaire en radian
[175]	29	c pu-----------input-R-vitesse dans la direction X (de O a E) en m/s (1ere couche)
	30	c pv-----------input-R-vitesse Y (de S a N) en m/s (1ere couche)
	31	c pzlev--------input-R-altitude pour chaque inter-couche (en m)
	32	c pzlay--------input-R-altitude pour chaque couche (en m)
	33	c ftsol--------input-R-temperature au sol (en K)
[3]	34	c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
	35	c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
	36	c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
[175]	37	c fclat--------output-R-flux de chaleur latente d'evaporation du reservoir CH4 (J/m2/s)
	38	c reservoir----outpur-R-un reservoir de surface !!! (m)
[3]	39	c======================================================================
[102]	40	USE infotrac
	41	use dimphy
[1056]	42	USE common_mod, only: rmcbar,xfbar,ncount,
	43	& flxesp_i,tau_drop,tau_aer,solesp,precip,
	44	& evapch4,occcld_m,occcld,satch4,satc2h6,satc2h2,rmcloud
	45	USE moyzon_mod
[1072]	46	USE write_field_phy
[102]	47	IMPLICIT none
[3]	48	#include "dimensions.h"
	49	#include "clesphys.h"
	50	#include "YOMCST.h"
[175]	51	#include "microtab.h"
	52	#include "varmuphy.h"
	53	#include "itemps.h"
[1056]	54	#include "logic.h"
[3]	55
	56	c======================================================================
	57	c Variables argument:
	58	c
	59	LOGICAL firstcall,lastcall
	60	INTEGER nqmax,nmicro,nlat,appkim
	61	REAL ptimestep,dtkim
	62	REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
	63	REAL ptemp(klon,klev)
	64	REAL pmu0(klon), pfract(klon), pdecli, lonsol
[175]	65	REAL pu(klon),pv(klon)
	66	REAL pzlev(klon,klev+1),pzlay(klon,klev)
	67	REAL ftsol(klon)
[3]	68	REAL tr_seri(klon,klev,nqmax)
[104]	69	REAL qaer(klon,klev,nqmax)
[3]	70	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
[175]	71	REAL fclat(klon)
	72	REAL reservoir(klon)
[3]	73
	74	c======================================================================
	75	c Local variables
[175]	76	REAL qaer0(klon,klev,nqmax)
[474]	77	REAL prec(klon,5)
[3]	78
[175]	79	REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)
	80
	81	REAL vcl,nuc,xgsn,xmsn,xesn,xasn
	82
	83
	84	ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)
	85
	86	REAL socccld
	87
[3]	88	c grandeurs en moyennes zonales
[1056]	89	REAL zplev(klon,klev+1),zplay(klon,klev)
	90	REAL zzlev(klon,klev+1),zzlay(klon,klev)
	91	REAL ztemp(klon,klev), delpbar(klon,klev)
[3]	92	real temp_eq(klev),press_eq(klev)
[1056]	93	REAL qaer0bar(klon,klev,nqmax) ! et non nmicro... Permet nmicro=0.
	94	REAL zdqmufi(klon,klev,nqmax)
	95	REAL ychim(klon,klev,nqmax-nmicro)
[3]	96	c La saturation n est calculee qu une seule fois: sauvegarde qysat
	97	c La chimie n est pas calculee tous les pas, il faut donc
	98	c sauvegarder les sorties de la chimie
[104]	99	REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:)
[3]	100
[104]	101	character*10 nomqy(nqmax-nmicro+1)
[175]	102	integer i,j,k,l,iq,ig0
[104]	103
[175]	104	c indice des esp chimiques utilisees dans la microfi
	105	integer icldch4,icldc2h6,icldc2h2
	106	save icldch4,icldc2h6,icldc2h2
	107
	108	real fte,ftm,Lvch4
	109
	110	REAL tmp,ex,kmin,kmax,dqsq
	111	REAL dqch4
	112
[3]	113	c======================================================================
	114	c======================================================================
	115
[104]	116	if (firstcall) then
[1056]	117	allocate(qysat(klev,nqmax-nmicro),pdyfi(klon,klev,nqmax-nmicro))
[175]	118
	119	c -------- Quelques verifications au demarrage sur les tailles des tableaux.
	120	IF (microfi.ge.1) then
	121	c Faire de la microphysique sans traceurs... bon courage !
	122	if (nmicro.le.0) then
	123	print*,"aLeRtE cRiTiQuE !!!"
	124	print*,"Vous faites de la microphysique sans traceurs"
	125	print*,"microphysique..."
	126	print*,"Je m'arrete et vous laisse reflechir !"
	127	stop
	128	endif
	129	c Nombre de traceurs incompatibles avec la microphysique.
[176]	130	if ((nmicro.ne.ntype*nrad).and.(clouds.eq.1)) then
[175]	131	print*,"aLeRtE cRiTiQuE !!!"
	132	print*,"Nb trac imcompatible avec la microphysique."
	133	print,nmicro,ntypenrad
	134	stop
	135	endif
[176]	136	if ((nmicro.ne.nrad).and.(clouds.eq.0)) then
	137	print*,"aLeRtE cRiTiQuE !!!"
	138	print*,"Nb trac imcompatible avec la microphysique."
	139	print*,nmicro,nrad
	140	stop
	141	endif
[175]	142	ENDIF
	143
[1056]	144	endif ! firstcall
[104]	145
[175]	146	c RAZ des sorties : les moyennes se font directement dans IOIPSL :
	147	c
	148	flxesp_i(:,:,:) = 0.
	149	tau_drop(:,:) = 0.
	150	tau_aer(:,:,:) = 0.
	151	solesp(:,:,:) = 0.
[474]	152	precip(:,:) = 0. ! c'est uniquement une sortie en um/s
	153	c
	154	prec(:,:) = 0. ! c'est la variable temporaire des precipitions de la microfi
	155	! prec est en m (metre precipitable)
[175]	156
[3]	157	c-----------------------------------------------------------------------
	158	c convertion moyennes zonales et changement d unites pour microphy
	159	c ---------------------------------
	160
	161	c print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
	162
[175]	163	c -------------------
	164	c Gestion de la temperature et de la pression :
[1056]	165	c Utilisation des moyennes zonales:
	166
[175]	167	c soit la chimie est active, soit la microphysique se fait en 2D.
	168	IF (chimi.or.microfi.eq.1) THEN
[1056]	169	zplev(:,:) = zplevbar(:,:)
	170	zplay(:,:) = zplaybar(:,:)
	171	zzlev(:,:) = zzlevbar(:,:)
	172	zzlay(:,:) = zzlaybar(:,:)
	173	ztemp(:,:) = ztfibar(:,:)
[175]	174	ychim = 0.0
[1056]	175	ENDIF
[3]	176
[1056]	177	c Si la microphysique est faite en 2D:
	178	IF (microfi.eq.1) THEN
[175]	179	DO l=1,llm
	180	DO i = 1, klon
[1056]	181	delpbar(i,l) = zplevbar(i,l) - zplevbar(i,l+1)
[175]	182	ENDDO
	183	ENDDO
[1056]	184	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
	185	DO iq=1,nmicro
	186	qaer(:,:,iq) = zqfibar(:,:,iq)*delpbar(:,:)/RG
[1072]	187	DO l=1,klev
	188	DO i = 1, klon
	189	if (qaer(i,l,iq).lt.0.) then
	190	print*,"NEGS ICI ICI !!!!",qaer(i,l,iq),i,l,iq
	191	qaer(i,l,iq)=0.
	192	c stop
	193	endif
	194	if (delpbar(i,l).lt.0.) then
	195	print*,"NEGS DELP ICI !!!!",i,l,iq,delpbar(i,l)
	196	stop
	197	endif
	198	ENDDO
	199	ENDDO
[1056]	200	qaer0(:,:,iq)= tr_seri(:,:,iq)*delp(:,:)/RG
	201	qaer0bar(:,:,iq) = qaer(:,:,iq)
	202	ENDDO
	203	ENDIF
	204
	205	c Si la microphysique est faite en 3D:
	206	IF (microfi.eq.2) THEN
	207	zplev(:,:) = pplev(:,:)
	208	zplay(:,:) = pplay(:,:)
	209	zzlev(:,:) = pzlev(:,:)
	210	zzlay(:,:) = pzlay(:,:)
	211	ztemp(:,:) = ptemp(:,:)
	212	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
	213	DO iq=1,nmicro
	214	qaer(:,:,iq) = tr_seri(:,:,iq)*delp(:,:)/RG
	215	qaer0(:,:,iq)= tr_seri(:,:,iq)*delp(:,:)/RG
	216	ENDDO
	217	ENDIF
	218
	219	do l=1,llm
	220	temp_eq = tmoy
	221	press_eq = playmoy/100. ! en mbar
	222	enddo
	223
[175]	224	c -------------------
	225	c Extraction des gaz pour les nuages
[1056]	226	IF ((microfi.ge.1).and.(clouds.eq.1)) THEN
	227
[175]	228	c recuperation des indices des gaz qui nous interesse
	229	if (firstcall) then
	230	icldch4=-1
	231	icldc2h6=-1
	232	icldc2h2=-1
	233	do i=1,nqmax
	234	if (tname(i).eq."CH4") then
	235	icldch4=i
	236	c ich4=i
	237	elseif (tname(i).eq."C2H6") then
	238	icldc2h6=i
	239	elseif (tname(i).eq."C2H2") then
	240	icldc2h2=i
	241	endif
	242	enddo
	243	if (icldch4 .eq.-1 .or.
	244	& icldc2h6.eq.-1 .or.
	245	& icldc2h2.eq.-1 ) then
	246	print*, "Sacrebleu !!!"
	247	print*, "Vous voulez faire des nuages sans gaz."
	248	print*, "Mais vous etes inconscient. Je vais m'arreter la"
	249	print*, "pour vous laisser reflechir au probleme"
	250	STOP
	251	endif
	252	endif ! firstcall
	253
	254	c Saturation et fraction molaire CLOUD
	255	c Calcul des saturations pour les esp chimique de la muphy des nuages.
	256	c On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine.
	257	c Elles passent ensuite dans un common pour passer dans les I/O.
[1056]	258
[175]	259	DO l=1,llm
	260	DO i = 1, klon
	261	call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
	262	satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)
	263
	264	call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
	265	satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)
	266
	267	call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
	268	satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
	269	ENDDO
	270	ENDDO
	271
	272	c Copie des gaz (en 3D) <== UNIQUEMENT SI ON FAIT DES NUAGES
[1056]	273	if (moyzon_mu) then
	274	gaz1(:,:) = zqfibar(:,:,icldch4)
	275	gaz2(:,:) = zqfibar(:,:,icldc2h6)
	276	gaz3(:,:) = zqfibar(:,:,icldc2h2)
	277	else
	278	gaz1(:,:) = tr_seri(:,:,icldch4)
	279	gaz2(:,:) = tr_seri(:,:,icldc2h6)
	280	gaz3(:,:) = tr_seri(:,:,icldc2h2)
	281	endif
[175]	282
[1056]	283	endif ! microfi.ge.1 + clouds.eq.1
	284	c -------------------
[175]	285
[3]	286	c AUTRES TRACEURS
	287
	288	if (nqmax.gt.nmicro) then
[1056]	289	do iq=nmicro+1,nqmax
	290	if (moyzon_ch) then
	291	ychim(:,:,iq-nmicro) = zqfibar(:,:,iq)
	292	else
	293	ychim(:,:,iq-nmicro) = tr_seri(:,:,iq)
	294	endif
	295	nomqy(iq-nmicro) = tname(iq)
	296	c print*,iq-nmicro,nomqy(iq-nmicro)
[3]	297	enddo
[1056]	298	nomqy(nqmax-nmicro+1) = "HV"
[3]	299	endif
	300
	301	c-----------------------------------------------------------------------
	302	c initialisation des qysat au premier appel:
	303	c ---------------------------------
	304
	305	c!! ATTENTION, qysat pris uniquement a l'equateur
	306	c!! justifie puisque dans cette region, les var de t et p sont faibles...
	307
	308	if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	309	call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
	310	endif
	311
	312	c-----------------------------------------------------------------------
[175]	313	c Appel de la microphysique en 2D/3D !!!!!!
[3]	314	c --------------------------
	315
	316	IF(firstcall) THEN
	317	print*,'MICROPHYSIQUE ',MICROFI
	318	ENDIF
	319
[175]	320	c call begintime(tt0)
	321	IF (MICROFI.eq.0) THEN
	322	c PAS DE MICROPHYSIQUE :
	323	IF (firstcall) THEN
	324	print*,'MICROPHYSIQUE OFF-LINE',MICROFI
[3]	325	ENDIF
[1056]	326	ELSE
[175]	327	zdqmufi = 0. ! ne sert que pour chimi pour condensation
[1056]	328	call muphys(klon,
[175]	329	& zplev,zplay,zzlev,zzlay,
[1056]	330	& ztemp,qaer,gaz1,gaz2,gaz3,
[175]	331	& nmicro,ptimestep,
[1056]	332	& pmu0,pfract,
[175]	333	c -------- sorties diagnostiques
	334	& flxesp_i,
	335	& tau_drop,tau_aer,
	336	& solesp,prec)
[1056]	337
[175]	338	c NOTES :
	339	c Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
	340	c On a rien a faire mis a part copier les dq dans les d_tr
[1056]	341
[175]	342	ENDIF
	343	c call endtime(tt0,tt1)
	344	c ttmuphys=ttmuphys+tt1
[3]	345
	346	c-----------------------------------------------------------------------
[175]	347	c Gestion des sources
	348	c -------------
	349	c
	350	IF (clouds.eq.1) THEN
[1056]	351	IF (microfi.eq.1) THEN
	352	c On repasse les gaz en 3D si on a fait de la microphysique en 2D
	353	gaz1(:,:)=gaz1(:,:)*tr_seri(:,:,icldch4)/zqfibar(:,:,icldch4)
	354	gaz2(:,:)=gaz2(:,:)*tr_seri(:,:,icldc2h6)/zqfibar(:,:,icldc2h6)
	355	gaz3(:,:)=gaz3(:,:)*tr_seri(:,:,icldc2h2)/zqfibar(:,:,icldc2h2)
	356	ENDIF
	357	c Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
	358	DO i=1,klon
	359	reservoir(i) = reservoir(i)+prec(i,1)
	360	ENDDO
[175]	361
[1056]	362	CALL sources(klon,klev,ptimestep,z0,
[175]	363	& pu,pv,pplev,pzlay,pzlev,
	364	& gaz1,gaz2,gaz3,
[474]	365	& ftsol,evapch4,reservoir)
[175]	366
	367	ENDIF
	368	c-----------------------------------------------------------------------
[3]	369	c Condensation
	370	c -------------
	371
[175]	372	IF ((chimi).and.(nqmax.gt.nmicro)) then
[3]	373
[175]	374	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
[3]	375	c print*,'Condensation'
	376
	377	do iq=1,nqmax-nmicro
	378	do l=1,llm
[1056]	379	do j=1,klon
[3]	380	if (ychim(j,l,iq).gt.qysat(l,iq)) then
[175]	381	zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
	382	. / ptimestep ! / dt
[3]	383	endif
	384	enddo
	385	enddo
	386	enddo
	387
	388	ENDIF
	389
	390	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	391	c eventuellement, modif initiale de la compo
	392	c
[175]	393	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
[3]	394	c
	395	c if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	396	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	397	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
	398	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	399	c do iq=1,nqmax-nmicro
	400	c if (nomqy(iq).eq."CH4") then
	401	c do l=1,llm
[1056]	402	c do j=1,klon
[3]	403	c if (ychim(j,l,iq).le.0.015) then
[175]	404	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
[3]	405	c . / ptimestep ! / dt
	406	c endif
	407	c enddo
	408	c enddo
	409	c endif
	410	c enddo
	411	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	412	c
	413	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	414	c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
	415	c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
	416	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	417	c do iq=1,nqmax-nmicro
	418	c if (nomqy(iq).eq."C2H2") then
	419	c do l=1,llm
[1056]	420	c do j=1,klon
[3]	421	c if (ychim(j,l,iq).gt.1.e-5) then
[175]	422	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
[3]	423	c . / ptimestep ! / dt
	424	c endif
	425	c enddo
	426	c enddo
	427	c endif
	428	c if (nomqy(iq).eq."C2H6") then
	429	c do l=1,llm
[1056]	430	c do j=1,klon
[3]	431	c if (ychim(j,l,iq).gt.3.e-5) then
[175]	432	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
[3]	433	c . / ptimestep ! / dt
	434	c endif
	435	c enddo
	436	c enddo
	437	c endif
	438	c enddo
	439	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	440	c endif
	441	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	442
[175]	443	c ----- commentaire de fin (mise a jour des profil de fraction molaire)
	444
[3]	445	c-----------------------------------------------------------------------
	446	c Appel de la chimie
	447	c --------------------------
	448
	449	if((appkim.eq.1).and.(chimi)) then
	450	print*,'On passe dans la CHIMIE'
	451
	452	c do iq=1,nqmax-nmicro
	453	c if (nomqy(iq).eq."C2H2") then
	454	c print*,"C2H2top=",ychim(:,klev,iq)
	455	c endif
	456	c enddo
	457
	458	c Appel Chimie
	459	c ------------
[1056]	460	CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
[1126]	461	. ztemp,zplay,zplev,zzlay,zzlev,
[3]	462	. pdyfi)
	463	c ychim ne doit pas etre modifie, pdyfi en /s
	464
	465	endif
	466
	467	c-----------------------------------------------------------------------
	468	c retour des tendances vers 3D
	469	c ---------------------------------
	470
[1072]	471	c call WriteField_phy('phytrac_qaer01',qaer(:,:,1),klev)
	472
	473	c===============================
[3]	474	c TRACEURS MICROPHYSIQUES
[175]	475	c
[1072]	476	c===============================
[175]	477	c ---> pas de microphysique
	478	IF (microfi.eq.0) THEN
	479	DO iq=1,nmicro
	480	d_tr_mph(:,:,iq)=0.
	481	ENDDO
	482	ENDIF
[1072]	483	c===============================
[175]	484	c ---> microphysique 2D
[1072]	485	c IF (microfi.eq.1) THEN
	486	c DO iq=1,nmicro
	487	c DO l=1,llm
	488	c DO i=1,klon
[1056]	489	c on repasse le champ de traceurs en 3D (pas les tendances)
	490	c qaer est ce qui entre dans muphy, donc la moyenne zonale
	491	c qaer0 est la valeur initiale du champ
[1071]	492	c qaer0bar est la moyenne zonale initiale
[1056]	493	c la variation relative pour une bande de latitude est donc (qaer/qaer0bar)
	494	c la nouvelle valeur en un point (3D) est donc qaer0*(qaer/qaer0bar)
	495	c et la tendance: qaer0*(qaer/qaer0bar)-qaer0
[175]	496	c un petit patch :
	497	c Si la moyenne zonale au depart est "nulle" :
	498	c On a quand meme le droit de produire des traceurs dans la cellule.
	499	c On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
[888]	500	c Cela permet aussi entre autre d eviter les NaN pour les traceurs des nuages !
[175]	501	c (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
[1072]	502	c IF (qaer0bar(i,l,iq).gt.1e-100) THEN
	503	c qaer(i,l,iq) = qaer0(i,l,iq) *
	504	c & qaer(i,l,iq)/qaer0bar(i,l,iq)
	505	c ENDIF
[1056]	506	c La tendance correspond a (qaer-qaer0)/ptimestep
[1072]	507	c d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
	508	c & ptimestep
	509	c ENDDO
	510	c ENDDO
	511	c ENDDO
	512	c ---> microphysique 3D
	513	c ELSEIF(microfi.gt.1) THEN
	514	c DO iq=1,nmicro
	515	c d_tr_mph(:,:,iq)=(qaer(:,:,iq)-qaer0(:,:,iq))/ptimestep
	516	c ENDDO
	517	c ENDIF ! microfi
	518
	519	c DO iq=1,nmicro
	520	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
	521	c d_tr_mph(:,:,iq) = d_tr_mph(:,:,iq)*RG/delp(:,:)
	522	c ENDDO
	523
	524	c===============================
	525	c TOUT CE QUI EST AU-DESSUS NE MARCHE PAS: PLEIN DE NEGS...
	526	c CA MARCHE EN 3D, MAIS PAS EN MOY ZON...
	527	c===============================
	528
	529	c ---> microphysique 2D
	530	IF (microfi.eq.1) THEN
	531	DO iq=1,nmicro
	532	DO l=1,llm
	533	DO i=1,klon
	534	c ici, qaer correspond a la moy zonale modifiee par la microphys.
	535	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
	536	c en mettant ceci:
	537	d_tr_mph(i,l,iq) = (qaer(i,l,iq)*RG/delpbar(i,l)
	538	& -qaer0(i,l,iq)*RG/delp(i,l))/
[175]	539	& ptimestep
[1072]	540	c on remplace le champ 3D initial par la valeur modifiee de sa moyenne zonale
	541	c => on remet un champ uniforme en zonal...
[175]	542	ENDDO
	543	ENDDO
	544	ENDDO
	545	c ---> microphysique 3D
	546	ELSEIF(microfi.gt.1) THEN
	547	DO iq=1,nmicro
[1056]	548	d_tr_mph(:,:,iq)=(qaer(:,:,iq)-qaer0(:,:,iq))/ptimestep
[1072]	549	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
	550	d_tr_mph(:,:,iq) = d_tr_mph(:,:,iq)*RG/delp(:,:)
[175]	551	ENDDO
	552	ENDIF ! microfi
[3]	553
[1072]	554	c===============================
[808]	555
[3]	556	c AUTRES TRACEURS
	557
	558	if ((chimi).and.(nqmax.gt.nmicro)) then
	559	c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
	560	c a d_tr_kim (tendance chimique 3D en /s, passee a physiq)
[175]	561	c et de zdqmufi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
[3]	562
	563	DO iq=nmicro+1,nqmax
[1056]	564	d_tr_kim(:,:,iq) = pdyfi(:,:,iq-nmicro)
	565	& *tr_seri(:,:,iq)/ychim(:,:,iq-nmicro)
	566	d_tr_mph(:,:,iq) = zdqmufi(:,:,iq)
	567	& *tr_seri(:,:,iq)/ychim(:,:,iq-nmicro)
[3]	568	ENDDO
	569
	570	endif ! chimi
	571
[175]	572	c--------------------------------------------------
	573	c CONDENSATION VIA MICROFI
	574	c----------------------
[474]	575	c La microphysique avec nuages doit se faire obligatoirement en 3D. (FAUX ACTUELLEMENT)
[888]	576	c Rien n empeche de faire la chimie en 2D. Cependant pour prendre en compte la
[175]	577	c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
	578	c les especes concernees (CH4, C2H6 pour le moment).
	579	IF (microfi.ge.1.and.clouds.eq.1) THEN
	580	c condensation CH4
[1056]	581	d_tr_mph(:,:,icldch4) =(gaz1(:,:)-tr_seri(:,:,icldch4))
[175]	582	& /ptimestep
	583	c condensation C2H6
[1056]	584	d_tr_mph(:,:,icldc2h6)=(gaz2(:,:)-tr_seri(:,:,icldc2h6))
	585	& /ptimestep
[175]	586	c condensation C2H2
[1056]	587	d_tr_mph(:,:,icldc2h2)=(gaz3(:,:)-tr_seri(:,:,icldc2h2))
	588	& /ptimestep
[175]	589	ENDIF
	590
	591	c--------------------------------------------------
	592	c MISE A JOUR CH4 : (pour refixer la fraction
	593	c molaire)
	594	c--------------------------------------------------
	595	c IF (firstcall) THEN
	596	c do i=1,klon
	597	c do j=1,llm
	598	c call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
	599	c tmp=0.950.85tmp*28./16.
	600	c if (pplay(i,j).lt.20000.) then
	601	c dqch4 = 1.4e-2
	602	c else
	603	c dqch4 = tmp
	604	c endif
	605	c d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/
	606	c & ptimestep
	607	c enddo
	608	c enddo
	609	c
	610	c ENDIF
	611
	612	c--------------------------------------------------
[474]	613	c CONVERSION PRECIPITATION :
	614	c en microns/secondes
	615	c--------------------------------------------------
	616	precip = prec * 1.e6 / ptimestep
	617
	618
	619	c--------------------------------------------------
[888]	620	c CALCUL DU FLUX DE CHALEUR LATENTE D EVAPORATION
[175]	621	c DU METHANE
	622	c--------------------------------------------------
	623	IF (clouds.eq.1) THEN
	624	DO i=1,klon
	625	fte= (1.-ftsol(i)/305.5)
	626	ftm= (1.-ftsol(i)/190.5)
	627	if(ftm.le.1.e-3) ftm=1.e-3
	628	if(fte.le.1.e-3) fte=1.e-3
	629	Lvch4 =8.314190.4
	630	& (7.08ftm0.354+10.951.1e-2ftm*0.456)
	631	& /mch4
[474]	632	! evapch4 en m3/m2 {ok}
[175]	633	! 425 en kg/m3
	634	! Lv en J/kg {ok}
	635	! ptimestep en s {ok}
[474]	636	fclat(i)=(evapch4(i)Lvch4rhoi_ch4) ! en J/m2/s
[175]	637	ENDDO
	638	ENDIF
	639
	640	c--------------------------------------------------
	641	c GESTION DES RAYONS DE GOUTTES POUR TR
	642	c--------------------------------------------------
	643	IF (clouds.eq.1) THEN
	644
	645	c Calcul du rayon des gouttes par bin ...
	646	c----------------------------------------
	647	DO i=1,klon
	648	DO j=1,klev
	649	DO iq=1,nrad
[1056]	650	* Rayon minimum selon la quantite de noyaux
	651	IF (qaer(i,j,iq+nrad) .le. 1.e-5) THEN
	652	rcloud(i,j,iq) = 1.e-10
[175]	653	ELSE
	654	rcloud(i,j,iq)=
	655	& ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
	656	& qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
[474]	657	& v_e(iq))0.75/RPI)*(1./3.)
[175]	658	ENDIF
	659	ENDDO
	660	ENDDO
	661	ENDDO
	662
	663	c .... et de leur rayon moyen total (tt bins confondu)
	664	c------------------------------------------------------
	665	DO i=1,klon
	666	socccld=0.
	667	DO j=klev,1,-1 !de haut en bas pour le calcul des opacites
	668	vcl=0.
	669	nuc=0.
	670	xgsn=0.
	671	xmsn=0.
	672	xesn=0.
	673	xasn=0.
	674	DO iq=1,nrad
	675	vcl=vcl+qaer(i,j,iq+2*nrad)+
	676	& qaer(i,j,iq+3*nrad)+
	677	& qaer(i,j,iq+4*nrad)+
	678	& v_e(iq)*qaer(i,j,iq+nrad) ! volume des gouttes
	679	nuc=nuc+qaer(i,j,iq+nrad) ! nombre de noyaux
	680	xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq) ! volume de noyaux
	681	xmsn=xmsn+qaer(i,j,iq+2*nrad) ! volume de methane
	682	xesn=xesn+qaer(i,j,iq+3*nrad) ! volume d' ethane
	683	xasn=xasn+qaer(i,j,iq+4*nrad) ! volume d' acethylene
	684	ENDDO
	685	IF (nuc .le. 1.e-5) THEN
	686	rmcloud(i,j)=1.0e-10
	687	xfrac(i,j,:)=0.
	688	ELSE
	689	IF(xgsn/vcl.lt.0. .or. xgsn/vcl.gt.1.001)
	690	& print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl
	691	rmcloud(i,j)= ! rayon moyen des gouttes
[474]	692	& (vcl/nuc0.75/RPI)*(1./3.)
[175]	693	xfrac(i,j,1)=xgsn/vcl ! fraction volumique noyau/goutte
	694	xfrac(i,j,2)=xmsn/vcl ! fraction volumique CH4/goutte
	695	xfrac(i,j,3)=xesn/vcl ! fraction volumique C2H6/goutte
	696	xfrac(i,j,4)=xasn/vcl ! fraction volumique C2H2/goutte
	697	c calcul du rayon moyen (moyenne temporelle)
	698	rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
	699	xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
	700	ncount(i,j) = ncount(i,j)+1
	701	ENDIF
[474]	702	socccld=socccld+RPI(rmcloud(i,j)2.)nuc
[175]	703	occcld(i,j)=socccld
	704	ENDDO
	705	ENDDO
	706	c
	707	c OCCCLD
[1056]	708	c Calcul le nombre d occurence d un nuage
	709	c d opacite comprise en kmin et kmax
[175]	710	c k kmin kmax
	711	c 1 0.0000000 0.10000000
	712	c 2 0.10000000 0.17782794
	713	c 3 0.17782794 0.31622776
	714	c 4 0.31622776 0.56234139
	715	c 5 0.56234139 1.0000000
	716	c 6 1.0000000 1.7782795
	717	c 7 1.7782795 3.1622777
	718	c 8 3.1622777 5.6234136
	719	c 9 5.6234136 10.000000
	720	c 10 10.000000 17.782795
	721	c 11 17.782795 31.622778
	722	c 12 31.622778 100000.00
	723	c
[474]	724	c mise a zero de occld_m
	725	occcld_m=0.
[175]	726	DO i=1,klon
	727	DO j=1,klev
	728	DO k=1,12
	729	ex=10.**(0.25)
	730	kmin=0.
	731	kmax=1.e5
	732	if(k.ne.1) kmin=0.1ex*(k-2)
	733	if(k.ne.12) kmax=0.1ex*(k-1)
	734	if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
	735	& occcld_m(i,j,k)=1.
	736	ENDDO
	737	ENDDO
	738	ENDDO
	739	ENDIF ! fin condition clouds => pas besoin de calculer des rayons
	740
[3]	741	RETURN
	742	END
[808]	743

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

Download in other formats:

Original Format