source: LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90 @ 1999

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

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