source: LMDZ5/trunk/libf/phylmd/cv3p_mixing.F @ 1622

Last change on this file since 1622 was 1573, checked in by jghattas, 13 years ago
  • Added missing attribute SAVE.
  • Changed FORALL into DO loop because some OpenMP compilers do not permet FORALL within the body of an OpenMP SINGLE directive. (Titane)

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