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