Context Navigation

muphys3D.F @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 17.3 KB

Rev	Line
[175]	1	subroutine muphys(ngrid,
	2	& plev,play,zlev,zlay,
	3	& tpt,tq,gaz1,gaz2,gaz3,
	4	& nmicro,ptimestep,
	5	& pmu0,pfract,
	6	* Sorties diagnostiques
	7	& flxesp_i,
	8	& tau_drop,tau_aer,
	9	& solesp,prec)
	10
	11	c
	12	c
	13	c CETTE NOUVELLE ROUTINE DE MICROPHYSIQUE GERE
	14	c LA MICROPHYSIQUE DES AEROSOLS ET CELLE DES NUAGES
	15	c EN UN SEUL APPEL...
	16	c
	17	c TOUS LES TRACEURS AEROSOL+NOYAUX+2*GLACES SONT CONTENUS DANS
	18	c LE TABLEAU TQ ET LEUR TENDANCES DANS TDQ. CES TABLEAUX
	19	c SONT COMPOSES DE TROIS PARTIES:
	20	c
	21	c TQ( 1, nmicro/4 pour les aerosols
	22	c + nmicro/4+1, 2*nmicro/4 pour les noyaux
	23	c + 2nmicro/4+1, 3nmicro/4 pour la glace 1
	24	c + 3*nmicro/4+1, nmicro ) pour la glace 2
	25	c
	26	c pour les aerosols, les noyaux, la glace 1 et la glace 2. la separation
	27	c puis la concatenation de fait juste avant l'appel aux routines
	28	c dans lesquels la separation est necessaire.
	29	c
	30	c _______
	31	c \| \| ____
	32	c EX: --------> QAER() ----> \| M \| \
	33	c / \| I \| \
	34	c / \| C \| \
	35	c TQ() ----------> QGLACE1() ---> \| R \| ------ TDQ()
	36	c \ \| O \| /
	37	c \ \| P \| /
	38	c --------> QGLACE2() ----> \| H \| ----/
	39	c \ \| Y \| /
	40	c \ \| \| /
	41	c -----> QNOYAUX() \| \| _/
	42	c \| \|
	43	c -------
	44	c
	45	c
	46	c
	47	c DANS LA MICROPHYSIQUE, SE TROUVENT IMBRIQUES LES PROCESSUS SUIVANTS
	48	c
	49	c
	50	c
	51	c
	52	c
	53	c NUCLEATION/SEDIMENTATION ---> n_ethane.F / n_methane.F
	54	c
	55	c MICROPHYSIQUE AEROSOLS ---> brume.F
	56	c
	57	c SEDIMENTATION DES GOUTTES ---> snuages.F
	58	c
	59	c
	60	c
	61	c
	62	c
	63	c
	64	c------------------------------------------------------
	65	use dimphy
[1056]	66	c use radcommon_h, only : volume,rayon,vrat,drayon,dvolume
	67	USE comgeomphy, only: rlatd
	68
[175]	69	IMPLICIT NONE
	70	#include "dimensions.h"
	71	#include "microtab.h"
	72	#include "varmuphy.h"
	73	#include "clesphys.h"
	74	#include "itemps.h"
	75
	76	integer ngrid
	77
	78	integer iq,nmicro
[1056]	79	real ptimestep
	80	real pdpsrf(ngrid)
[175]	81
[1056]	82	c a la place de radcommon_h:
[175]	83	common/part/vaer,raer,vrat,draer,dvaer
	84	real vaer(nrad),raer(nrad),vrat,
	85	& draer(nrad),dvaer(nrad)
	86
	87	c*************************************
	88	c declaration des variables internes *
	89	c*************************************
	90
	91	c sources *
	92	c------------------*
	93
	94	REAL plev(ngrid,klev+1)
	95	REAL play(ngrid,klev)
	96	REAL zlev(ngrid,klev+1)
	97	REAL zlay(ngrid,klev)
	98	REAL pu(ngrid),pv(ngrid)
	99	REAL pmu0(ngrid),pfract(ngrid)
	100	REAL tpt(ngrid,klev)
	101	REAL tq(ngrid,klev,nmicro),
	102	& gaz1(ngrid,klev),
	103	& gaz2(ngrid,klev),
	104	& gaz3(ngrid,klev)
	105
	106	c OUTPUT !!!!! c
	107	c note : gaz1,...,gazN sont aussi des outputs, ils ont modifié tout au long de muphys.
	108	c--------------------c
	109	REAL pdq(ngrid,klev,nmicro)
	110	REAL flxesp_i(ngrid,klev,3) ! flx esp GLACE
