Context Navigation

cv3p_mixing.F90 @ 354

Last change on this file since 354 was 265, checked in by lfita, 10 years ago
Adding control error values on convection, since with this it has ran!!!
File size: 22.8 KB

Line
1	SUBROUTINE cv3p_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb &
2	& ,ph,t,rr,rs,u,v,tra,h,lv,qnk &
3	& ,unk,vnk,hp,tv,tvp,ep,clw,sig &
4	& ,ment,qent,hent,uent,vent,nent &
5	& ,sigij,elij,supmax,ments,qents,traent)
6	!***************************************************************
7	!* *
8	!* CV3P_MIXING : compute mixed draught properties and, *
9	!* within a scaling factor, mixed draught *
10	!* mass fluxes. *
11	!* written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
12	!* modified by : *
13	!***************************************************************
14	!*
15	implicit none
16	!c
17	#include "cvthermo.h"
18	#include "cv3param.h"
19	#include "YOMCST2.h"
20
21	!c inputs:
22	integer ncum, nd, na, ntra, nloc
23	integer icb(nloc), inb(nloc), nk(nloc)
24	real sig(nloc,nd)
25	real qnk(nloc),unk(nloc),vnk(nloc)
26	real ph(nloc,nd+1)
27	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
28	real u(nloc,nd), v(nloc,nd)
29	real tra(nloc,nd,ntra) ! input of convect3
30	real lv(nloc,na)
31	real h(nloc,na) !liquid water static energy of environment
32	real hp(nloc,na) !liquid water static energy of air shed from adiab. asc.
33	real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
34
35	!c outputs:
36	real ment(nloc,na,na), qent(nloc,na,na)
37	real uent(nloc,na,na), vent(nloc,na,na)
38	real sigij(nloc,na,na), elij(nloc,na,na)
39	real supmax(nloc,na) ! Highest mixing fraction of mixed updraughts
40	! with the sign of (h-hp)
41	real traent(nloc,nd,nd,ntra)
42	real ments(nloc,nd,nd), qents(nloc,nd,nd)
43	real hent(nloc,nd,nd)
44	integer nent(nloc,nd)
45
46	!c local variables:
47	integer i, j, k, il, im, jm
48	integer num1, num2
49	real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
50	real alt, delp, delm
51	real Qmixmax(nloc), Rmixmax(nloc), SQmRmax(nloc)
52	real Qmixmin(nloc), Rmixmin(nloc), SQmRmin(nloc)
53	real signhpmh(nloc)
54	real Sx(nloc), Scrit2
55	real smid(nloc), sjmin(nloc), sjmax(nloc)
56	real Sbef(nloc), Sup(nloc), Smin(nloc)
57	real asij(nloc), smax(nloc), scrit(nloc)
58	real sij(nloc,nd,nd)
59	real csum(nloc,nd)
60	real awat
61	logical lwork(nloc)
62	!c
63	REAL amxupcrit, df, ff
64	INTEGER nstep
65
66	!C
67	!c -- Mixing probability distribution functions
68	!c
69	real Qcoef1,Qcoef2,QFF,QFFF,Qmix,Rmix,Qmix1,Rmix1,Qmix2,Rmix2,F
70	Qcoef1(F) = tanh(F/gammas)
71	Qcoef2(F) = ( tanh(F/gammas) + gammas * &
72	& log(cosh((1.- F)/gammas)/cosh(F/gammas)))
73	QFF(F) = Max(Min(F,1.),0.)
74	QFFF(F) = Min(QFF(F),scut)
75	Qmix1(F) = ( tanh((QFF(F) - Fmax)/gammas)+Qcoef1max )/ &
76	& Qcoef2max
77	Rmix1(F) = ( gammas*log(cosh((QFF(F)-Fmax)/gammas)) &
78	& +QFF(F)*Qcoef1max ) / Qcoef2max
79	Qmix2(F) = -Log(1.-QFFF(F))/scut
80	Rmix2(F) = (QFFF(F)+(1.-QFF(F))*Log(1.-QFFF(F)))/scut
81	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
82	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
83	!C
84	INTEGER, SAVE :: ifrst
85	DATA ifrst/0/
86	!$OMP THREADPRIVATE(ifrst)
87	!C
88
89	! L. Fita, LMD. February 2015.
90	INTEGER :: kl,kl2
91	INTEGER, DIMENSION(nloc) :: Nnodetr
92	CHARACTER(LEN=50) :: errmsg, fname
93
94	errmsg = 'ERROR -- error -- ERROR -- error'
95	fname = 'cv3p_mixing'
96
97	!c=====================================================================
98	!c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
99	!c=====================================================================
100	!c
101	!c -- Initialize mixing PDF coefficients
102	IF (ifrst .EQ. 0) THEN
103	ifrst = 1
104	Qcoef1max = Qcoef1(Fmax)
105	Qcoef2max = Qcoef2(Fmax)
106	!c
107	ENDIF
108	!c
109
110	!c ori do 360 i=1,ncum*nlp
111	do 361 j=1,nl
112	do 360 i=1,ncum
113	nent(i,j)=0
114	!c in convect3, m is computed in cv3_closure
115	!c ori m(i,1)=0.0
116	360 continue
117	361 continue
118
119	!c ori do 400 k=1,nlp
120	!c ori do 390 j=1,nlp
121	do 400 j=1,nl
122	do 390 k=1,nl
123	do 385 i=1,ncum
124	qent(i,k,j)=rr(i,j)
125	uent(i,k,j)=u(i,j)
126	vent(i,k,j)=v(i,j)
127	elij(i,k,j)=0.0
128	hent(i,k,j)=0.0
129	!AC! ment(i,k,j)=0.0
130	!AC! sij(i,k,j)=0.0
131	385 continue
132	390 continue
133	400 continue
134
135	!AC!
136	ment(1:ncum,1:nd,1:nd)=0.0
137	sij(1:ncum,1:nd,1:nd)=0.0
138	!AC!
139
140	do k=1,ntra
141	do j=1,nd ! instead nlp
142	do i=1,nd ! instead nlp
143	do il=1,ncum
144	traent(il,i,j,k)=tra(il,j,k)
145	enddo
146	enddo
147	enddo
148	enddo
149
150	!c=====================================================================
151	!c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
152	!c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
153	!c --- FRACTION (sij)
154	!c=====================================================================
155	! J.Y. Granpeix et L. Fita, LMD Feburary 2015. Counting the number of points without
156	! detrainment
157	Nnodetr = 0
158
159	do 750 i=minorig+1, nl
160
161	do 710 j=minorig,nl
162	do 700 il=1,ncum
163	if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
164	& (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
165
166	rti=qnk(il)-ep(il,i)*clw(il,i)
167	bf2=1.+lv(il,j)lv(il,j)rs(il,j)/(rrvt(il,j)t(il,j)*cpd)
168	anum=h(il,j)-hp(il,i)+(cpv-cpd)t(il,j)(rti-rr(il,j))
169	denom=h(il,i)-hp(il,i)+(cpd-cpv)(rr(il,i)-rti)t(il,j)
170	dei=denom
171	if(abs(dei).lt.0.01)dei=0.01
172	sij(il,i,j)=anum/dei
173	sij(il,i,i)=1.0
174	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
175	altem=altem/bf2
176	cwat=clw(il,j)*(1.-ep(il,j))
177	stemp=sij(il,i,j)
178	if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &
179	& .and.j.gt.i)then
180	anum=anum-lv(il,j)(rti-rs(il,j)-cwatbf2)
181	denom=denom+lv(il,j)*(rr(il,i)-rti)
182	if(abs(denom).lt.0.01)denom=0.01
183	sij(il,i,j)=anum/denom
184	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
185	altem=altem-(bf2-1.)*cwat
186	end if
187	if(sij(il,i,j).gt.0.0)then
188	!ccc ment(il,i,j)=m(il,i)
189	ment(il,i,j)=1.
190	elij(il,i,j)=altem
191	elij(il,i,j)=amax1(0.0,elij(il,i,j))
192	nent(il,i)=nent(il,i)+1
193	endif
194
195	sij(il,i,j)=amax1(0.0,sij(il,i,j))
196	sij(il,i,j)=amin1(1.0,sij(il,i,j))
197	endif ! new
198	700 continue
199	710 continue
200
201	!c
202	!c * if no air can entrain at level i assume that updraft detrains *
203	!c * at that level and calculate detrained air flux and properties *
204	!c
205
206	!c@ do 170 i=icb(il),inb(il)
207
208	do 740 il=1,ncum
209	if ((i.ge.icb(il)).and.(i.le.inb(il)) &
210	& .and.(nent(il,i).eq.0)) then
211	!c@ if(nent(il,i).eq.0)then
212	!ccc ment(il,i,i)=m(il,i)
213	ment(il,i,i)=1.
214	qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
215	uent(il,i,i)=unk(il)
216	vent(il,i,i)=vnk(il)
217	elij(il,i,i)=clw(il,i)*(1.-ep(il,i))
218	sij(il,i,i)=0.0
219	end if
220	740 continue
221	750 continue
222
223	do j=1,ntra
224	do i=minorig+1,nl
225	do il=1,ncum
226	if (i.ge.icb(il) .and. i.le.inb(il) &
227	& .and. nent(il,i).eq.0) then
228	traent(il,i,i,j)=tra(il,nk(il),j)
229	endif
230	enddo
231	enddo
232	enddo
233
234	do 100 j=minorig,nl
235	do 101 i=minorig,nl
236	do 102 il=1,ncum
237	if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
238	& .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
239	sigij(il,i,j)=sij(il,i,j)
240	endif
241	102 continue
242	101 continue
243	100 continue
244	!c@ enddo
245
246	!c@170 continue
247
248	!c=====================================================================
249	!c --- NORMALIZE ENTRAINED AIR MASS FLUXES
250	!c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
251	!c=====================================================================
252
253	call zilch(csum,nloc*nd)
254
255	do il=1,ncum
256	lwork(il) = .FALSE.
257	enddo
258
259	!c---------------------------------------------------------------
260	DO 789 i=minorig+1,nl !Loop on origin level "i"
261	!c---------------------------------------------------------------
262
263	num1=0
264	do il=1,ncum
265	if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
266	enddo
267	if (num1.le.0) goto 789
268
269	!c
270	!cjyg1 Find maximum of SIJ for J>I, if any.
271	!c
272	Sx(:) = 0.
273	!c
274	DO il = 1,ncum
275	IF ( i.ge.icb(il) .and. i.le.inb(il) ) THEN
276	Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
277	Sbef(il) = max(0.,signhpmh(il))
278	ENDIF
279	ENDDO
280
281	DO j = i+1,nl
282	DO il = 1,ncum
283	IF ( i.ge.icb(il) .and. i.le.inb(il) &
284	& .and. j.le.inb(il) ) THEN
285	IF (Sbef(il) .LT. Sij(il,i,j)) THEN
286	Sx(il) = max(Sij(il,i,j),Sx(il))
287	ENDIF
288	Sbef(il) = Sij(il,i,j)
289	ENDIF
290	ENDDO
291	ENDDO
292	!c
293
294	do 781 il=1,ncum
295	if ( i.ge.icb(il) .and. i.le.inb(il) ) then
296	lwork(il)=(nent(il,i).ne.0)
297	qp=qnk(il)-ep(il,i)*clw(il,i)
298	anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) &
299	& +(cpv-cpd)t(il,i)(qp-rr(il,i))
300	denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) &
301	& +(cpd-cpv)t(il,i)(rr(il,i)-qp)
302	if(abs(denom).lt.0.01)denom=0.01
303	scrit(il)=min(anum/denom,1.)
304	alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
305	!c
306	!cjyg1 Find new critical value Scrit2
307	!c such that : Sij > Scrit2 => mixed draught will detrain at J<I
308	!c Sij < Scrit2 => mixed draught will detrain at J>I
309	!c
310	Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il)) &
311	& +Scrit(il)*max(0.,signhpmh(il))
312	!c
313	Scrit(il) = Scrit2
314	!c
315	!cjyg Correction pour la nouvelle logique; la correction pour ALT
316	!c est un peu au hazard
317	if(scrit(il).le.0.0)scrit(il)=0.0
318	if(alt.le.0.0) scrit(il)=1.0
319	!C
320	smax(il)=0.0
321	asij(il)=0.0
322	Sup(il)=0. ! upper S-value reached by descending draughts
323	endif
324	781 continue
325
326	!c---------------------------------------------------------------
327	do 175 j=minorig,nl !Loop on destination level "j"
328	!c---------------------------------------------------------------
329
330	num2=0
331	do il=1,ncum
332	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
333	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
334	& .and. lwork(il) ) num2=num2+1
335	enddo
336	if (num2.le.0) goto 175
337
338	!c -----------------------------------------------
339	IF (j .GT. i) THEN
340	!c -----------------------------------------------
341	do 782 il=1,ncum
342	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
343	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
344	& .and. lwork(il) ) then
345	if(sij(il,i,j).gt.0.0)then
346	Smid(il)=min(Sij(il,i,j),Scrit(il))
347	Sjmax(il)=Smid(il)
348	Sjmin(il)=Smid(il)
349	IF (Smid(il) .LT. Smin(il) .AND. &
350	& Sij(il,i,j+1) .LT. Smid(il)) THEN
351	Smin(il)=Smid(il)
352	Sjmax(il)=min( (Sij(il,i,j+1)+Sij(il,i,j))/2. , &
353	& Sij(il,i,j) , &
354	& Scrit(il) )
355	Sjmin(il)=max( (Sbef(il)+Sij(il,i,j))/2. , &
356	& Sij(il,i,j) )
357	Sjmin(il)=min(Sjmin(il),Scrit(il))
358	Sbef(il) = Sij(il,i,j)
359	ENDIF
360	endif
361	endif
362	782 continue
363	!c -----------------------------------------------
364	ELSE IF (j .EQ. i) THEN
365	!c -----------------------------------------------
366	do 783 il=1,ncum
367	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
368	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
369	& .and. lwork(il) ) then
370	if(sij(il,i,j).gt.0.0)then
371	Smid(il) = 1.
372	Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2.,Scrit(il)) &
373	& *max(0.,-signhpmh(il)) &
374	& +min((Sij(il,i,j+1)+Smid(il))/2.,Scrit(il)) &
375	& *max(0., signhpmh(il))
376	Sjmin(il) = max(Sjmin(il),Sup(il))
377	Sjmax(il) = 1.
378	!c
379	!c- preparation des variables Scrit, Smin et Sbef pour la partie j>i
380	Scrit(il) = min(Sjmin(il),Sjmax(il),Scrit(il))
381
382	Smin(il) = 1.
383	Sbef(il) = max(0.,signhpmh(il))
384	Supmax(il,i) = sign(Scrit(il),-signhpmh(il))
385	endif
386	endif
387	783 continue
388	!c -----------------------------------------------
389	ELSE IF ( j .LT. i) THEN
390	!c -----------------------------------------------
391	do 784 il=1,ncum
392	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
393	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
394	& .and. lwork(il) ) then
395	if(sij(il,i,j).gt.0.0)then
396	Smid(il)=max(Sij(il,i,j),Scrit(il))
397	Sjmax(il) = Smid(il)
398	Sjmin(il) = Smid(il)
399	IF (Smid(il) .GT. Smax(il) .AND. &
400	& Sij(il,i,j+1) .GT. Smid(il)) THEN
401	Smax(il) = Smid(il)
402	Sjmax(il) = max( (Sij(il,i,j+1)+Sij(il,i,j))/2. , &
403	& Sij(il,i,j) )
404	Sjmax(il) = max(Sjmax(il),Scrit(il))
405	Sjmin(il) = min( (Sbef(il)+Sij(il,i,j))/2. , &
406	& Sij(il,i,j) )
407	Sjmin(il) = max(Sjmin(il),Scrit(il))
408	Sbef(il) = Sij(il,i,j)
409	ENDIF
410	IF (abs(Sjmin(il)-Sjmax(il)) .GT. 1.e-10) Sup(il)= &
411	& max(Sjmin(il),Sjmax(il),Sup(il))
412	endif
413	endif
414	784 continue
415	!c -----------------------------------------------
416	END IF
417	!c -----------------------------------------------
418	!c
419	!c
420	do il=1,ncum
421	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
422	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
423	& .and. lwork(il) ) then
424	if(sij(il,i,j).gt.0.0)then
425	rti=qnk(il)-ep(il,i)*clw(il,i)
426	Qmixmax(il)=Qmix(Sjmax(il))
427	Qmixmin(il)=Qmix(Sjmin(il))
428	Rmixmax(il)=Rmix(Sjmax(il))
429	Rmixmin(il)=Rmix(Sjmin(il))
430	SQmRmax(il)= Sjmax(il)*Qmix(Sjmax(il))-Rmix(Sjmax(il))
431	SQmRmin(il)= Sjmin(il)*Qmix(Sjmin(il))-Rmix(Sjmin(il))
432	!c
433	Ment(il,i,j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il,i,j)
434	!c
435	!c Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
436	IF (abs(Qmixmax(il)-Qmixmin(il)) .GT. 1.e-10) THEN
437	Sigij(il,i,j) = &
438	& (SQmRmax(il)-SQmRmin(il))/(Qmixmax(il)-Qmixmin(il))
439	ELSE
440	Sigij(il,i,j) = 0.
441	ENDIF
442	!c
443	!c -- Compute Qent, uent, vent according to the true mixing fraction
444	Qent(il,i,j) = (1.-Sigij(il,i,j))*rti &
445	& + Sigij(il,i,j)*rr(il,i)
446	uent(il,i,j) = (1.-Sigij(il,i,j))*unk(il) &
447	& + Sigij(il,i,j)*u(il,i)
448	vent(il,i,j) = (1.-Sigij(il,i,j))*vnk(il) &
449	& + Sigij(il,i,j)*v(il,i)
450	!c
451	!c-- Compute liquid water static energy of mixed draughts
452	!c IF (j .GT. i) THEN
453	!c awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
454	!c awat=amax1(awat,0.0)
455	!c ELSE
456	!c awat = 0.
457	!c ENDIF
458	!c Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
459	!c : + Sigij(il,i,j)*H(il,i)
460	!c : + (LV(il,j)+(cpd-cpv)t(il,j))awat
461	!IM 301008 beg
462	Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i) &
463	& + Sigij(il,i,j)*H(il,i)
464
465	Elij(il,i,j) = Qent(il,i,j)-rs(il,j)
466	Elij(il,i,j) = Elij(il,i,j) &
467	& + ((h(il,j)-Hent(il,i,j))rs(il,j)LV(il,j) &
468	& / ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv) &
469	& * rrvt(il,j)t(il,j)))
470	Elij(il,i,j) = Elij(il,i,j) &
471	& / (1.+LV(il,j)LV(il,j)rs(il,j) &
472	& / ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv) &
473	& * rrvt(il,j)t(il,j)))
474
475	Elij(il,i,j) = max(elij(il,i,j),0.)
476
477	Elij(il,i,j) = min(elij(il,i,j),Qent(il,i,j))
478
479	IF (j .GT. i) THEN
480	awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
481	awat=amax1(awat,0.0)
482	ELSE
483	awat = 0.
484	ENDIF
485
486	!c print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
487	!c : t(il,j))
488
489	Hent(il,i,j) = Hent(il,i,j)+(LV(il,j)+(cpd-cpv)*t(il,j)) &
490	& * awat
491	!IM 301008 end
492	!c
493	!c print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
494	!c : i,j,hent(il,i,j),sigij(il,i,j)
495	!c
496	!c -- ASij is the integral of P(F) over the relevant F interval
497	ASij(il) = ASij(il) &
498	& + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) &
499	& -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
500	!c
501	endif
502	endif
503	enddo
504	do k=1,ntra
505	do il=1,ncum
506	if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
507	& (j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
508	& .and. lwork(il) ) then
509	if(sij(il,i,j).gt.0.0)then
510	traent(il,i,j,k)=sigij(il,i,j)*tra(il,i,k) &
511	& +(1.-sigij(il,i,j))*tra(il,nk(il),k)
512	endif
513	endif
514	enddo
515	enddo
516	!c
517	!c -- If I=J (detrainement and entrainement at the same level), then only the
518	!c -- adiabatic ascent part of the mixture is considered
519	IF (I .EQ. J) THEN
520	do il=1,ncum
521	if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
522	& j.ge.(icb(il)-1) .and. j.le.inb(il) &
523	& .and. lwork(il) ) then
524	if(sij(il,i,j).gt.0.0)then
525	rti=qnk(il)-ep(il,i)*clw(il,i)
526	!ccc Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
527	Ment(il,i,i) = abs(Qmixmax(il)*(1.-Sjmax(il)) &
528	& +Rmixmax(il) &
529	& -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
530	Qent(il,i,i) = rti
531	uent(il,i,i) = unk(il)
532	vent(il,i,i) = vnk(il)
533	Hent(il,i,i) = hp(il,i)
534	Elij(il,i,i) = clw(il,i)*(1.-ep(il,i))
535	Sigij(il,i,i) = 0.
536	endif
537	endif
538	enddo
539	do k=1,ntra
540	do il=1,ncum
541	if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
542	& (j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
543	& .and. lwork(il) ) then
544	if(sij(il,i,j).gt.0.0)then
545	traent(il,i,i,k)=tra(il,nk(il),k)
546	endif
547	endif
548	enddo
549	enddo
550	!c
551	ENDIF
552	!c
553	175 continue
554
555	do il=1,ncum
556	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
557	asij(il)=amax1(1.0e-16,asij(il))
558	asij(il)=1.0/asij(il)
559	csum(il,i)=0.0
560	endif
561	enddo
562
563	do 180 j=minorig,nl
564	do il=1,ncum
565	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
566	& .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
567	! J.Y. Grandpeix, Janury 2015 LMD. No level of detrainment found. No mixing!
568	IF (asij(il) < 100.) THEN
569	ment(il,i,j)=ment(il,i,j)*asij(il)
570	ELSE
571	ment(il,i,j)=0.
572	Nnodetr(il) = Nnodetr(il) + 1
573	ENDIF
574	! L. Fita, LMD. Feburary 2015, Checkings...
575	IF (ment(il,i,j) > 10.) THEN
576	PRINT *,TRIM(errmsg)
577	PRINT ,' ' // TRIM(fname) // ' * after asij '
578	PRINT *,' ' // TRIM(fname) // ': wrong ment= ', ment(il,i,j), &
579	' value at ', il,', ',i,', ',j, '! (< 10.) !!'
580	PRINT *,' ment(i,j) asij _________________'
581	PRINT *,ment(il,i,j),asij(il)
582	END IF
583	endif
584	enddo
585	180 continue
586
587	do 197 j=minorig,nl
588	do il=1,ncum
589	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
590	& .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
591	csum(il,i)=csum(il,i)+ment(il,i,j)
592	endif
593	enddo
594	197 continue
595
596	do il=1,ncum
597	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
598	& .and. csum(il,i).lt.1. ) then
599	!ccc : .and. csum(il,i).lt.m(il,i) ) then
600	nent(il,i)=0
601	!ccc ment(il,i,i)=m(il,i)
602	ment(il,i,i)=1.
603	qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
604	uent(il,i,i)=unk(il)
605	vent(il,i,i)=vnk(il)
606	elij(il,i,i)=clw(il,i)*(1.-ep(il,i))
607	sij(il,i,i)=0.0
608	endif
609	enddo ! il
610
611	do j=1,ntra
612	do il=1,ncum
613	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
614	& .and. csum(il,i).lt.1. ) then
615	!ccc : .and. csum(il,i).lt.m(il,i) ) then
616	traent(il,i,i,j)=tra(il,nk(il),j)
617	endif
618	enddo
619	enddo
620	!c
621	789 continue
622	!c
623
624	! J.Y. Granpeix et L. Fita, LMD Feburary 2015. Counting the number of points without
625	! detrainment
626	DO il=1,ncum
627	IF (Nnodetr(il) > 0) THEN
628	PRINT *,' '//TRIM(fname)// ': No level of detrainment found. No mixing!'
629	PRINT *,' at point: ',il,' no detrainment found at ',Nnodetr(il), &
630	' levels!'
631	END IF
632	END DO
633
634	return
635	end
636

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

Download in other formats:

Original Format