Context Navigation

yamada4.F @ 1242

Last change on this file since 1242 was 1224, checked in by aslmd, 11 years ago

LMDZ.MARS. If not major, a quite important commit.

No more SAVE,ALLOCATABLE arrays outside modules.

This is important to solve the nesting conundrum in MESOSCALE.
And overall this is good for the harmony of the universe.
(Joke apart, this is good for any interfacing task. And compliant with a F90 spirit).
Note that bit-to-bit compatibility of results in debug mode was checked.

inifis is split in two : phys_state_var_init + conf_phys

This makes interfacing with MESOSCALE more transparent.
This is also clearer for LMDZ.MARS.
Before, inifis has two very different tasks to do.

a bit of cleaning as far as modules and saves are concerned

Point 1

Removed SAVE,ALLOCATABLE arrays from

physiq, aeropacity, updatereffrad, soil

and put those in

dimradmars_mod, surfdat_h, tracer_mod, comsoil_h

and changed accordingly the initialization subroutines associated to each module.
Allocating these arrays is thus done at initialization.

Point 2

Created a subroutine phys_state_var_init which does all the allocation / initialization work for modules. This was previously done in inifis.
Replaced inifis which was then (after the previous modification) just about reading callphys.def and setting a few constants by conf_phys. This mimics the new LMDZ terminology (cf. LMDZ.VENUS for instance)
Bye bye inifis.

Point 3

Removed comdiurn and put everything in comgeomfi
Created a turb_mod module for turbulence variables (e.g. l0 in yamada4)
dryness had nothing to do in tracer_h, put it in surfdat_h (like watercaptag)
topdust0 does not need to be SAVE in aeropacity. better use sinlat.
emisref does not need to be SAVE in newcondens. made it automatic array.
Removed useless co2ice argument in initracer.

File size: 23.6 KB

