Context Navigation

vdif_kc.F @ 706

Last change on this file since 706 was 325, checked in by acolaitis, 13 years ago
Included variation of mean molar mass in potential temperature definition: this modification realistically accounts for mixing in the polar night when the atm is enriched in Ar. This follows present computations in convadj for the polar night case.
File size: 21.5 KB

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

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

Download in other formats:

Original Format