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

Optimization of routines from the new physics

Optimisation de routines de la nouvelle physique

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.