GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/cv3p_mixing.F90 Lines: 238 258 92.2 %
Date: 2023-06-30 12:56:34 Branches: 274 316 86.7 %

Line Branch Exec Source
1
40068127
SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2
144
                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
3
144
                       unk, vnk, hp, tv, tvp, ep, clw, sig, &
4
144
                       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
  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
16
  USE ioipsl_getin_p_mod, ONLY: getin_p
17
  USE add_phys_tend_mod, ONLY: fl_cor_ebil
18
19
  IMPLICIT NONE
20
21
  include "cvthermo.h"
22
  include "cv3param.h"
23
  include "YOMCST2.h"
24
  include "cvflag.h"
25
26
!inputs:
27
  INTEGER, INTENT (IN)                               :: ncum, nd, na
28
  INTEGER, INTENT (IN)                               :: ntra, nloc
29
  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
30
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
31
  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
32
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
33
  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
34
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
35
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
36
  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
37
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
38
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
39
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
40
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
41
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
42
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
43
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
44
45
!outputs:
46
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
47
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
48
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
49
  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax           ! Highest mixing fraction of mixed
50
                                                                         ! updraughts with the sign of (h-hp)
51
  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
52
  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
53
  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
54
  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
55
56
!local variables:
57
  INTEGER i, j, k, il, im, jm
58
  INTEGER num1, num2
59
  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
60
  REAL                               :: alt, delp, delm
61
288
  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
62
288
  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
63
288
  REAL, DIMENSION (nloc)             :: signhpmh
64
288
  REAL, DIMENSION (nloc)             :: Sx
65
  REAL                               :: Scrit2
66
288
  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
67
288
  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
68
288
  REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
69
288
  REAL, DIMENSION (nloc, nd, nd)     :: Sij
70
288
  REAL, DIMENSION (nloc, nd)         :: csum
71
  REAL                               :: awat
72
  REAL                               :: cpm        !Mixed draught heat capacity
73
  REAL                               :: Tm         !Mixed draught temperature
74
288
  LOGICAL, DIMENSION (nloc)          :: lwork
75
76
  REAL amxupcrit, df, ff
77
  INTEGER nstep
78
79
  INTEGER,SAVE                                       :: igout=1
80
!$OMP THREADPRIVATE(igout)
81
82
! --   Mixing probability distribution functions
83
84
  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
85
86
  Qcoef1(F) = tanh(F/gammas)
87
  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
88
  QFF(F) = max(min(F,1.), 0.)
89
  QFFf(F) = min(QFF(F), scut)
90
  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
91
  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
92
  Qmix2(F) = -log(1.-QFFf(F))/scut
93
  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
94
  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
95
  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
96
97
  INTEGER, SAVE :: ifrst
98
  DATA ifrst/0/
99
!$OMP THREADPRIVATE(ifrst)
100
101
102
! =====================================================================
103
! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
104
! =====================================================================
105
106
! -- Initialize mixing PDF coefficients
107
144
  IF (ifrst==0) THEN
108
1
    ifrst = 1
109
1
    Qcoef1max = Qcoef1(Fmax)
110
1
    Qcoef2max = Qcoef2(Fmax)
111
!<jyg
112
1
   print*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
113
2
            fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
114
!>jyg
115
!
116
  END IF
117
118
119
! ori        do 360 i=1,ncum*nlp
120
4032
  DO j = 1, nl
121
1866114
    DO i = 1, ncum
122
1865970
      nent(i, j) = 0
123
! in convect3, m is computed in cv3_closure
124
! ori          m(i,1)=0.0
125
    END DO
126
  END DO
127
128
! ori      do 400 k=1,nlp
129
! ori       do 390 j=1,nlp
130
4032
  DO j = 1, nl
131
109008
    DO k = 1, nl
132
50385078
      DO i = 1, ncum
133
50276214
        Qent(i, k, j) = rr(i, j)
134
50276214
        uent(i, k, j) = u(i, j)
135
50276214
        vent(i, k, j) = v(i, j)
136
50276214
        elij(i, k, j) = 0.0
137
50381190
        hent(i, k, j) = 0.0
138
!AC!            Ment(i,k,j)=0.0
139
!AC!            Sij(i,k,j)=0.0
140
      END DO
141
    END DO
142
  END DO
143
144
!AC!
145

105122070
  Ment(1:ncum, 1:nd, 1:nd) = 0.0
146

105122070
  Sij(1:ncum, 1:nd, 1:nd) = 0.0
147
!AC!
148
!ym
149

105122070
  Sigij(1:ncum, 1:nd, 1:nd) = 0.0