[474]	111	REAL solesp(ngrid,klev,3) ! tx prod glace (puit/source)
[175]	112	REAL tau_drop(ngrid,klev)
	113	REAL tau_aer(ngrid,klev,nrad)
	114	REAL prec(ngrid,5)
	115
	116	c LOCAL
	117	c--------------------c
	118	real q(ngrid,klev,nmicro)
	119	REAL taused(klev,nrad)
[1056]	120	integer jsup,jinf,h,jalt,ihor,k,im1
[175]	121
	122	c microphysique *
	123	c---------------------*
	124	real c(klev,nrad), cni(klev,nrad)
	125	real c1i(klev,nrad),c2i(klev,nrad),c3i(klev,nrad)
	126	real gazc1(klev),gazc2(klev),gazc3(klev)
	127	real ddt
	128
	129	real vcl,nuc,r,xgsn,xmsn
	130	real zz,effg,xlog,rapport
	131
	132	integer IPREM,i,n,j,l
	133	integer ibid
	134	save IPREM
	135	data IPREM/0/
	136
	137
[474]	138	c real ttq(ngrid,klev,nmicro,2)
	139	c real tttq(ngrid,klev,nmicro,2)
[175]	140
	141
	142	c **************************
	143	c INITIALISATION DE TABLEAUX
	144	c **************************
	145	c A NE FAIRE QU'UNE FOIS
	146	c **************************
	147	c
	148	IF (IPREM.eq.0) THEN
	149
	150	IF (ngrid.ne.klon) THEN
	151	print*,"aLeRte :"
[1056]	152	print*,"microfi, mais ngrid.ne.klon"
[175]	153	print*,ngrid,klon
[1056]	154	stop "je m'arrete... (muphys3D)"
[175]	155	ENDIF
	156
	157	c initialisation des constantes de la microphysique :
	158	c ----------------------------------------------
	159	call inimphycst()
	160
	161
	162	c initialisation des c(z,r), c1i(z,r), c2i(z,r)
	163	c ----------------------------------------------
	164
	165	do i=1,nmicro
	166	do n=1,ngrid
	167	do j=klev,1,-1
	168	q(n,klev+1-j,i)=tq(n,j,i) ! glaces...
	169	if(i.le.2nrad) q(n,klev+1-j,i)=tq(n,j,i)tcorrect ! noyaux+aerosol
	170	pdq(n,j,i)=0.
	171	enddo
	172	enddo
	173	enddo
	174	c ici, les tableaux definissant la structure des aerosols sont
	175	c remplis: rf,df(nmicro),r(nmicro),v(nmicro)......
	176	call rdf()
	177	c ici on recopie la grille dans un common specifique a la microfi...
[1056]	178	c v_e = volume
	179	c r_e = rayon
	180	c vrat_e = vrat
	181	c dr_e = drayon
	182	c dv_e = dvolume
[175]	183	v_e = vaer
	184	r_e = raer
	185	vrat_e = vrat
	186	dr_e = draer
	187	dv_e = dvaer
	188	c
	189	ELSE
	190
	191	c les tq() doivent etre en nombre d'aerosols / cases
	192
	193	do j=1,klev ! j de 1 a 119
	194	do n=1,ngrid
	195	do i=1,nmicro
	196	q(n,j,i)=tq(n,klev-j+1,i)
	197	pdq(n,j,i)=0.
	198	enddo
	199	enddo
	200	enddo
	201
	202	ENDIF ! FIN IPREM
	203
	204
	205	c-----------------------------------------------------
	206	c !! La premiere fois, on ne passe pas par
	207	c !! q--->c et par pg3.F
	208	c !! on passe directement au remplissage c-->q
	209	IF (IPREM.eq.0) goto 102
	210	c-----------------------------------------------------
	211
	212	c****************************************
	213	c *
	214	c ADAPTATION GCM > micro *
	215	c *
	216	c****************************************
	217
	218
	219	c correpondance des couches / sens GCM > microphysique
	220	c-----------------------------------------------------
[1056]	221
	222	c***************************************************************
[175]	223	do IHOR=1,NGRID ! GRANDE BOUCLE HORIZONTALE / SEPARATION DES COLONNES
	224
