Context Navigation

phytrac.F @ 124

Last change on this file since 124 was 104, checked in by slebonnois, 14 years ago

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

Phytitan: modifications pour pouvoir compiler correctement.
Il ne manque plus que physiq.F a faire.

File size: 11.5 KB

Rev	Line
[3]	1	SUBROUTINE phytrac (firstcall,lastcall,
	2	. nqmax,nmicro,ptimestep,appkim,dtkim,
	3	. pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
[104]	4	. lonsol,tr_seri,qaer,d_tr_mph,d_tr_kim)
[3]	5
	6	c======================================================================
	7	c S. Lebonnois, mai 2008
	8	c
	9	c Arguments:
	10	c
	11	c firstcall----input-L-variable logique indiquant le premier passage
	12	c lastcall-----input-L-variable logique indiquant le dernier passage
	13	c nqmax--------input-I-nombre de traceurs (total)
	14	c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
	15	c ptimestep----input-R-pas d'integration pour la physique (seconde)
	16	c appkim-------input-I-appel a la chimie
	17	c dtkim--------input-R-pas de temps chimique (seconde)
	18	c pplev--------input-R-pression pour chaque inter-couche (en Pa)
	19	c pplay--------input-R-pression pour chaque couche (en Pa)
	20	c delp---------input-R-epaisseur d'une couche (en Pa)
	21	c ptemp--------input-R-temperature (K)
	22	c pmu0---------input-R-cos angle zenithal
	23	c pfract-------input-R-fractional day
	24	c pdecli-------input-R-declinaison en radian
	25	c lonsol-------input-R-longitude solaire en radian
	26	c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
	27	c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
	28	c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
	29	c======================================================================
[102]	30	USE infotrac
	31	use dimphy
	32	IMPLICIT none
[3]	33	#include "dimensions.h"
	34	#include "clesphys.h"
	35	#include "YOMCST.h"
	36
	37	c======================================================================
	38	c Variables argument:
	39	c
	40	LOGICAL firstcall,lastcall
	41	INTEGER nqmax,nmicro,nlat,appkim
	42	REAL ptimestep,dtkim
	43	REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
	44	REAL ptemp(klon,klev)
	45	REAL pmu0(klon), pfract(klon), pdecli, lonsol
	46	REAL tr_seri(klon,klev,nqmax)
[104]	47	REAL qaer(klon,klev,nqmax)
[3]	48	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
	49
	50	c======================================================================
	51	c Local variables
	52
	53	c grandeurs en moyennes zonales
	54	REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev)
	55	REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
	56	real temp_eq(klev),press_eq(klev)
	57	REAL zqaer(jjm+1,klev,nqmax) ! et non nmicro... Permet nmicro=0.
	58	REAL zqaer0(jjm+1,klev,nqmax)
	59	REAL pdqmfi(jjm+1,klev,nqmax)
	60	REAL ychim(jjm+1,klev,nqmax-nmicro)
	61	c La saturation n est calculee qu une seule fois: sauvegarde qysat
	62	c La chimie n est pas calculee tous les pas, il faut donc
	63	c sauvegarder les sorties de la chimie
[104]	64	REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:)
[3]	65
[104]	66	character*10 nomqy(nqmax-nmicro+1)
	67	integer i,j,l,iq,ig0
	68
[3]	69	c======================================================================
	70	c======================================================================
	71
[104]	72	if (firstcall) then
	73	allocate(qysat(klev,nqmax-nmicro),pdyfi(jjm+1,klev,nqmax-nmicro))
	74	endif
	75