150
!ym
151
152
!jyg!  DO k = 1, ntra
153
!jyg!    DO j = 1, nd ! instead nlp
154
!jyg!      DO i = 1, nd ! instead nlp
155
!jyg!        DO il = 1, ncum
156
!jyg!          traent(il, i, j, k) = tra(il, j, k)
157
!jyg!        END DO
158
!jyg!      END DO
159
!jyg!    END DO
160
!jyg!  END DO
161
162
! =====================================================================
163
! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
164
! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
165
! --- FRACTION (Sij)
166
! =====================================================================
167
168
3888
  DO i = minorig + 1, nl
169
170
3744
    IF (ok_entrain) THEN
171
104832
      DO j = minorig, nl
172
48518964
        DO il = 1, ncum
173
          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
174


48515220
                           .AND. (j<=inb(il))) THEN
175
176
!!            rti = qnk(il) - ep(il, i)*clw(il, i)
177
7095306
            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
178
7095306
            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
179
!jyg(from aj)<
180
7095306
            IF (cvflag_ice) THEN
181
! print*,cvflag_ice,'cvflag_ice dans do 700'
182
7095306
              IF (t(il,j)<=263.15) THEN
183
                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
184
2736968
                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
185
              END IF
186
            END IF
187
!>jyg
188
7095306
            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
189
7095306
            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
190
            dei = denom
191
7095306
            IF (abs(dei)<0.01) dei = 0.01
192
7095306
            Sij(il, i, j) = anum/dei
193
7095306
            Sij(il, i, i) = 1.0
194
7095306
            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
195
7095306
            altem = altem/bf2
196
7095306
            cwat = clw(il, j)*(1.-ep(il,j))
197
            stemp = Sij(il, i, j)
198


7095306
            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
199
!jyg(from aj)<
200
2959186
              IF (cvflag_ice) THEN
201
2959186
                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
202
2959186
                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
203
              ELSE
204
                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
205
                denom = denom + lv(il, j)*(rr(il,i)-rti)
206
              END IF
207
!>jyg
208
2959186
              IF (abs(denom)<0.01) denom = 0.01
209
2959186
              Sij(il, i, j) = anum/denom
210
2959186
              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
211
2959186
              altem = altem - (bf2-1.)*cwat
212
            END IF
213
7095306
            IF (Sij(il,i,j)>0.0) THEN
214
!!!                 Ment(il,i,j)=m(il,i)
215
6174981
              Ment(il, i, j) = 1.
216
              elij(il, i, j) = altem
217
6174981
              elij(il, i, j) = amax1(0.0, elij(il,i,j))
218
6174981
              nent(il, i) = nent(il, i) + 1
219
            END IF
220
221
7095306
            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
222
7095306
            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
223
41318826
          ELSE IF (j > i) THEN
224
19453004
            IF (prt_level >= 10) THEN
225
              print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
226
            ENDIF
227
          END IF ! new
228
        END DO
229
      END DO
230
    ELSE  ! (ok_entrain)
231
      DO il = 1,ncum
232
        nent(il,i) = 0
233
      ENDDO
234
    ENDIF ! (ok_entrain)
235
236
!jygdebug<
237
3744
    IF (prt_level >= 10) THEN
238
      print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
239
      IF (nent(igout,i) .gt. 0) THEN
240
        print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
241
      ENDIF
242
    ENDIF
243
!>jygdebug
244
245
! ***   if no air can entrain at level i assume that updraft detrains  ***
246
! ***   at that level and calculate detrained air flux and properties  ***
247
248
249
! @      do 170 i=icb(il),inb(il)
250
251
1797004
    DO il = 1, ncum
252

1796860
      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
253
! @      if(nent(il,i).eq.0)then
254
!!!       Ment(il,i,i)=m(il,i)
255
        Ment(il, i, i) = 1.
256
!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
257
        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
258
        uent(il, i, i) = unk(il)
259
        vent(il, i, i) = vnk(il)
260
        IF (fl_cor_ebil .GE. 2) THEN
261
          hent(il, i, i) = hp(il,i)
262
        ENDIF
263
        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
264
        Sij(il, i, i) = 0.0
265
      END IF
266
    END DO
267
  END DO ! i = minorig + 1, nl
268
269
!jyg!  DO j = 1, ntra
270
!jyg!    DO i = minorig + 1, nl
271
!jyg!      DO il = 1, ncum
272
!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
273
!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
274
!jyg!        END IF
275
!jyg!      END DO
276
!jyg!    END DO
277
!jyg!  END DO
278
279
4032
  DO j = minorig, nl