[1056]	225	if (IHOR.eq.1) then
	226	im1=1
	227	else
	228	im1=IHOR-1
	229	endif
	230
	231	c***************************************************************
	232	c On refait les calculs si on est au premier point
	233	c OU si on change de latitude
	234	c OU si on calcule la microfi en 3D
	235	c***************************************************************
	236	if((IHOR.eq.1)
	237	& .or.(rlatd(IHOR).ne.rlatd(im1))
	238	& .or.(microfi.eq.2)) then
	239	c***************************************************************
[175]	240
	241	c Ici, on initialise la grille verticale et les
	242	c variables communes aux routines de microphysique.
	243	c*******************************************************
	244	call inimuphy(ihor,plev(ihor,:),play(ihor,:),
	245	& zlev(ihor,:),zlay(ihor,:),
	246	& tpt(ihor,:))
	247
	248	c Ici, on scinde les tableaux aerosols et glaces
	249	c*******************************************************
	250
[176]	251	if (clouds.eq.0) then
	252	if(nrad .ne. nmicro) then
	253	print*,"aLeRte : nrad != nmicro"
	254	print*,'nmicro= ',nmicro
	255	print*,'nrad= ',nrad
	256	stop "je m'arrete..."
	257	endif
	258	else
[175]	259	if(nrad .ne. nmicro/ntype) then
	260	print*,"aLeRte : nrad != nmicro/ntype"
	261	print*,'nmicro= ',nmicro
	262	print*,'ntype=',ntype
	263	print*,'nmicro/ntype= ',nmicro/ntype
	264	print*,'nrad= ',nrad
	265	stop "je m'arrete..."
	266	endif
[176]	267	endif
[175]	268
	269	do i=1,nrad
	270	do j=1,klev
	271	c(j,i) =q(IHOR,j,i )/dzb(j) ! concentration aerosols/m^3
	272	if (clouds.eq.1) then
	273	cni(j,i)=q(IHOR,j,i+ nrad)/dzb(j) ! concentration noyaux /m^3
	274	c1i(j,i)=q(IHOR,j,i+2*nrad)/dzb(j) ! concentration volume glace/m^3
	275	c2i(j,i)=q(IHOR,j,i+3*nrad)/dzb(j) ! concentration volume glace /m^3
	276	c3i(j,i)=q(IHOR,j,i+4*nrad)/dzb(j) ! concentration volume glace /m^3
	277	endif
	278	enddo
	279	enddo
	280	if (clouds.eq.1) then
	281	do j=1,klev
	282	gazc1(j) =gaz1(IHOR,klev-j+1) ! fraction molaire CH4
	283	gazc2(j) =gaz2(IHOR,klev-j+1) ! fraction molaire C2H6
	284	gazc3(j) =gaz3(IHOR,klev-j+1) ! fraction molaire C2H2
	285	enddo
	286	endif
	287
	288	c****************************************
	289	c
	290	c FIN DE LA PREPARATION:
	291	c
	292	c
	293	c ON APPELLE LES MODELES MICROPHYSIQUES
	294	c
	295	c - brume (coagulation + sedimentation)
	296	c - nuages (nucleation + condensation)
	297	c - sedimentation nuages
	298	c
	299	c
	300	c****************************************
	301
	302
	303	do j=1,klev
	304	solesp(ihor,klev+1-j,:) = 0.
	305	do i=1,nrad
	306	tau_aer(ihor,klev+1-j,i)=0.
	307	enddo
	308	enddo
	309
	310	ddt=ptimestep
	311
	312	* concerne les aerosols (tableau c):
	313	c call begintime(tt0)
	314	call brume(ngrid,c,ddt,klev,nrad,taused,ihor,
	315	& pmu0(ihor),pfract(ihor),
	316	& prec)
	317	c call endtime(tt0,tt1)
	318	c tthaze=tthaze+tt1
	319
	320	* concerne aerosols + gouttes (tableaux c,cn,c1,c2):
	321
	322	do j=1,klev
	323	do i=1,nrad
	324	tau_aer(ihor,klev+1-j,i)=tau_aer(ihor,klev+1-j,i)
	325	& +taused(j,i)
	326	enddo
	327	enddo
	328
	329	IF (clouds.eq.1) THEN
	330
	331	c do j=1,klev
	332	c do i=1,nrad
	333	c ttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
	334	c ttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
	335	c ttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
	336	c ttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
	337	c ttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
	338	c enddo
	339	c enddo
	340
	341	c call begintime(tt0)
	342	call cnuages(c,c1i,c2i,c3i,cni,
	343	& gazc1,gazc2,gazc3,ddt)
	344	c call endtime(tt0,tt1)
	345	c ttcclds=ttcclds+tt1
	346
	347	c verification des valeurs de c,cni,c1i,c2i et c3i.
	348	c Lorsque l'on vide completement une case, on peut avoir des chiffres negatifs :s
	349	do j=1,klev
	350	do i=1,nrad
	351	c(j,i)= MAX(c(j,i),0.)
	352	cni(j,i)= MAX(cni(j,i),0.)
	353	c1i(j,i)= MAX(c1i(j,i),0.)
	354	c2i(j,i)= MAX(c2i(j,i),0.)
	355	c3i(j,i)= MAX(c3i(j,i),0.)
	356	c ttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
	357	c ttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
	358	c ttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
	359	c ttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
	360	c ttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
	361	enddo
	362	enddo
	363
	364
	365	do j=1,klev
	366	do i=1,nrad
