phylmd/CV30_routines/cv30_closure.f

module cv30_closure_m

  implicit none

contains

  SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
       sig, w0, cape, m)

    ! Vectorization: S. Bony

    use cv30_param_m, only: alpha, beta, dtcrit, minorig, nl
    use cv_thermo_m, only: rrd

    ! input:
    integer, intent(in):: ncum, nd, nloc
    integer, intent(in):: icb(nloc), inb(nloc)
    real pbase(nloc)
    real p(nloc, nd), ph(nloc, nd+1)
    real tv(nloc, nd), buoy(nloc, nd)

    ! input/output:
    real sig(nloc, nd), w0(nloc, nd)

    ! output:
    real cape(nloc)
    real m(nloc, nd)

    ! Local:
    integer i, j, k, icbmax
    real deltap, fac, w, amu
    real dtmin(nloc, nd), sigold(nloc, nd)

    !-------------------------------------------------------

    ! Initialization

    do k=1, nl
       do i=1, ncum
          m(i, k)=0.0
       enddo
    enddo

    ! Reset sig(i) and w0(i) for i>inb and i<icb

    ! update sig and w0 above LNB:

    do k=1, nl-1
       do i=1, ncum
          if ((inb(i) < (nl-1)).and.(k >= (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
       end do
    end do

    ! compute icbmax:

    icbmax=2
    do i=1, ncum
       icbmax=MAX(icbmax, icb(i))
    end do

    ! update sig and w0 below cloud base:

    do k=1, icbmax
       do i=1, ncum
          if (k <= 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
       end do
    end do

    ! Reset fractional areas of updrafts and w0 at initial time
    ! and after 10 time steps of no convection

    do k=1, nl-1
       do i=1, ncum
          if (sig(i, nd) < 1.5.or.sig(i, nd) > 12.0)then
             sig(i, k)=0.0
             w0(i, k)=0.0
          endif
       end do
    end do

    ! Calculate convective available potential energy (cape),
    ! vertical velocity (w), fractional area covered by
    ! undilute updraft (sig), and updraft mass flux (m)

    do i=1, ncum
       cape(i)=0.0
    end do

    ! compute dtmin (minimum buoyancy between ICB and given level k):

    do i=1, ncum
       do k=1, nl
          dtmin(i, k)=100.0
       enddo
    enddo

    do i=1, ncum
       do k=1, nl
          do j=minorig, nl
             if ((k >= (icb(i)+1)).and.(k <= inb(i)).and. &
                  (j >= icb(i)).and.(j <= (k-1)))then
                dtmin(i, k)=AMIN1(dtmin(i, k), buoy(i, j))
             endif
          end do
       end do
    end do

    ! the interval on which cape is computed starts at pbase :

    do k=1, nl
       do i=1, ncum

          if ((k >= (icb(i)+1)).and.(k <= 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)

             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

       end do
    end do

    do 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))
    end do

  end SUBROUTINE cv30_closure

end module cv30_closure_m
1	guez	187	module cv30_closure_m
2	guez	47
3	guez	187	implicit none
4	guez	47
5	guez	187	contains
6	guez	47
7	guez	187	SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
8			sig, w0, cape, m)
9	guez	47
10	guez	187	! Vectorization: S. Bony
11	guez	47
12	guez	187	use cv30_param_m, only: alpha, beta, dtcrit, minorig, nl
13	guez	190	use cv_thermo_m, only: rrd
14	guez	47
15	guez	187	! input:
16			integer, intent(in):: ncum, nd, nloc
17			integer, intent(in):: icb(nloc), inb(nloc)
18			real pbase(nloc)
19			real p(nloc, nd), ph(nloc, nd+1)
20			real tv(nloc, nd), buoy(nloc, nd)
21	guez	47
22	guez	187	! input/output:
23			real sig(nloc, nd), w0(nloc, nd)
24	guez	47
25	guez	187	! output:
26			real cape(nloc)
27			real m(nloc, nd)
28	guez	47
29	guez	187	! Local:
30			integer i, j, k, icbmax
31			real deltap, fac, w, amu
32			real dtmin(nloc, nd), sigold(nloc, nd)
33	guez	47
34	guez	187	!-------------------------------------------------------
35	guez	47
36	guez	187	! Initialization
37	guez	47
38	guez	187	do k=1, nl
39			do i=1, ncum
40			m(i, k)=0.0
41			enddo
42			enddo
43	guez	47
44	guez	187	! Reset sig(i) and w0(i) for i>inb and i<icb
45	guez	47
46	guez	187	! update sig and w0 above LNB:
47	guez	47
48	guez	187	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.alphabuoy(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	guez	47
59	guez	187	! compute icbmax:
60	guez	47
61	guez	187	icbmax=2
62			do i=1, ncum
63			icbmax=MAX(icbmax, icb(i))
64			end do
65	guez	47
66	guez	187	! update sig and w0 below cloud base:
67	guez	47
68	guez	187	do k=1, icbmax
69			do i=1, ncum
70			if (k <= icb(i))then
71			sig(i, k)=betasig(i, k)-2.alphabuoy(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	guez	47
78	guez	187	! Reset fractional areas of updrafts and w0 at initial time
79			! and after 10 time steps of no convection
80	guez	47
81	guez	187	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	guez	47
90	guez	187	! Calculate convective available potential energy (cape),
91			! vertical velocity (w), fractional area covered by
92			! undilute updraft (sig), and updraft mass flux (m)
93	guez	47
94	guez	187	do i=1, ncum
95	guez	47	cape(i)=0.0
96	guez	187	end do
97	guez	47
98	guez	187	! compute dtmin (minimum buoyancy between ICB and given level k):
99	guez	47
100	guez	187	do i=1, ncum
101			do k=1, nl
102			dtmin(i, k)=100.0
103	guez	47	enddo
104	guez	187	enddo
105	guez	47
106	guez	187	do i=1, ncum
107			do k=1, nl
108			do j=minorig, nl
109			if ((k >= (icb(i)+1)).and.(k <= inb(i)).and. &
110			(j >= icb(i)).and.(j <= (k-1)))then
111			dtmin(i, k)=AMIN1(dtmin(i, k), buoy(i, j))
112			endif
113			end do
114			end do
115			end do
116	guez	47
117	guez	187	! the interval on which cape is computed starts at pbase :
118	guez	47
119	guez	187	do k=1, nl
120			do i=1, ncum
121	guez	47
122	guez	187	if ((k >= (icb(i)+1)).and.(k <= inb(i))) then
123	guez	47
124	guez	187	deltap = MIN(pbase(i), ph(i, k-1))-MIN(pbase(i), ph(i, k))
125			cape(i)=cape(i)+rrdbuoy(i, k-1)deltap/p(i, k-1)
126			cape(i)=AMAX1(0.0, cape(i))
127			sigold(i, k)=sig(i, k)
128	guez	47
129	guez	187	sig(i, k)=betasig(i, k)+alphadtmin(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)facSQRT(cape(i))+beta*w0(i, k)
134			amu=0.5(sig(i, k)+sigold(i, k))w
135			m(i, k)=amu0.007p(i, k)*(ph(i, k)-ph(i, k+1))/tv(i, k)
136			w0(i, k)=w
137			endif
138	guez	47
139	guez	187	end do
140			end do
141	guez	47
142	guez	187	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	guez	47
151	guez	187	end SUBROUTINE cv30_closure
152	guez	47
153	guez	187	end module cv30_closure_m