280
109008
    DO i = minorig, nl
281
50385078
      DO il = 1, ncum
282
        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
283


50381190
            (i>=icb(il)) .AND. (i<=inb(il))) THEN
284
7095306
          Sigij(il, i, j) = Sij(il, i, j)
285
        END IF
286
      END DO
287
    END DO
288
  END DO
289
! @      enddo
290
291
! @170   continue
292
293
! =====================================================================
294
! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
295
! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
296
! =====================================================================
297
298
144
  CALL zilch(csum, nloc*nd)
299
300
69110
  DO il = 1, ncum
301
69110
    lwork(il) = .FALSE.
302
  END DO
303
304
! ---------------------------------------------------------------
305
3888
  DO i = minorig + 1, nl      !Loop on origin level "i"
306
! ---------------------------------------------------------------
307
308
    num1 = 0
309
1796860
    DO il = 1, ncum
310

1796860
      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
311
    END DO
312
3744
    IF (num1<=0) GO TO 789
313
314
315
!JYG1    Find maximum of SIJ for J>I, if any.
316
317
2725305
    Sx(:) = 0.
318
319
1314694
    DO il = 1, ncum
320

1314694
      IF (i>=icb(il) .AND. i<=inb(il)) THEN
321
586707
        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
322
586707
        Sbef(il) = max(0., signhpmh(il))
323
      END IF
324
    END DO
325
326
43794
    DO j = i + 1, nl
327
19707109
      DO il = 1, ncum
328

19704370
        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
329
2960946
          IF (Sbef(il)<Sij(il,i,j)) THEN
330
539847
            Sx(il) = max(Sij(il,i,j), Sx(il))
331
          END IF
332
2960946
          Sbef(il) = Sij(il, i, j)
333
        END IF
334
      END DO
335
    END DO
336
337
338
1314694
    DO il = 1, ncum
339

1314694
      IF (i>=icb(il) .AND. i<=inb(il)) THEN
340
586707
        lwork(il) = (nent(il,i)/=0)
341
!!        rti = qnk(il) - ep(il, i)*clw(il, i)
342
586707
        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
343
!jyg<
344
586707
        IF (cvflag_ice) THEN
345
346
          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
347
586707
                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
348
          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
349
586707
                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
350
        ELSE
351
352
          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
353
                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
354
          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
355
                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
356
        END IF
357
!>jyg
358
586707
        IF (abs(denom)<0.01) denom = 0.01
359
586707
        Scrit(il) = min(anum/denom, 1.)
360
586707
        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
361
362
!JYG1    Find new critical value Scrit2
363
!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
364
!                     Sij < Scrit2  => mixed draught will detrain at J>I
365
366
        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
367
586707
                 Scrit(il)*max(0., signhpmh(il))
368
369
586707
        Scrit(il) = Scrit2
370
371
!JYG    Correction pour la nouvelle logique; la correction pour ALT
372
! est un peu au hazard
373
586707
        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
374
586707
        IF (alt<=0.0) Scrit(il) = 1.0
375
376
586707
        smax(il) = 0.0
377
586707
        ASij(il) = 0.0
378
586707
        sup(il) = 0.      ! upper S-value reached by descending draughts
379
      END IF
380
    END DO
381
382
! ---------------------------------------------------------------
383
76692
    DO j = minorig, nl         !Loop on destination level "j"
384
! ---------------------------------------------------------------
385
386
      num2 = 0
387
35496738
      DO il = 1, ncum
388
        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
389


35422785
            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
390
7169259
            lwork(il)) num2 = num2 + 1
391
      END DO
392
73953
      IF (num2<=0) GO TO 175
393
394
! -----------------------------------------------
395
54206
      IF (j>i) THEN
396
! -----------------------------------------------
397
11696409
        DO il = 1, ncum
398
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
399


11672045
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
400
24364
              lwork(il)) THEN
401
2960946
            IF (Sij(il,i,j)>0.0) THEN
402
2520477
              Smid(il) = min(Sij(il,i,j), Scrit(il))
403
2520477
              Sjmax(il) = Smid(il)
404
2520477
              Sjmin(il) = Smid(il)
405

2520477
              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
406
1139716
                smin(il) = Smid(il)
407
1139716
                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
408
1139716
                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
409
1139716
                Sjmin(il) = min(Sjmin(il), Scrit(il))
410
1139716
                Sbef(il) = Sij(il, i, j)
411
              END IF
412
            END IF
413
          END IF
414
        END DO
415
! -----------------------------------------------
416
29842
      ELSE IF (j==i) THEN
