Context Navigation

vdif_kc.F @ 3523

Last change on this file since 3523 was 3236, checked in by gmilcareck, 9 months ago
Add virtual correction for convective adjustment and turbulent diffusion (turbdiff and vdifc) + correction of allocated tables in physiq_mod for varspec.
File size: 19.5 KB

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

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

Download in other formats:

Original Format