Ignore:
Timestamp:
Jul 29, 2024, 12:37:08 PM (3 months ago)
Author:
abarral
Message:

Put cvthermo.h, cv30param.h, cv3param.h into modules

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/cv_routines.F90

    r5117 r5141  
    1 
    21! $Id$
    32
    … …  
    3837  include "cvparam.h"
    3938  INTEGER nd
    40   CHARACTER (LEN=20) :: modname = 'cv_routines'
    41   CHARACTER (LEN=80) :: abort_message
     39  CHARACTER (LEN = 20) :: modname = 'cv_routines'
     40  CHARACTER (LEN = 80) :: abort_message
    4241
    4342  ! noff: integer limit for convection (nd-noff)
    … …  
    7170  delta = 0.01 ! cld
    7271
    73 
    7472END SUBROUTINE cv_param
    7573
    7674SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
     75  USE lmdz_cvthermo
     76
    7777  IMPLICIT NONE
    7878
    … …  
    9393  REAL cpx(len, nd)
    9494
    95   include "cvthermo.h"
    9695  include "cvparam.h"
    97 
    9896
    9997  DO k = 1, nlp
    10098    DO i = 1, len
    101       lv(i, k) = lv0 - clmcpv*(t(i,k)-t0)
    102       cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
    103       cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
    104       tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
     99      lv(i, k) = lv0 - clmcpv * (t(i, k) - t0)
     100      cpn(i, k) = cpd * (1.0 - q(i, k)) + cpv * q(i, k)
     101      cpx(i, k) = cpd * (1.0 - q(i, k)) + cl * q(i, k)
     102      tv(i, k) = t(i, k) * (1.0 + q(i, k) * epsim1)
    105103    END DO
    106104  END DO
    … …  
    113111  DO k = 2, nlp
    114112    DO i = 1, len
    115       gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
    116         k)
     113      gz(i, k) = gz(i, k - 1) + hrd * (tv(i, k - 1) + tv(i, k)) * (p(i, k - 1) - p(i, k)) / ph(i, &
     114              k)
    117115    END DO
    118116  END DO
    … …  
    123121  DO k = 1, nlp
    124122    DO i = 1, len
    125       h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
    126       hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
    127     END DO
    128   END DO
    129 
     123      h(i, k) = gz(i, k) + cpn(i, k) * t(i, k)
     124      hm(i, k) = gz(i, k) + cpx(i, k) * (t(i, k) - t(i, 1)) + lv(i, k) * q(i, k)
     125    END DO
     126  END DO
    130127
    131128END SUBROUTINE cv_prelim
    132129
    133130SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
    134     qnk, gznk, plcl)
     131        qnk, gznk, plcl)
    135132  IMPLICIT NONE
    136133
    … …  
    169166  DO k = 2, nlp
    170167    DO i = 1, len
    171       IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
     168      IF ((hm(i, k)<work(i)) .AND. (hm(i, k)<hm(i, k - 1))) THEN
    172169        work(i) = hm(i, k)
    173170        ihmin(i) = k
    … …  
    193190  DO k = minorig + 1, nl
    194191    DO i = 1, len
    195       IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
     192      IF ((hm(i, k)>work(i)) .AND. (k<=ihmin(i))) THEN
    196193        work(i) = hm(i, k)
    197194        nk(i) = k
    … …  
    204201  ! -------------------------------------------------------------------
    205202  DO i = 1, len
    206     IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
    207       400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
     203    IF (((t(i, nk(i))<250.0) .OR. (q(i, nk(i))<=0.0) .OR. (p(i, ihmin(i))< &
     204            400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
    208205  END DO
    209206  ! -------------------------------------------------------------------
    … …  
    218215    qsnk(i) = qs(i, nk(i))
    219216
    220     rh(i) = qnk(i)/qsnk(i)
     217    rh(i) = qnk(i) / qsnk(i)
    221218    rh(i) = min(1.0, rh(i))
    222     chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
    223     plcl(i) = pnk(i)*(rh(i)**chi(i))
     219    chi(i) = tnk(i) / (1669.0 - 122.0 * rh(i) - tnk(i))
     220    plcl(i) = pnk(i) * (rh(i)**chi(i))
    224221    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
    225       ) = 8
     222            ) = 8
    226223  END DO
    227224  ! -------------------------------------------------------------------
    … …  
    234231  DO k = minorig, nl
    235232    DO i = 1, len
    236       IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
     233      IF ((k>=(nk(i) + 1)) .AND. (p(i, k)<plcl(i))) icb(i) = min(icb(i), k)
    237234    END DO
    238235  END DO
    … …  
    249246  END DO
    250247
    251 
    252248END SUBROUTINE cv_feed
    253249
    254250SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
    255     clw)
     251        clw)
     252  USE lmdz_cvthermo
     253
    256254  IMPLICIT NONE
    257255
    258   include "cvthermo.h"
    259256  include "cvparam.h"
    260257
    … …  
    292289
    293290  DO i = 1, len
    294     ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
    295       273.15)) + gznk(i)
    296     cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
     291    ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
     292            273.15)) + gznk(i)
     293    cpp(i) = cpd * (1. - qnk(i)) + qnk(i) * cpv
    297294  END DO
    298295
    … …  
    301298  DO k = minorig, icbmax - 1
    302299    DO i = 1, len
    303       tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
    304       tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
     300      tp(i, k) = tnk(i) - (gz(i, k) - gznk(i)) / cpp(i)
     301      tvp(i, k) = tp(i, k) * (1. + qnk(i) * epsi)
    305302    END DO
    306303  END DO
    … …  
    311308    tg = ticb(i)
    312309    qg = qs(i, icb(i))
    313     alv = lv0 - clmcpv*(ticb(i)-t0)
     310    alv = lv0 - clmcpv * (ticb(i) - t0)
    314311
    315312    ! First iteration.
    316313
    317     s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
    318     s = 1./s
    319     ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
    320     tg = tg + s*(ah0(i)-ahg)
     314    s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
     315    s = 1. / s
     316    ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
     317    tg = tg + s * (ah0(i) - ahg)
    321318    tg = max(tg, 35.0)
    322319    tc = tg - t0
    323320    denom = 243.5 + tc
    324321    IF (tc>=0.0) THEN
    325       es = 6.112*exp(17.67*tc/denom)
     322      es = 6.112 * exp(17.67 * tc / denom)
    326323    ELSE
    327       es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     324      es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
    328325    END IF
    329     qg = eps*es/(p(i,icb(i))-es*(1.-eps))
     326    qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
    330327
    331328    ! Second iteration.
    332329
    333     s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
    334     s = 1./s
    335     ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
    336     tg = tg + s*(ah0(i)-ahg)
     330    s = cpd + alv * alv * qg / (rrv * ticb(i) * ticb(i))
     331    s = 1. / s
     332    ahg = cpd * tg + (cl - cpd) * qnk(i) * ticb(i) + alv * qg + gzicb(i)
     333    tg = tg + s * (ah0(i) - ahg)
    337334    tg = max(tg, 35.0)
    338335    tc = tg - t0
    339336    denom = 243.5 + tc
    340337    IF (tc>=0.0) THEN
    341       es = 6.112*exp(17.67*tc/denom)
     338      es = 6.112 * exp(17.67 * tc / denom)
    342339    ELSE
    343       es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     340      es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
    344341    END IF
    345     qg = eps*es/(p(i,icb(i))-es*(1.-eps))
    346 
    347     alv = lv0 - clmcpv*(ticb(i)-273.15)
    348     tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
     342    qg = eps * es / (p(i, icb(i)) - es * (1. - eps))
     343
     344    alv = lv0 - clmcpv * (ticb(i) - 273.15)
     345    tp(i, icb(i)) = (ah0(i) - (cl - cpd) * qnk(i) * ticb(i) - gz(i, icb(i)) - alv * qg) / cpd
    349346    clw(i, icb(i)) = qnk(i) - qg
    350     clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
    351     rg = qg/(1.-qnk(i))
    352     tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
     347    clw(i, icb(i)) = max(0.0, clw(i, icb(i)))
     348    rg = qg / (1. - qnk(i))
     349    tvp(i, icb(i)) = tp(i, icb(i)) * (1. + rg * epsi)
    353350  END DO
    354351
    355352  DO k = minorig, icbmax
    356353    DO i = 1, len
    357       tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
    358     END DO
    359   END DO
    360 
     354      tvp(i, k) = tvp(i, k) - tp(i, k) * qnk(i)
     355    END DO
     356  END DO
    361357
    362358END SUBROUTINE cv_undilute1
    … …  
    383379  INTEGER i
    384380
    385 
    386381  DO i = 1, len
    387382    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
    388       icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
    389   END DO
    390 
     383            icb(i))<=(tv(i, icb(i)) - dtmax))) iflag(i) = 4
     384  END DO
    391385
    392386END SUBROUTINE cv_trigger
    393387
    394388SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
    395     tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
    396     tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
    397     v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
     389        tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
     390        tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
     391        v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
    398392  USE lmdz_print_control, ONLY: lunout
    399393  USE lmdz_abort_physic, ONLY: abort_physic
    … …  
    408402  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
    409403  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
    410   REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
     404  REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd)
    411405  REAL tvp1(len, nd), clw1(len, nd)
    412406
    … …  
    416410  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
    417411  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
    418   REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
     412  REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd)
    419413  REAL tvp(nloc, nd), clw(nloc, nd)
    420414  REAL dph(nloc, nd)
    … …  
    422416  ! local variables:
    423417  INTEGER i, k, nn
    424   CHARACTER (LEN=20) :: modname = 'cv_compress'
    425   CHARACTER (LEN=80) :: abort_message
    426 
     418  CHARACTER (LEN = 20) :: modname = 'cv_compress'
     419  CHARACTER (LEN = 80) :: abort_message
    427420
    428421  DO k = 1, nl + 1
    … …  
    473466  DO k = 1, nl
    474467    DO i = 1, ncum
    475       dph(i, k) = ph(i, k) - ph(i, k+1)
    476     END DO
    477   END DO
    478 
     468      dph(i, k) = ph(i, k) - ph(i, k + 1)
     469    END DO
     470  END DO
    479471
    480472END SUBROUTINE cv_compress
    481473
    482474SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
    483     gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
     475        gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
     476  USE lmdz_cvthermo
     477
    484478  IMPLICIT NONE
    485479
    … …  
    494488  ! ---------------------------------------------------------------------
    495489
    496   include "cvthermo.h"
    497490  include "cvparam.h"
    498491
    … …  
    538531  ! ***  Calculate certain parcel quantities, including static energy   ***
    539532
    540 
    541   DO i = 1, ncum
    542     ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
    543       t0)) + gznk(i)
     533  DO i = 1, ncum
     534    ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - &
     535            t0)) + gznk(i)
    544536  END DO
    545537
    … …  
    547539  ! ***  Find lifted parcel quantities above cloud base    ***
    548540
    549 
    550541  DO k = minorig + 1, nl
    551542    DO i = 1, ncum
    552       IF (k>=(icb(i)+1)) THEN
     543      IF (k>=(icb(i) + 1)) THEN
    553544        tg = t(i, k)
    554545        qg = qs(i, k)
    555         alv = lv0 - clmcpv*(t(i,k)-t0)
     546        alv = lv0 - clmcpv * (t(i, k) - t0)
    556547
    557548        ! First iteration.
    558549
    559         s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
    560         s = 1./s
    561         ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
    562         tg = tg + s*(ah0(i)-ahg)
     550        s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
     551        s = 1. / s
     552        ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
     553        tg = tg + s * (ah0(i) - ahg)
    563554        tg = max(tg, 35.0)
    564555        tc = tg - t0
    565556        denom = 243.5 + tc
    566557        IF (tc>=0.0) THEN
    567           es = 6.112*exp(17.67*tc/denom)
     558          es = 6.112 * exp(17.67 * tc / denom)
    568559        ELSE
    569           es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     560          es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
    570561        END IF
    571         qg = eps*es/(p(i,k)-es*(1.-eps))
     562        qg = eps * es / (p(i, k) - es * (1. - eps))
    572563
    573564        ! Second iteration.
    574565
    575         s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
    576         s = 1./s
    577         ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
    578         tg = tg + s*(ah0(i)-ahg)
     566        s = cpd + alv * alv * qg / (rrv * t(i, k) * t(i, k))
     567        s = 1. / s
     568        ahg = cpd * tg + (cl - cpd) * qnk(i) * t(i, k) + alv * qg + gz(i, k)
     569        tg = tg + s * (ah0(i) - ahg)
    579570        tg = max(tg, 35.0)
    580571        tc = tg - t0
    581572        denom = 243.5 + tc
    582573        IF (tc>=0.0) THEN
    583           es = 6.112*exp(17.67*tc/denom)
     574          es = 6.112 * exp(17.67 * tc / denom)
    584575        ELSE
    585           es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     576          es = exp(23.33086 - 6111.72784 / tg + 0.15215 * log(tg))
    586577        END IF
    587         qg = eps*es/(p(i,k)-es*(1.-eps))
    588 
    589         alv = lv0 - clmcpv*(t(i,k)-t0)
     578        qg = eps * es / (p(i, k) - es * (1. - eps))
     579
     580        alv = lv0 - clmcpv * (t(i, k) - t0)
    590581        ! PRINT*,'cpd dans convect2 ',cpd
    591582        ! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
    592583        ! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
    593         tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
     584        tp(i, k) = (ah0(i) - (cl - cpd) * qnk(i) * t(i, k) - gz(i, k) - alv * qg) / cpd
    594585        ! if (.NOT.cpd.gt.1000.) THEN
    595586        ! PRINT*,'CPD=',cpd
    … …  
    597588        ! END IF
    598589        clw(i, k) = qnk(i) - qg
    599         clw(i, k) = max(0.0, clw(i,k))
    600         rg = qg/(1.-qnk(i))
    601         tvp(i, k) = tp(i, k)*(1.+rg*epsi)
     590        clw(i, k) = max(0.0, clw(i, k))
     591        rg = qg / (1. - qnk(i))
     592        tvp(i, k) = tp(i, k) * (1. + rg * epsi)
    602593      END IF
    603594    END DO
    … …  
    612603  DO k = minorig + 1, nl
    613604    DO i = 1, ncum
    614       IF (k>=(nk(i)+1)) THEN
     605      IF (k>=(nk(i) + 1)) THEN
    615606        tca = tp(i, k) - t0
    616607        IF (tca>=0.0) THEN
    617608          elacrit = elcrit
    618609        ELSE
    619           elacrit = elcrit*(1.0-tca/tlcrit)
     610          elacrit = elcrit * (1.0 - tca / tlcrit)
    620611        END IF
    621612        elacrit = max(elacrit, 0.0)
    622         ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
    623         ep(i, k) = max(ep(i,k), 0.0)
    624         ep(i, k) = min(ep(i,k), 1.0)
     613        ep(i, k) = 1.0 - elacrit / max(clw(i, k), 1.0E-8)
     614        ep(i, k) = max(ep(i, k), 0.0)
     615        ep(i, k) = min(ep(i, k), 1.0)
    625616        sigp(i, k) = sigs
    626617      END IF
    … …  
    635626  DO k = minorig + 1, nl
    636627    DO i = 1, ncum
    637       IF (k>=(icb(i)+1)) THEN
    638         tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
     628      IF (k>=(icb(i) + 1)) THEN
     629        tvp(i, k) = tvp(i, k) * (1.0 - qnk(i) + ep(i, k) * clw(i, k))
    639630        ! PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
    640631        ! PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
    … …  
    643634  END DO
    644635  DO i = 1, ncum
    645     tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
     636    tvp(i, nlp) = tvp(i, nl) - (gz(i, nlp) - gz(i, nl)) / cpd
    646637  END DO
    647638
    … …  
    721712    DO i = 1, ncum
    722713      IF (cape(i)<0.0) lcape(i) = .FALSE.
    723       IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
    724         by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
    725         byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
     714      IF ((k>=(icb(i) + 1)) .AND. lcape(i)) THEN
     715        by = (tvp(i, k) - tv(i, k)) * dph(i, k) / p(i, k)
     716        byp(i) = (tvp(i, k + 1) - tv(i, k + 1)) * dph(i, k + 1) / p(i, k + 1)
    726717        cape(i) = cape(i) + by
    727718        IF (by>=0.0) inb1(i) = k + 1
    … …  
    737728    defrac = capem(i) - cape(i)
    738729    defrac = max(defrac, 0.001)
    739     frac(i) = -cape(i)/defrac
     730    frac(i) = -cape(i) / defrac
    740731    frac(i) = min(frac(i), 1.0)
    741732    frac(i) = max(frac(i), 0.0)
    … …  
    747738
    748739  ! initialization:
    749   DO i = 1, ncum*nlp
     740  DO i = 1, ncum * nlp
    750741    hp(i, 1) = h(i, 1)
    751742  END DO
    … …  
    754745    DO i = 1, ncum
    755746      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
    756         hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
    757           )
    758       END IF
    759     END DO
    760   END DO
    761 
     747        hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k &
     748                )
     749      END IF
     750    END DO
     751  END DO
    762752
    763753END SUBROUTINE cv_undilute2
    764754
    765755SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
    766     cpn, iflag, cbmf)
     756        cpn, iflag, cbmf)
     757  USE lmdz_cvthermo
     758
    767759  IMPLICIT NONE
    768760
    … …  
    771763  INTEGER nk(nloc), icb(nloc)
    772764  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
    773   REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
     765  REAL ph(nloc, nd + 1) ! caution nd instead ndp1 to be consistent...
    774766  REAL plcl(nloc), cpn(nloc, nd)
    775767
    … …  
    783775  REAL work(nloc)
    784776
    785   include "cvthermo.h"
    786777  include "cvparam.h"
    787778
    … …  
    805796  DO i = 1, ncum
    806797    dtpbl(i) = 0.0
    807     tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
    808       i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
    809     tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
    810       ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
     798    tvpplcl(i) = tvp(i, icb(i) - 1) - rrd * tvp(i, icb(i) - 1) * (p(i, icb(i) - 1) - plcl(&
     799            i)) / (cpn(i, icb(i) - 1) * p(i, icb(i) - 1))
     800    tvaplcl(i) = tv(i, icb(i)) + (tvp(i, icb(i)) - tvp(i, icb(i) + 1)) * (plcl(i) - p(i &
     801            , icb(i))) / (p(i, icb(i)) - p(i, icb(i) + 1))
    811802  END DO
    812803
    … …  
    820811  DO k = minorig, icbmax
    821812    DO i = 1, ncum
    822       IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
    823         dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
    824       END IF
    825     END DO
    826   END DO
    827   DO i = 1, ncum
    828     dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
     813      IF ((k>=nk(i)) .AND. (k<=(icb(i) - 1))) THEN
     814        dtpbl(i) = dtpbl(i) + (tvp(i, k) - tv(i, k)) * dph(i, k)
     815      END IF
     816    END DO
     817  END DO
     818  DO i = 1, ncum
     819    dtpbl(i) = dtpbl(i) / (ph(i, nk(i)) - ph(i, icb(i)))
    829820    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
    830821  END DO
    … …  
    836827  DO i = 1, ncum
    837828    work(i) = cbmf(i)
    838     cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
     829    cbmf(i) = max(0.0, (1.0 - damp) * cbmf(i) + 0.1 * alpha * dtmin(i))
    839830    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
    840831      iflag(i) = 3
    … …  
    842833  END DO
    843834
    844 
    845835END SUBROUTINE cv_closure
    846836
    847837SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
    848     h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
    849     sij, elij)
     838        h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
     839        sij, elij)
     840  USE lmdz_cvthermo
     841
    850842  IMPLICIT NONE
    851843
    852   include "cvthermo.h"
    853844  include "cvparam.h"
    854845
    … …  
    857848  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
    858849  REAL cbmf(nloc), qnk(nloc)
    859   REAL ph(nloc, nd+1)
     850  REAL ph(nloc, nd + 1)
    860851  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
    861852  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
    … …  
    881872  ! =====================================================================
    882873
    883   DO i = 1, ncum*nlp
     874  DO i = 1, ncum * nlp
    884875    nent(i, 1) = 0
    885876    m(i, 1) = 0.0
    … …  
    907898  DO j = minorig + 1, nl
    908899    DO i = 1, ncum
    909       IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
     900      IF ((j>=(icb(i) + 1)) .AND. (j<=inb(i))) THEN
    910901        k = min(j, inb1(i))
    911         dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
    912           entp*0.04*(ph(i,k)-ph(i,k+1))
     902        dbo = abs(tv(i, k + 1) - tvp(i, k + 1) - tv(i, k - 1) + tvp(i, k - 1)) + &
     903                entp * 0.04 * (ph(i, k) - ph(i, k + 1))
    913904        work(i) = work(i) + dbo
    914         m(i, j) = cbmf(i)*dbo
     905        m(i, j) = cbmf(i) * dbo
    915906      END IF
    916907    END DO
    … …  
    918909  DO k = minorig + 1, nl
    919910    DO i = 1, ncum
    920       IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
    921         m(i, k) = m(i, k)/work(i)
     911      IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN
     912        m(i, k) = m(i, k) / work(i)
    922913      END IF
    923914    END DO
    … …  
    931922  ! =====================================================================
    932923
    933 
    934924  DO i = minorig + 1, nl
    935925    DO j = minorig + 1, nl
    936926      DO ij = 1, ncum
    937         IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
    938             inb(ij))) THEN
    939           qti = qnk(ij) - ep(ij, i)*clw(ij, i)
    940           bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
    941           anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
    942           denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
     927        IF ((i>=(icb(ij) + 1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
     928                inb(ij))) THEN
     929          qti = qnk(ij) - ep(ij, i) * clw(ij, i)
     930          bf2 = 1. + lv(ij, j) * lv(ij, j) * qs(ij, j) / (rrv * t(ij, j) * t(ij, j) * cpd)
     931          anum = h(ij, j) - hp(ij, i) + (cpv - cpd) * t(ij, j) * (qti - q(ij, j))
     932          denom = h(ij, i) - hp(ij, i) + (cpd - cpv) * (q(ij, i) - qti) * t(ij, j)
    943933          dei = denom
    944934          IF (abs(dei)<0.01) dei = 0.01
    945           sij(ij, i, j) = anum/dei
     935          sij(ij, i, j) = anum / dei
    946936          sij(ij, i, i) = 1.0
    947           altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
    948           altem = altem/bf2
    949           cwat = clw(ij, j)*(1.-ep(ij,j))
     937          altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
     938          altem = altem / bf2
     939          cwat = clw(ij, j) * (1. - ep(ij, j))
    950940          stemp = sij(ij, i, j)
    951941          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
    952             anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
    953             denom = denom + lv(ij, j)*(q(ij,i)-qti)
     942            anum = anum - lv(ij, j) * (qti - qs(ij, j) - cwat * bf2)
     943            denom = denom + lv(ij, j) * (q(ij, i) - qti)
    954944            IF (abs(denom)<0.01) denom = 0.01
    955             sij(ij, i, j) = anum/denom
    956             altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
    957             altem = altem - (bf2-1.)*cwat
     945            sij(ij, i, j) = anum / denom
     946            altem = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti - qs(ij, j)
     947            altem = altem - (bf2 - 1.) * cwat
    958948          END IF
    959           IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
    960             qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
    961             uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
    962               (1.-sij(ij,i,j))*u(ij, nk(ij))
    963             vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
    964               (1.-sij(ij,i,j))*v(ij, nk(ij))
     949          IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
     950            qent(ij, i, j) = sij(ij, i, j) * q(ij, i) + (1. - sij(ij, i, j)) * qti
     951            uent(ij, i, j) = sij(ij, i, j) * u(ij, i) + &
     952                    (1. - sij(ij, i, j)) * u(ij, nk(ij))
     953            vent(ij, i, j) = sij(ij, i, j) * v(ij, i) + &
     954                    (1. - sij(ij, i, j)) * v(ij, nk(ij))
    965955            elij(ij, i, j) = altem
    966             elij(ij, i, j) = max(0.0, elij(ij,i,j))
    967             ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
     956            elij(ij, i, j) = max(0.0, elij(ij, i, j))
     957            ment(ij, i, j) = m(ij, i) / (1. - sij(ij, i, j))
    968958            nent(ij, i) = nent(ij, i) + 1
    969959          END IF
    970           sij(ij, i, j) = max(0.0, sij(ij,i,j))
    971           sij(ij, i, j) = min(1.0, sij(ij,i,j))
     960          sij(ij, i, j) = max(0.0, sij(ij, i, j))
     961          sij(ij, i, j) = min(1.0, sij(ij, i, j))
    972962        END IF
    973963      END DO
    … …  
    980970
    981971    DO ij = 1, ncum
    982       IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
     972      IF ((i>=(icb(ij) + 1)) .AND. (i<=inb(ij)) .AND. (nent(ij, i)==0)) THEN
    983973        ment(ij, i, i) = m(ij, i)
    984         qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
     974        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
    985975        uent(ij, i, i) = u(ij, nk(ij))
    986976        vent(ij, i, i) = v(ij, nk(ij))
    … …  
    1000990  ! =====================================================================
    1001991
    1002   CALL zilch(bsum, ncum*nlp)
     992  CALL zilch(bsum, ncum * nlp)
    1003993  DO ij = 1, ncum
    1004994    lwork(ij) = .FALSE.
    … …  
    1008998    num1 = 0
    1009999    DO ij = 1, ncum
    1010       IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
     1000      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) num1 = num1 + 1
    10111001    END DO
    10121002    IF (num1<=0) GO TO 789
    10131003
    10141004    DO ij = 1, ncum
    1015       IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
    1016         lwork(ij) = (nent(ij,i)/=0)
    1017         qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
    1018         anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
    1019         denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
     1005      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij))) THEN
     1006        lwork(ij) = (nent(ij, i)/=0)
     1007        qp1 = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
     1008        anum = h(ij, i) - hp(ij, i) - lv(ij, i) * (qp1 - qs(ij, i))
     1009        denom = h(ij, i) - hp(ij, i) + lv(ij, i) * (q(ij, i) - qp1)
    10201010        IF (abs(denom)<0.01) denom = 0.01
    1021         scrit(ij) = anum/denom
    1022         alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
     1011        scrit(ij) = anum / denom
     1012        alt = qp1 - qs(ij, i) + scrit(ij) * (q(ij, i) - qp1)
    10231013        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
    10241014        asij(ij) = 0.0
    … …  
    10301020      num2 = 0
    10311021      DO ij = 1, ncum
    1032         IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
    1033           ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
     1022        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
     1023                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
    10341024      END DO
    10351025      IF (num2<=0) GO TO 783
    10361026
    10371027      DO ij = 1, ncum
    1038         IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
    1039             ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
    1040           IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
     1028        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
     1029                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
     1030          IF (sij(ij, i, j)>0.0 .AND. sij(ij, i, j)<0.9) THEN
    10411031            IF (j>i) THEN
    1042               smid = min(sij(ij,i,j), scrit(ij))
     1032              smid = min(sij(ij, i, j), scrit(ij))
    10431033              sjmax = smid
    10441034              sjmin = smid
    1045               IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
     1035              IF (smid<smin(ij) .AND. sij(ij, i, j + 1)<smid) THEN
    10461036                smin(ij) = smid
    1047                 sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
    1048                 sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
     1037                sjmax = min(sij(ij, i, j + 1), sij(ij, i, j), scrit(ij))
     1038                sjmin = max(sij(ij, i, j - 1), sij(ij, i, j))
    10491039                sjmin = min(sjmin, scrit(ij))
    10501040              END IF
    10511041            ELSE
    1052               sjmax = max(sij(ij,i,j+1), scrit(ij))
    1053               smid = max(sij(ij,i,j), scrit(ij))
     1042              sjmax = max(sij(ij, i, j + 1), scrit(ij))
     1043              smid = max(sij(ij, i, j), scrit(ij))
    10541044              sjmin = 0.0
    1055               IF (j>1) sjmin = sij(ij, i, j-1)
     1045              IF (j>1) sjmin = sij(ij, i, j - 1)
    10561046              sjmin = max(sjmin, scrit(ij))
    10571047            END IF
    1058             delp = abs(sjmax-smid)
    1059             delm = abs(sjmin-smid)
    1060             asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
    1061             ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
     1048            delp = abs(sjmax - smid)
     1049            delm = abs(sjmin - smid)
     1050            asij(ij) = asij(ij) + (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
     1051            ment(ij, i, j) = ment(ij, i, j) * (delp + delm) * (ph(ij, j) - ph(ij, j + 1))
    10621052          END IF
    10631053        END IF
    10641054      END DO
    1065 783 END DO
     1055    783 END DO
    10661056    DO ij = 1, ncum
    1067       IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
     1057      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
    10681058        asij(ij) = max(1.0E-21, asij(ij))
    1069         asij(ij) = 1.0/asij(ij)
     1059        asij(ij) = 1.0 / asij(ij)
    10701060        bsum(ij, i) = 0.0
    10711061      END IF
    … …  
    10731063    DO j = minorig, nl + 1
    10741064      DO ij = 1, ncum
    1075         IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
    1076             ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
    1077           ment(ij, i, j) = ment(ij, i, j)*asij(ij)
     1065        IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (j>=icb(&
     1066                ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
     1067          ment(ij, i, j) = ment(ij, i, j) * asij(ij)
    10781068          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
    10791069        END IF
    … …  
    10811071    END DO
    10821072    DO ij = 1, ncum
    1083       IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
    1084           i)<1.0E-18) .AND. lwork(ij)) THEN
     1073      IF ((i>=icb(ij) + 1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
     1074              i)<1.0E-18) .AND. lwork(ij)) THEN
    10851075        nent(ij, i) = 0
    10861076        ment(ij, i, i) = m(ij, i)
    1087         qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
     1077        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i) * clw(ij, i)
    10881078        uent(ij, i, i) = u(ij, nk(ij))
    10891079        vent(ij, i, i) = v(ij, nk(ij))
    … …  
    10921082      END IF
    10931083    END DO
    1094 789 END DO
    1095 
     1084  789 END DO
    10961085
    10971086END SUBROUTINE cv_mixing
    10981087
    10991088SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
    1100     ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
     1089        ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
     1090  USE lmdz_cvthermo
     1091
    11011092  IMPLICIT NONE
    11021093
    1103 
    1104   include "cvthermo.h"
    11051094  include "cvparam.h"
    11061095
    … …  
    11101099  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
    11111100  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
    1112   REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
     1101  REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
    11131102  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
    11141103  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
    … …  
    11501139  DO k = 2, nl + 1
    11511140    DO i = 1, ncum
    1152       qp(i, k) = q(i, k-1)
    1153       up(i, k) = u(i, k-1)
    1154       vp(i, k) = v(i, k-1)
     1141      qp(i, k) = q(i, k - 1)
     1142      up(i, k) = u(i, k - 1)
     1143      vp(i, k) = v(i, k - 1)
    11551144    END DO
    11561145  END DO
    … …  
    11641153  ! ***                and condensed water flux                    ***
    11651154
    1166 
    11671155  DO i = 1, ncum
    11681156    jtt(i) = 2
    1169     IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
     1157    IF (ep(i, inb(i))<=0.0001) iflag(i) = 2
    11701158    IF (iflag(i)==0) THEN
    11711159      lwork(i) = .TRUE.
    … …  
    11771165  ! ***                    Begin downdraft loop                    ***
    11781166
    1179 
    11801167  CALL zilch(wdtrain, ncum)
    11811168  DO i = nl + 1, 1, -1
    … …  
    11921179    DO ij = 1, ncum
    11931180      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
    1194         wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
     1181        wdtrain(ij) = g * ep(ij, i) * m(ij, i) * clw(ij, i)
    11951182      END IF
    11961183    END DO
    … …  
    12001187        DO ij = 1, ncum
    12011188          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
    1202             awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
     1189            awat = elij(ij, j, i) - (1. - ep(ij, i)) * clw(ij, i)
    12031190            awat = max(0.0, awat)
    1204             wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
     1191            wdtrain(ij) = wdtrain(ij) + g * awat * ment(ij, j, i)
    12051192          END IF
    12061193        END DO
    … …  
    12231210        ! rain   ***
    12241211
    1225         IF (t(ij,i)>273.0) THEN
     1212        IF (t(ij, i)>273.0) THEN
    12261213          coeff = coeffr
    12271214          wt(ij, i) = omtrain
    12281215        END IF
    1229         qsm = 0.5*(q(ij,i)+qp(ij,i+1))
    1230         afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
     1216        qsm = 0.5 * (q(ij, i) + qp(ij, i + 1))
     1217        afac = coeff * ph(ij, i) * (qs(ij, i) - qsm) / (1.0E4 + 2.0E3 * ph(ij, i) * qs(ij, i))
    12311218        afac = max(afac, 0.0)
    12321219        sigt = sigp(ij, i)
    12331220        sigt = max(0.0, sigt)
    12341221        sigt = min(1.0, sigt)
    1235         b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
    1236         c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
    1237         revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
    1238         evap(ij, i) = sigt*afac*revap
    1239         water(ij, i) = revap*revap
     1222        b6 = 100. * (ph(ij, i) - ph(ij, i + 1)) * sigt * afac / wt(ij, i)
     1223        c6 = (water(ij, i + 1) * wt(ij, i + 1) + wdtrain(ij) / sigd) / wt(ij, i)
     1224        revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6))
     1225        evap(ij, i) = sigt * afac * revap
     1226        water(ij, i) = revap * revap
    12401227
    12411228        ! ***  Calculate precipitating downdraft mass flux under     ***
    … …  
    12431230
    12441231        IF (i>1) THEN
    1245           dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
     1232          dhdp = (h(ij, i) - h(ij, i - 1)) / (p(ij, i - 1) - p(ij, i))
    12461233          dhdp = max(dhdp, 10.0)
    1247           mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
    1248           mp(ij, i) = max(mp(ij,i), 0.0)
     1234          mp(ij, i) = 100. * ginv * lv(ij, i) * sigd * evap(ij, i) / dhdp
     1235          mp(ij, i) = max(mp(ij, i), 0.0)
    12491236
    12501237          ! ***   Add small amount of inertia to downdraft              ***
    12511238
    1252           fac = 20.0/(ph(ij,i-1)-ph(ij,i))
    1253           mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
     1239          fac = 20.0 / (ph(ij, i - 1) - ph(ij, i))
     1240          mp(ij, i) = (fac * mp(ij, i + 1) + mp(ij, i)) / (1. + fac)
    12541241
    12551242          ! ***      Force mp to decrease linearly to zero
    … …  
    12581245          ! ***
    12591246
    1260           IF (p(ij,i)>(0.949*p(ij,1))) THEN
     1247          IF (p(ij, i)>(0.949 * p(ij, 1))) THEN
    12611248            jtt(ij) = max(jtt(ij), i)
    1262             mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
    1263               (p(ij,1)-p(ij,jtt(ij)))
     1249            mp(ij, i) = mp(ij, jtt(ij)) * (p(ij, 1) - p(ij, i)) / &
     1250                    (p(ij, 1) - p(ij, jtt(ij)))
    12641251          END IF
    12651252        END IF
    … …  
    12711258            qstm = qs(ij, 1)
    12721259          ELSE
    1273             qstm = qs(ij, i-1)
     1260            qstm = qs(ij, i - 1)
    12741261          END IF
    1275           IF (mp(ij,i)>mp(ij,i+1)) THEN
    1276             rat = mp(ij, i+1)/mp(ij, i)
    1277             qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
    1278               100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
    1279             up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
    1280             vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
     1262          IF (mp(ij, i)>mp(ij, i + 1)) THEN
     1263            rat = mp(ij, i + 1) / mp(ij, i)
     1264            qp(ij, i) = qp(ij, i + 1) * rat + q(ij, i) * (1.0 - rat) + &
     1265                    100. * ginv * sigd * (ph(ij, i) - ph(ij, i + 1)) * (evap(ij, i) / mp(ij, i))
     1266            up(ij, i) = up(ij, i + 1) * rat + u(ij, i) * (1. - rat)
     1267            vp(ij, i) = vp(ij, i + 1) * rat + v(ij, i) * (1. - rat)
    12811268          ELSE
    1282             IF (mp(ij,i+1)>0.0) THEN
    1283               qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
    1284                 i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
    1285                 i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
    1286               up(ij, i) = up(ij, i+1)
    1287               vp(ij, i) = vp(ij, i+1)
     1269            IF (mp(ij, i + 1)>0.0) THEN
     1270              qp(ij, i) = (gz(ij, i + 1) - gz(ij, i) + qp(ij, i + 1) * (lv(ij, i + 1) + t(ij, &
     1271                      i + 1) * (cl - cpd)) + cpd * (t(ij, i + 1) - t(ij, &
     1272                      i))) / (lv(ij, i) + t(ij, i) * (cl - cpd))
     1273              up(ij, i) = up(ij, i + 1)
     1274              vp(ij, i) = vp(ij, i + 1)
    12881275            END IF
    12891276          END IF
    1290           qp(ij, i) = min(qp(ij,i), qstm)
    1291           qp(ij, i) = max(qp(ij,i), 0.0)
     1277          qp(ij, i) = min(qp(ij, i), qstm)
     1278          qp(ij, i) = max(qp(ij, i), 0.0)
    12921279        END IF
    12931280      END IF
    12941281    END DO
    1295 899 END DO
    1296 
     1282  899 END DO
    12971283
    12981284END SUBROUTINE cv_unsat
    12991285
    13001286SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
    1301     ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
    1302     ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
    1303     precip, cbmf, ft, fq, fu, fv, ma, qcondc)
     1287        ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
     1288        ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
     1289        precip, cbmf, ft, fq, fu, fv, ma, qcondc)
     1290  USE lmdz_cvthermo
     1291
    13041292  IMPLICIT NONE
    13051293
    1306   include "cvthermo.h"
    13071294  include "cvparam.h"
    13081295
    … …  
    13141301  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
    13151302  REAL gz(nloc, nd)
    1316   REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
     1303  REAL p(nloc, nd), ph(nloc, nd + 1), h(nloc, nd)
    13171304  REAL hp(nloc, nd), lv(nloc, nd)
    13181305  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
    … …  
    13441331  ! -- initializations:
    13451332
    1346   delti = 1.0/delt
     1333  delti = 1.0 / delt
    13471334
    13481335  DO i = 1, ncum
    … …  
    13561343      fv(i, k) = 0.0
    13571344      fq(i, k) = 0.0
    1358       lvcp(i, k) = lv(i, k)/cpn(i, k)
     1345      lvcp(i, k) = lv(i, k) / cpn(i, k)
    13591346      qcondc(i, k) = 0.0 ! cld
    13601347      qcond(i, k) = 0.0 ! cld
    … …  
    13711358      ! c     &                /(rowl*g)
    13721359      ! c            precip(i)=precip(i)*delt/86400.
    1373       precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
     1360      precip(i) = wt(i, 1) * sigd * water(i, 1) * 86400 / g
    13741361    END IF
    13751362  END DO
    … …  
    13801367
    13811368  DO i = 1, ncum
    1382     wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
    1383     qprime(i) = 0.5*(qp(i,1)-q(i,1))
    1384     tprime(i) = lv0*qprime(i)/cpd
     1369    wd(i) = betad * abs(mp(i, icb(i))) * 0.01 * rrd * t(i, icb(i)) / (sigd * p(i, icb(i)))
     1370    qprime(i) = 0.5 * (qp(i, 1) - q(i, 1))
     1371    tprime(i) = lv0 * qprime(i) / cpd
    13851372  END DO
    13861373
    … …  
    13891376
    13901377  DO i = 1, ncum
    1391     work(i) = 0.01/(ph(i,1)-ph(i,2))
     1378    work(i) = 0.01 / (ph(i, 1) - ph(i, 2))
    13921379    am(i) = 0.0
    13931380  END DO
    … …  
    14001387  END DO
    14011388  DO i = 1, ncum
    1402     IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
    1403     ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
    1404       1))/cpn(i,1))
    1405     ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
    1406     ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
    1407       work(i)/cpn(i, 1)
    1408     fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
    1409       sigd*evap(i, 1)
    1410     fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
    1411     fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
    1412       2)-u(i,1)))
    1413     fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
    1414       2)-v(i,1)))
     1389    IF ((g * work(i) * am(i))>=delti) iflag(i) = 1
     1390    ft(i, 1) = ft(i, 1) + g * work(i) * am(i) * (t(i, 2) - t(i, 1) + (gz(i, 2) - gz(i, &
     1391            1)) / cpn(i, 1))
     1392    ft(i, 1) = ft(i, 1) - lvcp(i, 1) * sigd * evap(i, 1)
     1393    ft(i, 1) = ft(i, 1) + sigd * wt(i, 2) * (cl - cpd) * water(i, 2) * (t(i, 2) - t(i, 1)) * &
     1394            work(i) / cpn(i, 1)
     1395    fq(i, 1) = fq(i, 1) + g * mp(i, 2) * (qp(i, 2) - q(i, 1)) * work(i) + &
     1396            sigd * evap(i, 1)
     1397    fq(i, 1) = fq(i, 1) + g * am(i) * (q(i, 2) - q(i, 1)) * work(i)
     1398    fu(i, 1) = fu(i, 1) + g * work(i) * (mp(i, 2) * (up(i, 2) - u(i, 1)) + am(i) * (u(i, &
     1399            2) - u(i, 1)))
     1400    fv(i, 1) = fv(i, 1) + g * work(i) * (mp(i, 2) * (vp(i, 2) - v(i, 1)) + am(i) * (v(i, &
     1401            2) - v(i, 1)))
    14151402  END DO
    14161403  DO j = 2, nl
    14171404    DO i = 1, ncum
    14181405      IF (j<=inb(i)) THEN
    1419         fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
    1420         fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
    1421         fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
     1406        fq(i, 1) = fq(i, 1) + g * work(i) * ment(i, j, 1) * (qent(i, j, 1) - q(i, 1))
     1407        fu(i, 1) = fu(i, 1) + g * work(i) * ment(i, j, 1) * (uent(i, j, 1) - u(i, 1))
     1408        fv(i, 1) = fv(i, 1) + g * work(i) * ment(i, j, 1) * (vent(i, j, 1) - v(i, 1))
    14221409      END IF
    14231410    END DO
    … …  
    14431430    DO k = i + 1, nl + 1
    14441431      DO ij = 1, ncum
    1445         IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
     1432        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij) + 1))) THEN
    14461433          amp1(ij) = amp1(ij) + m(ij, k)
    14471434        END IF
    … …  
    14521439      DO j = i + 1, nl + 1
    14531440        DO ij = 1, ncum
    1454           IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
     1441          IF ((j<=(inb(ij) + 1)) .AND. (i<=inb(ij))) THEN
    14551442            amp1(ij) = amp1(ij) + ment(ij, k, j)
    14561443          END IF
    … …  
    14701457    DO ij = 1, ncum
    14711458      IF (i<=inb(ij)) THEN
    1472         dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
    1473         cpinv = 1.0/cpn(ij, i)
    1474 
    1475         ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
    1476           i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
    1477           i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
    1478         ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
    1479           ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
    1480         ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
    1481           ij,i+1)-t(ij,i))*dpinv*cpinv
    1482         fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
    1483           i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
    1484         fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
    1485           i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
    1486         fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
    1487           i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
     1459        dpinv = 0.01 / (ph(ij, i) - ph(ij, i + 1))
     1460        cpinv = 1.0 / cpn(ij, i)
     1461
     1462        ft(ij, i) = ft(ij, i) + g * dpinv * (amp1(ij) * (t(ij, i + 1) - t(ij, &
     1463                i) + (gz(ij, i + 1) - gz(ij, i)) * cpinv) - ad(ij) * (t(ij, i) - t(ij, &
     1464                i - 1) + (gz(ij, i) - gz(ij, i - 1)) * cpinv)) - sigd * lvcp(ij, i) * evap(ij, i)
     1465        ft(ij, i) = ft(ij, i) + g * dpinv * ment(ij, i, i) * (hp(ij, i) - h(ij, i) + t(ij &
     1466                , i) * (cpv - cpd) * (q(ij, i) - qent(ij, i, i))) * cpinv
     1467        ft(ij, i) = ft(ij, i) + sigd * wt(ij, i + 1) * (cl - cpd) * water(ij, i + 1) * (t(&
     1468                ij, i + 1) - t(ij, i)) * dpinv * cpinv
     1469        fq(ij, i) = fq(ij, i) + g * dpinv * (amp1(ij) * (q(ij, i + 1) - q(ij, &
     1470                i)) - ad(ij) * (q(ij, i) - q(ij, i - 1)))
     1471        fu(ij, i) = fu(ij, i) + g * dpinv * (amp1(ij) * (u(ij, i + 1) - u(ij, &
     1472                i)) - ad(ij) * (u(ij, i) - u(ij, i - 1)))
     1473        fv(ij, i) = fv(ij, i) + g * dpinv * (amp1(ij) * (v(ij, i + 1) - v(ij, &
     1474                i)) - ad(ij) * (v(ij, i) - v(ij, i - 1)))
    14881475      END IF
    14891476    END DO
    … …  
    14911478      DO ij = 1, ncum
    14921479        IF (i<=inb(ij)) THEN
    1493           awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
     1480          awat = elij(ij, k, i) - (1. - ep(ij, i)) * clw(ij, i)
    14941481          awat = max(awat, 0.0)
    1495           fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
    1496             (ij,i))
    1497           fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
    1498             ))
    1499           fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
    1500             ))
     1482          fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - awat - q &
     1483                  (ij, i))
     1484          fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
     1485                  ))
     1486          fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
     1487                  ))
    15011488          ! (saturated updrafts resulting from mixing)               ! cld
    1502           qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
     1489          qcond(ij, i) = qcond(ij, i) + (elij(ij, k, i) - awat) ! cld
    15031490          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
    15041491        END IF
    … …  
    15081495      DO ij = 1, ncum
    15091496        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
    1510           fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
    1511             ))
    1512           fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
    1513             ))
    1514           fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
    1515             ))
     1497          fq(ij, i) = fq(ij, i) + g * dpinv * ment(ij, k, i) * (qent(ij, k, i) - q(ij, i &
     1498                  ))
     1499          fu(ij, i) = fu(ij, i) + g * dpinv * ment(ij, k, i) * (uent(ij, k, i) - u(ij, i &
     1500                  ))
     1501          fv(ij, i) = fv(ij, i) + g * dpinv * ment(ij, k, i) * (vent(ij, k, i) - v(ij, i &
     1502                  ))
    15161503        END IF
    15171504      END DO
    … …  
    15191506    DO ij = 1, ncum
    15201507      IF (i<=inb(ij)) THEN
    1521         fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
    1522           i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
    1523         fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
    1524           i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
    1525         fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
    1526           i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
     1508        fq(ij, i) = fq(ij, i) + sigd * evap(ij, i) + g * (mp(ij, i + 1) * (qp(ij, &
     1509                i + 1) - q(ij, i)) - mp(ij, i) * (qp(ij, i) - q(ij, i - 1))) * dpinv
     1510        fu(ij, i) = fu(ij, i) + g * (mp(ij, i + 1) * (up(ij, i + 1) - u(ij, &
     1511                i)) - mp(ij, i) * (up(ij, i) - u(ij, i - 1))) * dpinv
     1512        fv(ij, i) = fv(ij, i) + g * (mp(ij, i + 1) * (vp(ij, i + 1) - v(ij, &
     1513                i)) - mp(ij, i) * (vp(ij, i) - v(ij, i - 1))) * dpinv
    15271514        ! (saturated downdrafts resulting from mixing)               ! cld
    15281515        DO k = i + 1, inb(ij) ! cld
    … …  
    15311518        END DO ! cld
    15321519        ! (particular case: no detraining level is found)            ! cld
    1533         IF (nent(ij,i)==0) THEN ! cld
    1534           qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
     1520        IF (nent(ij, i)==0) THEN ! cld
     1521          qcond(ij, i) = qcond(ij, i) + (1. - ep(ij, i)) * clw(ij, i) ! cld
    15351522          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
    15361523        END IF ! cld
    1537         IF (nqcond(ij,i)/=0.) THEN ! cld
    1538           qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
     1524        IF (nqcond(ij, i)/=0.) THEN ! cld
     1525          qcond(ij, i) = qcond(ij, i) / nqcond(ij, i) ! cld
    15391526        END IF ! cld
    15401527      END IF
    15411528    END DO
    1542 1500 END DO
     1529  1500 END DO
    15431530
    15441531  ! *** Adjust tendencies at top of convection layer to reflect  ***
    … …  
    15471534  DO ij = 1, ncum
    15481535    fqold = fq(ij, inb(ij))
    1549     fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
    1550     fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
    1551       inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
    1552       inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
     1536    fq(ij, inb(ij)) = fq(ij, inb(ij)) * (1. - frac(ij))
     1537    fq(ij, inb(ij) - 1) = fq(ij, inb(ij) - 1) + frac(ij) * fqold * ((ph(ij, &
     1538            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
     1539            inb(ij)))) * lv(ij, inb(ij)) / lv(ij, inb(ij) - 1)
    15531540    ftold = ft(ij, inb(ij))
    1554     ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
    1555     ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
    1556       inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
    1557       inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
     1541    ft(ij, inb(ij)) = ft(ij, inb(ij)) * (1. - frac(ij))
     1542    ft(ij, inb(ij) - 1) = ft(ij, inb(ij) - 1) + frac(ij) * ftold * ((ph(ij, &
     1543            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, &
     1544            inb(ij)))) * cpn(ij, inb(ij)) / cpn(ij, inb(ij) - 1)
    15581545    fuold = fu(ij, inb(ij))
    1559     fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
    1560     fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
    1561       inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
     1546    fu(ij, inb(ij)) = fu(ij, inb(ij)) * (1. - frac(ij))
     1547    fu(ij, inb(ij) - 1) = fu(ij, inb(ij) - 1) + frac(ij) * fuold * ((ph(ij, &
     1548            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
    15621549    fvold = fv(ij, inb(ij))
    1563     fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
    1564     fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
    1565       inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
     1550    fv(ij, inb(ij)) = fv(ij, inb(ij)) * (1. - frac(ij))
     1551    fv(ij, inb(ij) - 1) = fv(ij, inb(ij) - 1) + frac(ij) * fvold * ((ph(ij, &
     1552            inb(ij)) - ph(ij, inb(ij) + 1)) / (ph(ij, inb(ij) - 1) - ph(ij, inb(ij))))
    15661553  END DO
    15671554
    … …  
    15741561    vav(ij) = 0.0
    15751562    DO i = 1, inb(ij)
    1576       ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
    1577         ph(ij,i+1))
    1578       uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
    1579       vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
     1563      ents(ij) = ents(ij) + (cpn(ij, i) * ft(ij, i) + lv(ij, i) * fq(ij, i)) * (ph(ij, i) - &
     1564              ph(ij, i + 1))
     1565      uav(ij) = uav(ij) + fu(ij, i) * (ph(ij, i) - ph(ij, i + 1))
     1566      vav(ij) = vav(ij) + fv(ij, i) * (ph(ij, i) - ph(ij, i + 1))
    15801567    END DO
    15811568  END DO
    15821569  DO ij = 1, ncum
    1583     ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
    1584     uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
    1585     vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
     1570    ents(ij) = ents(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
     1571    uav(ij) = uav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
     1572    vav(ij) = vav(ij) / (ph(ij, 1) - ph(ij, inb(ij) + 1))
    15861573  END DO
    15871574  DO ij = 1, ncum
    15881575    DO i = 1, inb(ij)
    1589       ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
    1590       fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
    1591       fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
     1576      ft(ij, i) = ft(ij, i) - ents(ij) / cpn(ij, i)
     1577      fu(ij, i) = (1. - cu) * (fu(ij, i) - uav(ij))
     1578      fv(ij, i) = (1. - cu) * (fv(ij, i) - vav(ij))
    15921579    END DO
    15931580  END DO
    … …  
    15951582  DO k = 1, nl + 1
    15961583    DO i = 1, ncum
    1597       IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
    1598     END DO
    1599   END DO
    1600 
     1584      IF ((q(i, k) + delt * fq(i, k))<0.0) iflag(i) = 10
     1585    END DO
     1586  END DO
    16011587
    16021588  DO i = 1, ncum
    … …  
    16251611  DO k = nl, 1, -1
    16261612    DO i = 1, ncum
    1627       ma(i, k) = ma(i, k+1) + m(i, k)
     1613      ma(i, k) = ma(i, k + 1) + m(i, k)
    16281614    END DO
    16291615  END DO
    … …  
    16471633      ax(ij, i) = 0. ! cld
    16481634      DO j = icb(ij), i ! cld
    1649         ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
    1650           *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
     1635        ax(ij, i) = ax(ij, i) + rrd * (tvp(ij, j) - tv(ij, j)) & ! cld
     1636                * (ph(ij, j) - ph(ij, j + 1)) / p(ij, j) ! cld
    16511637      END DO ! cld
    1652       IF (ax(ij,i)>0.0) THEN ! cld
    1653         wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
     1638      IF (ax(ij, i)>0.0) THEN ! cld
     1639        wa(ij, i) = sqrt(2. * ax(ij, i)) ! cld
    16541640      END IF ! cld
    16551641    END DO ! cld
    16561642    DO i = 1, nl ! cld
    1657       IF (wa(ij,i)>0.0) &          ! cld
    1658         siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
    1659         *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
    1660       siga(ij, i) = min(siga(ij,i), 1.0) ! cld
    1661       qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
    1662         +(1.-siga(ij,i))*qcond(ij, i) ! cld
     1643      IF (wa(ij, i)>0.0) &          ! cld
     1644              siga(ij, i) = mac(ij, i) / wa(ij, i) & ! cld
     1645                      * rrd * tvp(ij, i) / p(ij, i) / 100. / delta ! cld
     1646      siga(ij, i) = min(siga(ij, i), 1.0) ! cld
     1647      qcondc(ij, i) = siga(ij, i) * clw(ij, i) * (1. - ep(ij, i)) & ! cld
     1648              + (1. - siga(ij, i)) * qcond(ij, i) ! cld
    16631649    END DO ! cld
    16641650  END DO ! cld
    16651651
    1666 
    16671652END SUBROUTINE cv_yield
    16681653
    16691654SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
    1670     fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
    1671     qcondc1)
     1655        fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
     1656        qcondc1)
    16721657  IMPLICIT NONE
    16731658
    … …  
    17101695  END DO
    17111696
    1712 
    17131697END SUBROUTINE cv_uncompress
    17141698
Note: See TracChangeset for help on using the changeset viewer.