417
! -----------------------------------------------
418
1314694
        DO il = 1, ncum
419
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
420


1311955
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
421
2739
              lwork(il)) THEN
422
586707
            IF (Sij(il,i,j)>0.0) THEN
423
586707
              Smid(il) = 1.
424
              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
425
586707
                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
426
586707
              Sjmin(il) = max(Sjmin(il), sup(il))
427
586707
              Sjmax(il) = 1.
428
429
! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
430
586707
              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
431
432
586707
              smin(il) = 1.
433
586707
              Sbef(il) = max(0., signhpmh(il))
434
586707
              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
435
            END IF
436
          END IF
437
        END DO
438
! -----------------------------------------------
439
27103
      ELSE IF (j<i) THEN
440
! -----------------------------------------------
441
13011103
        DO il = 1, ncum
442
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
443


12984000
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
444
27103
              lwork(il)) THEN
445
3547653
            IF (Sij(il,i,j)>0.0) THEN
446
3067797
              Smid(il) = max(Sij(il,i,j), Scrit(il))
447
3067797
              Sjmax(il) = Smid(il)
448
3067797
              Sjmin(il) = Smid(il)
449

3067797
              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
450
36752
                smax(il) = Smid(il)
451
36752
                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
452
36752
                Sjmax(il) = max(Sjmax(il), Scrit(il))
453
36752
                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
454
36752
                Sjmin(il) = max(Sjmin(il), Scrit(il))
455
36752
                Sbef(il) = Sij(il, i, j)
456
              END IF
457
3067797
              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
458
34108
                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
459
            END IF
460
          END IF
461
        END DO
462
! -----------------------------------------------
463
      END IF
464
! -----------------------------------------------
465
466
467
26022206
      DO il = 1, ncum
468
        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
469


25968000
            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
470
54206
            lwork(il)) THEN
471
7095306
          IF (Sij(il,i,j)>0.0) THEN
472
!!            rti = qnk(il) - ep(il, i)*clw(il, i)
473
6174981
            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
474
6174981
            Qmixmax(il) = Qmix(Sjmax(il))
475
6174981
            Qmixmin(il) = Qmix(Sjmin(il))
476
6174981
            Rmixmax(il) = Rmix(Sjmax(il))
477
6174981
            Rmixmin(il) = Rmix(Sjmin(il))
478
6174981
            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
479
6174981
            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
480
481
6174981
            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
482
483
! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
484
6174981
            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
485
1653182
              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
486
            ELSE
487
4521799
              Sigij(il, i, j) = 0.
488
            END IF
489
490
! --    Compute Qent, uent, vent according to the true mixing fraction
491
6174981
            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
492
6174981
            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
493
6174981
            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
494
495
! --     Compute liquid water static energy of mixed draughts
496
!    IF (j .GT. i) THEN
497
!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
498
!      awat=amax1(awat,0.0)
499
!    ELSE
500
!      awat = 0.
501
!    ENDIF
502
!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
503
!    :         + Sigij(il,i,j)*H(il,i)
504
!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
505
!IM 301008 beg
506
6174981
            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
507
508
!jyg<
509
!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
510
!            elij(il, i, j) = elij(il, i, j) + &
511
!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
512
!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
513
!            elij(il, i, j) = elij(il, i, j) / &
514
!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
515
!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
516
!
517
!       Computation of condensate amount Elij, taking into account the ice fraction frac
518
!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
519
!                 should be corrected.
520
!
521
!  Heat capacity of mixed draught
522
6174981
    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
523
!
524

6174981
    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
525
2422208
            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
526
            elij(il, i, j) = elij(il, i, j) + &
527
                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
528
2422208
                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
529
            elij(il, i, j) = elij(il, i, j) / &
530
                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
531
2422208
                              (cpm*rrv*t(il,j)*t(il,j)))
532
    ELSE
533
3752773
            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
534
            elij(il, i, j) = elij(il, i, j) + &
535
                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
536
3752773
                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
537
            elij(il, i, j) = elij(il, i, j) / &
538
                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
539
3752773
                              (cpm*rrv*t(il,j)*t(il,j)))
540
    ENDIF
541
!>jyg
542
6174981
            elij(il, i, j) = max(elij(il,i,j), 0.)
543
544
6174981
            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
545
546
6174981
            IF (j>i) THEN
547
2520477
              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
548
2520477
              awat = amax1(awat, 0.0)
549
            ELSE
550
              awat = 0.
551
            END IF
552
553
! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
554
! :         t(il,j))
555
556
!jyg<
557
!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
558
! Mixed draught temperature at level j
559