[474]	367	! solesp en m3/m3 pour passer en m3/m2 il faut faire :
	368	! (c1i(j,i)dzb(j) -q(IHOR,j,i+2nrad))
[175]	369	solesp(ihor,klev+1-j,1)=solesp(ihor,klev+1-j,1) +
	370	& (c1i(j,i)-q(IHOR,j,i+2*nrad)/dzb(j))
	371	solesp(ihor,klev+1-j,2)=solesp(ihor,klev+1-j,2) +
	372	& (c2i(j,i)-q(IHOR,j,i+3*nrad)/dzb(j))
	373	solesp(ihor,klev+1-j,3)=solesp(ihor,klev+1-j,3) +
	374	& (c3i(j,i)-q(IHOR,j,i+4*nrad)/dzb(j))
	375	enddo
	376	enddo
	377
	378	* concerne les gouttes (tableaux cn,c1,c2):
	379
	380	c do j=1,klev
	381	c do i=1,nrad
	382	c tttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
	383	c tttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
	384	c tttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
	385	c tttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
	386	c tttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
	387	c enddo
	388	c enddo
	389
	390	c call begintime(tt0)
	391	call snuages(ngrid,cni,c1i,c2i,c3i,c,ddt,
	392	& klev,nrad,ihor,
	393	& flxesp_i,tau_drop,prec)
	394	c call endtime(tt0,tt1)
	395	c ttsclds=ttsclds+tt1
	396
	397	c do j=1,klev
	398	c do i=1,nrad
	399	c tttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
	400	c tttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
	401	c tttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
	402	c tttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
	403	c tttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
	404	c enddo
	405	c enddo
	406	ENDIF ! flag nuages :)
	407
	408
	409	* on recompose le tableau de traceurs ici.
	410	*--------------------------------------------------
	411
	412	do i=1,nrad
	413	do j=1,klev
	414	q(IHOR,j,i) = c(j,i)*dzb(j) ! nombre aerosols /m^2
	415	if (clouds.eq.1) then
	416	q(IHOR,j,i+ nrad) = cni(j,i)*dzb(j) ! nombre noyaux /m^2
	417	q(IHOR,j,i+2nrad) = c1i(j,i)dzb(j) ! concentration volume glace/m^2
	418	q(IHOR,j,i+3nrad) = c2i(j,i)dzb(j) ! concentration volume glace/m^2
	419	q(IHOR,j,i+4nrad) = c3i(j,i)dzb(j) ! concentration volume glace/m^2
	420	endif
	421	enddo
	422	enddo
	423	if (clouds.eq.1) then
	424	do j=1,klev
	425	gaz1(IHOR,klev-j+1) = gazc1(j) ! fraction molaire
	426	gaz2(IHOR,klev-j+1) = gazc2(j) ! fraction molaire
	427	gaz3(IHOR,klev-j+1) = gazc3(j) ! fraction molaire
	428	enddo
	429	endif
	430
[1056]	431	c***************************************************************
	432	else ! same latitude, we don't do calculations again
	433	q(ihor,:,:) = q(im1,:,:)
	434	tau_aer(ihor,:,:) = tau_aer(im1,:,:)
	435	prec(ihor,:) = prec(im1,:)
	436	if (clouds.eq.1) then
	437	solesp(ihor,:,:) = solesp(im1,:,:)
	438	flxesp_i(ihor,:,:) = flxesp_i(im1,:,:)
	439	tau_drop(ihor,:) = tau_drop(im1,:)
	440	gaz1(ihor,:) = gaz1(im1,:)
	441	gaz2(ihor,:) = gaz2(im1,:)
	442	gaz3(ihor,:) = gaz3(im1,:)
	443	endif
	444	endif
	445
