Changeset 1494 for LMDZ5/trunk


Ignore:
Timestamp:
Mar 9, 2011, 11:05:02 AM (14 years ago)
Author:
Laurent Fairhead
Message:

Optimization of routines from the new physics

Optimisation de routines de la nouvelle physique

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

Legend:

Unmodified
Added
Removed

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

    r1479 r1494  
    3838      CHARACTER (LEN=20) :: modname='cv3_param'
    3939      CHARACTER (LEN=80) :: abort_message
    40       LOGICAL, SAVE :: FIRST=.TRUE.
    41 !$OMP THREADPRIVATE(FIRST)
    42 
    4340
    4441c noff: integer limit for convection (nd-noff)
    … …  
    103100      CLOSE(99)
    1041019999  Continue
    105       if (first) then
    106         WRITE(*,*)'dpbase=',dpbase
    107         WRITE(*,*)'pbcrit=',pbcrit
    108         WRITE(*,*)'ptcrit=',ptcrit
    109         WRITE(*,*)'sigdz=',sigdz
    110         WRITE(*,*)'spfac=',spfac
    111       endif
     102      WRITE(*,*)'dpbase=',dpbase
     103      WRITE(*,*)'pbcrit=',pbcrit
     104      WRITE(*,*)'ptcrit=',ptcrit
     105      WRITE(*,*)'sigdz=',sigdz
     106      WRITE(*,*)'spfac=',spfac
    112107
    113108      return
    … …  
    10551050           pden=ptcrit-pbcrit
    10561051           ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax
    1057            ep(i,k)=amax1(ep(i,k),0.0)
     1052           ep(i,k)=max(ep(i,k),0.0)
    10581053           ep(i,k)=amin1(ep(i,k),epmax)
    10591054           sigp(i,k)=spfac
    … …  
    11571152      enddo
    11581153
    1159       do 530 i=1,ncum
    1160        do 535 k=1,nl-1
     1154      do 535 k=1,nl-1
     1155       do 530 i=1,ncum
    11611156        if ((k.ge.iposit(i)).and.(buoy(i,k).lt.dtovsh)) then
    11621157         inb(i)=MIN(inb(i),k)
    11631158        endif
    1164  535   continue
    1165  530  continue
    1166 
     1159 530   continue
     1160 535  continue
     1161c
    11671162c-- end convect3
    11681163
    … …  
    13541349        if (k.le.icb(i))then
    13551350         sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i))
    1356          sig(i,k)=amax1(sig(i,k),0.0)
     1351         sig(i,k)=max(sig(i,k),0.0)
    13571352         w0(i,k)=beta*w0(i,k)
    13581353        endif
    … …  
    13641359c!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
    13651360c!     1              abs(buoy(inb))
    1366 c!            sig(i)=amax1(sig(i),0.0)
     1361c!            sig(i)=max(sig(i),0.0)
    13671362c!            w0(i)=beta*w0(i)
    13681363c!   85    continue
    … …  
    13711366c!      do 87 i=1,icb
    13721367c!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
    1373 c!         sig(i)=amax1(sig(i),0.0)
     1368c!         sig(i)=max(sig(i),0.0)
    13741369c!         w0(i)=beta*w0(i)
    13751370c!   87 continue
    … …  
    14361431
    14371432         sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*ABS(dtmin(i,k))
    1438          sig(i,k)=amax1(sig(i,k),0.0)
     1433         sig(i,k)=max(sig(i,k),0.0)
    14391434         sig(i,k)=amin1(sig(i,k),0.01)
    14401435         fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
    … …  
    14961491c!         dcape=rrd*buoy(i-1)*deltap/p(i-1)
    14971492c!         dlnp=deltap/p(i-1)
    1498 c!         cape=amax1(0.0,cape)
     1493c!         cape=max(0.0,cape)
    14991494c!         sigold=sig(i)
    15001495
    … …  
    15051500
    15061501c!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
    1507 c!         sig(i)=amax1(sig(i),0.0)
     1502c!         sig(i)=max(sig(i),0.0)
    15081503c!         sig(i)=amin1(sig(i),0.01)
    15091504c!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
    … …  
    16541649c!!!      end do
    16551650          elij(il,i,j)=altem
    1656           elij(il,i,j)=amax1(0.0,elij(il,i,j))
     1651          elij(il,i,j)=max(0.0,elij(il,i,j))
    16571652          ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
    16581653          nent(il,i)=nent(il,i)+1
    16591654         end if
    1660          sij(il,i,j)=amax1(0.0,sij(il,i,j))
     1655         sij(il,i,j)=max(0.0,sij(il,i,j))
    16611656         sij(il,i,j)=amin1(1.0,sij(il,i,j))
    16621657         endif ! new
    … …  
    17781773        wgh=1.0
    17791774        if(j.gt.i)then
    1780          sjmax=amax1(sij(il,i,j+1),smax(il))
     1775         sjmax=max(sij(il,i,j+1),smax(il))
    17811776         sjmax=amin1(sjmax,scrit(il))
    1782          smax(il)=amax1(sij(il,i,j),smax(il))
    1783          sjmin=amax1(sij(il,i,j-1),smax(il))
     1777         smax(il)=max(sij(il,i,j),smax(il))
     1778         sjmin=max(sij(il,i,j-1),smax(il))
    17841779         sjmin=amin1(sjmin,scrit(il))
    17851780         if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
    17861781         smid=amin1(sij(il,i,j),scrit(il))
    17871782        else
    1788          sjmax=amax1(sij(il,i,j+1),scrit(il))
    1789          smid=amax1(sij(il,i,j),scrit(il))
     1783         sjmax=max(sij(il,i,j+1),scrit(il))
     1784         smid=max(sij(il,i,j),scrit(il))
    17901785         sjmin=0.0
    17911786         if(j.gt.1)sjmin=sij(il,i,j-1)
    1792          sjmin=amax1(sjmin,scrit(il))
     1787         sjmin=max(sjmin,scrit(il))
    17931788        endif
    17941789        delp=abs(sjmax-smid)
    … …  
    18041799      do il=1,ncum
    18051800       if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
    1806         asij(il)=amax1(1.0e-16,asij(il))
     1801        asij(il)=max(1.0e-16,asij(il))
    18071802        asij(il)=1.0/asij(il)
    18081803        asum(il,i)=0.0
    … …  
    18341829      do il=1,ncum
    18351830       if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
    1836         bsum(il,i)=amax1(bsum(il,i),1.0e-16)
     1831        bsum(il,i)=max(bsum(il,i),1.0e-16)
    18371832        bsum(il,i)=1.0/bsum(il,i)
    18381833       endif
    … …  
    19511946      real tinv, delti
    19521947      real awat, afac, afac1, afac2, bfac
    1953       real pr1, pr2, sigt, b6, c6, revap, tevap, delth
     1948      real pr1, pr2, sigt, b6, c6, revap, delth
    19541949      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
    19551950      real ampmax
     1951      real tevap(nloc)
    19561952      real lvcp(nloc,na)
    19571953      real h(nloc,na),hm(nloc,na)
    19581954      real wdtrain(nloc)
    1959       logical lwork(nloc)
     1955      logical lwork(nloc),mplus(nloc)
    19601956
    19611957
    … …  
    19931989
    19941990        do il=1,ncum
    1995           lwork(il)=.TRUE.
    1996           if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
     1991!!          lwork(il)=.TRUE.
     1992!!          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
     1993          lwork(il)= ep(il,inb(il)) .ge. 0.0001
    19971994        enddo
    1998 
    1999         call zilch(wdtrain,ncum)
    20001995
    20011996c   ***  Set the fractionnal area sigd of precipitating downdraughts
    … …  
    20041999        enddo
    20052000
     2001
     2002c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2003c
     2004c    ***                    begin downdraft loop                    ***
     2005c
     2006c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2007c
    20062008        DO 400 i=nl+1,1,-1
    20072009
    … …  
    20122014        if (num1.le.0) goto 400
    20132015
     2016        call zilch(wdtrain,ncum)
     2017
    20142018c
    20152019c   ***  integrate liquid water equation to find condensed water   ***
    20162020c   ***                and condensed water flux                    ***
    20172021c
    2018 
    2019 c
    2020 c    ***                    begin downdraft loop                    ***
    2021 c
    2022 
    20232022c
    20242023c    ***              calculate detrained precipitation             ***
    … …  
    20392038          if (i.le.inb(il) .and. lwork(il)) then
    20402039           awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
    2041            awat=amax1(awat,0.0)
     2040           awat=max(awat,0.0)
    20422041           if (cvflag_grav) then
    20432042            wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)
    … …  
    20552054c
    20562055
    2057       do 999 il=1,ncum
    2058 
     2056      do 995 il=1,ncum
    20592057       if (i.le.inb(il) .and. lwork(il)) then
    20602058
    … …  
    20662064       rp(il,i)=0.5*(rp(il,i)+rr(il,i))
    20672065      endif
    2068       rp(il,i)=amax1(rp(il,i),0.0)
     2066      rp(il,i)=max(rp(il,i),0.0)
    20692067      rp(il,i)=amin1(rp(il,i),rs(il,i))
    20702068      rp(il,inb(il))=rr(il,inb(il))
    … …  
    20772075       rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
    20782076       rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
    2079        rp(il,i-1)=amax1(rp(il,i-1),0.0)
     2077       rp(il,i-1)=max(rp(il,i-1),0.0)
    20802078       afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i))
    20812079       afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
    … …  
    20842082      endif
    20852083      if(i.eq.inb(il))afac=0.0
    2086       afac=amax1(afac,0.0)
     2084      afac=max(afac,0.0)
    20872085      bfac=1./(sigd(il)*wt(il,i))
    20882086c
    … …  
    21482146     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    21492147c
     2148       endif !(i.le.inb(il) .and. lwork(il))
     2149995   Continue
     2150c----------------------------------------------------------------
     2151c
    21502152ccc
    21512153c    ***  calculate precipitating downdraft mass flux under     ***
    21522154c    ***              hydrostatic approximation                 ***
    21532155c
    2154       if (i.ne.1) then
    2155 
    2156       tevap=amax1(0.0,evap(il,i))
    2157       delth=amax1(0.001,(th(il,i)-th(il,i-1)))
     2156      Do 996 il = 1,ncum
     2157       if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
     2158c
     2159      tevap(il)=max(0.0,evap(il,i))
     2160      delth=max(0.001,(th(il,i)-th(il,i-1)))
    21582161      if (cvflag_grav) then
    2159        mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap
     2162       mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il)
    21602163     :              *(p(il,i-1)-p(il,i))/delth
    21612164      else
    2162        mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap
     2165       mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap(il)
    21632166     :         *(p(il,i-1)-p(il,i))/delth
    21642167      endif
     2168c
     2169       endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
     2170996   Continue
     2171c----------------------------------------------------------------
    21652172c
    21662173c    ***           if hydrostatic assumption fails,             ***
    … …  
    21682175c    ***  and mass flux from two simultaneous differential eqns ***
    21692176c
     2177      Do 997 il = 1,ncum
     2178       if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
     2179c
    21702180      amfac=sigd(il)*sigd(il)*70.0*ph(il,i)*(p(il,i-1)-p(il,i))
    21712181     :          *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
    21722182      amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
     2183c
    21732184      if(amp2.gt.(0.1*amfac))then
    21742185       xf=100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
    … …  
    21772188       af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
    21782189       bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
    2179      :            +50.*(p(il,i-1)-p(il,i))*xf*tevap
     2190     :            +50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
    21802191       fac2=1.0
    21812192       if(bf.lt.0.0)fac2=-1.0
    … …  
    21932204        mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv)
    21942205       endif
    2195        mp(il,i)=amax1(0.0,mp(il,i))
     2206       mp(il,i)=max(0.0,mp(il,i))
    21962207
    21972208       if (cvflag_grav) then
    … …  
    21992210C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
    22002211C Et il faut bien revoir les facteurs 100.
    2201         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
     2212        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
    22022213     2   /(mp(il,i)+sigd(il)*0.1)
    22032214     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
    22042215     : *sigd(il)*th(il,i))
    22052216       else
    2206         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
     2217        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
    22072218     2   /(mp(il,i)+sigd(il)*0.1)
    22082219     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
    22092220     : *sigd(il)*th(il,i))
    22102221       endif
    2211        b(il,i-1)=amax1(b(il,i-1),0.0)
    2212       endif
     2222       b(il,i-1)=max(b(il,i-1),0.0)
     2223c
     2224      endif !(amp2.gt.(0.1*amfac))
    22132225c
    22142226c   ***         limit magnitude of mp(i) to meet cfl condition      ***
    … …  
    22162228      ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti
    22172229      amp2=2.0*(ph(il,i-1)-ph(il,i))*delti
    2218       ampmax=amin1(ampmax,amp2)
    2219       mp(il,i)=amin1(mp(il,i),ampmax)
     2230      ampmax=min(ampmax,amp2)
     2231      mp(il,i)=min(mp(il,i),ampmax)
    22202232c
    22212233c    ***      force mp to decrease linearly to zero                 ***
    … …  
    22312243      endif
    22322244
    2233 360   continue
    2234       endif ! i.eq.1
     2245       endif ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
     2246997   Continue
     2247c----------------------------------------------------------------
    22352248c
    22362249c    ***       find mixing ratio of precipitating downdraft     ***
    22372250c
    2238 
    2239       if (i.ne.inb(il)) then
    2240 
     2251      Do il = 1,ncum
     2252       if (i.lt.inb(il) .and. lwork(il)) then
     2253        mplus(il) = mp(il,i).gt.mp(il,i+1)
     2254       endif ! (i.lt.inb(il) .and. lwork(il))
     2255      enddo
     2256c
     2257      Do 999 il = 1,ncum
     2258       if (i.lt.inb(il) .and. lwork(il)) then
     2259c
    22412260      rp(il,i)=rr(il,i)
    22422261
    2243       if(mp(il,i).gt.mp(il,i+1))then
     2262      if(mplus(il))then
    22442263
    22452264       if (cvflag_grav) then
    … …  
    22582277      vp(il,i)=vp(il,i)/mp(il,i)
    22592278
    2260       do j=1,ntra
    2261       trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
    2262      :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
    2263       trap(il,i,j)=trap(il,i,j)/mp(il,i)
    2264       end do
    2265 
    2266       else
     2279      else ! if (mplus(il))
    22672280
    22682281       if(mp(il,i+1).gt.1.0e-16)then
    … …  
    22782291       up(il,i)=up(il,i+1)
    22792292       vp(il,i)=vp(il,i+1)
    2280 
    2281        do j=1,ntra
    2282        trap(il,i,j)=trap(il,i+1,j)
    2283        end do
    2284 
    2285        endif
    2286       endif
     2293       endif ! (mp(il,i+1).gt.1.0e-16)
     2294      endif ! (mplus(il)) else if (.not.mplus(il))
     2295c
    22872296      rp(il,i)=amin1(rp(il,i),rs(il,i))
    2288       rp(il,i)=amax1(rp(il,i),0.0)
    2289 
    2290       endif
    2291       endif
     2297      rp(il,i)=max(rp(il,i),0.0)
     2298
     2299      endif ! (i.lt.inb(il) .and. lwork(il))
    22922300999   continue
     2301c----------------------------------------------------------------
     2302c
     2303c    ***       find tracer concentrations in precipitating downdraft     ***
     2304c
     2305      do j=1,ntra
     2306       do il = 1,ncum
     2307       if (i.lt.inb(il) .and. lwork(il)) then
     2308c
     2309         if(mplus(il))then
     2310          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2311     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
     2312          trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2313         else ! if (mplus(il))
     2314          if(mp(il,i+1).gt.1.0e-16)then
     2315           trap(il,i,j)=trap(il,i+1,j)
     2316          endif
     2317         endif ! (mplus(il)) else if (.not.mplus(il))
     2318c
     2319        endif ! (i.lt.inb(il) .and. lwork(il))
     2320       enddo
     2321      end do
    22932322
    22942323400   continue
     2324c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2325c
     2326c    ***                    end of downdraft loop                    ***
     2327c
     2328c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2329c
    22952330
    22962331       return
    … …  
    28142849       do il = 1,ncum
    28152850        awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
    2816         awat(il)=amax1(awat(il),0.0)
     2851        awat(il)=max(awat(il),0.0)
    28172852       enddo
    28182853c

  • LMDZ5/trunk/libf/phylmd/cv3p_mixing.F

    r1050 r1494  
    5353      real Qmixmin(nloc), Rmixmin(nloc), SQmRmin(nloc)
    5454      real signhpmh(nloc)
    55       real Sx, Scrit2
    56       integer Jx
     55      real Sx(nloc), Scrit2
    5756      real smid(nloc), sjmin(nloc), sjmax(nloc)
    5857      real Sbef(nloc), Sup(nloc), Smin(nloc)
    … …  
    240239      enddo
    241240
    242       DO 789 i=minorig+1,nl
     241c---------------------------------------------------------------
     242      DO 789 i=minorig+1,nl     !Loop on origin level "i"
     243c---------------------------------------------------------------
    243244
    244245      num1=0
    … …  
    248249      if (num1.le.0) goto 789
    249250
     251c
     252cjyg1    Find maximum of SIJ for J>I, if any.
     253c
     254       Sx(:) = 0.
     255c
     256       DO il = 1,ncum
     257        IF ( i.ge.icb(il) .and. i.le.inb(il) ) THEN
     258         Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
     259         Sbef(il) = max(0.,signhpmh(il))
     260        ENDIF
     261       ENDDO
     262
     263       DO j = i+1,nl
     264        DO il = 1,ncum
     265         IF ( i.ge.icb(il) .and. i.le.inb(il)
     266     :         .and. j.le.inb(il) ) THEN
     267          IF (Sbef(il) .LT. Sij(il,i,j)) THEN
     268            Sx(il) = max(Sij(il,i,j),Sx(il))
     269          ENDIF
     270          Sbef(il) = Sij(il,i,j)
     271         ENDIF
     272        ENDDO
     273       ENDDO
     274c
    250275
    251276      do 781 il=1,ncum
    252277       if ( i.ge.icb(il) .and. i.le.inb(il) ) then
    253278        lwork(il)=(nent(il,i).ne.0)
    254         Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
    255279        qp=qnk(il)-ep(il,i)*clw(il,i)
    256280        anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
    … …  
    262286        alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
    263287c
    264 cjyg1    Find maximum of SIJ for J>I, if any, and new critical value Scrit2
     288cjyg1    Find new critical value Scrit2
    265289c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
    266290c                    Sij < Scrit2  => mixed draught will detrain at J>I
    267291c
    268        Sx = 0.
    269        Jx = 0.
    270        Sbef(il) = max(0.,signhpmh(il))
    271        DO j = i+1,inb(il)
    272          IF (Sbef(il) .LT. Sij(il,i,j)) THEN
    273            Sx = max(Sij(il,i,j),Sx)
    274            Jx = J
    275          ENDIF
    276          Sbef(il) = Sij(il,i,j)
    277        ENDDO
    278 c
    279        Scrit2 = min(Scrit(il),Sx)*max(0.,-signhpmh(il))
     292       Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il))
    280293     :         +Scrit(il)*max(0.,signhpmh(il))
    281294c
    … …  
    293306781   continue
    294307
    295       do 175 j=minorig,nl
     308c---------------------------------------------------------------
     309      do 175 j=minorig,nl   !Loop on destination level "j"
     310c---------------------------------------------------------------
    296311
    297312      num2=0

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

    r1403 r1494  
    8080!$OMP THREADPRIVATE(lev_out)
    8181
    82       INTEGER ig,k,l,ll
     82      INTEGER ig,k,l,ll,ierr
    8383      real zsortie1d(klon)
    8484      INTEGER lmax(klon),lmin(klon),lalim(klon)
    … …  
    233233!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
    234234     do ig=1,klon
    235       if (prt_level.ge.20) then
    236        print*,'th_main ig f0',ig,f0(ig)
    237       endif
    238235         f0(ig)=max(f0(ig),1.e-2)
    239236         zmax0(ig)=max(zmax0(ig),40.)
    … …  
    241238     enddo
    242239
     240      if (prt_level.ge.20) then
     241       do ig=1,ngrid
     242          print*,'th_main ig f0',ig,f0(ig)
     243       enddo
     244      endif
    243245!-----------------------------------------------------------------------
    244246! Calcul de T,q,ql a partir de Tl et qT dans l environnement
    … …  
    290292!-----------------------------------------------------------------------
    291293
    292       do l=1,nlay
    293          rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
    294       enddo
    295 
    296 !IM
     294     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
     295
    297296     if (prt_level.ge.10)write(lunout,*)                                &
    298297    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
    … …  
    619618      enddo
    620619!IM
     620      ierr=0
    621621      do ig=1,ngrid
    622622        if (pcon(ig).le.pplay(ig,nlay)) then
    623623           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
     624           ierr=1
     625        endif
     626      enddo
     627      if (ierr==1) then
    624628           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
    625629           CALL abort_gcm (modname,abort_message,1)
    626         endif
    627       enddo
     630      endif
     631
    628632      if (prt_level.ge.1) print*,'14b OK convect8'
    629633      do k=nlay,1,-1
    … …  
    655659            zf2=zf/(1.-zf)
    656660!
    657       if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
    658 !
    659       if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
    660661            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
    661662            if(zw2(ig,l).gt.1.e-10) then
    … …  
    664665             wth2(ig,l)=0.
    665666            endif
    666 !           print*,'wth2=',wth2(ig,l)
    667667            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
    668668     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
    669       if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
    670669            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    671670!test: on calcul q2/po=ratqsc
    … …  
    682681      enddo
    683682!
     683      if (prt_level.ge.10) then
     684         ig=igout
     685         do l=1,nlay
     686            print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
     687            print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
     688         enddo
     689      endif
     690
    684691!      print*,'avant calcul ale et alp'
    685692!calcul de ALE et ALP pour la convection
    … …  
    705712!initialisations
    706713!      print*,'ponderation'
    707       do ig=1,ngrid
    708            fm_tot(ig)=0.
    709       enddo
    710        do ig=1,ngrid
    711         do k=1,klev
    712            wght_th(ig,k)=1.
    713         enddo
    714        enddo
    715        do ig=1,ngrid
    716 !         lalim_conv(ig)=lmix_bis(ig)
    717 !la hauteur de la couche alim_conv = hauteur couche alim_therm
    718          lalim_conv(ig)=lalim(ig)
    719 !         zentr(ig)=zlev(ig,lalim(ig))
    720       enddo
    721       do ig=1,ngrid
    722         do k=1,lalim_conv(ig)
    723            fm_tot(ig)=fm_tot(ig)+fm(ig,k)
    724         enddo
    725       enddo
    726       do ig=1,ngrid
    727         do k=1,lalim_conv(ig)
    728            if (fm_tot(ig).gt.1.e-10) then
    729 !           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
    730            endif
    731 !on pondere chaque couche par a*
    732              if (alim_star(ig,k).gt.1.e-10) then
    733              wght_th(ig,k)=alim_star(ig,k)
    734              else
    735              wght_th(ig,k)=1.
    736              endif
    737         enddo
    738       enddo
     714
     715      fm_tot(:)=0.
     716      wght_th(:,:)=1.
     717      lalim_conv(:)=lalim(:)
     718
     719      do k=1,klev
     720         do ig=1,ngrid
     721            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
     722         enddo
     723      enddo
     724
     725! assez bizarre car, si on est dans la couche d'alim et que alim_star et
     726! plus petit que 1.e-10, on prend wght_th=1.
     727      do k=1,klev
     728         do ig=1,ngrid
     729            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
     730               wght_th(ig,k)=alim_star(ig,k)
     731            endif
     732         enddo
     733      enddo
     734
    739735!      print*,'apres wght_th'
    740736!test pour prolonger la convection
    … …  
    748744      enddo
    749745
     746
    750747!calcul du ratqscdiff
    751748      if (prt_level.ge.1) print*,'14e OK convect8'
    … …  
    753750      vardiff=0.
    754751      ratqsdiff(:,:)=0.
    755       do ig=1,ngrid
    756          do l=1,lalim(ig)
     752
     753      do l=1,klev
     754         do ig=1,ngrid
     755            if (l<=lalim(ig)) then
    757756            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
    758          enddo
    759       enddo
     757            endif
     758         enddo
     759      enddo
     760
    760761      if (prt_level.ge.1) print*,'14f OK convect8'
    761       do ig=1,ngrid
    762           do l=1,lalim(ig)
    763           zf=fraca(ig,l)
    764           zf2=zf/(1.-zf)
    765           vardiff=vardiff+alim_star(ig,l)  &
    766      &           *(zqta(ig,l)*1000.-var)**2
    767 !         ratqsdiff=ratqsdiff+alim_star(ig,l)*
    768 !     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    769           enddo
    770       enddo
     762
     763      do l=1,klev
     764         do ig=1,ngrid
     765            if (l<=lalim(ig)) then
     766               zf=fraca(ig,l)
     767               zf2=zf/(1.-zf)
     768               vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
     769            endif
     770         enddo
     771      enddo
     772
    771773      if (prt_level.ge.1) print*,'14g OK convect8'
    772774      do l=1,nlay
    … …  
    779781!
    780782!ecriture des fichiers sortie
    781 !     print*,'15 OK convect8'
     783!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
    782784
    783785#ifdef wrgrads_thermcell
Note: See TracChangeset for help on using the changeset viewer.