Ignore:
Timestamp:
Mar 4, 2004, 4:11:16 PM (21 years ago)
Author:
lmdzadmin
Message:

Optimisation de differentes routines, IM, MAF, FH
LF

Location:
LMDZ.3.3/branches/rel-LF/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F

    r418 r495  
    4848      REAL xx(klon), aux(klon), coeff(klon), block(klon)
    4949      REAL  dist(klon), fprime(klon), det(klon)
    50       REAL pi, u(klon), v(klon), erfu(klon), erfv(klon)
     50      REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon)
    5151      REAL  xx1(klon), xx2(klon)
    5252      real erf,kkk
     53      real sqrtpi,sqrt2,zx1,zx2,exdel
    5354c lconv = true si le calcul a converge (entre autre si qsub < min_q)
    5455       LOGICAL lconv(klon)
    … …  
    5657
    5758      pi = ACOS(-1.)
     59      sqrtpi=sqrt(pi)
     60      sqrt2=sqrt(2.)
    5861
    5962      ptconv=.false.
    … …  
    123126          xx(i) = -0.0001
    124127        else
    125          xx1(i)=-SQRT(2.)*vmax(i)*(1.0-SQRT(1.0+delta(i)/(vmax(i)**2.)))
    126          xx2(i)=-SQRT(2.)*vmax(i)*(1.0+SQRT(1.0+delta(i)/(vmax(i)**2.)))
     128         zx1=-sqrt2*vmax(i)
     129         zx2=SQRT(1.0+delta(i)/(vmax(i)**2.))
     130         xx1(i)=zx1*(1.0-zx2)
     131         xx2(i)=zx1*(1.0+zx2)
    127132         xx(i) = 1.01 * xx1(i)
    128133         if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i)
    … …  
    143148        if (.not.lconv(i)) then
    144149
    145           u(i) = delta(i)/(xx(i)*sqrt(2.)) + xx(i)/(2.*sqrt(2.))
    146           v(i) = delta(i)/(xx(i)*sqrt(2.)) - xx(i)/(2.*sqrt(2.))
     150          u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
     151          v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)
    147152
    148153          IF ( v(i) .GT. vmax(i) ) THEN
    … …  
    153158c -- use asymptotic expression of erf for u and v large:
    154159c ( -> analytic solution for xx )
    155 
    156              aux(i) = 2.0*delta(i)*(1.-beta(i)*EXP(delta(i)))
    157      :                       /(1.+beta(i)*EXP(delta(i)))
     160             exdel=beta(i)*EXP(delta(i))
     161             aux(i) = 2.0*delta(i)*(1.-exdel)
     162     :                       /(1.+exdel)
    158163             if (aux(i).lt.0.) then
    159                 print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)
     164c                print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)
    160165                aux(i)=0.
    161166             endif
    162167             xx(i) = -SQRT(aux(i))
    163              block(i) = EXP(-v(i)*v(i)) / v(i) / SQRT(pi)
     168             block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi
    164169             dist(i) = 0.0
    165170             fprime(i) = 1.0
    … …  
    169174c -- erfv -> 1.0, use an asymptotic expression of erfv for v large:
    170175
    171              erfu(i) = ERF(u(i))
     176             erfcu(i) = 1.0-ERF(u(i))
    172177c  !!! ATTENTION : rajout d'un seuil pour l'exponentiel
    173              aux(i) = SQRT(pi)*(1.0-erfu(i))*EXP(min(v(i)*v(i),100.))
     178             aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i),100.))
    174179             coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.)
    175              block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / SQRT(pi)
     180             block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi
    176181             dist(i) = v(i) * aux(i) / coeff(i) - beta(i)
    177182             fprime(i) = 2.0 / xx(i) * (v(i)**2.)
    … …  
    185190c -- general case:
    186191
    187            erfu(i) = ERF(u(i))
    188            erfv(i) = ERF(v(i))
    189            block(i) = 1.0-erfv(i)
    190            dist(i) = (1.0 - erfu(i)) / (1.0 - erfv(i)) - beta(i)
     192           erfcu(i) = 1.0-ERF(u(i))
     193           erfcv(i) = 1.0-ERF(v(i))
     194           block(i) = erfcv(i)
     195           dist(i) = erfcu(i) / erfcv(i) - beta(i)
    191196           zu2(i)=u(i)*u(i)
    192197           zv2(i)=v(i)*v(i)
    193198           if(zu2(i).gt.20..or. zv2(i).gt.20.) then
    194               print*,'ATTENTION !!! xx(',i,') =', xx(i)
    195            print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
    196      .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
    197      .CLDF(i,k)
    198               print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i)
     199c              print*,'ATTENTION !!! xx(',i,') =', xx(i)
     200c           print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
     201c     .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
     202c     .CLDF(i,k)
     203c              print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i)
    199204              zu2(i)=20.
    200205              zv2(i)=20.
    201206             fprime(i) = 0.
    202207           else
    203              fprime(i) = 2. /SQRT(pi) /xx(i) /(1.0-erfv(i))**2.
    204      :           * (   (1.0-erfv(i))*v(i)*EXP(-zu2(i))
    205      :               - (1.0-erfu(i))*u(i)*EXP(-zv2(i)) )
     208             fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2.
     209     :           * (   erfcv(i)*v(i)*EXP(-zu2(i))
     210     :               - erfcu(i)*u(i)*EXP(-zv2(i)) )
    206211           endif
    207212          ENDIF ! x
    … …  
    219224            ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.)
    220225            CLDF(i,K) = 0.5 * block(i)
    221 cccccccccccccccccccccccc
    222 c           kkk=-sqrt(log(1.+ratqsc(i,k)**2))
    223 c           u(i)=delta(i)/(kkk*sqrt(2.))-kkk/(2.*sqrt(2.))
    224 c           v(i)=delta(i)/(kkk*sqrt(2.))+kkk/(2.*sqrt(2.))
    225 c           erfu(i)=erf(u(i))
    226 c           erfv(i)=erf(v(i))
    227 c           print*,'SIG ',k,qsub(i,k)
    228 c    s      ,mu(i)*((1.-erfv(i))/(1.-erfu(i))-qsat(i)/mu(i))
    229 c    s      ,0.5*erfu(i)
    230 cccccccccccccccccccccccc
    231226          else
    232227            xx(i) = xx(i) - dist(i)/fprime(i)

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F

    r486 r495  
    804804 110  continue
    805805
    806       do 121 j=1,ntra
    807 ccccc      do 111 k=1,nl+1
    808       do 111 k=1,nd
    809        nn=0
    810       do 101 i=1,len
    811       if(iflag1(i).eq.0)then
    812        nn=nn+1
    813        tra(nn,k,j)=tra1(i,k,j)
    814       endif
    815  101  continue
    816  111  continue
    817  121  continue
     806c      do 121 j=1,ntra
     807c      do 111 k=1,nd
     808c       nn=0
     809c      do 101 i=1,len
     810c      if(iflag1(i).eq.0)then
     811c       nn=nn+1
     812c       tra(nn,k,j)=tra1(i,k,j)
     813c      endif
     814c 101  continue
     815c 111  continue
     816c 121  continue
    818817
    819818      if (nn.ne.ncum) then
    … …  
    14931492 400  continue
    14941493
    1495       do k=1,ntra
    1496        do j=1,nd  ! instead nlp
    1497         do i=1,nd ! instead nlp
    1498          do il=1,ncum
    1499             traent(il,i,j,k)=tra(il,j,k)
    1500          enddo
    1501         enddo
    1502        enddo
    1503       enddo
     1494c      do k=1,ntra
     1495c       do j=1,nd  ! instead nlp
     1496c        do i=1,nd ! instead nlp
     1497c         do il=1,ncum
     1498c            traent(il,i,j,k)=tra(il,j,k)
     1499c         enddo
     1500c        enddo
     1501c       enddo
     1502c      enddo
    15041503      zm(:,:)=0.
    15051504
    … …  
    15571556 710  continue
    15581557
    1559        do k=1,ntra
    1560         do j=minorig,nl
    1561          do il=1,ncum
    1562           if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
    1563      :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
    1564             traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    1565      :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
    1566           endif
    1567          enddo
    1568         enddo
    1569        enddo
     1558c       do k=1,ntra
     1559c        do j=minorig,nl
     1560c         do il=1,ncum
     1561c          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
     1562c     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
     1563c            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1564c     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1565c          endif
     1566c         enddo
     1567c        enddo
     1568c       enddo
    15701569
    15711570c
    … …  
    15901589 750  continue
    15911590 
    1592       do j=1,ntra
    1593        do i=minorig+1,nl
    1594         do il=1,ncum
    1595          if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
    1596           traent(il,i,i,j)=tra(il,nk(il),j)
    1597          endif
    1598         enddo
    1599        enddo
    1600       enddo
     1591c      do j=1,ntra
     1592c       do i=minorig+1,nl
     1593c        do il=1,ncum
     1594c         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
     1595c          traent(il,i,i,j)=tra(il,nk(il),j)
     1596c         endif
     1597c        enddo
     1598c       enddo
     1599c      enddo
    16011600
    16021601      do 100 j=minorig,nl
    … …  
    17641763      enddo ! il
    17651764
    1766       do j=1,ntra
    1767        do il=1,ncum
    1768         if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
    1769      :     .and. csum(il,i).lt.m(il,i) ) then
    1770          traent(il,i,i,j)=tra(il,nk(il),j)
    1771         endif
    1772        enddo
    1773       enddo
    1774 
     1765c      do j=1,ntra
     1766c       do il=1,ncum
     1767c        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
     1768c     :     .and. csum(il,i).lt.m(il,i) ) then
     1769c         traent(il,i,i,j)=tra(il,nk(il),j)
     1770c        endif
     1771c       enddo
     1772c      enddo
    17751773789   continue
    17761774c     
    … …  
    18691867        enddo
    18701868
    1871         do k=1,ntra
    1872          do i=1,nd
    1873           do il=1,ncum
    1874            trap(il,i,k)=tra(il,i,k)
    1875           enddo
    1876          enddo
    1877         enddo
     1869c        do k=1,ntra
     1870c         do i=1,nd
     1871c          do il=1,ncum
     1872c           trap(il,i,k)=tra(il,i,k)
     1873c          enddo
     1874c         enddo
     1875c        enddo
    18781876
    18791877c
    … …  
    21032101      vp(il,i)=vp(il,i)/mp(il,i)
    21042102
    2105       do j=1,ntra
    2106       trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2103c      do j=1,ntra
     2104c      trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
    21072105ctestmaf     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
    2108      :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
    2109       trap(il,i,j)=trap(il,i,j)/mp(il,i)
    2110       end do
     2106c     :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
     2107c      trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2108c      end do
    21112109
    21122110      else
    … …  
    21252123       vp(il,i)=vp(il,i+1)
    21262124
    2127        do j=1,ntra
    2128        trap(il,i,j)=trap(il,i+1,j)
    2129        end do
     2125c       do j=1,ntra
     2126c       trap(il,i,j)=trap(il,i+1,j)
     2127c       end do
    21302128
    21312129       endif
    … …  
    22262224      enddo
    22272225
    2228       do j=1,ntra
    2229        do i=1,nd
    2230         do il=1,ncum
    2231           ftra(il,i,j)=0.0
    2232         enddo
    2233        enddo
    2234       enddo
     2226c      do j=1,ntra
     2227c       do i=1,nd
     2228c        do il=1,ncum
     2229c          ftra(il,i,j)=0.0
     2230c        enddo
     2231c       enddo
     2232c      enddo
    22352233
    22362234      do i=1,nl
    … …  
    23302328      enddo ! il
    23312329
    2332       do j=1,ntra
    2333        do il=1,ncum
    2334         if (cvflag_grav) then
    2335          ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
    2336      :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2337      :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2338         else
    2339          ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
    2340      :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2341      :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2342         endif
    2343        enddo
    2344       enddo
     2330c      do j=1,ntra
     2331c       do il=1,ncum
     2332c        if (cvflag_grav) then
     2333c         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
     2334c     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2335c     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2336c        else
     2337c         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
     2338c     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2339c     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2340c        endif
     2341c       enddo
     2342c      enddo
    23452343
    23462344      do j=2,nl
    … …  
    23662364      enddo
    23672365
    2368       do k=1,ntra
    2369        do j=2,nl
    2370         do il=1,ncum
    2371          if (j.le.inb(il)) then
    2372 
    2373           if (cvflag_grav) then
    2374            ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
    2375      :                *(traent(il,j,1,k)-tra(il,1,k))
    2376           else
    2377            ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
    2378      :                *(traent(il,j,1,k)-tra(il,1,k))
    2379           endif
    2380 
    2381          endif
    2382         enddo
    2383        enddo
    2384       enddo
     2366c      do k=1,ntra
     2367c       do j=2,nl
     2368c        do il=1,ncum
     2369c         if (j.le.inb(il)) then
     2370
     2371c          if (cvflag_grav) then
     2372c           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
     2373c     :                *(traent(il,j,1,k)-tra(il,1,k))
     2374c          else
     2375c           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
     2376c     :                *(traent(il,j,1,k)-tra(il,1,k))
     2377c          endif
     2378
     2379c         endif
     2380c        enddo
     2381c       enddo
     2382c      enddo
    23852383
    23862384c
    … …  
    248824861350  continue
    24892487
    2490       do k=1,ntra
    2491        do il=1,ncum
    2492         if (i.le.inb(il)) then
    2493          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2494          cpinv=1.0/cpn(il,i)
    2495          if (cvflag_grav) then
    2496            ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
    2497      :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2498      :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2499          else
    2500            ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
    2501      :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2502      :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2503          endif
    2504         endif
    2505        enddo
    2506       enddo
     2488c      do k=1,ntra
     2489c       do il=1,ncum
     2490c        if (i.le.inb(il)) then
     2491c         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2492c         cpinv=1.0/cpn(il,i)
     2493c         if (cvflag_grav) then
     2494c           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
     2495c     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2496c     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2497c         else
     2498c           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
     2499c     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2500c     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2501c         endif
     2502c        endif
     2503c       enddo
     2504c      enddo
    25072505
    25082506      do 480 k=1,i-1
    … …  
    25382536480   continue
    25392537
    2540       do j=1,ntra
    2541        do k=1,i-1
    2542         do il=1,ncum
    2543          if (i.le.inb(il)) then
    2544           dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2545           cpinv=1.0/cpn(il,i)
    2546           if (cvflag_grav) then
    2547            ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    2548      :        *(traent(il,k,i,j)-tra(il,i,j))
    2549           else
    2550            ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    2551      :        *(traent(il,k,i,j)-tra(il,i,j))
    2552           endif
    2553          endif
    2554         enddo
    2555        enddo
    2556       enddo
     2538c      do j=1,ntra
     2539c       do k=1,i-1
     2540c        do il=1,ncum
     2541c         if (i.le.inb(il)) then
     2542c          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2543c          cpinv=1.0/cpn(il,i)
     2544c          if (cvflag_grav) then
     2545c           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     2546c     :        *(traent(il,k,i,j)-tra(il,i,j))
     2547c          else
     2548c           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     2549c     :        *(traent(il,k,i,j)-tra(il,i,j))
     2550c          endif
     2551c         endif
     2552c        enddo
     2553c       enddo
     2554c      enddo
    25572555
    25582556      do 490 k=i,nl+1
    … …  
    25812579490   continue
    25822580
    2583       do j=1,ntra
    2584        do k=i,nl+1
    2585         do il=1,ncum
    2586          if (i.le.inb(il) .and. k.le.inb(il)) then
    2587           dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2588           cpinv=1.0/cpn(il,i)
    2589           if (cvflag_grav) then
    2590            ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    2591      :         *(traent(il,k,i,j)-tra(il,i,j))
    2592           else
    2593            ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    2594      :             *(traent(il,k,i,j)-tra(il,i,j))
    2595           endif
    2596          endif ! i and k
    2597         enddo
    2598        enddo
    2599       enddo
     2581c      do j=1,ntra
     2582c       do k=i,nl+1
     2583c        do il=1,ncum
     2584c         if (i.le.inb(il) .and. k.le.inb(il)) then
     2585c          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2586c          cpinv=1.0/cpn(il,i)
     2587c          if (cvflag_grav) then
     2588c           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     2589c     :         *(traent(il,k,i,j)-tra(il,i,j))
     2590c          else
     2591c           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     2592c     :             *(traent(il,k,i,j)-tra(il,i,j))
     2593c          endif
     2594c         endif ! i and k
     2595c        enddo
     2596c       enddo
     2597c      enddo
    26002598
    26012599      do 1400 il=1,ncum
    … …  
    26542652      enddo
    26552653
    2656       do j=1,ntra
    2657        do il=1,ncum
    2658         if (i.le.inb(il)) then
    2659          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2660          cpinv=1.0/cpn(il,i)
    2661 
    2662          if (cvflag_grav) then
    2663           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
    2664      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    2665      :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
    2666          else
    2667           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
    2668      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    2669      :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
    2670          endif
    2671         endif ! i
    2672        enddo
    2673       enddo
    2674 
     2654c      do j=1,ntra
     2655c       do il=1,ncum
     2656c        if (i.le.inb(il)) then
     2657c         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2658c         cpinv=1.0/cpn(il,i)
     2659
     2660c         if (cvflag_grav) then
     2661c          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
     2662c     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     2663c     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
     2664c         else
     2665c          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
     2666c     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     2667c     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
     2668c         endif
     2669c        endif ! i
     2670c       enddo
     2671c      enddo
    26752672
    26762673500   continue
    … …  
    27152712503   continue
    27162713
    2717       do j=1,ntra
    2718        do il=1,ncum
    2719         ex=0.1*ment(il,inb(il),inb(il))
    2720      :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    2721      :      /(ph(il,inb(il))-ph(il,inb(il)+1))
    2722         ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    2723         ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    2724      :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    2725      :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    2726        enddo
    2727       enddo
     2714c      do j=1,ntra
     2715c       do il=1,ncum
     2716c        ex=0.1*ment(il,inb(il),inb(il))
     2717c     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     2718c     :      /(ph(il,inb(il))-ph(il,inb(il)+1))
     2719c        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     2720c        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     2721c     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     2722c     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     2723c       enddo
     2724c      enddo
    27282725
    27292726c
    … …  
    29812978        end
    29822979
     2980      SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
     2981     &                        ment,sij,da,phi)
     2982        implicit none
     2983c inputs:
     2984        integer ncum, nd, na, nloc,len
     2985        real ment(nloc,na,na),sij(nloc,na,na)
     2986c ouputs:
     2987        real da(nloc,na),phi(nloc,na,na)
     2988c local variables:
     2989        integer i,j,k
     2990c       
     2991        da(:,:)=0.
     2992c
     2993        do j=1,na
     2994          do k=1,na
     2995            do i=1,ncum
     2996            da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
     2997            phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
     2998c            print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
     2999            end do
     3000          end do
     3001        end do
     3002   
     3003        return
     3004        end
     3005
    29833006
    29843007      SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
    … …  
    30513074
    30523075
    3053         do 2100 j=1,ntra
    3054          do 2110 k=1,nd ! oct3
    3055           do 2120 i=1,ncum
    3056             ftra1(idcum(i),k,j)=ftra(i,k,j)
    3057  2120     continue
    3058  2110    continue
    3059  2100   continue
     3076c        do 2100 j=1,ntra
     3077c         do 2110 k=1,nd ! oct3
     3078c          do 2120 i=1,ncum
     3079c            ftra1(idcum(i),k,j)=ftra(i,k,j)
     3080c 2120     continue
     3081c 2110    continue
     3082c 2100   continue
    30603083
    30613084        return

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F

    r493 r495  
    315315      call gwprofil
    316316     *       (  nlon , nlev
    317      *       , kgwd   , kdx
     317     *       , kgwd   , kdx , ktest
    318318     *       , ikcrith, icrit
    319319     *       , paphm1, zrho   , zstab ,  zvph
    … …  
    343343c
    344344c
    345       do 523 jl=1,kgwd
    346       ji=kdx(jl)
     345c     do 523 jl=1,kgwd
     346c     ji=kdx(jl)
     347c  Modif vectorisation 02/04/2004
     348      do 523 ji=kidia,kfdia
     349      if(ktest(ji).eq.1) then
     350
    347351      zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
    348352      ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
    … …  
    401405      pte(ji,jk)=0.0
    402406
     407      endif
    403408  523 continue
    404409
    … …  
    10071012      SUBROUTINE GWPROFIL
    10081013     *         ( NLON, NLEV
    1009      *         , kgwd, kdx
     1014     *         , kgwd, kdx , ktest
    10101015     *         , KKCRITH, KCRIT
    10111016     *         , PAPHM1, PRHO   , PSTAB  , PVPH , PRI , PTAU
    … …  
    10751080      integer nlon,nlev
    10761081      INTEGER KKCRITH(NLON),KCRIT(NLON)
    1077      *       ,kdx(nlon)
     1082     *       ,kdx(nlon) , ktest(nlon)
     1083
    10781084C
    10791085      REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1),
    … …  
    11091115      ilevh=KLEV/3
    11101116C
    1111       DO 400 ji=1,kgwd
    1112       jl=kdx(ji)
     1117c     DO 400 ji=1,kgwd
     1118c     jl=kdx(ji)
     1119c  Modif vectorisation 02/04/2004
     1120      DO 400 jl=kidia,kfdia
     1121      if (ktest(jl).eq.1) then
    11131122      Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
    11141123      ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
     1124      endif
    11151125  400 CONTINUE
    11161126 
    … …  
    11231133  410 CONTINUE
    11241134C
    1125       DO 411 ji=1,kgwd
    1126       jl=kdx(ji)
     1135c     DO 411 ji=1,kgwd
     1136c     jl=kdx(ji)
     1137c  Modif vectorisation 02/04/2004
     1138      do 411 jl=kidia,kfdia
     1139      if (ktest(jl).eq.1) then
    11271140           IF(JK.GT.KKCRITH(JL)) THEN
    11281141           PTAU(JL,JK)=ZTAU(JL,KLEV+1)
    … …  
    11321145           PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1)
    11331146           ENDIF
     1147      endif
    11341148 411  CONTINUE             
    11351149C
    … …  
    11431157  420 CONTINUE
    11441158C
    1145       DO 421 ji=1,kgwd
    1146       jl=kdx(ji)
     1159c     DO 421 ji=1,kgwd
     1160c     jl=kdx(ji)
     1161c  Modif vectorisation 02/04/2004
     1162      do 421 jl=kidia,kfdia
     1163      if(ktest(jl).eq.1) then
    11471164      IF(JK.LT.KKCRITH(JL)) THEN
    11481165      ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
    … …  
    11501167      ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
    11511168      ENDIF
     1169      endif
    11521170  421 CONTINUE
    11531171C
    … …  
    11571175C
    11581176                         
    1159       DO 431 ji=1,kgwd
    1160       jl=kdx(ji)
     1177c     DO 431 ji=1,kgwd
     1178c     jl=Kdx(ji)
     1179c  Modif vectorisation 02/04/2004
     1180      do 431 jl=kidia,kfdia
     1181      if(ktest(jl).eq.1) then
     1182
    11611183          IF(JK.LT.KKCRITH(JL)) THEN
    11621184          IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
    … …  
    11781200          ENDIF
    11791201          ENDIF
     1202      endif
    11801203  431 CONTINUE
    11811204 
    … …  
    11851208C  REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
    11861209
    1187       DO 530 ji=1,kgwd
    1188       jl=kdx(ji)
     1210c     DO 530 ji=1,kgwd
     1211c     jl=kdx(ji)
     1212c  Modif vectorisation 02/04/2004
     1213      do 530 jl=kidia,kfdia
     1214      if(ktest(jl).eq.1) then
    11891215      ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
    11901216      ZTAU(JL,NSTRA)=PTAU(JL,NSTRA)
     1217      endif
    11911218 530  CONTINUE     
    11921219
    11931220      DO 531 JK=1,KLEV
    11941221     
    1195       DO 532 ji=1,kgwd
    1196       jl=kdx(ji)
     1222c     DO 532 ji=1,kgwd
     1223c     jl=kdx(ji)
     1224c  Modif vectorisation 02/04/2004
     1225      do 532 jl=kidia,kfdia
     1226      if(ktest(jl).eq.1) then
     1227
    11971228               
    11981229         IF(JK.GT.KKCRITH(JL))THEN
    … …  
    12061237        ENDIF
    12071238           
     1239      endif
    12081240 532  CONTINUE   
    12091241 
    12101242C  REORGANISATION IN THE STRATOSPHERE
    12111243
    1212       DO 533 ji=1,kgwd
    1213       jl=kdx(ji)
     1244c     DO 533 ji=1,kgwd
     1245c     jl=kdx(ji)
     1246c  Modif vectorisation 02/04/2004
     1247      do 533 jl=kidia,kfdia
     1248      if(ktest(jl).eq.1) then
     1249
    12141250
    12151251         IF(JK.LT.NSTRA)THEN
    … …  
    12211257        ENDIF
    12221258
     1259      endif
    12231260 533  CONTINUE
    12241261
    12251262C REORGANISATION IN THE TROPOSPHERE
    12261263
    1227        DO 534 ji=1,kgwd
    1228        jl=kdx(ji)
     1264c      DO 534 ji=1,kgwd
     1265c      jl=kdx(ji)
     1266c  Modif vectorisation 02/04/2004
     1267      do 534 jl=kidia,kfdia
     1268      if(ktest(jl).eq.1) then
     1269
    12291270
    12301271         IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN
    … …  
    12361277
    12371278       ENDIF
     1279      endif
    12381280 534   CONTINUE
    12391281
Note: See TracChangeset for help on using the changeset viewer.