/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_closure.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_closure.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/Sources/phylmd/CV3_routines/cv3_closure.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC trunk/Sources/phylmd/CV30_routines/cv30_closure.f revision 187 by guez, Mon Mar 21 18:01:02 2016 UTC
# Line 1  Line 1 
1    module cv30_closure_m
2    
3        SUBROUTINE cv3_closure(nloc,ncum,nd,icb,inb &    implicit none
4                              ,pbase,p,ph,tv,buoy &  
5                              ,sig,w0,cape,m)  contains
6              use cv3_param_m  
7              use cvthermo    SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
8        implicit none         sig, w0, cape, m)
9    
10  !===================================================================      ! Vectorization: S. Bony
11  ! ---  CLOSURE OF CONVECT3  
12  !      use cv30_param_m, only: alpha, beta, dtcrit, minorig, nl
13  ! vectorization: S. Bony      use cvthermo, only: rrd
14  !===================================================================  
15        ! input:
16        integer, intent(in):: ncum, nd, nloc
17  ! input:      integer, intent(in):: icb(nloc), inb(nloc)
18        integer, intent(in):: ncum, nd, nloc      real pbase(nloc)
19        integer icb(nloc), inb(nloc)      real p(nloc, nd), ph(nloc, nd+1)
20        real pbase(nloc)      real tv(nloc, nd), buoy(nloc, nd)
21        real p(nloc,nd), ph(nloc,nd+1)  
22        real tv(nloc,nd), buoy(nloc,nd)      ! input/output:
23        real sig(nloc, nd), w0(nloc, nd)
24  ! input/output:  
25        real sig(nloc,nd), w0(nloc,nd)      ! output:
26        real cape(nloc)
27  ! output:      real m(nloc, nd)
28        real cape(nloc)  
29        real m(nloc,nd)      ! Local:
30        integer i, j, k, icbmax
31  ! local variables:      real deltap, fac, w, amu
32        integer i, j, k, icbmax      real dtmin(nloc, nd), sigold(nloc, nd)
33        real deltap, fac, w, amu  
34        real dtmin(nloc,nd), sigold(nloc,nd)      !-------------------------------------------------------
35    
36        ! Initialization
37  ! -------------------------------------------------------  
38  ! -- Initialization      do k=1, nl
39  ! -------------------------------------------------------         do i=1, ncum
40              m(i, k)=0.0
       do k=1,nl  
        do i=1,ncum  
         m(i,k)=0.0  
41         enddo         enddo
42        enddo      enddo
43    
44        ! Reset sig(i) and w0(i) for i>inb and i<icb
45    
46  ! -------------------------------------------------------      ! update sig and w0 above LNB:
 ! -- Reset sig(i) and w0(i) for i>inb and i<icb  
 ! -------------------------------------------------------  
   
 ! update sig and w0 above LNB:  
   
       do 100 k=1,nl-1  
        do 110 i=1,ncum  
         if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then  
          sig(i,k)=beta*sig(i,k) &  
                   +2.*alpha*buoy(i,inb(i))*ABS(buoy(i,inb(i)))  
          sig(i,k)=AMAX1(sig(i,k),0.0)  
          w0(i,k)=beta*w0(i,k)  
         endif  
  110   continue  
  100  continue  
   
 ! compute icbmax:  
   
       icbmax=2  
       do 200 i=1,ncum  
         icbmax=MAX(icbmax,icb(i))  
  200  continue  
   
 ! update sig and w0 below cloud base:  
   
       do 300 k=1,icbmax  
        do 310 i=1,ncum  
         if (k.le.icb(i))then  
          sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i))  
          sig(i,k)=amax1(sig(i,k),0.0)  
          w0(i,k)=beta*w0(i,k)  
         endif  
 310    continue  
 300    continue  
   
 !!      if(inb.lt.(nl-1))then  
 !!         do 85 i=inb+1,nl-1  
 !!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*  
 !!     1              abs(buoy(inb))  
 !!            sig(i)=amax1(sig(i),0.0)  
 !!            w0(i)=beta*w0(i)  
 !!   85    continue  
 !!      end if  
   
 !!      do 87 i=1,icb  
 !!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)  
 !!         sig(i)=amax1(sig(i),0.0)  
 !!         w0(i)=beta*w0(i)  
 !!   87 continue  
   
 ! -------------------------------------------------------------  
 ! -- Reset fractional areas of updrafts and w0 at initial time  
 ! -- and after 10 time steps of no convection  
 ! -------------------------------------------------------------  
   
       do 400 k=1,nl-1  
        do 410 i=1,ncum  
         if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then  
          sig(i,k)=0.0  
          w0(i,k)=0.0  
         endif  
  410   continue  
  400  continue  
   
 ! -------------------------------------------------------------  
 ! -- Calculate convective available potential energy (cape),  
 ! -- vertical velocity (w), fractional area covered by  
 ! -- undilute updraft (sig), and updraft mass flux (m)  
 ! -------------------------------------------------------------  