[3]	76	c-----------------------------------------------------------------------
	77	c convertion moyennes zonales et changement d unites pour microphy
	78	c ---------------------------------
	79
	80	c print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
	81
	82	zplev = 0.0
	83	zplay = 0.0
	84	ztemp = 0.0
	85	zqaer = 0.0
	86	ychim = 0.0
	87	zmu0 = 0.0
	88	zfract= 0.0
	89
	90	do l=1,llm+1
	91	zplev(1,l) = pplev(1,l)
	92	do j=2,jjm
	93	ig0=1+(j-2)*iim
	94	do i=1,iim
	95	zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
	96	enddo
	97	enddo
	98	zplev(jjm+1,l) = pplev(klon,l)
	99	enddo
	100
	101	do l=1,llm
	102	ztemp(1,l) = ptemp(1,l)
	103	zplay(1,l) = pplay(1,l)
	104	do j=2,jjm
	105	ig0=1+(j-2)*iim
	106	do i=1,iim
	107	ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
	108	zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
	109	enddo
	110	enddo
	111	ztemp(jjm+1,l) = ptemp(klon,l)
	112	zplay(jjm+1,l) = pplay(klon,l)
	113	temp_eq = ztemp((jjm+1)/2,:)
	114	press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
	115	enddo
	116
	117	zmu0(1) = pmu0(1)
	118	zfract(1) = pfract(1)
	119	do j=2,jjm
	120	ig0=1+(j-2)*iim
	121	do i=1,iim
	122	zmu0(j) = zmu0(j) + pmu0(ig0+i)/iim
	123	zfract(j) = zfract(j) + pfract(ig0+i)/iim
	124	enddo
	125	enddo
	126	zmu0(jjm+1) = pmu0(klon)
	127	zfract(jjm+1) = pfract(klon)
	128
	129	c TRACEURS MICROPHYSIQUES
	130
	131	if (microfi.eq.1) then
	132
	133	do iq=1,nmicro
	134	c print*,tname(iq)
	135
	136	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
	137	DO l=1,llm
	138	DO i = 1, klon
	139	qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
	140	ENDDO
	141	ENDDO
	142
	143	do l=1,llm
	144	zqaer(1,l,iq) = qaer(1,l,iq)
	145	do j=2,jjm
	146	ig0=1+(j-2)*iim
	147	do i=1,iim
	148	zqaer(j,l,iq) = zqaer(j,l,iq) + qaer(ig0+i,l,iq)/iim
	149	enddo
	150	enddo
	151	zqaer(jjm+1,l,iq) = qaer(klon,l,iq)
	152	enddo
	153	enddo
	154
	155	endif ! microfi
	156
	157	c AUTRES TRACEURS
	158
	159	if (nqmax.gt.nmicro) then
	160	do iq=nmicro+1,nqmax
	161	do l=1,llm
	162	ychim(1,l,iq-nmicro) = tr_seri(1,l,iq)
	163	do j=2,jjm
	164	ig0=1+(j-2)*iim
	165	do i=1,iim
	166	ychim(j,l,iq-nmicro) = ychim(j,l,iq-nmicro)
	167	. + tr_seri(ig0+i,l,iq)/iim
	168	enddo
	169	enddo
	170	ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
	171	enddo
	172	nomqy(iq-nmicro) = tname(iq)
	173
	174	c print*,iq-nmicro,nomqy(iq-nmicro)
	175
	176	enddo
	177	nomqy(nqmax-nmicro+1) = "HV"
	178	endif
	179
	180	c-----------------------------------------------------------------------
	181	c initialisation des qysat au premier appel:
	182	c ---------------------------------
	183
	184	c!! ATTENTION, qysat pris uniquement a l'equateur
	185	c!! justifie puisque dans cette region, les var de t et p sont faibles...
	186
	187	if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	188	call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
	189	endif
	190
	191	c-----------------------------------------------------------------------
	192	c Appel de la microphysique
	193	c --------------------------
	194
	195	pdqmfi = 0.
	196	zqaer0 = zqaer
	197
	198	IF(firstcall) THEN
	199	print*,'MICROPHYSIQUE ',MICROFI
	200	ENDIF
	201
	202	IF(MICROFI.eq.1) THEN
	203
	204	IF(firstcall) THEN
	205	print*,'OH! On passe dans la microphysique'
	206	ENDIF
	207
	208	CALL pg2(zplev,ztemp,zqaer,pdqmfi ! tendances 2D, /s
	209	& ,nmicro,ptimestep,zmu0,zfract,lastcall)
	210
	211	ELSE
	212
	213	IF(firstcall) THEN
	214	print*,'MICROPHYSIQUE OFF-LINE',MICROFI
	215	if (nmicro.gt.0) then
	216	CALL pg2(zplev,ztemp,zqaer,pdqmfi
	217	& ,nmicro,ptimestep,zmu0,zfract,lastcall)
	218	endif
	219	ENDIF
	220
	221	ENDIF
	222
	223	zqaer = zqaer+pdqmfi*ptimestep
	224	pdqmfi = (zqaer-zqaer0)/ptimestep
	225
	226	c-----------------------------------------------------------------------
	227	c Condensation
	228	c -------------
	229
	230	if ((chimi).and.(nqmax.gt.nmicro)) then
	231
	232	c tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
	233	c print*,'Condensation'
	234
	235	do iq=1,nqmax-nmicro
	236	do l=1,llm
	237	do j=1,jjm+1
	238	if (ychim(j,l,iq).gt.qysat(l,iq)) then
	239	pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
	240	. / ptimestep ! / dt
	241	endif
	242	enddo
	243	enddo
	244	enddo
	245
	246	ENDIF
	247
	248	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	249	c eventuellement, modif initiale de la compo
	250	c
	251	c tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
	252	c
	253	c if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	254	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	255	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
	256	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	257	c do iq=1,nqmax-nmicro
	258	c if (nomqy(iq).eq."CH4") then
	259	c do l=1,llm
	260	c do j=1,jjm+1
	261	c if (ychim(j,l,iq).le.0.015) then
	262	c pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
	263	c . / ptimestep ! / dt
	264	c endif
	265	c enddo
	266	c enddo
	267	c endif
	268	c enddo
	269	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	270	c
	271	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	272	c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
	273	c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
	274	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	275	c do iq=1,nqmax-nmicro
	276	c if (nomqy(iq).eq."C2H2") then
	277	c do l=1,llm
	278	c do j=1,jjm+1
	279	c if (ychim(j,l,iq).gt.1.e-5) then
	280	c pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
	281	c . / ptimestep ! / dt
	282	c endif
	283	c enddo
	284	c enddo
	285	c endif
	286	c if (nomqy(iq).eq."C2H6") then
	287	c do l=1,llm
	288	c do j=1,jjm+1
	289	c if (ychim(j,l,iq).gt.3.e-5) then
	290	c pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
	291	c . / ptimestep ! / dt
	292	c endif
	293	c enddo
	294	c enddo
	295	c endif
	296	c enddo
	297	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	298	c endif
	299	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	300
	301	c-----------------------------------------------------------------------
	302	c Appel de la chimie
	303	c --------------------------
	304
	305	if((appkim.eq.1).and.(chimi)) then
	306	print*,'On passe dans la CHIMIE'
	307
	308	c do iq=1,nqmax-nmicro
	309	c if (nomqy(iq).eq."C2H2") then
	310	c print*,"C2H2top=",ychim(:,klev,iq)
	311	c endif
	312	c enddo
	313
	314	c Appel Chimie
	315	c ------------
	316	CALL calchim(nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
	317	. ztemp,zplay,zplev,
	318	. pdyfi)
	319	c ychim ne doit pas etre modifie, pdyfi en /s
	320
	321	endif
	322
	323	c-----------------------------------------------------------------------
	324	c retour des tendances vers 3D
	325	c ---------------------------------
	326
	327	c TRACEURS MICROPHYSIQUES
	328
	329	if (microfi.eq.1) then
	330
	331	DO iq=1,nmicro
	332	DO l=1,llm
	333	d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
	334	ig0 = 2
	335	DO j=2,jjm
	336	DO i = 1, iim
	337	d_tr_mph(ig0,l,iq) = pdqmfi(j,l,iq)
	338	& *qaer(ig0,l,iq)/zqaer0(j,l,iq)
	339	ig0 = ig0 + 1
	340	ENDDO
	341	ENDDO
	342	d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
	343	ENDDO
	344	ENDDO
	345
	346	do iq=1,nmicro
	347	DO l=1,llm
	348	DO i = 1, klon
	349	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
	350	d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
	351	ENDDO
	352	ENDDO
	353	enddo
	354
	355	endif ! microfi
	356
	357	c AUTRES TRACEURS
	358
	359	if ((chimi).and.(nqmax.gt.nmicro)) then
	360	c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
	361	c a d_tr_kim (tendance chimique 3D en /s, passee a physiq)
	362	c et de pdqmfi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
	363
	364	DO iq=nmicro+1,nqmax
	365	DO l=1,llm
	366	d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
	367	d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
	368	ig0 = 2
	369	DO j=2,jjm
	370	DO i = 1, iim
	371	d_tr_kim(ig0,l,iq) = pdyfi(j,l,iq-nmicro)
	372	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
	373	d_tr_mph(ig0,l,iq) = pdqmfi(j,l,iq)
	374	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
	375	ig0 = ig0 + 1
	376	ENDDO
	377	ENDDO
	378	d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
	379	d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
	380	ENDDO
	381	ENDDO
	382
	383	endif ! chimi
	384
	385	RETURN
	386	END

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

Download in other formats:

Original Format