Changeset 2420


Ignore:
Timestamp:
Jan 5, 2016, 4:37:49 PM (9 years ago)
Author:
crio
Message:

Nouvelle option d'epluchage de l'ascendance adiabatique dans le schema d'Emanuel: epluchage fonction de B/w2 au lieu de w. S'active avec iflag_mix_adiab=1 (valeur par defaut iflag_mix_adiab=0). Fonctionne avec iflag_mix=0 et iflag_mix=1.
Correction de bugs dans le schema de convection pour le calcul de inb, cape et buoy (sous le meme flag pour l'instant).
New option for the erosion of the adiabatic ascent in the Emanuel scheme: erosion function of B/w2 instead of w. Activated by iflag_mix_adiab=1 (standard value iflag_mix_adiab=0). Should work with iflag_mix=0 and iflag_mix=1.
Various bug corrections in the convection scheme for the computation of inb, cape, buoy (protected by the same flag for now).

Location:
LMDZ5/trunk/libf/phylmd
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/YOMCST2.h

    r2201 r2420  
    11
    2       INTEGER choice, iflag_mix
     2      INTEGER choice, iflag_mix, iflag_mix_adiab
    33      REAL  gammas, alphas, betas, Fmax, qqa1, qqa2, qqa3, scut
    44      REAL  Qcoef1max,Qcoef2max,Supcrit1,Supcrit2
    … …  
    99     &               Qcoef1max,Qcoef2max,                               &
    1010     &               Supcrit1, Supcrit2,                                &
    11      &               choice,iflag_mix,coef_clos_ls
     11     &               choice,iflag_mix,coef_clos_ls,iflag_mix_adiab
    1212!$OMP THREADPRIVATE(/YOMCST2/)
    1313!    --------------------------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/conf_phys_m.F90

    r2415 r2420  
    143143    REAL, SAVE :: supcrit1_omp, supcrit2_omp
    144144    INTEGER, SAVE :: iflag_mix_omp
     145    INTEGER, SAVE :: iflag_mix_adiab_omp
    145146    real, save :: scut_omp, qqa1_omp, qqa2_omp, gammas_omp, Fmax_omp, alphas_omp
    146147    REAL, SAVE :: tmax_fonte_cv_omp
    … …  
    17551756    iflag_mix_omp = 1
    17561757    call getin('iflag_mix',iflag_mix_omp)
     1758
     1759!
     1760    ! PARAMETERS FOR THE EROSION OF THE ADIABATIC ASCENTS
     1761    ! iflag_mix_adiab: 0=OLD,
     1762    !                  1=NEW (CR),           
     1763    !           
     1764    !
     1765    !Config Key  = iflag_mix_adiab
     1766    !Config Desc =
     1767    !Config Def  = 1
     1768    !Config Help =
     1769    !
     1770    iflag_mix_adiab_omp = 0
     1771    call getin('iflag_mix_adiab',iflag_mix_adiab_omp)
    17571772
    17581773    !
    … …  
    21332148    supcrit2 = supcrit2_omp
    21342149    iflag_mix = iflag_mix_omp
     2150    iflag_mix_adiab = iflag_mix_adiab_omp
    21352151    scut = scut_omp
    21362152    qqa1 = qqa1_omp
    … …  
    23572373    write(lunout,*)' supcrit2 = ', supcrit2
    23582374    write(lunout,*)' iflag_mix = ', iflag_mix
     2375    write(lunout,*)' iflag_mix_adiab = ', iflag_mix_adiab
    23592376    write(lunout,*)' scut = ', scut
    23602377    write(lunout,*)' qqa1 = ', qqa1

  • LMDZ5/trunk/libf/phylmd/cv3_buoy.F90

    r1992 r2420  
    1414  include "cvthermo.h"
    1515  include "cv3param.h"
     16  include "YOMCST2.h"
    1617
    1718  ! input:
    … …  
    139140  END DO
    140141
    141 
     142!CR:Correction of buoy for what comes next
     143!keep flag or to modify in all cases?
     144  IF (iflag_mix_adiab.eq.1) THEN
     145  DO k = 1, nl
     146    DO il = 1, ncum
     147       IF ((k>=kmx(il)) .AND. (k<=inb(il)) .AND. (buoy(il,k).lt.0.)) THEN
     148          buoy(il,k)=buoy(il,k-1)
     149       END IF
     150    ENDDO
     151  ENDDO
     152  ENDIF
    142153
    143154  RETURN

  • LMDZ5/trunk/libf/phylmd/cv3_routines.F90

    r2398 r2420  
    10871087SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
    10881088                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
    1089                          p, h, tv, lv, lf, pbase, buoybase, plcl, &
     1089                         p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
    10901090                         inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
    10911091  IMPLICIT NONE
    … …  
    11131113  include "conema3.h"
    11141114  include "cvflag.h"
     1115  include "YOMCST2.h"
    11151116
    11161117!inputs:
    … …  
    11191120  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, q, qs, gz
    11201121  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
     1122  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
    11211123  REAL, DIMENSION (nloc), INTENT (IN)                :: tnk, qnk, gznk
    11221124  REAL, DIMENSION (nloc), INTENT (IN)                :: hnk
    … …  
    11441146  INTEGER iposit(nloc)
    11451147  REAL fracg
     1148  REAL deltap
    11461149
    11471150! =====================================================================
    … …  
    14761479    END DO
    14771480  END DO
     1481
     1482!CR fix computation of inb
     1483!keep flag or modify in all cases?
     1484  IF (iflag_mix_adiab.eq.1) THEN
     1485  DO i = 1, ncum
     1486     cape(i)=0.
     1487     inb(i)=icb(i)+1
     1488  ENDDO
     1489 
     1490  DO k = 2, nl
     1491    DO i = 1, ncum
     1492       IF ((k>=iposit(i))) THEN
     1493       deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
     1494       cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
     1495       IF (cape(i).gt.0.) THEN
     1496        inb(i) = max(inb(i), k)
     1497       END IF
     1498       ENDIF
     1499    ENDDO
     1500  ENDDO
     1501
     1502!  DO i = 1, ncum
     1503!     print*,"inb",inb(i)
     1504!  ENDDO
     1505
     1506  endif
    14781507
    14791508! -- end convect3

  • LMDZ5/trunk/libf/phylmd/cv3p1_closure.F90

    r2311 r2420  
    5858  INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
    5959  REAL deltap, fac, w, amu
    60   REAL rhodp
     60  REAL rhodp, dz
    6161  REAL pbmxup
    6262  REAL dtmin(nloc, nd), sigold(nloc, nd)
    … …  
    7979  REAL term1, term2, term3
    8080  REAL alp2(nloc) ! Alp with offset
     81  !CR: variables for new erosion of adiabiatic ascent
     82  REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd)
     83  REAL med(nloc, nd), md(nloc,nd)
     84  REAL coef_peel
     85  PARAMETER (coef_peel=0.25)
    8186
    8287  REAL sigmax
    … …  
    110115    END DO
    111116  END DO
     117
     118!CR: initializations for erosion of adiabatic ascent
     119  DO k = 1,nl
     120    DO il = 1, ncum
     121        mad(il,k)=0.
     122        me(il,k)=0.
     123        betalim(il,k)=1.
     124        wlim(il,k)=0.
     125    ENDDO
     126  ENDDO
    112127
    113128  ! -------------------------------------------------------
    … …  
    431446
    432447      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
    433 
     448        IF (iflag_mix_adiab.eq.1) THEN
     449!CR:computation of cape from LCL: keep flag or to modify in all cases?
     450        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
     451        ELSE
    434452        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
     453        ENDIF
    435454        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
    436455        cape(il) = amax1(0.0, cape(il))
    … …  
    601620  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
    602621
     622!CR: new erosion of adiabatic ascent: modification of m
     623!computation of the sum of ascending fluxes
     624  IF (iflag_mix_adiab.eq.1) THEN
     625
     626!Verification sum(me)=sum(m)
     627  DO k = 1,nl+1
     628    DO il = 1, ncum
     629       md(il,k)=0.
     630       med(il,k)=0.
     631    ENDDO
     632  ENDDO
     633
     634  DO k = nl,1,-1
     635    DO il = 1, ncum
     636           md(il,k)=md(il,k+1)+m(il,k+1)
     637    ENDDO
     638  ENDDO
     639
     640  DO k = nl,1,-1
     641    DO il = 1, ncum
     642        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
     643           mad(il,k)=mad(il,k+1)+m(il,k+1)
     644        ENDIF
     645!        print*,"mad",il,k,mad(il,k)
     646    ENDDO
     647  ENDDO
     648
     649!CR: erosion of each adiabatic ascent during its ascent
     650
     651!Computation of erosion coefficient beta_coef
     652  DO k = 1, nl
     653    DO il = 1, ncum
     654       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
     655!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
     656          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
     657       ELSE
     658          beta_coef(il,k)=0.
     659       ENDIF
     660    ENDDO
     661  ENDDO
     662
     663!  print*,"apres beta_coef"
     664
     665  DO k = 1, nl
     666    DO il = 1, ncum
     667
     668      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
     669
     670!        print*,"dz",il,k,tv(il, k-1)
     671        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
     672        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
     673!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
     674!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
     675        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
     676!        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))
     677        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
     678!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
     679     
     680      END IF
     681       
     682!Modification of m
     683      m(il,k)=me(il,k)
     684    END DO
     685  END DO
     686 
     687!  DO il = 1, ncum
     688!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
     689!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
     690!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
     691!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
     692!  ENDDO
     693
     694!Verification sum(me)=sum(m)
     695  DO k = nl,1,-1
     696    DO il = 1, ncum
     697           med(il,k)=med(il,k+1)+m(il,k+1)
     698!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
     699    ENDDO
     700  ENDDO
     701
     702
     703  ENDIF !(iflag_mix_adiab)
     704!RC
     705
     706
     707
    603708  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
    604709  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if

  • LMDZ5/trunk/libf/phylmd/cv3p2_closure.F90

    r2398 r2420  
    5959  REAL                                               :: deltap, fac, w, amu
    6060  REAL, DIMENSION (nloc, nd)                         :: rhodp               ! Factor such that m=rhodp*sig*w
     61  REAL                                               :: dz
    6162  REAL                                               :: pbmxup
    6263  REAL, DIMENSION (nloc, nd)                         :: dtmin, sigold
    … …  
    8586  REAL, DIMENSION (nloc)                             :: alp2                  ! Alp with offset
    8687
     88!CR: variables for new erosion of adiabiatic ascent
     89  REAL, DIMENSION (nloc, nd)                         :: mad, me, betalim, beta_coef
     90  REAL, DIMENSION (nloc, nd)                         :: med, md
     91  REAL                                               :: coef_peel
     92  PARAMETER (coef_peel=0.25)
     93
    8794  REAL                                               :: sigmax
    8895  PARAMETER (sigmax=0.1)
    … …  
    120127    END DO
    121128  END DO
     129
     130!CR: initializations for erosion of adiabatic ascent
     131  DO k = 1,nl
     132    DO il = 1, ncum
     133        mad(il,k)=0.
     134        me(il,k)=0.
     135        betalim(il,k)=1.
     136        wlim(il,k)=0.
     137    ENDDO
     138  ENDDO
    122139
    123140  ! -------------------------------------------------------
    … …  
    496513      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
    497514
     515        IF (iflag_mix_adiab.eq.1) THEN
     516!CR:computation of cape from LCL: keep flag or to modify in all cases?
     517        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
     518        ELSE
    498519        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
     520        ENDIF
    499521        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
    500522        cape(il) = amax1(0.0, cape(il))
    … …  
    681703                         (k,w0(igout,k),sig(igout,k), k=icb(igout),inb(igout))
    682704
     705!CR: new erosion of adiabatic ascent: modification of m
     706!computation of the sum of ascending fluxes
     707  IF (iflag_mix_adiab.eq.1) THEN
     708
     709!Verification sum(me)=sum(m)
     710  DO k = 1,nl+1
     711    DO il = 1, ncum
     712       md(il,k)=0.
     713       med(il,k)=0.
     714    ENDDO
     715  ENDDO
     716
     717  DO k = nl,1,-1
     718    DO il = 1, ncum
     719           md(il,k)=md(il,k+1)+m(il,k+1)
     720    ENDDO
     721  ENDDO
     722
     723  DO k = nl,1,-1
     724    DO il = 1, ncum
     725        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
     726           mad(il,k)=mad(il,k+1)+m(il,k+1)
     727        ENDIF
     728!        print*,"mad",il,k,mad(il,k)
     729    ENDDO
     730  ENDDO
     731
     732!CR: erosion of each adiabatic ascent during its ascent
     733
     734!Computation of erosion coefficient beta_coef
     735  DO k = 1, nl
     736    DO il = 1, ncum
     737       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
     738!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
     739          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
     740       ELSE
     741          beta_coef(il,k)=0.
     742       ENDIF
     743    ENDDO
     744  ENDDO
     745
     746!  print*,"apres beta_coef"
     747
     748  DO k = 1, nl
     749    DO il = 1, ncum
     750
     751      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
     752
     753!        print*,"dz",il,k,tv(il, k-1)
     754        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
     755        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
     756!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
     757!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
     758        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
     759!        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))
     760        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
     761!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
     762     
     763      END IF
     764       
     765!Modification of m
     766      m(il,k)=me(il,k)
     767    END DO
     768  END DO
     769 
     770!  DO il = 1, ncum
     771!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
     772!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
     773!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
     774!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
     775!  ENDDO
     776
     777!Verification sum(me)=sum(m)
     778  DO k = nl,1,-1
     779    DO il = 1, ncum
     780           med(il,k)=med(il,k+1)+m(il,k+1)
     781!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
     782    ENDDO
     783  ENDDO
     784
     785
     786  ENDIF !(iflag_mix_adiab)
     787!RC
     788
    683789  ! c 3. Compute final cloud base mass flux;
    684790  ! c    set iflag to 3 if cloud base mass flux is exceedingly small and is

  • LMDZ5/trunk/libf/phylmd/cva_driver.F90

    r2407 r2420  
    880880      CALL cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &              !na->nd
    881881                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
    882                          p, h, tv, lv, lf, pbase, buoybase, plcl, &
     882                         p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
    883883                         inb, tp, tvp, clw, hp, ep, sigp, buoy, &
    884884                         frac)
Note: See TracChangeset for help on using the changeset viewer.