source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/cv3p_mixing.F90 @ 3069

Last change on this file since 3069 was 2056, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1997:2055 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.1 KB
Line 
1SUBROUTINE 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
17  include "cvthermo.h"
18  include "cv3param.h"
19  include "YOMCST2.h"
20
21!inputs:
22  INTEGER, INTENT (IN)                               :: ncum, nd, na
23  INTEGER, INTENT (IN)                               :: ntra, nloc
24  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
25  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
26  REAL, DIMENSION (nloc), INTENT (IN)                :: qnk, unk, vnk
27  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
28  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
29  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
30  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
31  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
32  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environMent
33  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
34  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
35  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
36
37!outputs:
38  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
39  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
40  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
41  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax(nloc, na) ! Highest mixing fraction of mixed
42                                                                         ! updraughts with the sign of (h-hp)
43  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
44  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
45  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
46  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
47
48!local variables:
49  INTEGER i, j, k, il, im, jm
50  INTEGER num1, num2
51  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
52  REAL                               :: alt, delp, delm
53  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
54  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
55  REAL, DIMENSION (nloc)             :: signhpmh
56  REAL, DIMENSION (nloc)             :: Sx
57  REAL                               :: Scrit2
58  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
59  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
60  REAL, DIMENSION (nloc)             :: ASij, smax, Scrit
61  REAL, DIMENSION (nloc, nd, nd)     :: Sij
62  REAL, DIMENSION (nloc, nd)         :: csum
63  REAL                               :: awat
64  LOGICAL, DIMENSION (nloc)          :: lwork
65
66  REAL amxupcrit, df, ff
67  INTEGER nstep
68
69! --   Mixing probability distribution functions
70
71  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
72
73  Qcoef1(F) = tanh(F/gammas)
74  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
75  QFF(F) = max(min(F,1.), 0.)
76  QFFf(F) = min(QFF(F), scut)
77  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
78  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+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) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
82  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
83
84  INTEGER, SAVE :: ifrst
85  DATA ifrst/0/
86!$OMP THREADPRIVATE(ifrst)
87
88
89! =====================================================================
90! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
91! =====================================================================
92
93! -- Initialize mixing PDF coefficients
94  IF (ifrst==0) THEN
95    ifrst = 1
96    Qcoef1max = Qcoef1(Fmax)
97    Qcoef2max = Qcoef2(Fmax)
98
99  END IF
100
101
102! ori        do 360 i=1,ncum*nlp
103  DO j = 1, nl
104    DO i = 1, ncum
105      nent(i, j) = 0
106! in convect3, m is computed in cv3_closure
107! ori          m(i,1)=0.0
108    END DO
109  END DO
110
111! ori      do 400 k=1,nlp
112! ori       do 390 j=1,nlp
113  DO j = 1, nl
114    DO k = 1, nl
115      DO i = 1, ncum
116        Qent(i, k, j) = rr(i, j)
117        uent(i, k, j) = u(i, j)
118        vent(i, k, j) = v(i, j)
119        elij(i, k, j) = 0.0
120        hent(i, k, j) = 0.0
121!AC!            Ment(i,k,j)=0.0
122!AC!            Sij(i,k,j)=0.0
123      END DO
124    END DO
125  END DO
126
127!AC!
128  Ment(1:ncum, 1:nd, 1:nd) = 0.0
129  Sij(1:ncum, 1:nd, 1:nd) = 0.0
130!AC!
131
132  DO k = 1, ntra
133    DO j = 1, nd ! instead nlp
134      DO i = 1, nd ! instead nlp
135        DO il = 1, ncum
136          traent(il, i, j, k) = tra(il, j, k)
137        END DO
138      END DO
139    END DO
140  END DO
141
142! =====================================================================
143! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
144! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
145! --- FRACTION (Sij)
146! =====================================================================
147
148  DO i = minorig + 1, nl
149
150    DO j = minorig, nl
151      DO il = 1, ncum
152        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
153                         .AND. (j<=inb(il))) THEN
154
155          rti = qnk(il) - ep(il, i)*clw(il, i)
156          bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
157          anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
158          denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
159          dei = denom
160          IF (abs(dei)<0.01) dei = 0.01
161          Sij(il, i, j) = anum/dei
162          Sij(il, i, i) = 1.0
163          altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
164          altem = altem/bf2
165          cwat = clw(il, j)*(1.-ep(il,j))
166          stemp = Sij(il, i, j)
167          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
168            anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
169            denom = denom + lv(il, j)*(rr(il,i)-rti)
170            IF (abs(denom)<0.01) denom = 0.01
171            Sij(il, i, j) = anum/denom
172            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
173            altem = altem - (bf2-1.)*cwat
174          END IF
175          IF (Sij(il,i,j)>0.0) THEN
176!!!                 Ment(il,i,j)=m(il,i)
177            Ment(il, i, j) = 1.
178            elij(il, i, j) = altem
179            elij(il, i, j) = amax1(0.0, elij(il,i,j))
180            nent(il, i) = nent(il, i) + 1
181          END IF
182
183          Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
184          Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
185        END IF ! new
186      END DO
187    END DO
188
189
190! ***   if no air can entrain at level i assume that updraft detrains  ***
191! ***   at that level and calculate detrained air flux and properties  ***
192
193
194! @      do 170 i=icb(il),inb(il)
195
196    DO il = 1, ncum
197      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
198! @      if(nent(il,i).eq.0)then
199!!!       Ment(il,i,i)=m(il,i)
200        Ment(il, i, i) = 1.
201        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
202        uent(il, i, i) = unk(il)
203        vent(il, i, i) = vnk(il)
204        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
205        Sij(il, i, i) = 0.0
206      END IF
207    END DO
208  END DO
209
210  DO j = 1, ntra
211    DO i = minorig + 1, nl
212      DO il = 1, ncum
213        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
214          traent(il, i, i, j) = tra(il, nk(il), j)
215        END IF
216      END DO
217    END DO
218  END DO
219
220  DO j = minorig, nl
221    DO i = minorig, nl
222      DO il = 1, ncum
223        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
224            (i>=icb(il)) .AND. (i<=inb(il))) THEN
225          Sigij(il, i, j) = Sij(il, i, j)
226        END IF
227      END DO
228    END DO
229  END DO
230! @      enddo
231
232! @170   continue
233
234! =====================================================================
235! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
236! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
237! =====================================================================
238
239  CALL zilch(csum, nloc*nd)
240
241  DO il = 1, ncum
242    lwork(il) = .FALSE.
243  END DO
244
245! ---------------------------------------------------------------
246  DO i = minorig + 1, nl      !Loop on origin level "i"
247! ---------------------------------------------------------------
248
249    num1 = 0
250    DO il = 1, ncum
251      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
252    END DO
253    IF (num1<=0) GO TO 789
254
255
256!JYG1    Find maximum of SIJ for J>I, if any.
257
258    Sx(:) = 0.
259
260    DO il = 1, ncum
261      IF (i>=icb(il) .AND. i<=inb(il)) THEN
262        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
263        Sbef(il) = max(0., signhpmh(il))
264      END IF
265    END DO
266
267    DO j = i + 1, nl
268      DO il = 1, ncum
269        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
270          IF (Sbef(il)<Sij(il,i,j)) THEN
271            Sx(il) = max(Sij(il,i,j), Sx(il))
272          END IF
273          Sbef(il) = Sij(il, i, j)
274        END IF
275      END DO
276    END DO
277
278
279    DO il = 1, ncum
280      IF (i>=icb(il) .AND. i<=inb(il)) THEN
281        lwork(il) = (nent(il,i)/=0)
282        rti = qnk(il) - ep(il, i)*clw(il, i)
283        anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
284               (cpv-cpd)*t(il, i)*(rti-rr(il,i))
285        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
286                (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
287        IF (abs(denom)<0.01) denom = 0.01
288        Scrit(il) = min(anum/denom, 1.)
289        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
290
291!JYG1    Find new critical value Scrit2
292!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
293!                     Sij < Scrit2  => mixed draught will detrain at J>I
294
295        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
296                 Scrit(il)*max(0., signhpmh(il))
297
298        Scrit(il) = Scrit2
299
300!JYG    Correction pour la nouvelle logique; la correction pour ALT
301! est un peu au hazard
302        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
303        IF (alt<=0.0) Scrit(il) = 1.0
304
305        smax(il) = 0.0
306        ASij(il) = 0.0
307        sup(il) = 0.      ! upper S-value reached by descending draughts
308      END IF
309    END DO
310
311! ---------------------------------------------------------------
312    DO j = minorig, nl         !Loop on destination level "j"
313! ---------------------------------------------------------------
314
315      num2 = 0
316      DO il = 1, ncum
317        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
318            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
319            lwork(il)) num2 = num2 + 1
320      END DO
321      IF (num2<=0) GO TO 175
322
323! -----------------------------------------------
324      IF (j>i) THEN
325! -----------------------------------------------
326        DO il = 1, ncum
327          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
328              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
329              lwork(il)) THEN
330            IF (Sij(il,i,j)>0.0) THEN
331              Smid(il) = min(Sij(il,i,j), Scrit(il))
332              Sjmax(il) = Smid(il)
333              Sjmin(il) = Smid(il)
334              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
335                smin(il) = Smid(il)
336                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
337                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
338                Sjmin(il) = min(Sjmin(il), Scrit(il))
339                Sbef(il) = Sij(il, i, j)
340              END IF
341            END IF
342          END IF
343        END DO
344! -----------------------------------------------
345      ELSE IF (j==i) THEN
346! -----------------------------------------------
347        DO il = 1, ncum
348          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
349              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
350              lwork(il)) THEN
351            IF (Sij(il,i,j)>0.0) THEN
352              Smid(il) = 1.
353              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
354                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
355              Sjmin(il) = max(Sjmin(il), sup(il))
356              Sjmax(il) = 1.
357
358! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
359              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
360
361              smin(il) = 1.
362              Sbef(il) = max(0., signhpmh(il))
363              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
364            END IF
365          END IF
366        END DO
367! -----------------------------------------------
368      ELSE IF (j<i) THEN
369! -----------------------------------------------
370        DO il = 1, ncum
371          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
372              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
373              lwork(il)) THEN
374            IF (Sij(il,i,j)>0.0) THEN
375              Smid(il) = max(Sij(il,i,j), Scrit(il))
376              Sjmax(il) = Smid(il)
377              Sjmin(il) = Smid(il)
378              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
379                smax(il) = Smid(il)
380                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
381                Sjmax(il) = max(Sjmax(il), Scrit(il))
382                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
383                Sjmin(il) = max(Sjmin(il), Scrit(il))
384                Sbef(il) = Sij(il, i, j)
385              END IF
386              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
387                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
388            END IF
389          END IF
390        END DO
391! -----------------------------------------------
392      END IF
393! -----------------------------------------------
394
395
396      DO il = 1, ncum
397        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
398            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
399            lwork(il)) THEN
400          IF (Sij(il,i,j)>0.0) THEN
401            rti = qnk(il) - ep(il, i)*clw(il, i)
402            Qmixmax(il) = Qmix(Sjmax(il))
403            Qmixmin(il) = Qmix(Sjmin(il))
404            Rmixmax(il) = Rmix(Sjmax(il))
405            Rmixmin(il) = Rmix(Sjmin(il))
406            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
407            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
408
409            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
410
411! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
412            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
413              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
414            ELSE
415              Sigij(il, i, j) = 0.
416            END IF
417
418! --    Compute Qent, uent, vent according to the true mixing fraction
419            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
420            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
421            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
422
423! --     Compute liquid water static energy of mixed draughts
424!    IF (j .GT. i) THEN
425!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
426!      awat=amax1(awat,0.0)
427!    ELSE
428!      awat = 0.
429!    ENDIF
430!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
431!    :         + Sigij(il,i,j)*H(il,i)
432!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
433!IM 301008 beg
434            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
435
436            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
437            elij(il, i, j) = elij(il, i, j) + &
438                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
439                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
440            elij(il, i, j) = elij(il, i, j) / &
441                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
442                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
443
444            elij(il, i, j) = max(elij(il,i,j), 0.)
445
446            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
447
448            IF (j>i) THEN
449              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
450              awat = amax1(awat, 0.0)
451            ELSE
452              awat = 0.
453            END IF
454
455! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
456! :         t(il,j))
457
458            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
459!IM 301008 end
460
461! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
462! :               i,j,hent(il,i,j),Sigij(il,i,j)
463
464! --      ASij is the integral of P(F) over the relevant F interval
465            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
466                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
467
468          END IF
469        END IF
470      END DO
471      DO k = 1, ntra
472        DO il = 1, ncum
473          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
474              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
475              lwork(il)) THEN
476            IF (Sij(il,i,j)>0.0) THEN
477              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
478                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
479            END IF
480          END IF
481        END DO
482      END DO
483
484! --    If I=J (detrainement and entrainement at the same level), then only the
485! --    adiabatic ascent part of the mixture is considered
486      IF (i==j) THEN
487        DO il = 1, ncum
488          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
489              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
490              lwork(il)) THEN
491            IF (Sij(il,i,j)>0.0) THEN
492              rti = qnk(il) - ep(il, i)*clw(il, i)
493!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
494              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
495                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
496              Qent(il, i, i) = rti
497              uent(il, i, i) = unk(il)
498              vent(il, i, i) = vnk(il)
499              hent(il, i, i) = hp(il, i)
500              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
501              Sigij(il, i, i) = 0.
502            END IF
503          END IF
504        END DO
505        DO k = 1, ntra
506          DO il = 1, ncum
507            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
508                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
509                lwork(il)) THEN
510              IF (Sij(il,i,j)>0.0) THEN
511                traent(il, i, i, k) = tra(il, nk(il), k)
512              END IF
513            END IF
514          END DO
515        END DO
516
517      END IF
518
519! ---------------------------------------------------------------
520175 END DO        ! End loop on destination level "j"
521! ---------------------------------------------------------------
522
523    DO il = 1, ncum
524      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
525        ASij(il) = amax1(1.0E-16, ASij(il))
526        ASij(il) = 1.0/ASij(il)
527        csum(il, i) = 0.0
528      END IF
529    END DO
530
531    DO j = minorig, nl
532      DO il = 1, ncum
533        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
534            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
535          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
536        END IF
537      END DO
538    END DO
539
540    DO j = minorig, nl
541      DO il = 1, ncum
542        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
543            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
544          csum(il, i) = csum(il, i) + Ment(il, i, j)
545        END IF
546      END DO
547    END DO
548
549    DO il = 1, ncum
550      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
551! cc     :     .and. csum(il,i).lt.m(il,i) ) then
552        nent(il, i) = 0
553! cc        Ment(il,i,i)=m(il,i)
554        Ment(il, i, i) = 1.
555        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
556        uent(il, i, i) = unk(il)
557        vent(il, i, i) = vnk(il)
558        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
559        Sij(il, i, i) = 0.0
560      END IF
561    END DO ! il
562
563    DO j = 1, ntra
564      DO il = 1, ncum
565        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
566! cc     :     .and. csum(il,i).lt.m(il,i) ) then
567          traent(il, i, i, j) = tra(il, nk(il), j)
568        END IF
569      END DO
570    END DO
571
572! ---------------------------------------------------------------
573789 END DO              ! End loop on origin level "i"
574! ---------------------------------------------------------------
575
576
577  RETURN
578END SUBROUTINE cv3p_mixing
579
Note: See TracBrowser for help on using the repository browser.