[175]	446	ENDDO ! Fin de la boucle IHOR
[1056]	447	c***************************************************************
[175]	448
	449	102 CONTINUE ! la premiere fois, c'est une boucle vide!
	450
	451
	452	c***************************************************************
	453	c FIN: on renvoie les nouvelles valeurs q(t+dt)=q(t) + dq(t)
	454	c
	455	c Pour les aerosols, les noyaux, les glaces et les vapeurs modifiees...
	456	c
	457	c***************************************************************
	458
	459	do n=1,ngrid
	460	do i=1,nmicro
	461	do j=1,klev ! j de 1 a 54
	462	tq(n,j,i) = q(n,klev+1-j,i)
	463	enddo
	464	enddo
	465	enddo
	466
	467
	468	c do j=1,15
	469	cc CH4 -- cnuages
	470	c write(210,'(I4,3(ES24.16,1X))') j,
	471	c & sum(ttq(40,j,2nrad+1:3nrad,1)),
	472	c & sum(ttq(40,j,2nrad+1:3nrad,2)),
	473	c & sum(ttq(40,j,2nrad+1:3nrad,2))-
	474	c & sum(ttq(40,j,2nrad+1:3nrad,1))
	475	cc C2H6 -- cnuages
	476	c write(211,'(I4,3(ES24.16,1X))') j,
	477	c & sum(ttq(40,j,3nrad+1:4nrad,1)),
	478	c & sum(ttq(40,j,3nrad+1:4nrad,2)),
	479	c & sum(ttq(40,j,3nrad+1:4nrad,2))-
	480	c & sum(ttq(40,j,3nrad+1:4nrad,1))
	481	cc C2H2 -- cnuages
	482	c write(212,'(I4,3(ES24.16,1X))') j,
	483	c & sum(ttq(40,j,4nrad+1:5nrad,1)),
	484	c & sum(ttq(40,j,4nrad+1:5nrad,2)),
	485	c & sum(ttq(40,j,4nrad+1:5nrad,2))-
	486	c & sum(ttq(40,j,4nrad+1:5nrad,1))
	487	c
	488	cc CH4 -- snuages
	489	c write(310,'(I4,3(ES24.16,1X))') j,
	490	c & sum(tttq(40,j,2nrad+1:3nrad,1)),
	491	c & sum(tttq(40,j,2nrad+1:3nrad,2)),
	492	c & sum(tttq(40,j,2nrad+1:3nrad,2))-
	493	c & sum(tttq(40,j,2nrad+1:3nrad,1))
	494	cc C2H6 -- snuages
	495	c write(311,'(I4,3(ES24.16,1X))') j,
	496	c & sum(tttq(40,j,3nrad+1:4nrad,1)),
	497	c & sum(tttq(40,j,3nrad+1:4nrad,2)),
	498	c & sum(tttq(40,j,3nrad+1:4nrad,2))-
	499	c & sum(tttq(40,j,3nrad+1:4nrad,1))
	500	cc C2H2 -- snuages
	501	c write(312,'(I4,3(ES24.16,1X))') j,
	502	c & sum(tttq(40,j,4nrad+1:5nrad,1)),
	503	c & sum(tttq(40,j,4nrad+1:5nrad,2)),
	504	c & sum(tttq(40,j,4nrad+1:5nrad,2))-
	505	c & sum(tttq(40,j,4nrad+1:5nrad,1))
	506	c enddo
	507	c write(210,*) "NEWLINE"
	508	c write(211,*) "NEWLINE"
	509	c write(212,*) "NEWLINE"
	510	c write(310,*) "NEWLINE"
	511	c write(311,*) "NEWLINE"
	512	c write(312,*) "NEWLINE"
	513
	514
	515	c do j=1,20
	516	c write(410,'(I4,3(ES24.16,1X))') j,
	517	c & flxesp_i(40,j,1),flxesp_i(40,j,2),flxesp_i(40,j,3)
	518	c write(510,'(I4,3(ES24.16,1X))') j,
	519	c & solesp(40,j,1),solesp(40,j,2),solesp(40,j,3)
	520	c enddo
	521	c write(410,*) "NEWLINE"
	522	c write(510,*) "NEWLINE"
	523
	524	IPREM=1 ! LA PROCHAINE FOIS NE SERA PLUS LA 1ERE
	525
	526
	527	16 return
	528
	529	end
	530
	531
	532	c---------------------------------------------------------------------

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

Download in other formats:

Original Format