Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

vdif_kc.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 19.8 KB

Rev	Line
[3175]	1	SUBROUTINE vdif_kc(dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
	2	IMPLICIT NONE
	3	c.......................................................................
	4	#include "dimensions.h"
	5	#include "dimphys.h"
	6	c.......................................................................
	7	c
	8	c dt : pas de temps
	9	c g : g
	10	c zlev : altitude a chaque niveau (interface inferieure de la couche
	11	c de meme indice)
	12	c zlay : altitude au centre de chaque couche
	13	c u,v : vitesse au centre de chaque couche
	14	c (en entree : la valeur au debut du pas de temps)
	15	c teta : temperature potentielle au centre de chaque couche
	16	c (en entree : la valeur au debut du pas de temps)
	17	c cd : cdrag
	18	c (en entree : la valeur au debut du pas de temps)
	19	c q2 : $q^2$ au bas de chaque couche
	20	c (en entree : la valeur au debut du pas de temps)
	21	c (en sortie : la valeur a la fin du pas de temps)
	22	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	23	c couche)
	24	c (en sortie : la valeur a la fin du pas de temps)
	25	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	26	c (en sortie : la valeur a la fin du pas de temps)
	27	c
	28	c.......................................................................
	29	REAL dt,g
	30	REAL zlev(ngridmx,nlayermx+1)
	31	REAL zlay(ngridmx,nlayermx)
	32	REAL u(ngridmx,nlayermx)
	33	REAL v(ngridmx,nlayermx)
	34	REAL teta(ngridmx,nlayermx)
	35	REAL cd(ngridmx)
	36	REAL q2(ngridmx,nlayermx+1)
	37	REAL km(ngridmx,nlayermx+1)
	38	REAL kn(ngridmx,nlayermx+1)
	39	c.......................................................................
	40	c
	41	c nlay : nombre de couches
	42	c nlev : nombre de niveaux
	43	c ngrid : nombre de points de grille
	44	c unsdz : 1 sur l'epaisseur de couche
	45	c unsdzdec : 1 sur la distance entre le centre de la couche et le
	46	c centre de la couche inferieure
	47	c q : echelle de vitesse au bas de chaque couche
	48	c (valeur a la fin du pas de temps)
	49	c
	50	c.......................................................................
	51	INTEGER nlay,nlev,ngrid
	52	REAL unsdz(ngridmx,nlayermx)
	53	REAL unsdzdec(ngridmx,nlayermx+1)
	54	REAL q(ngridmx,nlayermx+1)
	55	c.......................................................................
	56	c
	57	c kmpre : km au debut du pas de temps
	58	c qcstat : q : solution stationnaire du probleme couple
	59	c (valeur a la fin du pas de temps)
	60	c q2cstat : q2 : solution stationnaire du probleme couple
	61	c (valeur a la fin du pas de temps)
	62	c
	63	c.......................................................................
	64	REAL kmpre(ngridmx,nlayermx+1)
	65	REAL qcstat
	66	REAL q2cstat
	67	c.......................................................................
	68	c
	69	c long : longueur de melange calculee selon Blackadar
	70	c
	71	c.......................................................................
	72	REAL long(ngridmx,nlayermx+1)
	73	c.......................................................................
	74	c
	75	c kmq3 : terme en q^3 dans le developpement de km
	76	c (valeur au debut du pas de temps)
	77	c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
	78	c (valeur a la fin du pas de temps)
	79	c knq3 : terme en q^3 dans le developpement de kn
	80	c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
	81	c (valeur a la fin du pas de temps)
	82	c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
	83	c (valeur a la fin du pas de temps)
	84	c m : valeur a la fin du pas de temps
	85	c mpre : valeur au debut du pas de temps
	86	c m2 : valeur a la fin du pas de temps
	87	c n2 : valeur a la fin du pas de temps
	88	c
	89	c.......................................................................
	90	REAL kmq3
	91	REAL kmcstat
	92	REAL knq3
	93	REAL mcstat
	94	REAL m2cstat
	95	REAL m(ngridmx,nlayermx+1)
	96	REAL mpre(ngridmx,nlayermx+1)
	97	REAL m2(ngridmx,nlayermx+1)
	98	REAL n2(ngridmx,nlayermx+1)
	99	c.......................................................................
	100	c
	101	c gn : intermediaire pour les coefficients de stabilite
	102	c gnmin : borne inferieure de gn (-0.23 ou -0.28)
	103	c gnmax : borne superieure de gn (0.0233)
	104	c gninf : vrai si gn est en dessous de sa borne inferieure
	105	c gnsup : vrai si gn est en dessus de sa borne superieure
	106	c gm : drole d'objet bien utile
	107	c ri : nombre de Richardson
	108	c sn : coefficient de stabilite pour n
	109	c snq2 : premier terme du developement limite de sn en q2
	110	c sm : coefficient de stabilite pour m
	111	c smq2 : premier terme du developement limite de sm en q2
	112	c
	113	c.......................................................................
	114	REAL gn
	115	REAL gnmin
	116	REAL gnmax
	117	LOGICAL gninf
	118	LOGICAL gnsup
	119	REAL gm
	120	c REAL ri(ngridmx,nlayermx+1)
	121	REAL sn(ngridmx,nlayermx+1)
	122	REAL snq2(ngridmx,nlayermx+1)
	123	REAL sm(ngridmx,nlayermx+1)
	124	REAL smq2(ngridmx,nlayermx+1)
	125	c.......................................................................
	126	c
	127	c kappa : consatnte de Von Karman (0.4)
	128	c long0 : longueur de reference pour le calcul de long (160)
	129	c a1,a2,b1,b2,c1 : constantes d'origine pour les coefficients
	130	c de stabilite (0.92/0.74/16.6/10.1/0.08)
	131	c cn1,cn2 : constantes pour sn
	132	c cm1,cm2,cm3,cm4 : constantes pour sm
	133	c
	134	c.......................................................................
	135	REAL kappa
	136	REAL long0
	137	REAL a1,a2,b1,b2,c1
	138	REAL cn1,cn2
	139	REAL cm1,cm2,cm3,cm4
	140	c.......................................................................
	141	c
	142	c termq : termes en $q$ dans l'equation de q2
	143	c termq3 : termes en $q^3$ dans l'equation de q2
	144	c termqm2 : termes en $q*m^2$ dans l'equation de q2
	145	c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
	146	c
	147	c.......................................................................
	148	REAL termq
	149	REAL termq3
	150	REAL termqm2
	151	REAL termq3m2
	152	c.......................................................................
	153	c
	154	c q2min : borne inferieure de q2
	155	c q2max : borne superieure de q2
	156	c
	157	c.......................................................................
	158	REAL q2min
	159	REAL q2max
	160	c.......................................................................
	161	c knmin : borne inferieure de kn
	162	c kmmin : borne inferieure de km
	163	c.......................................................................
	164	REAL knmin
	165	REAL kmmin
	166	c.......................................................................
	167	INTEGER ilay,ilev,igrid
	168	REAL tmp1,tmp2
	169	c.......................................................................
	170	PARAMETER (kappa=0.4E+0)
	171	PARAMETER (long0=160.E+0)
	172	PARAMETER (gnmin=-10.E+0)
	173	PARAMETER (gnmax=0.0233E+0)
	174	PARAMETER (a1=0.92E+0)
	175	PARAMETER (a2=0.74E+0)
	176	PARAMETER (b1=16.6E+0)
	177	PARAMETER (b2=10.1E+0)
	178	PARAMETER (c1=0.08E+0)
	179	PARAMETER (knmin=1.E-5)
	180	PARAMETER (kmmin=1.E-5)
	181	PARAMETER (q2min=1.E-3)
	182	PARAMETER (q2max=1.E+2)
	183	PARAMETER (nlay=nlayermx)
	184	PARAMETER (nlev=nlayermx+1)
	185	PARAMETER (ngrid=ngridmx)
	186	c
	187	PARAMETER (
	188	& cn1=a2(1.E+0 -6.E+0 a1/b1)
	189	& )
	190	PARAMETER (
	191	& cn2=-3.E+0 a2(6.E+0 *a1+b2)
	192	& )
	193	PARAMETER (
	194	& cm1=a1(1.E+0 -3.E+0 c1-6.E+0 *a1/b1)
	195	& )
	196	PARAMETER (
	197	& cm2=a1(-3.E+0 a2((b2-3.E+0 a2)(1.E+0 -6.E+0 a1/b1)
	198	& -3.E+0 c1(b2+6.E+0 *a1)))
	199	& )
	200	PARAMETER (
	201	& cm3=-3.E+0 a2(6.E+0 *a1+b2)
	202	& )
	203	PARAMETER (
	204	& cm4=-9.E+0 a1a2
	205	& )
	206	c.......................................................................
	207	c traitment des valeur de q2 en entree
	208	c.......................................................................
	209	c
	210	DO ilev=1,nlev
	211	DO igrid=1,ngrid
	212	q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
	213	q(igrid,ilev)=sqrt(q2(igrid,ilev))
	214	ENDDO
	215	ENDDO
	216	c
	217	DO igrid=1,ngrid
	218	tmp1=cd(igrid)(u(igrid,1)2+v(igrid,1)*2)
	219	q2(igrid,1)=b1*(2.E+0/3.E+0)tmp1
	220	q2(igrid,1)=amax1(q2(igrid,1),q2min)
	221	q(igrid,1)=sqrt(q2(igrid,1))
	222	ENDDO
	223	c
	224	c.......................................................................
	225	c les increments verticaux
	226	c.......................................................................
	227	c
	228	c!!!!! allerte !!!!!c
	229	c!!!!! zlev n'est pas declare a nlev !!!!!c
	230	c!!!!! ---->
	231	DO igrid=1,ngrid
	232	zlev(igrid,nlev)=zlay(igrid,nlay)
	233	& +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
	234	ENDDO
	235	c!!!!! <----
	236	c!!!!! allerte !!!!!c
	237	c
	238	DO ilay=1,nlay
	239	DO igrid=1,ngrid
	240	unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
	241	ENDDO
	242	ENDDO
	243	DO igrid=1,ngrid
	244	unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
	245	ENDDO
	246	DO ilay=2,nlay
	247	DO igrid=1,ngrid
	248	unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
	249	ENDDO
	250	ENDDO
	251	DO igrid=1,ngrid
	252	unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
	253	ENDDO
	254	c
	255	c.......................................................................
	256	c le cisaillement et le gradient de temperature
	257	c.......................................................................
	258	c
	259	DO igrid=1,ngrid
	260	m2(igrid,1)=(unsdzdec(igrid,1)
	261	& u(igrid,1))*2
	262	& +(unsdzdec(igrid,1)
	263	& v(igrid,1))*2
	264	m(igrid,1)=sqrt(m2(igrid,1))
	265	mpre(igrid,1)=m(igrid,1)
	266	ENDDO
	267	c
	268	c-----------------------------------------------------------------------
	269	DO ilev=2,nlev-1
	270	DO igrid=1,ngrid
	271	c-----------------------------------------------------------------------
	272	c
	273	n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
	274	& *(teta(igrid,ilev)-teta(igrid,ilev-1))
	275	& /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
	276	c
	277	c --->
	278	c on ne sais traiter que les cas stratifies. et l'ajustement
	279	c convectif est cense faire en sorte que seul des configurations
	280	c stratifiees soient rencontrees en entree de cette routine.
	281	c mais, bon ... on sait jamais (meme on sait que n2 prends
	282	c quelques valeurs negatives ... parfois) alors :
	283	c<---
	284	c
	285	IF (n2(igrid,ilev).lt.0.E+0) THEN
	286	n2(igrid,ilev)=0.E+0
	287	ENDIF
	288	c
	289	m2(igrid,ilev)=(unsdzdec(igrid,ilev)
	290	& (u(igrid,ilev)-u(igrid,ilev-1)))*2
	291	& +(unsdzdec(igrid,ilev)
	292	& (v(igrid,ilev)-v(igrid,ilev-1)))*2
	293	m(igrid,ilev)=sqrt(m2(igrid,ilev))
	294	mpre(igrid,ilev)=m(igrid,ilev)
	295	c
	296	c-----------------------------------------------------------------------
	297	ENDDO
	298	ENDDO
	299	c-----------------------------------------------------------------------
	300	c
	301	DO igrid=1,ngrid
	302	m2(igrid,nlev)=m2(igrid,nlev-1)
	303	m(igrid,nlev)=m(igrid,nlev-1)
	304	mpre(igrid,nlev)=m(igrid,nlev)
	305	ENDDO
	306	c
	307	c.......................................................................
	308	c calcul des fonctions de stabilite
	309	c.......................................................................
	310	c
	311	c-----------------------------------------------------------------------
	312	DO ilev=2,nlev-1
	313	DO igrid=1,ngrid
	314	c-----------------------------------------------------------------------
	315	c
	316	tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
	317	long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0)
	318	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
	319	& * n2(igrid,ilev)
	320	gm=long(igrid,ilev)**2 / q2(igrid,ilev)
	321	& * m2(igrid,ilev)
	322	c
	323	gninf=.false.
	324	gnsup=.false.
	325	long(igrid,ilev)=long(igrid,ilev)
	326	long(igrid,ilev)=long(igrid,ilev)
	327	c
	328	IF (gn.lt.gnmin) THEN
	329	gninf=.true.
	330	gn=gnmin
	331	ENDIF
	332	c
	333	IF (gn.gt.gnmax) THEN
	334	gnsup=.true.
	335	gn=gnmax
	336	ENDIF
	337	c
	338	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
	339	sm(igrid,ilev)=
	340	& (cm1+cm2*gn)
	341	& /( (1.E+0 +cm3*gn)
	342	& (1.E+0 +cm4gn) )
	343	c
	344	IF ((gninf).or.(gnsup)) THEN
	345	snq2(igrid,ilev)=0.E+0
	346	smq2(igrid,ilev)=0.E+0
	347	ELSE
	348	snq2(igrid,ilev)=
	349	& -gn
	350	& (-cn1cn2/(1.E+0 +cn2gn)*2 )
	351	smq2(igrid,ilev)=
	352	& -gn
	353	& ( cm2(1.E+0 +cm3*gn)
	354	& (1.E+0 +cm4gn)
	355	& -( cm3(1.E+0 +cm4gn)
	356	& +cm4(1.E+0 +cm3gn) )
	357	& (cm1+cm2gn) )
	358	& /( (1.E+0 +cm3*gn)
	359	& (1.E+0 +cm4gn) )**2
	360	ENDIF
	361	c
	362	c --->
	363	c la decomposition de Taylor en q2 n'a de sens que
	364	c dans les cas stratifies ou sn et sm sont quasi
	365	c proportionnels a q2. ailleurs on laisse le meme
	366	c algorithme car l'ajustement convectif fait le travail.
	367	c mais c'est delirant quand sn et snq2 n'ont pas le meme
	368	c signe : dans ces cas, on ne fait pas la decomposition.
	369	c<---
	370	c
	371	IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
	372	& snq2(igrid,ilev)=0.E+0
	373	IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
	374	& smq2(igrid,ilev)=0.E+0
	375	c
	376	c-----------------------------------------------------------------------
	377	ENDDO
	378	ENDDO
	379	c-----------------------------------------------------------------------
	380	c
	381	c.......................................................................
	382	c calcul de km et kn au debut du pas de temps
	383	c.......................................................................
	384	c
	385	DO igrid=1,ngrid
	386	kn(igrid,1)=knmin
	387	km(igrid,1)=kmmin
	388	kmpre(igrid,1)=km(igrid,1)
	389	ENDDO
	390	c
	391	c-----------------------------------------------------------------------
	392	DO ilev=2,nlev-1
	393	DO igrid=1,ngrid
	394	c-----------------------------------------------------------------------
	395	c
	396	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	397	& *sn(igrid,ilev)
	398	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	399	& *sm(igrid,ilev)
	400	kmpre(igrid,ilev)=km(igrid,ilev)
	401	c
	402	c-----------------------------------------------------------------------
	403	ENDDO
	404	ENDDO
	405	c-----------------------------------------------------------------------
	406	c
	407	DO igrid=1,ngrid
	408	kn(igrid,nlev)=kn(igrid,nlev-1)
	409	km(igrid,nlev)=km(igrid,nlev-1)
	410	kmpre(igrid,nlev)=km(igrid,nlev)
	411	ENDDO
	412	c
	413	c.......................................................................
	414	c boucle sur les niveaux 2 a nlev-1
	415	c.......................................................................
	416	c
	417	c---->
	418	DO 10001 ilev=2,nlev-1
	419	c---->
	420	DO 10002 igrid=1,ngrid
	421	c
	422	c.......................................................................
	423	c
	424	c calcul des termes sources et puits de l'equation de q2
	425	c ------------------------------------------------------
	426	c
	427	knq3=kn(igrid,ilev)*snq2(igrid,ilev)
	428	& /sn(igrid,ilev)
	429	kmq3=km(igrid,ilev)*smq2(igrid,ilev)
	430	& /sm(igrid,ilev)
	431	c
	432	termq=0.E+0
	433	termq3=0.E+0
	434	termqm2=0.E+0
	435	termq3m2=0.E+0
	436	c
	437	tmp1=dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
	438	tmp2=dt2.E+0 kmq3*m2(igrid,ilev)
	439	termqm2=termqm2
	440	& +dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
	441	& -dt2.E+0 kmq3*m2(igrid,ilev)
	442	termq3m2=termq3m2
	443	& +dt2.E+0 kmq3*m2(igrid,ilev)
	444	c
	445	termq=termq
	446	& -dt2.E+0 kn(igrid,ilev)*n2(igrid,ilev)
	447	& +dt2.E+0 knq3*n2(igrid,ilev)
	448	termq3=termq3
	449	& -dt2.E+0 knq3*n2(igrid,ilev)
	450	c
	451	termq3=termq3
	452	& -dt2.E+0 q(igrid,ilev)*3 / (b1long(igrid,ilev))
	453	c
	454	c.......................................................................
	455	c
	456	c resolution stationnaire couplee avec le gradient de vitesse local
	457	c -----------------------------------------------------------------
	458	c
	459	c -----{on cherche le cisaillement qui annule l'equation de q^2
	460	c supposee en q3}
	461	c
	462	tmp1=termq+termq3
	463	tmp2=termqm2+termq3m2
	464	m2cstat=m2(igrid,ilev)
	465	& -(tmp1+tmp2)/(dt2.E+0km(igrid,ilev))
	466	mcstat=sqrt(m2cstat)
	467	c
	468	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	469	c supposee en q3}
	470	c
	471	IF (ilev.eq.2) THEN
	472	kmcstat=1.E+0 / mcstat
	473	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
	474	& *mpre(igrid,ilev+1)
	475	& +unsdz(igrid,ilev-1)
	476	& *cd(igrid)
	477	& ( sqrt(u(igrid,3)2+v(igrid,3)*2)
	478	& -mcstat/unsdzdec(igrid,ilev)
	479	& -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
	480	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
	481	ELSE
	482	kmcstat=1.E+0 / mcstat
	483	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
	484	& *mpre(igrid,ilev+1)
	485	& +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
	486	& *mpre(igrid,ilev-1) )
	487	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
	488	ENDIF
	489	tmp2=kmcstat
	490	& /( sm(igrid,ilev)/q2(igrid,ilev) )
	491	& /long(igrid,ilev)
	492	qcstat=tmp2**(1.E+0/3.E+0)
	493	q2cstat=qcstat**2
	494	c
	495	c.......................................................................
	496	c
	497	c choix de la solution finale
	498	c ---------------------------
	499	c
	500	q(igrid,ilev)=qcstat
	501	q2(igrid,ilev)=q2cstat
	502	m(igrid,ilev)=mcstat
	503	m2(igrid,ilev)=m2cstat
	504	c
	505	c --->
	506	c pour des raisons simples q2 est minore
	507	c<---
	508	c
	509	IF (q2(igrid,ilev).lt.q2min) THEN
	510	q2(igrid,ilev)=q2min
	511	q(igrid,ilev)=sqrt(q2min)
	512	ENDIF
	513	c
	514	c.......................................................................
	515	c
	516	c calcul final de kn et km
	517	c ------------------------
	518	c
	519	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
	520	& * n2(igrid,ilev)
	521	IF (gn.lt.gnmin) gn=gnmin
	522	IF (gn.gt.gnmax) gn=gnmax
	523	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
	524	sm(igrid,ilev)=
	525	& (cm1+cm2*gn)
	526	& /( (1.E+0 +cm3gn)(1.E+0 +cm4*gn) )
	527	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	528	& *sn(igrid,ilev)
	529	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	530	& *sm(igrid,ilev)
	531	c
	532	c.......................................................................
	533	c
	534	10002 CONTINUE
	535	c
	536	10001 CONTINUE
	537	c
	538	c.......................................................................
	539	c
	540	c
	541	DO igrid=1,ngrid
	542	kn(igrid,1)=knmin
	543	km(igrid,1)=kmmin
	544	q2(igrid,nlev)=q2(igrid,nlev-1)
	545	q(igrid,nlev)=q(igrid,nlev-1)
	546	kn(igrid,nlev)=kn(igrid,nlev-1)
	547	km(igrid,nlev)=km(igrid,nlev-1)
	548	ENDDO
	549
	550	c
	551	c.......................................................................
	552	c
	553	RETURN
	554	END

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

Download in other formats:

Original Format