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

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

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

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