47    
48        do 500 i=1,ncum      do k=1, nl-1
49           do i=1, ncum
50              if ((inb(i) < (nl-1)).and.(k >= (inb(i)+1)))then
51                 sig(i, k)=beta*sig(i, k) &
52                      +2.*alpha*buoy(i, inb(i))*ABS(buoy(i, inb(i)))
53                 sig(i, k)=AMAX1(sig(i, k), 0.0)
54                 w0(i, k)=beta*w0(i, k)
55              endif
56           end do
57        end do
58    
59        ! compute icbmax:
60    
61        icbmax=2
62        do i=1, ncum
63           icbmax=MAX(icbmax, icb(i))
64        end do
65    
66        ! update sig and w0 below cloud base:
67    
68        do k=1, icbmax
69           do i=1, ncum
70              if (k <= icb(i))then
71                 sig(i, k)=beta*sig(i, k)-2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
72                 sig(i, k)=amax1(sig(i, k), 0.0)
73                 w0(i, k)=beta*w0(i, k)
74              endif
75           end do
76        end do
77    
78        ! Reset fractional areas of updrafts and w0 at initial time
79        ! and after 10 time steps of no convection
80    
81        do k=1, nl-1
82           do i=1, ncum
83              if (sig(i, nd) < 1.5.or.sig(i, nd) > 12.0)then
84                 sig(i, k)=0.0
85                 w0(i, k)=0.0
86              endif
87           end do
88        end do
89    
90        ! Calculate convective available potential energy (cape),
91        ! vertical velocity (w), fractional area covered by
92        ! undilute updraft (sig), and updraft mass flux (m)
93    
94        do i=1, ncum
95         cape(i)=0.0         cape(i)=0.0
96   500  continue      end do
97    
98  ! compute dtmin (minimum buoyancy between ICB and given level k):      ! compute dtmin (minimum buoyancy between ICB and given level k):
99    
100        do i=1,ncum      do i=1, ncum
101         do k=1,nl         do k=1, nl
102           dtmin(i,k)=100.0            dtmin(i, k)=100.0
103         enddo         enddo
104        enddo      enddo
105    
106        do 550 i=1,ncum      do i=1, ncum
107         do 560 k=1,nl         do k=1, nl
108           do 570 j=minorig,nl            do j=minorig, nl
109            if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and. &               if ((k >= (icb(i)+1)).and.(k <= inb(i)).and. &
110                 (j.ge.icb(i)).and.(j.le.(k-1)) )then                    (j >= icb(i)).and.(j <= (k-1)))then
111             dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))                  dtmin(i, k)=AMIN1(dtmin(i, k), buoy(i, j))
112                 endif
113              end do
114           end do
115        end do
116    
117        ! the interval on which cape is computed starts at pbase :
118    
119        do k=1, nl
120           do i=1, ncum
121    
122              if ((k >= (icb(i)+1)).and.(k <= inb(i))) then
123    
124                 deltap = MIN(pbase(i), ph(i, k-1))-MIN(pbase(i), ph(i, k))
125                 cape(i)=cape(i)+rrd*buoy(i, k-1)*deltap/p(i, k-1)
126                 cape(i)=AMAX1(0.0, cape(i))
127                 sigold(i, k)=sig(i, k)
128    
129                 sig(i, k)=beta*sig(i, k)+alpha*dtmin(i, k)*ABS(dtmin(i, k))
130                 sig(i, k)=amax1(sig(i, k), 0.0)
131                 sig(i, k)=amin1(sig(i, k), 0.01)
132                 fac=AMIN1(((dtcrit-dtmin(i, k))/dtcrit), 1.0)
133                 w=(1.-beta)*fac*SQRT(cape(i))+beta*w0(i, k)
134                 amu=0.5*(sig(i, k)+sigold(i, k))*w
135                 m(i, k)=amu*0.007*p(i, k)*(ph(i, k)-ph(i, k+1))/tv(i, k)
136                 w0(i, k)=w
137            endif            endif
  570     continue  
  560   continue  
  550  continue  
   
 ! the interval on which cape is computed starts at pbase :  
   
       do 600 k=1,nl  
        do 610 i=1,ncum  
   
         if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then  
   
          deltap = MIN(pbase(i),ph(i,k-1))-MIN(pbase(i),ph(i,k))  
          cape(i)=cape(i)+rrd*buoy(i,k-1)*deltap/p(i,k-1)  
          cape(i)=AMAX1(0.0,cape(i))  
          sigold(i,k)=sig(i,k)  
   
 !         dtmin(i,k)=100.0  
 !         do 97 j=icb(i),k-1 ! mauvaise vectorisation  
 !          dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))  
 !  97     continue  
   
          sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*ABS(dtmin(i,k))  
          sig(i,k)=amax1(sig(i,k),0.0)  
          sig(i,k)=amin1(sig(i,k),0.01)  
          fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)  
          w=(1.-beta)*fac*SQRT(cape(i))+beta*w0(i,k)  
          amu=0.5*(sig(i,k)+sigold(i,k))*w  
          m(i,k)=amu*0.007*p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k)  
          w0(i,k)=w  
         endif  
   
  610   continue  
  600  continue  
   
       do 700 i=1,ncum  
        w0(i,icb(i))=0.5*w0(i,icb(i)+1)  
        m(i,icb(i))=0.5*m(i,icb(i)+1) &  
                    *(ph(i,icb(i))-ph(i,icb(i)+1)) &  
                    /(ph(i,icb(i)+1)-ph(i,icb(i)+2))  
        sig(i,icb(i))=sig(i,icb(i)+1)  
        sig(i,icb(i)-1)=sig(i,icb(i))  
  700  continue  
   
   
 !!      cape=0.0  
 !!      do 98 i=icb+1,inb  
 !!         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))  
 !!         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)  
 !!         dcape=rrd*buoy(i-1)*deltap/p(i-1)  
 !!         dlnp=deltap/p(i-1)  
 !!         cape=amax1(0.0,cape)  
 !!         sigold=sig(i)  
   
 !!         dtmin=100.0  
 !!         do 97 j=icb,i-1  
 !!            dtmin=amin1(dtmin,buoy(j))  
 !!   97    continue  
   
 !!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)  
 !!         sig(i)=amax1(sig(i),0.0)  
 !!         sig(i)=amin1(sig(i),0.01)  
 !!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)  
 !!         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)  
 !!         amu=0.5*(sig(i)+sigold)*w  
 !!         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)  
 !!         w0(i)=w  
 !!   98 continue  
 !!      w0(icb)=0.5*w0(icb+1)  
 !!      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))  
 !!      sig(icb)=sig(icb+1)  
 !!      sig(icb-1)=sig(icb)  
138    
139         return         end do
140         end      end do
141    
142        do i=1, ncum
143           w0(i, icb(i))=0.5*w0(i, icb(i)+1)
144           m(i, icb(i))=0.5*m(i, icb(i)+1) &
145                *(ph(i, icb(i))-ph(i, icb(i)+1)) &
146                /(ph(i, icb(i)+1)-ph(i, icb(i)+2))
147           sig(i, icb(i))=sig(i, icb(i)+1)
148           sig(i, icb(i)-1)=sig(i, icb(i))
149        end do
150    
151      end SUBROUTINE cv30_closure
152    
153    end module cv30_closure_m
Legend:
Removed from v.134  
changed lines
  Added in v.187
  ViewVC Help
Powered by ViewVC 1.1.21