6174981
    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
560
2422208
          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))
561
2422208
          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
562
    ELSE
563
3752773
          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))
564
3752773
          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
565
    ENDIF
566
!>jyg
567
568
!IM 301008 end
569
570
! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
571
! :               i,j,hent(il,i,j),Sigij(il,i,j)
572
573
! --      ASij is the integral of P(F) over the relevant F interval
574
            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
575
6174981
                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
576
577
          END IF
578
        END IF
579
      END DO
580
!jyg!      DO k = 1, ntra
581
!jyg!        DO il = 1, ncum
582
!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
583
!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
584
!jyg!              lwork(il)) THEN
585
!jyg!            IF (Sij(il,i,j)>0.0) THEN
586
!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
587
!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
588
!jyg!            END IF
589
!jyg!          END IF
590
!jyg!        END DO
591
!jyg!      END DO
592
593
! --    If I=J (detrainement and entrainement at the same level), then only the
594
! --    adiabatic ascent part of the mixture is considered
595
54206
      IF (i==j) THEN
596
1314694
        DO il = 1, ncum
597
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
598


1311955
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
599
2739
              lwork(il)) THEN
600
586707
            IF (Sij(il,i,j)>0.0) THEN
601
!!              rti = qnk(il) - ep(il, i)*clw(il, i)
602
586707
              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
603
!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
604
              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
605
586707
                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
606
586707
              Qent(il, i, i) = rti
607
586707
              uent(il, i, i) = unk(il)
608
586707
              vent(il, i, i) = vnk(il)
609
586707
              hent(il, i, i) = hp(il, i)
610
586707
              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
611
586707
              Sigij(il, i, i) = 0.
612
            END IF
613
          END IF
614
        END DO
615
!jyg!        DO k = 1, ntra
616
!jyg!          DO il = 1, ncum
617
!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
618
!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
619
!jyg!                lwork(il)) THEN
620
!jyg!              IF (Sij(il,i,j)>0.0) THEN
621
!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
622
!jyg!              END IF
623
!jyg!            END IF
624
!jyg!          END DO
625
!jyg!        END DO
626
627
      END IF
628
629
! ---------------------------------------------------------------
630
2739
175 END DO        ! End loop on destination level "j"
631
! ---------------------------------------------------------------
632
633
1314694
    DO il = 1, ncum
634

1314694
      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
635
586707
        ASij(il) = amax1(1.0E-16, ASij(il))
636
!jyg+lluis<
637
!!        ASij(il) = 1.0/ASij(il)
638
586707
        ASij_inv(il) = 1.0/ASij(il)
639
!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
640
586707
        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
641
!>jyg+lluis
642
586707
        csum(il, i) = 0.0
643
      END IF
644
    END DO
645
646
76692
    DO j = minorig, nl
647
35499477
      DO il = 1, ncum
648
        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
649


35496738
            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
650
!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
651
7095306
          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
652
        END IF
653
      END DO
654
    END DO
655
656
76692
    DO j = minorig, nl
657
35499477
      DO il = 1, ncum
658
        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
659


35496738
            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
660
7095306
          csum(il, i) = csum(il, i) + Ment(il, i, j)
661
        END IF
662
      END DO
663
    END DO
664
665
1314694
    DO il = 1, ncum
666


1315699
      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
667
! cc     :     .and. csum(il,i).lt.m(il,i) ) then
668
41487
        nent(il, i) = 0
669
! cc        Ment(il,i,i)=m(il,i)
670
41487
        Ment(il, i, i) = 1.
671
!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
672
41487
        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
673
41487
        uent(il, i, i) = unk(il)
674
41487
        vent(il, i, i) = vnk(il)
675
41487
        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
676
41487
        IF (fl_cor_ebil .GE. 2) THEN
677
          hent(il, i, i) = hp(il,i)
678
          Sigij(il, i, i) = 0.0
679
        ELSE
680
41487
          Sij(il, i, i) = 0.0
681
        ENDIF
682
      END IF
683
    END DO ! il
684
685
!jyg!    DO j = 1, ntra
686
!jyg!      DO il = 1, ncum
687
!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
688
!jyg!! cc     :     .and. csum(il,i).lt.m(il,i) ) then
689
!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
690
!jyg!        END IF
691
!jyg!      END DO
692
!jyg!    END DO
693
694
! ---------------------------------------------------------------
695
144
789 END DO              ! End loop on origin level "i"
696
! ---------------------------------------------------------------
697
698
699
144
  RETURN
700
END SUBROUTINE cv3p_mixing
701