source: LMDZ6/branches/Amaury_dev/libf/phylmd/cv3p_mixing.F90 @ 5409

Last change on this file since 5409 was 5160, checked in by abarral, 5 months ago

Put .h into modules

  • 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: 25.3 KB
RevLine 
[2007]1SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
[3496]2                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
[2007]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! **************************************************************
[879]14
[5112]15  USE lmdz_print_control, ONLY: mydebug=>debug , lunout, prt_level
16  USE lmdz_ioipsl_getin_p, ONLY: getin_p
[2900]17  USE add_phys_tend_mod, ONLY: fl_cor_ebil
[5140]18  USE lmdz_cvflag
[5141]19  USE lmdz_cvthermo
20  USE lmdz_cv3param
[5158]21  USE lmdz_yomcst2
[2393]22
[1992]23  IMPLICIT NONE
[879]24
[2007]25!inputs:
26  INTEGER, INTENT (IN)                               :: ncum, nd, na
27  INTEGER, INTENT (IN)                               :: ntra, nloc
28  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
29  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
[3496]30  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
31  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
[2007]32  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
33  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
34  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
35  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
36  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
[2397]37  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
38  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
[2393]39  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
[2007]40  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
41  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
42  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
[879]43
[2007]44!outputs:
45  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
46  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
47  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
[2393]48  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax           ! Highest mixing fraction of mixed
[2007]49                                                                         ! updraughts with the sign of (h-hp)
50  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
51  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
52  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
53  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
[879]54
[2007]55!local variables:
[1992]56  INTEGER i, j, k, il, im, jm
57  INTEGER num1, num2
[2397]58  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
[2007]59  REAL                               :: alt, delp, delm
60  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
61  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
62  REAL, DIMENSION (nloc)             :: signhpmh
63  REAL, DIMENSION (nloc)             :: Sx
64  REAL                               :: Scrit2
65  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
66  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
[2226]67  REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
[2007]68  REAL, DIMENSION (nloc, nd, nd)     :: Sij
69  REAL, DIMENSION (nloc, nd)         :: csum
70  REAL                               :: awat
[2397]71  REAL                               :: cpm        !Mixed draught heat capacity
72  REAL                               :: Tm         !Mixed draught temperature
[2007]73  LOGICAL, DIMENSION (nloc)          :: lwork
[879]74
[1992]75  REAL amxupcrit, df, ff
76  INTEGER nstep
[879]77
[2393]78  INTEGER,SAVE                                       :: igout=1
79!$OMP THREADPRIVATE(igout)
80
[2007]81! --   Mixing probability distribution functions
[1650]82
[2007]83  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
[879]84
[2007]85  Qcoef1(F) = tanh(F/gammas)
86  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
87  QFF(F) = max(min(F,1.), 0.)
88  QFFf(F) = min(QFF(F), scut)
89  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
90  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
91  Qmix2(F) = -log(1.-QFFf(F))/scut
92  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
93  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
94  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
[879]95
[1992]96  INTEGER, SAVE :: ifrst
97  DATA ifrst/0/
[2007]98!$OMP THREADPRIVATE(ifrst)
[879]99
100
[2007]101! =====================================================================
102! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
103! =====================================================================
[879]104
[2007]105! -- Initialize mixing PDF coefficients
[1992]106  IF (ifrst==0) THEN
107    ifrst = 1
[2007]108    Qcoef1max = Qcoef1(Fmax)
109    Qcoef2max = Qcoef2(Fmax)
[2393]110!<jyg
[5103]111   PRINT*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
[2393]112            fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
113!>jyg
[5099]114
[1992]115  END IF
[879]116
117
[2007]118! ori        do 360 i=1,ncum*nlp
[1992]119  DO j = 1, nl
120    DO i = 1, ncum
121      nent(i, j) = 0
[2007]122! in convect3, m is computed in cv3_closure
123! ori          m(i,1)=0.0
[1992]124    END DO
125  END DO
[879]126
[2007]127! ori      do 400 k=1,nlp
128! ori       do 390 j=1,nlp
[1992]129  DO j = 1, nl
130    DO k = 1, nl
131      DO i = 1, ncum
[2007]132        Qent(i, k, j) = rr(i, j)
[1992]133        uent(i, k, j) = u(i, j)
[2478]134        vent(i, k, j) = v(i, j)
[1992]135        elij(i, k, j) = 0.0
136        hent(i, k, j) = 0.0
[2007]137!AC!            Ment(i,k,j)=0.0
138!AC!            Sij(i,k,j)=0.0
[1992]139      END DO
140    END DO
141  END DO
[879]142
[2007]143!AC!
144  Ment(1:ncum, 1:nd, 1:nd) = 0.0
145  Sij(1:ncum, 1:nd, 1:nd) = 0.0
146!AC!
[2478]147!ym
148  Sigij(1:ncum, 1:nd, 1:nd) = 0.0
149!ym
[879]150
[2393]151!jyg!  DO k = 1, ntra
152!jyg!    DO j = 1, nd ! instead nlp
153!jyg!      DO i = 1, nd ! instead nlp
154!jyg!        DO il = 1, ncum
155!jyg!          traent(il, i, j, k) = tra(il, j, k)
156!jyg!        END DO
157!jyg!      END DO
158!jyg!    END DO
159!jyg!  END DO
[879]160
[2007]161! =====================================================================
162! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
163! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
164! --- FRACTION (Sij)
165! =====================================================================
[879]166
[1992]167  DO i = minorig + 1, nl
[879]168
[2900]169    IF (ok_entrain) THEN
170      DO j = minorig, nl
171        DO il = 1, ncum
172          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
173                           .AND. (j<=inb(il))) THEN
[879]174
[3496]175!!            rti = qnk(il) - ep(il, i)*clw(il, i)
176            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
[2900]177            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
[2397]178!jyg(from aj)<
[2900]179            IF (cvflag_ice) THEN
[5103]180! PRINT*,cvflag_ice,'cvflag_ice dans do 700'
[2900]181              IF (t(il,j)<=263.15) THEN
182                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
183                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
184              END IF
[2397]185            END IF
186!>jyg
[2900]187            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
188            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
189            dei = denom
190            IF (abs(dei)<0.01) dei = 0.01
191            Sij(il, i, j) = anum/dei
192            Sij(il, i, i) = 1.0
193            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
194            altem = altem/bf2
195            cwat = clw(il, j)*(1.-ep(il,j))
196            stemp = Sij(il, i, j)
197            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
[2397]198!jyg(from aj)<
[2900]199              IF (cvflag_ice) THEN
200                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
201                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
202              ELSE
203                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
204                denom = denom + lv(il, j)*(rr(il,i)-rti)
205              END IF
206!>jyg
207              IF (abs(denom)<0.01) denom = 0.01
208              Sij(il, i, j) = anum/denom
209              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
210              altem = altem - (bf2-1.)*cwat
[2397]211            END IF
[2900]212            IF (Sij(il,i,j)>0.0) THEN
[2007]213!!!                 Ment(il,i,j)=m(il,i)
[2900]214              Ment(il, i, j) = 1.
215              elij(il, i, j) = altem
216              elij(il, i, j) = amax1(0.0, elij(il,i,j))
217              nent(il, i) = nent(il, i) + 1
218            END IF
[879]219
[2900]220            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
221            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
[3496]222          ELSE IF (j > i) THEN
223            IF (prt_level >= 10) THEN
[5160]224              PRINT *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
[3496]225            ENDIF
[2900]226          END IF ! new
227        END DO
[1992]228      END DO
[2900]229    ELSE  ! (ok_entrain)
230      DO il = 1,ncum
231        nent(il,i) = 0
232      ENDDO
233    ENDIF ! (ok_entrain)
[879]234
[2393]235!jygdebug<
236    IF (prt_level >= 10) THEN
[5160]237      PRINT *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
[5082]238      IF (nent(igout,i) > 0) THEN
[5160]239        PRINT *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
[2393]240      ENDIF
241    ENDIF
242!>jygdebug
[879]243
[2007]244! ***   if no air can entrain at level i assume that updraft detrains  ***
245! ***   at that level and calculate detrained air flux and properties  ***
[1494]246
[879]247
[2007]248! @      do 170 i=icb(il),inb(il)
[879]249
[1992]250    DO il = 1, ncum
251      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
[5117]252! @      IF(nent(il,i).EQ.0)THEN
[2007]253!!!       Ment(il,i,i)=m(il,i)
254        Ment(il, i, i) = 1.
[3496]255!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
256        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
[1992]257        uent(il, i, i) = unk(il)
258        vent(il, i, i) = vnk(il)
[5082]259        IF (fl_cor_ebil >= 2) THEN
[2902]260          hent(il, i, i) = hp(il,i)
261        ENDIF
[1992]262        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
[2007]263        Sij(il, i, i) = 0.0
[1992]264      END IF
265    END DO
[2902]266  END DO ! i = minorig + 1, nl
[879]267
[2393]268!jyg!  DO j = 1, ntra
269!jyg!    DO i = minorig + 1, nl
270!jyg!      DO il = 1, ncum
271!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
272!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
273!jyg!        END IF
274!jyg!      END DO
275!jyg!    END DO
276!jyg!  END DO
[879]277
[1992]278  DO j = minorig, nl
279    DO i = minorig, nl
280      DO il = 1, ncum
[2007]281        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
282            (i>=icb(il)) .AND. (i<=inb(il))) THEN
283          Sigij(il, i, j) = Sij(il, i, j)
[1992]284        END IF
285      END DO
286    END DO
287  END DO
[2007]288! @      enddo
[1041]289
[2007]290! @170   continue
[1041]291
[2007]292! =====================================================================
293! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
294! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
295! =====================================================================
[1041]296
[1992]297  CALL zilch(csum, nloc*nd)
[1041]298
[1992]299  DO il = 1, ncum
300    lwork(il) = .FALSE.
301  END DO
[1041]302
[2007]303! ---------------------------------------------------------------
304  DO i = minorig + 1, nl      !Loop on origin level "i"
305! ---------------------------------------------------------------
[1041]306
[1992]307    num1 = 0
308    DO il = 1, ncum
309      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
310    END DO
311    IF (num1<=0) GO TO 789
[879]312
313
[2007]314!JYG1    Find maximum of SIJ for J>I, if any.
[879]315
[2007]316    Sx(:) = 0.
[879]317
[1992]318    DO il = 1, ncum
319      IF (i>=icb(il) .AND. i<=inb(il)) THEN
320        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
[2007]321        Sbef(il) = max(0., signhpmh(il))
[1992]322      END IF
323    END DO
[879]324
[1992]325    DO j = i + 1, nl
326      DO il = 1, ncum
327        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
[2007]328          IF (Sbef(il)<Sij(il,i,j)) THEN
329            Sx(il) = max(Sij(il,i,j), Sx(il))
[1992]330          END IF
[2007]331          Sbef(il) = Sij(il, i, j)
[1992]332        END IF
333      END DO
334    END DO
[879]335
[1992]336
337    DO il = 1, ncum
338      IF (i>=icb(il) .AND. i<=inb(il)) THEN
339        lwork(il) = (nent(il,i)/=0)
[3496]340!!        rti = qnk(il) - ep(il, i)*clw(il, i)
341        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
[2397]342!jyg<
343        IF (cvflag_ice) THEN
344
345          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
346                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
347          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
348                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
349        ELSE
350
351          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
352                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
353          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
354                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
355        END IF
356!>jyg
[1992]357        IF (abs(denom)<0.01) denom = 0.01
[2007]358        Scrit(il) = min(anum/denom, 1.)
359        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
[1992]360
[2007]361!JYG1    Find new critical value Scrit2
362!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
363!                     Sij < Scrit2  => mixed draught will detrain at J>I
[1992]364
[2007]365        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
366                 Scrit(il)*max(0., signhpmh(il))
[1992]367
[2007]368        Scrit(il) = Scrit2
[1992]369
[2007]370!JYG    Correction pour la nouvelle logique; la correction pour ALT
371! est un peu au hazard
372        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
373        IF (alt<=0.0) Scrit(il) = 1.0
[1992]374
375        smax(il) = 0.0
[2007]376        ASij(il) = 0.0
377        sup(il) = 0.      ! upper S-value reached by descending draughts
[1992]378      END IF
379    END DO
380
[2007]381! ---------------------------------------------------------------
382    DO j = minorig, nl         !Loop on destination level "j"
383! ---------------------------------------------------------------
[1992]384
385      num2 = 0
386      DO il = 1, ncum
[2007]387        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
388            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
389            lwork(il)) num2 = num2 + 1
[1992]390      END DO
391      IF (num2<=0) GO TO 175
392
[2007]393! -----------------------------------------------
[1992]394      IF (j>i) THEN
[2007]395! -----------------------------------------------
[1992]396        DO il = 1, ncum
[2007]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              Smid(il) = min(Sij(il,i,j), Scrit(il))
402              Sjmax(il) = Smid(il)
403              Sjmin(il) = Smid(il)
404              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
405                smin(il) = Smid(il)
406                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
407                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
408                Sjmin(il) = min(Sjmin(il), Scrit(il))
409                Sbef(il) = Sij(il, i, j)
[1992]410              END IF
411            END IF
412          END IF
413        END DO
[2007]414! -----------------------------------------------
[1992]415      ELSE IF (j==i) THEN
[2007]416! -----------------------------------------------
[1992]417        DO il = 1, ncum
[2007]418          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
419              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
420              lwork(il)) THEN
421            IF (Sij(il,i,j)>0.0) THEN
422              Smid(il) = 1.
423              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
424                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
425              Sjmin(il) = max(Sjmin(il), sup(il))
426              Sjmax(il) = 1.
[1992]427
[2007]428! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
429              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
[1992]430
431              smin(il) = 1.
[2007]432              Sbef(il) = max(0., signhpmh(il))
433              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
[1992]434            END IF
435          END IF
436        END DO
[2007]437! -----------------------------------------------
[1992]438      ELSE IF (j<i) THEN
[2007]439! -----------------------------------------------
[1992]440        DO il = 1, ncum
[2007]441          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
442              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
443              lwork(il)) THEN
444            IF (Sij(il,i,j)>0.0) THEN
445              Smid(il) = max(Sij(il,i,j), Scrit(il))
446              Sjmax(il) = Smid(il)
447              Sjmin(il) = Smid(il)
448              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
449                smax(il) = Smid(il)
450                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
451                Sjmax(il) = max(Sjmax(il), Scrit(il))
452                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
453                Sjmin(il) = max(Sjmin(il), Scrit(il))
454                Sbef(il) = Sij(il, i, j)
[1992]455              END IF
[2007]456              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
457                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
[1992]458            END IF
459          END IF
460        END DO
[2007]461! -----------------------------------------------
[1992]462      END IF
[2007]463! -----------------------------------------------
[1992]464
465
466      DO il = 1, ncum
[2007]467        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
468            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
469            lwork(il)) THEN
470          IF (Sij(il,i,j)>0.0) THEN
[3496]471!!            rti = qnk(il) - ep(il, i)*clw(il, i)
472            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
[2007]473            Qmixmax(il) = Qmix(Sjmax(il))
474            Qmixmin(il) = Qmix(Sjmin(il))
475            Rmixmax(il) = Rmix(Sjmax(il))
476            Rmixmin(il) = Rmix(Sjmin(il))
477            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
478            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
[1992]479
[2007]480            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
[1992]481
[2007]482! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
483            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
484              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
[1992]485            ELSE
[2007]486              Sigij(il, i, j) = 0.
[1992]487            END IF
488
[2007]489! --    Compute Qent, uent, vent according to the true mixing fraction
490            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
491            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
492            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
[1992]493
[2007]494! --     Compute liquid water static energy of mixed draughts
495!    IF (j .GT. i) THEN
496!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
497!      awat=amax1(awat,0.0)
498!    ELSE
499!      awat = 0.
500!    ENDIF
501!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
502!    :         + Sigij(il,i,j)*H(il,i)
503!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
504!IM 301008 beg
505            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
[1992]506
[2397]507!jyg<
508!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
509!            elij(il, i, j) = elij(il, i, j) + &
510!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
511!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
512!            elij(il, i, j) = elij(il, i, j) / &
513!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
514!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
[5099]515
[2397]516!       Computation of condensate amount Elij, taking into account the ice fraction frac
517!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
518!                 should be corrected.
[5099]519
[2397]520!  Heat capacity of mixed draught
521    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
[5099]522
[5117]523    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
[2007]524            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
525            elij(il, i, j) = elij(il, i, j) + &
[2397]526                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
527                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
[2007]528            elij(il, i, j) = elij(il, i, j) / &
[2397]529                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
530                              (cpm*rrv*t(il,j)*t(il,j)))
531    ELSE
532            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
533            elij(il, i, j) = elij(il, i, j) + &
534                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
535                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
536            elij(il, i, j) = elij(il, i, j) / &
[2007]537                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
[2397]538                              (cpm*rrv*t(il,j)*t(il,j)))
539    ENDIF
540!>jyg
[1992]541            elij(il, i, j) = max(elij(il,i,j), 0.)
542
[2007]543            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
[1992]544
545            IF (j>i) THEN
546              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
547              awat = amax1(awat, 0.0)
548            ELSE
549              awat = 0.
550            END IF
551
[5160]552! PRINT *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
[2007]553! :         t(il,j))
[1992]554
[2397]555!jyg<
556!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
557! Mixed draught temperature at level j
[5117]558    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
[2397]559          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
560          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
561    ELSE
562          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
563          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
564    ENDIF
565!>jyg
566
[2007]567!IM 301008 end
[1992]568
[5160]569! PRINT *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
[2007]570! :               i,j,hent(il,i,j),Sigij(il,i,j)
[1992]571
[2007]572! --      ASij is the integral of P(F) over the relevant F interval
573            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
574                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
[1992]575
576          END IF
577        END IF
578      END DO
[2393]579!jyg!      DO k = 1, ntra
580!jyg!        DO il = 1, ncum
581!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
582!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
583!jyg!              lwork(il)) THEN
584!jyg!            IF (Sij(il,i,j)>0.0) THEN
585!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
586!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
587!jyg!            END IF
588!jyg!          END IF
589!jyg!        END DO
590!jyg!      END DO
[1992]591
[2007]592! --    If I=J (detrainement and entrainement at the same level), then only the
593! --    adiabatic ascent part of the mixture is considered
[1992]594      IF (i==j) THEN
595        DO il = 1, ncum
[2007]596          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
597              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
598              lwork(il)) THEN
599            IF (Sij(il,i,j)>0.0) THEN
[3496]600!!              rti = qnk(il) - ep(il, i)*clw(il, i)
601              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
[2007]602!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
603              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
604                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
605              Qent(il, i, i) = rti
[1992]606              uent(il, i, i) = unk(il)
607              vent(il, i, i) = vnk(il)
608              hent(il, i, i) = hp(il, i)
609              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
[2007]610              Sigij(il, i, i) = 0.
[1992]611            END IF
612          END IF
613        END DO
[2393]614!jyg!        DO k = 1, ntra
615!jyg!          DO il = 1, ncum
616!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
617!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
618!jyg!                lwork(il)) THEN
619!jyg!              IF (Sij(il,i,j)>0.0) THEN
620!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
621!jyg!              END IF
622!jyg!            END IF
623!jyg!          END DO
624!jyg!        END DO
[1992]625
626      END IF
627
[2007]628! ---------------------------------------------------------------
629175 END DO        ! End loop on destination level "j"
630! ---------------------------------------------------------------
[1992]631
632    DO il = 1, ncum
633      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
[2007]634        ASij(il) = amax1(1.0E-16, ASij(il))
[2226]635!jyg+lluis<
636!!        ASij(il) = 1.0/ASij(il)
637        ASij_inv(il) = 1.0/ASij(il)
638!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
639        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
640!>jyg+lluis
[1992]641        csum(il, i) = 0.0
642      END IF
643    END DO
644
645    DO j = minorig, nl
646      DO il = 1, ncum
[2007]647        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
648            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
[2226]649!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
650          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
[1992]651        END IF
652      END DO
653    END DO
654
655    DO j = minorig, nl
656      DO il = 1, ncum
[2007]657        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
658            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
659          csum(il, i) = csum(il, i) + Ment(il, i, j)
[1992]660        END IF
661      END DO
662    END DO
663
664    DO il = 1, ncum
[2007]665      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
[5117]666! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
[1992]667        nent(il, i) = 0
[2007]668! cc        Ment(il,i,i)=m(il,i)
669        Ment(il, i, i) = 1.
[3496]670!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
671        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
[1992]672        uent(il, i, i) = unk(il)
673        vent(il, i, i) = vnk(il)
674        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
[5082]675        IF (fl_cor_ebil >= 2) THEN
[2902]676          hent(il, i, i) = hp(il,i)
677          Sigij(il, i, i) = 0.0
678        ELSE
679          Sij(il, i, i) = 0.0
680        ENDIF
[1992]681      END IF
682    END DO ! il
683
[2393]684!jyg!    DO j = 1, ntra
685!jyg!      DO il = 1, ncum
686!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
[5117]687!jyg!! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
[2393]688!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
689!jyg!        END IF
690!jyg!      END DO
691!jyg!    END DO
[1992]692
[2007]693! ---------------------------------------------------------------
694789 END DO              ! End loop on origin level "i"
695! ---------------------------------------------------------------
[1992]696
[2007]697
[5105]698
[1992]699END SUBROUTINE cv3p_mixing
700
Note: See TracBrowser for help on using the repository browser.