source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p2_closure.F90 @ 5017

Last change on this file since 5017 was 3571, checked in by fhourdin, 5 years ago

iflag_wb devient un real pour pouvoir lui donner une valeur continue
quand on l'utilise comme valeur de wbmax en surface.
Pas très joli mais important pour le tuning automatique pour lequel on
genere des valeurs continues de wb_max.

File size: 26.3 KB
Line 
1
2
3SUBROUTINE cv3p2_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
4    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
5    iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmflast, plfc, &
6    wbeff)
7
8
9  ! **************************************************************
10  ! *
11  ! CV3P2_CLOSURE                                               *
12  ! Ale & Alp Closure of Convect3              *
13  ! *
14  ! written by   :   Kerry Emanuel                              *
15  ! vectorization:   S. Bony                                    *
16  ! modified by :    Jean-Yves Grandpeix, 18/06/2003, 19.32.10  *
17  ! Julie Frohwirth,     14/10/2005  17.44.22  *
18  ! **************************************************************
19
20  USE print_control_mod, ONLY: prt_level, lunout
21  IMPLICIT NONE
22
23  include "cvthermo.h"
24  include "cv3param.h"
25  include "cvflag.h"
26  include "YOMCST2.h"
27  include "YOMCST.h"
28  include "conema3.h"
29
30  ! input:
31  INTEGER, INTENT (IN)                               :: ncum, nd, nloc
32  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
33  REAL, DIMENSION (nloc), INTENT (IN)                :: pbase, plcl
34  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
35  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
36  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, buoy
37  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: supmax
38  LOGICAL, INTENT (IN)                               :: ok_inhib ! enable convection inhibition by dryness
39  REAL, DIMENSION (nloc), INTENT (IN)                :: ale, alp
40  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: omega
41
42  ! input/output:
43  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
44  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2
45
46  ! output:
47  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
48  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
49  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
50  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
51  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
52  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
53  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmflast, plfc
54  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff
55  INTEGER, DIMENSION (nloc), INTENT (OUT)            :: iflag
56
57  ! local variables:
58  INTEGER                                            :: il, i, j, k, icbmax
59  INTEGER, DIMENSION (nloc)                          :: i0, klfc
60  REAL                                               :: deltap, fac, w, amu
61  REAL, DIMENSION (nloc, nd)                         :: rhodp               ! Factor such that m=rhodp*sig*w
62  REAL                                               :: dz
63  REAL                                               :: pbmxup
64  REAL, DIMENSION (nloc, nd)                         :: dtmin, sigold
65  REAL, DIMENSION (nloc, nd)                         :: coefmix
66  REAL, DIMENSION (nloc)                             :: dtminmax
67  REAL, DIMENSION (nloc)                             :: pzero, ptop2old
68  REAL, DIMENSION (nloc)                             :: cina, cinb
69  INTEGER, DIMENSION (nloc)                          :: ibeg
70  INTEGER, DIMENSION (nloc)                          :: nsupmax
71  REAL                                               :: supcrit
72  REAL, DIMENSION (nloc, nd)                         :: temp
73  REAL, DIMENSION (nloc)                             :: p1, pmin
74  REAL, DIMENSION (nloc)                             :: asupmax0
75  LOGICAL, DIMENSION (nloc)                          :: ok
76  REAL, DIMENSION (nloc, nd)                         :: siglim, wlim, mlim
77  REAL, DIMENSION (nloc)                             :: wb2
78  REAL, DIMENSION (nloc)                             :: cbmf0        ! initial cloud base mass flux
79  REAL, DIMENSION (nloc)                             :: cbmflim      ! cbmf given by Cape closure
80  REAL, DIMENSION (nloc)                             :: cbmfalp      ! cbmf given by Alp closure
81  REAL, DIMENSION (nloc)                             :: cbmfalpb     ! bounded cbmf given by Alp closure
82  REAL, DIMENSION (nloc)                             :: cbmfmax      ! upper bound on cbmf
83  REAL, DIMENSION (nloc)                             :: coef
84  REAL, DIMENSION (nloc)                             :: xp, xq, xr, discr, b3, b4
85  REAL, DIMENSION (nloc)                             :: theta, bb
86  REAL                                               :: term1, term2, term3
87  REAL, DIMENSION (nloc)                             :: alp2                  ! Alp with offset
88
89!CR: variables for new erosion of adiabiatic ascent
90  REAL, DIMENSION (nloc, nd)                         :: mad, me, betalim, beta_coef
91  REAL, DIMENSION (nloc, nd)                         :: med, md
92!jyg<
93! coef_peel is now in the common cv3_param
94!!  REAL                                               :: coef_peel
95!!  PARAMETER (coef_peel=0.25)
96!>jyg
97
98  REAL                                               :: sigmax
99  PARAMETER (sigmax=0.1)
100!!  PARAMETER (sigmax=10.)
101
102  CHARACTER (LEN=20)                                 :: modname = 'cv3p2_closure'
103  CHARACTER (LEN=80)                                 :: abort_message
104
105  INTEGER,SAVE                                       :: igout=1
106!$OMP THREADPRIVATE(igout)
107
108 IF (prt_level>=20) print *,' -> cv3p2_closure, Ale ',ale(igout)
109
110
111  ! -------------------------------------------------------
112  ! -- Initialization
113  ! -------------------------------------------------------
114
115
116  DO il = 1, ncum
117    alp2(il) = max(alp(il), 1.E-5)
118    ! IM
119    alp2(il) = max(alp(il), 1.E-12)
120  END DO
121
122  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
123  ! exist (if any)
124
125  IF (prt_level>=20) PRINT *, 'cv3p2_closure nloc ncum nd icb inb nl', nloc, &
126    ncum, nd, icb(nloc), inb(nloc), nl
127  DO k = 1, nl
128    DO il = 1, ncum
129      rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
130    END DO
131  END DO
132
133!CR+jyg: initializations (up to nd) for erosion of adiabatic ascent and of m and wlim
134  DO k = 1,nd
135    DO il = 1, ncum
136        mad(il,k)=0.
137        me(il,k)=0.
138        betalim(il,k)=1.
139        wlim(il,k)=0.
140        m(il, k) = 0.0
141    ENDDO
142  ENDDO
143
144  ! -------------------------------------------------------
145  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
146  ! -------------------------------------------------------
147
148  ! update sig and w0 above LNB:
149
150  DO k = 1, nl - 1
151    DO il = 1, ncum
152      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
153        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il,inb(il)))
154        sig(il, k) = amax1(sig(il,k), 0.0)
155        w0(il, k) = beta*w0(il, k)
156      END IF
157    END DO
158  END DO
159
160  ! if(prt.level.GE.20) print*,'cv3p2_closure apres 100'
161  ! compute icbmax:
162
163  icbmax = 2
164  DO il = 1, ncum
165    icbmax = max(icbmax, icb(il))
166  END DO
167  ! if(prt.level.GE.20) print*,'cv3p2_closure apres 200'
168
169  ! update sig and w0 below cloud base:
170
171  DO k = 1, icbmax
172    DO il = 1, ncum
173      IF (k<=icb(il)) THEN
174        sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il,icb(il))
175        sig(il, k) = amax1(sig(il,k), 0.0)
176        w0(il, k) = beta*w0(il, k)
177      END IF
178    END DO
179  END DO
180  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 300'
181
182  ! -------------------------------------------------------------
183  ! -- Reset fractional areas of updrafts and w0 at initial time
184  ! -- and after 10 time steps of no convection
185  ! -------------------------------------------------------------
186
187!jyg<
188  IF (ok_convstop) THEN
189    DO k = 1, nl - 1
190      DO il = 1, ncum
191        IF (sig(il,nd)<1.5 .OR. sig(il,nd)>noconv_stop) THEN
192          sig(il, k) = 0.0
193          w0(il, k) = 0.0
194        END IF
195      END DO
196    END DO
197  ELSE
198  DO k = 1, nl - 1
199    DO il = 1, ncum
200      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
201        sig(il, k) = 0.0
202        w0(il, k) = 0.0
203      END IF
204    END DO
205  END DO
206  ENDIF  ! (ok_convstop)
207!>jyg
208  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 400'
209
210  ! -------------------------------------------------------
211  ! -- Compute initial cloud base mass flux (Cbmf0)
212  ! -------------------------------------------------------
213  DO il = 1, ncum
214    cbmf0(il) = 0.0
215  END DO
216
217  DO k = 1, nl
218    DO il = 1, ncum
219      IF (k>=icb(il) .AND. k<=inb(il) &
220          .AND. icb(il)+1<=inb(il)) THEN
221        cbmf0(il) = cbmf0(il) + sig(il, k)*w0(il,k)*rhodp(il,k)
222      END IF
223    END DO
224  END DO
225
226  ! -------------------------------------------------------------
227  ! jyg1
228  ! --  Calculate adiabatic ascent top pressure (ptop)
229  ! -------------------------------------------------------------
230
231
232  ! c 1. Start at first level where precipitations form
233  DO il = 1, ncum
234    pzero(il) = plcl(il) - pbcrit
235  END DO
236
237  ! c 2. Add offset
238  DO il = 1, ncum
239    pzero(il) = pzero(il) - pbmxup
240  END DO
241  DO il = 1, ncum
242    ptop2old(il) = ptop2(il)
243  END DO
244
245  DO il = 1, ncum
246    ! CR:c est quoi ce 300??
247    p1(il) = pzero(il) - 300.
248  END DO
249
250  ! compute asupmax=abs(supmax) up to lnm+1
251
252  DO il = 1, ncum
253    ok(il) = .TRUE.
254    nsupmax(il) = inb(il)
255  END DO
256
257  DO i = 1, nl
258    DO il = 1, ncum
259      IF (i>icb(il) .AND. i<=inb(il)) THEN
260        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
261          nsupmax(il) = i
262          ok(il) = .FALSE.
263        END IF ! end IF (P(i) ...  )
264      END IF ! end IF (icb+1 le i le inb)
265    END DO
266  END DO
267
268  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 2.'
269  DO i = 1, nl
270    DO il = 1, ncum
271      asupmax(il, i) = abs(supmax(il,i))
272    END DO
273  END DO
274
275
276  DO il = 1, ncum
277    asupmaxmin(il) = 10.
278    pmin(il) = 100.
279    ! IM ??
280    asupmax0(il) = 0.
281  END DO
282
283  ! c 3.  Compute in which level is Pzero
284
285  ! IM bug      i0 = 18
286  DO il = 1, ncum
287    i0(il) = nl
288  END DO
289
290  DO i = 1, nl
291    DO il = 1, ncum
292      IF (i>icb(il) .AND. i<=inb(il)) THEN
293        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
294          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
295            i0(il) = i
296          END IF
297        END IF
298      END IF
299    END DO
300  END DO
301  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 3.'
302
303  ! c 4.  Compute asupmax at Pzero
304
305  DO i = 1, nl
306    DO il = 1, ncum
307      IF (i>icb(il) .AND. i<=inb(il)) THEN
308        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
309          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))- &
310            (pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il,i0(il)-1))
311        END IF
312      END IF
313    END DO
314  END DO
315
316
317  DO i = 1, nl
318    DO il = 1, ncum
319      IF (p(il,i)==pzero(il)) THEN
320        asupmax(i, il) = asupmax0(il)
321      END IF
322    END DO
323  END DO
324  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 4.'
325
326  ! c 5. Compute asupmaxmin, minimum of asupmax
327
328  DO i = 1, nl
329    DO il = 1, ncum
330      IF (i>icb(il) .AND. i<=inb(il)) THEN
331        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
332          IF (asupmax(il,i)<asupmaxmin(il)) THEN
333            asupmaxmin(il) = asupmax(il, i)
334            pmin(il) = p(il, i)
335          END IF
336        END IF
337      END IF
338    END DO
339  END DO
340
341  DO il = 1, ncum
342    ! IM
343    IF (prt_level>=20) THEN
344      PRINT *, 'cv3p2_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
345        asupmaxmin(il), pzero(il), pmin(il)
346    END IF
347    IF (asupmax0(il)<asupmaxmin(il)) THEN
348      asupmaxmin(il) = asupmax0(il)
349      pmin(il) = pzero(il)
350    END IF
351  END DO
352  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 5.'
353
354
355  ! Compute Supmax at Pzero
356
357  DO i = 1, nl
358    DO il = 1, ncum
359      IF (i>icb(il) .AND. i<=inb(il)) THEN
360        IF (p(il,i)<=pzero(il)) THEN
361          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)- &
362            (p(il,i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
363          GO TO 425
364        END IF ! end IF (P(i) ... )
365      END IF ! end IF (icb+1 le i le inb)
366    END DO
367  END DO
368
369425 CONTINUE
370  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 425.'
371
372  ! c 6. Calculate ptop2
373
374  DO il = 1, ncum
375    IF (asupmaxmin(il)<supcrit1) THEN
376      ptop2(il) = pmin(il)
377    END IF
378
379    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
380      ptop2(il) = ptop2old(il)
381    END IF
382
383    IF (asupmaxmin(il)>supcrit2) THEN
384      ptop2(il) = ph(il, inb(il))
385    END IF
386  END DO
387
388  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 6.'
389
390  ! c 7. Compute multiplying factor for adiabatic updraught mass flux
391
392
393  IF (ok_inhib) THEN
394
395    DO i = 1, nl
396      DO il = 1, ncum
397        IF (i<=nl) THEN
398          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph(il,i))
399          coefmix(il, i) = min(coefmix(il,i), 1.)
400        END IF
401      END DO
402    END DO
403
404
405  ELSE ! when inhibition is not taken into account, coefmix=1
406
407
408
409    DO i = 1, nl
410      DO il = 1, ncum
411        IF (i<=nl) THEN
412          coefmix(il, i) = 1.
413        END IF
414      END DO
415    END DO
416
417  END IF ! ok_inhib
418  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 7.'
419  ! -------------------------------------------------------------------
420  ! -------------------------------------------------------------------
421
422
423  ! jyg2
424
425  ! ==========================================================================
426
427
428  ! -------------------------------------------------------------
429  ! -- Calculate convective inhibition (CIN)
430  ! -------------------------------------------------------------
431
432  ! do i=1,nloc
433  ! print*,'avant cine p',pbase(i),plcl(i)
434  ! enddo
435  ! do j=1,nd
436  ! do i=1,nloc
437  ! print*,'avant cine t',tv(i),tvp(i)
438  ! enddo
439  ! enddo
440  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
441    cinb, plfc)
442
443  DO il = 1, ncum
444    cin(il) = cina(il) + cinb(il)
445  END DO
446  IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_cine: cina, cinb, cin ', &
447                              cina(igout), cinb(igout), cin(igout)
448  ! -------------------------------------------------------------
449  ! --Update buoyancies to account for Ale
450  ! -------------------------------------------------------------
451
452  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
453    tvp, buoy)
454  IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_buoy'
455
456  ! -------------------------------------------------------------
457  ! -- Calculate convective available potential energy (cape),
458  ! -- vertical velocity (w), fractional area covered by
459  ! -- undilute updraft (sig), and updraft mass flux (m)
460  ! -------------------------------------------------------------
461
462  DO il = 1, ncum
463    cape(il) = 0.0
464    dtminmax(il) = -100.
465  END DO
466
467  ! compute dtmin (minimum buoyancy between ICB and given level k):
468
469  DO k = 1, nl
470    DO il = 1, ncum
471      dtmin(il, k) = 100.0
472    END DO
473  END DO
474
475  DO k = 1, nl
476    DO j = minorig, nl
477      DO il = 1, ncum
478        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) &
479                             .AND. (j<=(k-1))) THEN
480          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
481        END IF
482      END DO
483    END DO
484  END DO
485!jyg<
486!  Store maximum of dtmin
487!  C est pas terrible d avoir ce test sur Ale+Cin encore une fois ici.
488!                      A REVOIR !
489  DO k = 1, nl
490    DO il = 1, ncum
491      IF (k>=(icb(il)+1) .AND. k<=inb(il) .AND. ale(il)+cin(il)>0.) THEN
492        dtminmax(il) = max(dtmin(il,k), dtminmax(il))
493      ENDIF
494    END DO
495  END DO
496!
497!    prevent convection when ale+cin <= 0
498  DO k = 1, nl
499    DO il = 1, ncum
500      IF (k>=(icb(il)+1) .AND. k<=inb(il)) THEN
501        dtmin(il,k) = min(dtmin(il,k), dtminmax(il))
502      ENDIF
503    END DO
504  END DO
505!>jyg
506!
507  IF (prt_level >= 20) THEN
508    print *,'cv3p2_closure: dtmin ', (k, dtmin(igout,k), k=1,nl)
509    print *,'cv3p2_closure: dtminmax ', dtminmax(igout)
510  ENDIF
511!
512  ! the interval on which cape is computed starts at pbase :
513
514  DO k = 1, nl
515    DO il = 1, ncum
516
517      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
518
519        IF (iflag_mix_adiab.eq.1) THEN
520!CR:computation of cape from LCL: keep flag or to modify in all cases?
521        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
522        ELSE
523        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
524        ENDIF
525        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
526        cape(il) = amax1(0.0, cape(il))
527        sigold(il, k) = sig(il, k)
528
529
530        ! jyg       Coefficient coefmix limits convection to levels where a
531        ! sufficient
532        ! fraction of mixed draughts are ascending.
533        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
534        siglim(il, k) = amax1(siglim(il,k), 0.0)
535        siglim(il, k) = amin1(siglim(il,k), 0.01)
536        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
537        fac = 1.
538        wlim(il, k) = fac*sqrt(cape(il))
539        amu = siglim(il, k)*wlim(il, k)
540!!        rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) !cor jyg : computed earlier
541        mlim(il, k) = amu*rhodp(il,k)
542        ! print*, 'siglim ', k,siglim(1,k)
543      END IF
544
545    END DO
546  END DO
547  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 600'
548
549  DO il = 1, ncum
550    ! IM beg
551    IF (prt_level>=20) THEN
552      PRINT *, 'cv3p2_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
553        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
554        ph(il, icb(il)+2)
555    END IF
556
557    IF (icb(il)+1<=inb(il)) THEN
558      ! IM end
559      mlim(il, icb(il)) = 0.5*mlim(il,icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
560                                               (ph(il,icb(il)+1)-ph(il,icb(il)+2))
561      ! IM beg
562    END IF !(icb(il.le.inb(il))) then
563    ! IM end
564  END DO
565  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 700'
566
567  !
568  ! ------------------------------------------------------------------------
569  ! c     Compute Cloud base mass flux given by Cape closure (cbmflim = cbmf of
570  ! c     elementary systems), cbmf given by Alp closure (cbmfalp), cbmf given by Alp
571  ! c     closure with an upper bound imposed (cbmfalpb) and cbmf resulting from
572  ! c     time integration (cbmflast).
573  ! ------------------------------------------------------------------------
574
575  DO il = 1, ncum
576    cbmflim(il) = 0.
577    cbmfalp(il) = 0.
578    cbmfalpb(il) = 0.
579    cbmflast(il) = 0.
580  END DO
581
582  ! c 1. Compute cloud base mass flux of elementary system (Cbmflim)
583
584  DO k = 1, nl
585    DO il = 1, ncum
586      ! old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
587      ! IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
588      IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
589          .AND. icb(il)+1<=inb(il)) THEN !cor jyg
590        cbmflim(il) = cbmflim(il) + mlim(il, k)
591      END IF
592    END DO
593  END DO
594  IF (prt_level>=20) PRINT *, 'cv3p2_closure after cbmflim: cbmflim ', cbmflim(igout)
595
596  ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmfalp), maximum
597  !     allowed mass flux (Cbmfmax) and bounded mass flux (Cbmfalpb)
598  !     Cbmfalpb is set to zero if Cbmflim (the mass flux of elementary cloud)
599  !     is exceedingly small.
600
601  DO il = 1, ncum
602    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
603  END DO
604
605  DO il = 1, ncum
606    IF (plfc(il)<100.) THEN
607      ! This is an irealistic value for plfc => no calculation of wbeff
608      wbeff(il) = 100.1
609    ELSE
610      ! Calculate wbeff
611      IF (NINT(flag_wb)==0) THEN
612        wbeff(il) = wbmax
613      ELSE IF (NINT(flag_wb)==1) THEN
614        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
615      ELSE IF (NINT(flag_wb)==2) THEN
616        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
617      END IF
618    END IF
619  END DO
620
621!CR:Compute k at plfc
622  DO il=1,ncum
623           klfc(il)=nl
624  ENDDO
625  DO k=1,nl
626     DO il=1,ncum
627        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
628           klfc(il)=k
629        endif
630     ENDDO
631  ENDDO
632!RC
633
634  DO il = 1, ncum
635    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
636    ! c       cbmfalp(il) = alp2(il)/(0.5*wb*wb-Cin(il))
637    cbmfalp(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
638!CR: Add large-scale component to the mass-flux
639!encore connu sous le nom "Experience du tube de dentifrice"
640    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
641       cbmfalp(il) = cbmfalp(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
642    endif
643!RC
644    IF (cbmfalp(il)==0 .AND. alp2(il)/=0.) THEN
645      WRITE (lunout, *) 'cv3p2_closure cbmfalp=0 and alp NE 0 il alp2 alp cin ' , &
646                         il, alp2(il), alp(il), cin(il)
647      abort_message = ''
648      CALL abort_physic(modname, abort_message, 1)
649    END IF
650    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
651  END DO
652
653!jyg<
654  IF (OK_intermittent) THEN
655    DO il = 1, ncum
656      IF (cbmflim(il)>1.E-6) THEN
657        cbmfalpb(il) = min(cbmfalp(il), (cbmfmax(il)-beta*cbmf0(il))/(1.-beta))
658        ! print*,'cbmfalpb',cbmfalpb(il),cbmfmax(il)
659      END IF
660    END DO
661  ELSE
662!>jyg
663  DO il = 1, ncum
664    IF (cbmflim(il)>1.E-6) THEN
665      ! ATTENTION TEST CR
666      ! if (cbmfmax(il).lt.1.e-12) then
667      cbmfalpb(il) = min(cbmfalp(il), cbmfmax(il))
668      ! else
669      ! cbmfalpb(il) = cbmfalp(il)
670      ! endif
671      ! print*,'cbmfalpb',cbmfalp(il),cbmfmax(il)
672    END IF
673  END DO
674  ENDIF  !(OK_intermittent)
675  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmfalpb: cbmfalpb ',cbmfalpb(igout)
676
677  ! c 2. Compute coefficient and apply correction
678
679  DO il = 1, ncum
680    coef(il) = (cbmfalpb(il)+1.E-10)/(cbmflim(il)+1.E-10)
681  END DO
682  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres coef_plantePLUS'
683
684     DO k = 1, nl
685       DO il = 1, ncum
686         IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
687           amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)*wlim(il, k)
688           w0(il, k) = wlim(il, k)
689           w0(il, k) = max(w0(il,k), 1.E-10)
690           sig(il, k) = amu/w0(il, k)
691           sig(il, k) = min(sig(il,k), 1.)
692           ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
693           !jyg m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
694           m(il, k) = amu*rhodp(il,k)
695         END IF
696       END DO
697     END DO
698  ! jyg2
699  DO il = 1, ncum
700    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
701    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
702                                         (ph(il,icb(il)+1)-ph(il,icb(il)+2))
703    sig(il, icb(il)) = sig(il, icb(il)+1)
704    sig(il, icb(il)-1) = sig(il, icb(il))
705  END DO
706  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres w0_sig_M: w0, sig ', &
707                         (k,w0(igout,k),sig(igout,k), k=icb(igout),inb(igout))
708
709!CR: new erosion of adiabatic ascent: modification of m
710!computation of the sum of ascending fluxes
711  IF (iflag_mix_adiab.eq.1) THEN
712
713!Verification sum(me)=sum(m)
714  DO k = 1,nd
715    DO il = 1, ncum
716       md(il,k)=0.
717       med(il,k)=0.
718    ENDDO
719  ENDDO
720
721  DO k = nl,1,-1
722    DO il = 1, ncum
723           md(il,k)=md(il,k+1)+m(il,k+1)
724    ENDDO
725  ENDDO
726
727  DO k = nl,1,-1
728    DO il = 1, ncum
729        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
730           mad(il,k)=mad(il,k+1)+m(il,k+1)
731        ENDIF
732!        print*,"mad",il,k,mad(il,k)
733    ENDDO
734  ENDDO
735
736!CR: erosion of each adiabatic ascent during its ascent
737
738!Computation of erosion coefficient beta_coef
739  DO k = 1, nl
740    DO il = 1, ncum
741       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
742!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
743          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
744       ELSE
745          beta_coef(il,k)=0.
746       ENDIF
747    ENDDO
748  ENDDO
749
750!  print*,"apres beta_coef"
751
752  DO k = 1, nl
753    DO il = 1, ncum
754
755      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
756
757!        print*,"dz",il,k,tv(il, k-1)
758        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
759        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
760!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
761!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
762        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
763!        me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k))
764        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
765!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
766     
767      END IF
768       
769!Modification of m
770      m(il,k)=me(il,k)
771    END DO
772  END DO
773 
774!  DO il = 1, ncum
775!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
776!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
777!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
778!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
779!  ENDDO
780
781!Verification sum(me)=sum(m)
782  DO k = nl,1,-1
783    DO il = 1, ncum
784           med(il,k)=med(il,k+1)+m(il,k+1)
785!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
786    ENDDO
787  ENDDO
788
789
790  ENDIF !(iflag_mix_adiab)
791!RC
792
793  ! c 3. Compute final cloud base mass flux;
794  ! c    set iflag to 3 if cloud base mass flux is exceedingly small and is
795  ! c     decreasing (i.e. if the final mass flux (cbmflast) is greater than
796  ! c     the target mass flux (cbmfalpb)).
797  ! c    If(ok_convstop): set iflag to 4 if no positive buoyancy has been met
798
799!jyg  DO il = 1, ncum
800!jyg    cbmflast(il) = 0.
801!jyg  END DO
802
803  DO k = 1, nl
804    DO il = 1, ncum
805      IF (k>=icb(il) .AND. k<=inb(il)) THEN
806          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
807        cbmflast(il) = cbmflast(il) + m(il, k)
808      END IF
809    END DO
810  END DO
811  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmflast: cbmflast ',cbmflast(igout)
812
813  DO il = 1, ncum
814    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmfalpb(il)) THEN
815      iflag(il) = 3
816    END IF
817  END DO
818
819!jyg<
820  IF (ok_convstop) THEN
821    DO il = 1, ncum
822      IF (dtminmax(il) .LE. 0.) THEN
823        iflag(il) = 4
824      END IF
825    END DO
826  ELSE
827!>jyg
828  DO k = 1, nl
829    DO il = 1, ncum
830      IF (iflag(il)>=3) THEN
831        m(il, k) = 0.
832        sig(il, k) = 0.
833        w0(il, k) = 0.
834      END IF
835    END DO
836  END DO
837  ENDIF ! (ok_convstop)
838!
839  IF (prt_level >= 10) THEN
840   print *,'cv3p2_closure: iflag ',iflag(igout)
841  ENDIF
842!
843
844  ! c 4. Introduce a correcting factor for coef, in order to obtain an
845  ! effective
846  ! c    sigdz larger in the present case (using cv3p2_closure) than in the
847  ! old
848  ! c    closure (using cv3_closure).
849  IF (1==0) THEN
850    DO il = 1, ncum
851      ! c      coef(il) = 2.*coef(il)
852      coef(il) = 5.*coef(il)
853    END DO
854    ! version CVS du ..2008
855  ELSE
856    IF (iflag_cvl_sigd==0) THEN
857      ! test pour verifier qu on fait la meme chose qu avant: sid constant
858      coef(1:ncum) = 1.
859    ELSE
860      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
861      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
862    END IF
863  END IF
864
865  IF (prt_level>=20) PRINT *, 'cv3p2_closure FIN'
866  RETURN
867END SUBROUTINE cv3p2_closure
868
869
Note: See TracBrowser for help on using the repository browser.