Line
1	!************************************************************
2	!************************************************************
3	!
4	! YAMADA4 EARTH =>>> MARS VERSION
5	! Modifications by: A.C. 02-03-2012 (marked by 'MARS')
6	! Original version given by F.H. 01-03-2012
7	!
8	!************************************************************
9	!************************************************************
10	SUBROUTINE yamada4(ngrid,nlay,nq,dt,g,rconst,plev,temp
11	s ,zlev,zlay,u,v,phc,pq,cd,q2,km,kn,kq,ustar
12	s ,iflag_pbl)
13	use tracer_mod, only: noms
14	use turb_mod, only: l0
15	IMPLICIT NONE
16	!.......................................................................
17	! MARS
18	#include "dimensions.h"
19	#include "dimphys.h"
20	!#include "tracer.h"
21	#include "callkeys.h"
22	!.......................................................................
23	!
24	! dt : pas de temps
25	! g : g
26	! zlev : altitude a chaque niveau (interface inferieure de la couche
27	! de meme indice)
28	! zlay : altitude au centre de chaque couche
29	! u,v : vitesse au centre de chaque couche
30	! (en entree : la valeur au debut du pas de temps)
31	! phc : temperature potentielle au centre de chaque couche
32	! (en entree : la valeur au debut du pas de temps)
33	! cd : cdrag
34	! (en entree : la valeur au debut du pas de temps)
35	! q2 : $q^2$ au bas de chaque couche
36	! (en entree : la valeur au debut du pas de temps)
37	! (en sortie : la valeur a la fin du pas de temps)
38	! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
39	! couche)
40	! (en sortie : la valeur a la fin du pas de temps)
41	! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
42	! (en sortie : la valeur a la fin du pas de temps)
43	!
44	! iflag_pbl doit valoir entre 6 et 9
45	! l=6, on prend systematiquement une longueur d equilibre
46	! iflag_pbl=6 : MY 2.0
47	! iflag_pbl=7 : MY 2.0.Fournier
48	! iflag_pbl=8 : MY 2.5
49	! iflag_pbl>=9 : MY 2.5 avec diffusion verticale
50	!
51	!.......................................................................
52
53	REAL, INTENT(IN) :: dt,g,rconst
54	REAL, INTENT(IN) :: u(ngrid,nlay)
55	REAL, INTENT(IN) :: v(ngrid,nlay)
56	REAL, INTENT(IN) :: phc(ngrid,nlay)
57	REAL, INTENT(IN) :: cd(ngrid)
58	REAL, INTENT(IN) :: temp(ngrid,nlay)
59	REAL, INTENT(IN) :: plev(ngrid,nlay+1)
60	REAL, INTENT(IN) :: ustar(ngrid)
61	REAL, INTENT(IN) :: zlev(ngrid,nlay+1)
62	REAL, INTENT(IN) :: zlay(ngrid,nlay)
63	INTEGER, INTENT(IN) :: iflag_pbl,ngrid
64	INTEGER, INTENT(IN) :: nlay
65	INTEGER, INTENT(IN) :: nq
66	REAL, INTENT(INOUT) :: q2(ngrid,nlay+1)
67	REAL, INTENT(OUT) :: km(ngrid,nlay+1)
68	REAL, INTENT(OUT) :: kn(ngrid,nlay+1)
69	REAL, INTENT(OUT) :: kq(ngrid,nlay+1)
70
71	REAL kmin,qmin,pblhmin(ngrid),coriol(ngrid)
72	REAL unsdz(ngrid,nlay)
73	REAL unsdzdec(ngrid,nlay+1)
74	REAL kmpre(ngrid,nlay+1),tmp2
75	REAL mpre(ngrid,nlay+1)
76	REAL ff(ngrid,nlay+1),delta(ngrid,nlay+1)
77	REAL aa(ngrid,nlay+1),aa0,aa1,qpre
78
79	LOGICAL first
80	INTEGER ipas,nlev
81	SAVE first,ipas
82	!FH/IM DATA first,ipas/.true.,0/
83	DATA first,ipas/.false.,0/
84	INTEGER ig,k
85
86
87	REAL ri,zrif,zalpha,zsm,zsn
88	REAL rif(ngrid,nlay+1),sm(ngrid,nlay+1),alpha(ngrid,nlay)
89
90	REAL m2(ngrid,nlay+1),dz(ngrid,nlay+1),zq,n2(ngrid,nlay+1)
91	REAL dtetadz(ngrid,nlay+1)
92	REAL m2cstat,mcstat,kmcstat
93	REAL l(ngrid,nlay+1)
94	REAL sq(ngrid),sqz(ngrid),zz(ngrid,nlay+1)
95	INTEGER iter
96
97	REAL ric,rifc,b1,kap
98	SAVE ric,rifc,b1,kap
99	DATA ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
100	REAL frif,falpha,fsm
101	REAL fl,zzz,zl0,zq2,zn2
102
103	REAL rino(ngrid,nlay+1),smyam(ngrid,nlay),styam(ngrid,nlay)
104	s ,lyam(ngrid,nlay),knyam(ngrid,nlay)
105	s ,w2yam(ngrid,nlay),t2yam(ngrid,nlay)
106	LOGICAL,SAVE :: firstcall=.true.
107	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
108	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
109	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
110	fl(zzz,zl0,zq2,zn2)=
111	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
112	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
113
114
115	! MARS
116	REAL,SAVE :: q2min,q2max,knmin,kmmin
117	DATA q2min,q2max,knmin,kmmin/1.E-10,1.E+2,1.E-5,1.E-5/
118	INTEGER ico2,iq
119	SAVE ico2
120	REAL m_co2, m_noco2, A , B
121	SAVE A, B
122	REAL teta(ngrid,nlay)
123	REAL pq(ngrid,nlay,nq)
124	REAL kminfact
125	INTEGER i
126	REAL ztimestep
127	INTEGER, SAVE :: ndt
128
129	nlev=nlay+1
130
131	c.......................................................................
132	c Initialization
133	c.......................................................................
134
135	if(firstcall) then
136	ico2=0
137	if (tracer) then
138	! Prepare Special treatment if one of the tracers is CO2 gas
139	do iq=1,nq
140	if (noms(iq).eq."co2") then
141	ico2=iq
142	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
143	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
144	! Compute A and B coefficient use to compute
145	! mean molecular mass Mair defined by
146	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
147	! 1/Mair = A*q(ico2) + B
148	A =(1/m_co2 - 1/m_noco2)
149	B=1/m_noco2
150	end if
151	enddo
152	endif
153	ndt=ceiling(3840./(3699.*24./dt))
154	firstcall=.false.
155	endif !of if firstcall
156
157	c.......................................................................
158	c Special treatment for co2
159	c.......................................................................
160
161	! if (ico2.ne.0) then
162	!! Special case if one of the tracers is CO2 gas
163	! DO k=1,nlay
164	! DO ig=1,ngrid
165	! teta(ig,k) = phc(ig,k)(Apq(ig,k,ico2)+B)
166	! ENDDO
167	! ENDDO
168	! else
169	teta(:,:)=phc(:,:)
170	! end if
171
172	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
173	stop'probleme de coherence dans appel a MY'
174	endif
175
176	ipas=ipas+1
177	! MARS
178	! if (0.eq.1.and.first) then
179	! do ig=1,1000
180	! ri=(ig-800.)/500.
181	! if (ri.lt.ric) then
182	! zrif=frif(ri)
183	! else
184	! zrif=rifc
185	! endif
186	! if(zrif.lt.0.16) then
187	! zalpha=falpha(zrif)
188	! zsm=fsm(zrif)
189	! else
190	! zalpha=1.12
191	! zsm=0.085
192	! endif
193	! print*,ri,rif,zalpha,zsm
194	! enddo
195	! endif
196
197	!.......................................................................
198	! les increments verticaux
199	!.......................................................................
200	!
201	!!!!!! allerte !!!!!
202	!!!!!! zlev n'est pas declare a nlev !!!!!
203	!!!!!! ---->
204	! MARS
205	!
206	! DO ig=1,ngrid
207	! zlev(ig,nlev)=zlay(ig,nlay)
208	! & +( zlay(ig,nlay) - zlev(ig,nlev-1) )
209	! ENDDO
210	!!!!! <----
211	!!!!! allerte !!!!!
212
213	DO k=1,nlay
214	DO ig=1,ngrid
215	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
216	ENDDO
217	ENDDO
218	DO ig=1,ngrid
219	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
220	ENDDO
221	DO k=2,nlay
222	DO ig=1,ngrid
223	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
224	ENDDO
225	ENDDO
226	DO ig=1,ngrid
227	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
228	ENDDO
229	!
230	!.......................................................................
231
232	do k=2,nlay
233	do ig=1,ngrid
234	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
235	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
236	s /(dz(ig,k)*dz(ig,k))
237	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
238	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
239	! n2(ig,k)=0.
240	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
241	if (ri.lt.ric) then
242	rif(ig,k)=frif(ri)
243	else
244	rif(ig,k)=rifc
245	endif
246	if(rif(ig,k).lt.0.16) then
247	alpha(ig,k)=falpha(rif(ig,k))
248	sm(ig,k)=fsm(rif(ig,k))
249	else
250	alpha(ig,k)=1.12
251	sm(ig,k)=0.085
252	endif
253	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
254	! print,'RIF L=',k,rif(ig,k),rialpha(ig,k)
255
256
257	enddo
258	enddo
259
260
261	!====================================================================
262	! Au premier appel, on determine l et q2 de facon iterative.
263	! iterration pour determiner la longueur de melange
264
265
266	if (first.or.iflag_pbl.eq.6) then
267	do ig=1,ngrid
268	! MARS
269	! l0(ig)=10.
270	l0(ig)=160.
271	enddo
272	do k=2,nlay-1
273	do ig=1,ngrid
274	l(ig,k)=l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
275	enddo
276	enddo
277
278	do iter=1,10
279	do ig=1,ngrid
280	sq(ig)=1.e-10
281	sqz(ig)=1.e-10
282	enddo
283	do k=2,nlay-1
284	do ig=1,ngrid
285	q2(ig,k)=l(ig,k)*2zz(ig,k)
286	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
287	zq=sqrt(q2(ig,k))
288	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
289	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
290	enddo
291	enddo
292	do ig=1,ngrid
293	l0(ig)=0.2*sqz(ig)/sq(ig)
294	! l0(ig)=30.
295	enddo
296	! print*,'ITER=',iter,' L0=',l0
297
298	enddo
299
300	! print*,'Fin de l initialisation de q2 et l0'
301
302	endif ! first
303
304	!====================================================================
305	! Calcul de la longueur de melange.
306	!====================================================================
307
308	! Mise a jour de l0
309	do ig=1,ngrid
310	sq(ig)=1.e-10
311	sqz(ig)=1.e-10
312	enddo
313	do k=2,nlay-1
314	do ig=1,ngrid
315	zq=sqrt(q2(ig,k))
316	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
317	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
318	enddo
319	enddo
320	do ig=1,ngrid
321	l0(ig)=0.2*sqz(ig)/sq(ig)
322	! l0(ig)=30.
323	enddo
324	! print*,'ITER=',iter,' L0=',l0
325	! calcul de l(z)
326	do k=2,nlay
327	do ig=1,ngrid
328	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
329	if(first) then
330	q2(ig,k)=l(ig,k)*2zz(ig,k)
331	endif
332	enddo
333	enddo
334
335	!====================================================================
336	! Yamada 2.0
337	!====================================================================
338	if (iflag_pbl.eq.6) then
339
340	do k=2,nlay
341	do ig=1,ngrid
342	q2(ig,k)=l(ig,k)*2zz(ig,k)
343	enddo
344	enddo
345
346
347	else if (iflag_pbl.eq.7) then
348	!====================================================================
349	! Yamada 2.Fournier
350	!====================================================================
351
352	! Calcul de l, km, au pas precedent
353	do k=2,nlay
354	do ig=1,ngrid
355	c print*,'SMML=',sm(ig,k),l(ig,k)
356	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
357	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
358	mpre(ig,k)=sqrt(m2(ig,k))
359	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
360	enddo
361	enddo
362
363	do k=2,nlay-1
364	do ig=1,ngrid
365	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
366	mcstat=sqrt(m2cstat)
367
368	! print*,'M2 L=',k,mpre(ig,k),mcstat
369	!
370	! -----{puis on ecrit la valeur de q qui annule l'equation de m
371	! supposee en q3}
372	!
373	IF (k.eq.2) THEN
374	kmcstat=1.E+0 / mcstat
375	& ( unsdz(ig,k)kmpre(ig,k+1)
376	& *mpre(ig,k+1)
377	& +unsdz(ig,k-1)
378	& *cd(ig)
379	& ( sqrt(u(ig,3)2+v(ig,3)*2)
380	& -mcstat/unsdzdec(ig,k)
381	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
382	& /( unsdz(ig,k)+unsdz(ig,k-1) )
383	ELSE
384	kmcstat=1.E+0 / mcstat
385	& ( unsdz(ig,k)kmpre(ig,k+1)
386	& *mpre(ig,k+1)
387	& +unsdz(ig,k-1)*kmpre(ig,k-1)
388	& *mpre(ig,k-1) )
389	& /( unsdz(ig,k)+unsdz(ig,k-1) )
390	ENDIF
391	! print*,'T2 L=',k,tmp2
392	tmp2=kmcstat
393	& /( sm(ig,k)/q2(ig,k) )
394	& /l(ig,k)
395
396	! MARS
397	! q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
398	q2(ig,k)=max(q2min,max(tmp2,1.e-12)**(2./3.))
399
400	! print*,'Q2 L=',k,q2(ig,k)
401	!
402	enddo
403	enddo
404
405	else if (iflag_pbl.ge.8) then
406	!====================================================================
407	! Yamada 2.5 a la Didi
408	!====================================================================
409
410	ztimestep=dt/real(ndt)
411	do i=1,ndt
412
413	! Calcul de l, km, au pas precedent
414	do k=2,nlay
415	do ig=1,ngrid
416	! print*,'SMML=',sm(ig,k),l(ig,k)
417	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
418	if (delta(ig,k).lt.1.e-20) then
419	! print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
420	delta(ig,k)=1.e-20
421	endif
422	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
423	aa0=
424	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
425	aa1=
426	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
427	! abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
428	aa(ig,k)=aa1ztimestep/(delta(ig,k)l(ig,k))
429	! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
430	qpre=sqrt(q2(ig,k))
431	if (iflag_pbl.eq.8 ) then
432	if (aa(ig,k).gt.0.) then
433	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
434	else
435	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
436	endif
437	else ! iflag_pbl=9
438	if (aa(ig,k)*qpre.gt.0.9) then
439	q2(ig,k)=(qpre10.)*2
440	else
441	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
442	endif
443	endif
444
445	! MARS
446	q2(ig,k)=min(max(q2(ig,k),q2min),q2max)
447	! q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
448
449	! print,'Q2 L=',k,q2(ig,k),qpreqpre
450	enddo
451	enddo
452
453	! MARS
454	q2(:,nlay+1)=q2(:,nlay)
455
456	if (iflag_pbl .eq. 9) then
457	do k=2,nlay
458	do ig=1,ngrid
459	zq=sqrt(q2(ig,k))
460	km(ig,k)=l(ig,k)zqsm(ig,k)
461	kn(ig,k)=km(ig,k)*alpha(ig,k)
462	kq(ig,k)=l(ig,k)zq0.2
463	enddo
464	enddo
465	! boundary conditions for km
466	km(:,nlay+1)=0
467	km(:,1)=km(:,2) ! km(:,1)=0
468	! boundary conditions for kn
469	kn(:,nlay+1)=0
470	kn(:,1)=kn(:,2) ! kn(:,1)=0
471	! boundary conditions for kq
472	kq(:,nlay+1)=0 ! zero at top of atmosphere
473	kq(:,1)=kq(:,2) ! no gradient at surface
474
475	q2(:,1)=q2(:,2)
476	call vdif_q2(ztimestep,g,rconst,ngrid,nlay,plev,temp,kq,q2)
477
478	endif ! of if iflag_pbl eq 9
479
480	enddo !of i=1,ndt
481
482	endif ! Fin du cas 8
483
484	! print*,'OK8'
485
486	!====================================================================
487	! Calcul des coefficients de melange
488	!====================================================================
489	if (iflag_pbl .ne. 9) then
490	do k=2,nlay
491	! print*,'k=',k
492	do ig=1,ngrid
493	!abde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
494	zq=sqrt(q2(ig,k))
495	km(ig,k)=l(ig,k)zqsm(ig,k)
496	kn(ig,k)=km(ig,k)*alpha(ig,k)
497	kq(ig,k)=l(ig,k)zq0.2
498	! print*,'KML=',km(ig,k),kn(ig,k)
499	enddo
500	enddo
501
502	! MARS
503	km(:,nlay+1)=km(:,nlay)
504	kn(:,nlay+1)=kn(:,nlay)
505	kq(:,nlay+1)=kq(:,nlay)
506
507	! Transport diffusif vertical de la TKE.
508	! if (iflag_pbl.ge.9) then
509	!! print*,'YAMADA VDIF'
510	! q2(:,1)=q2(:,2)
511	! call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
512	! endif
513
514	endif
515
516	! Traitement des cas noctrunes avec l'introduction d'une longueur
517	! minilale.
518	!
519	!====================================================================
520	! Traitement particulier pour les cas tres stables.
521	! D'apres Holtslag Boville.
522
523	! MARS
524	! callkmin=.true.
525	! call getin("callkmin",callkmin)
526	! IF (callkmin) THEN
527	do ig=1,ngrid
528	! coriol(ig)=1.e-4
529	! pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
530
531	if (ngrid .eq. 1) then
532	kminfact=0.3
533	else
534	kminfact=0.45
535	endif
536
537	pblhmin(ig)=kminfact0.07MAX(ustar(ig),1.e-3)/1.e-4
538	enddo
539	! print*,'pblhmin ',pblhmin
540	!CTest a remettre 21 11 02
541	! test abd 13 05 02 if(0.eq.1) then
542	! if(0.eq.1) then
543	do k=2,nlay
544	do ig=1,ngrid
545	if (teta(ig,2).gt.teta(ig,1)) then
546	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
547	! kmin=kapzlev(ig,k)qmin
548	kmin=fl(zlev(ig,k),l0(ig),qmin*2,n2(ig,k))qmin
549	else
550	kmin=-1. ! kmin n'est utilise que pour les SL stables.
551	endif
552	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
553	! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
554	! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
555	! kn(ig,k)=kmin
556	! km(ig,k)=kmin
557	! kq(ig,k)=kmin
558
559	kn(ig,k)=kmin*alpha(ig,k)
560	km(ig,k)=kmin
561	kq(ig,k)=kmin*0.2
562	! la longueur de melange est suposee etre l= kap z
563	! K=l q Sm d'ou q2=(K/l Sm)**2
564	! q2(ig,k)=(qmin/sm(ig,k))**2
565	q2(ig,k)=(kmin/
566	& (fl(zlev(ig,k),l0(ig),qmin*2,n2(ig,k))sm(ig,k)))**2
567	endif
568	enddo
569	enddo
570	! endif
571
572	! ENDIF
573
574	! Diagnostique pour stokage
575
576	if(1.eq.0)then
577	rino=rif
578	smyam(1:ngrid,1)=0.
579	styam(1:ngrid,1)=0.
580	lyam(1:ngrid,1)=0.
581	knyam(1:ngrid,1)=0.
582	w2yam(1:ngrid,1)=0.
583	t2yam(1:ngrid,1)=0.
584
585	smyam(1:ngrid,2:nlay)=sm(1:ngrid,2:nlay)
586	styam(1:ngrid,2:nlay)=sm(1:ngrid,2:nlay)*alpha(1:ngrid,2:nlay)
587	lyam(1:ngrid,2:nlay)=l(1:ngrid,2:nlay)
588	knyam(1:ngrid,2:nlay)=kn(1:ngrid,2:nlay)
589
590	! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
591
592	w2yam(1:ngrid,2:nlay)=q2(1:ngrid,2:nlay)*0.24
593	s +lyam(1:ngrid,2:nlay)5.17kn(1:ngrid,2:nlay)
594	s *n2(1:ngrid,2:nlay)/sqrt(q2(1:ngrid,2:nlay))
595
596	t2yam(1:ngrid,2:nlay)=9.1*kn(1:ngrid,2:nlay)
597	s dtetadz(1:ngrid,2:nlay)*2
598	s /sqrt(q2(1:ngrid,2:nlay))*lyam(1:ngrid,2:nlay)
599	endif
600
601	! print*,'OKFIN'
602	first=.false.
603	return
604	end
605	SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,nlay
606	& ,plev,temp,kmy,q2)
607	IMPLICIT NONE
608	!.......................................................................
609	! MARS
610	#include "dimensions.h"
611	#include "dimphys.h"
612	!#include "tracer.h"
613	#include "callkeys.h"
614	!.......................................................................
615	!
616	! dt : pas de temps
617	!
618	REAL plev(ngrid,nlay+1)
619	REAL temp(ngrid,nlay)
620	REAL timestep
621	REAL gravity,rconst
622	REAL kstar(ngrid,nlay+1),zz
623	REAL kmy(ngrid,nlay+1)
624	REAL q2(ngrid,nlay+1)
625	REAL deltap(ngrid,nlay+1)
626	REAL denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
627	INTEGER ngrid,nlay
628
629	INTEGER i,k
630
631	! print*,'RD=',rconst
632	do k=1,nlay
633	do i=1,ngrid
634	! test
635	! print*,'i,k',i,k
636	! print*,'temp(i,k)=',temp(i,k)
637	! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
638	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
639	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
640	s /(plev(i,k)-plev(i,k+1))*timestep
641	enddo
642	enddo
643
644	do k=2,nlay
645	do i=1,ngrid
646	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
647	enddo
648	enddo
649	do i=1,ngrid
650	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
651	deltap(i,nlay+1)=0.5*(plev(i,nlay)-plev(i,nlay+1))
652	denom(i,nlay+1)=deltap(i,nlay+1)+kstar(i,nlay)
653	alpha(i,nlay+1)=deltap(i,nlay+1)*q2(i,nlay+1)/denom(i,nlay+1)
654	beta(i,nlay+1)=kstar(i,nlay)/denom(i,nlay+1)
655	enddo
656
657	do k=nlay,2,-1
658	do i=1,ngrid
659	denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
660	s kstar(i,k)+kstar(i,k-1)
661	! correction d'un bug 10 01 2001
662	alpha(i,k)=(q2(i,k)*deltap(i,k)
663	s +kstar(i,k)*alpha(i,k+1))/denom(i,k)
664	beta(i,k)=kstar(i,k-1)/denom(i,k)
665	enddo
666	enddo
667
668	! Si on recalcule q2(1)
669	if(1.eq.0) then
670	do i=1,ngrid
671	denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
672	q2(i,1)=(q2(i,1)*deltap(i,1)
673	s +kstar(i,1)*alpha(i,2))/denom(i,1)
674	enddo
675	endif
676	! sinon, on peut sauter cette boucle...
677
678	do k=2,nlay+1
679	do i=1,ngrid
680	q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
681	enddo
682	enddo
683
684	return
685	end
686	SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,nlay,
687	& plev,temp,kmy,q2)
688	IMPLICIT NONE
689	!.......................................................................
690	! MARS
691	#include "dimensions.h"
692	#include "dimphys.h"
693	!#include "tracer.h"
694	#include "callkeys.h"
695	!.......................................................................
696	!
697	! dt : pas de temps
698
699	REAL plev(ngrid,nlay+1)
700	REAL temp(ngrid,nlay)
701	REAL timestep
702	REAL gravity,rconst
703	REAL kstar(ngrid,nlay+1),zz
704	REAL kmy(ngrid,nlay+1)
705	REAL q2(ngrid,nlay+1)
706	REAL deltap(ngrid,nlay+1)
707	REAL denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
708	INTEGER ngrid,nlay
709
710	INTEGER i,k
711
712	do k=1,nlay
713	do i=1,ngrid
714	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
715	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
716	s /(plev(i,k)-plev(i,k+1))*timestep
717	enddo
718	enddo
719
720	do k=2,nlay
721	do i=1,ngrid
722	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
723	enddo
724	enddo
725	do i=1,ngrid
726	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
727	deltap(i,nlay+1)=0.5*(plev(i,nlay)-plev(i,nlay+1))
728	enddo
729
730	do k=nlay,2,-1
731	do i=1,ngrid
732	q2(i,k)=q2(i,k)+
733	s ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
734	s -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
735	s /deltap(i,k)
736	enddo
737	enddo
738
739	do i=1,ngrid
740	q2(i,1)=q2(i,1)+
741	s ( kstar(i,1)*(q2(i,2)-q2(i,1))
742	s )
743	s /deltap(i,1)
744	q2(i,nlay+1)=q2(i,nlay+1)+
745	s (
746	s -kstar(i,nlay)*(q2(i,nlay+1)-q2(i,nlay)) )
747	s /deltap(i,nlay+1)
748	enddo
749
750	return
751	end

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

Download in other formats:

Original Format