phylmd/CV30_routines/cv30_closure.f

module cv30_closure_m

  implicit none

contains

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

    ! CLOSURE
    ! Vectorization: S. Bony

    use cv30_param_m, only: alpha, beta, dtcrit, minorig, nl
    USE dimphy, ONLY: klev, klon
    use suphec_m, only: rd

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

    ! input/output:
    real sig(klon, klev), w0(klon, klev)

    ! output:
    real cape(klon)
    real m(klon, klev)

    ! Local:
    integer ncum
    integer i, j, k, icbmax
    real deltap, fac, w, amu
    real dtmin(klon, klev), sigold(klon, klev)

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

    ncum = size(icb)

    ! 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, klev) < 1.5.or.sig(i, klev) > 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) + rd*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	module cv30_closure_m
2
3	implicit none
4
5	contains
6
7	SUBROUTINE cv30_closure(icb, inb, pbase, p, ph, tv, buoy, sig, w0, cape, m)
8
9	! CLOSURE
10	! Vectorization: S. Bony
11
12	use cv30_param_m, only: alpha, beta, dtcrit, minorig, nl
13	USE dimphy, ONLY: klev, klon
14	use suphec_m, only: rd
15
16	! input:
17	integer, intent(in):: icb(:), inb(:) ! (ncum)
18	real pbase(klon)
19	real p(:, :) ! (klon, klev)
20	real, intent(in):: ph(:, :) ! (ncum, klev + 1)
21	real tv(klon, klev), buoy(klon, klev)
22
23	! input/output:
24	real sig(klon, klev), w0(klon, klev)
25
26	! output:
27	real cape(klon)
28	real m(klon, klev)
29
30	! Local:
31	integer ncum
32	integer i, j, k, icbmax
33	real deltap, fac, w, amu
34	real dtmin(klon, klev), sigold(klon, klev)
35
36	!-------------------------------------------------------
37
38	ncum = size(icb)
39
40	! Initialization
41
42	do k=1, nl
43	do i=1, ncum
44	m(i, k)=0.0
45	enddo
46	enddo
47
48	! Reset sig(i) and w0(i) for i>inb and i<icb
49
50	! update sig and w0 above LNB:
51
52	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.alphabuoy(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)=betasig(i, k)-2.alphabuoy(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
100	end do
101
102	! compute dtmin (minimum buoyancy between ICB and given level k):
103
104	do i=1, ncum
105	do k=1, nl
106	dtmin(i, k)=100.0
107	enddo
108	enddo
109
110	do i=1, ncum
111	do k=1, nl
112	do j=minorig, nl
113	if ((k >= (icb(i) + 1)).and.(k <= inb(i)).and. &
114	(j >= icb(i)).and.(j <= (k-1)))then
115	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) + rdbuoy(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)=betasig(i, k) + alphadtmin(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)facSQRT(cape(i)) + beta*w0(i, k)
136	amu=0.5(sig(i, k) + sigold(i, k))w
137	m(i, k)=amu0.007p(i, k)*(ph(i, k)-ph(i, k + 1))/tv(i, k)
138	w0(i, k)=w
139	endif
140
141	end do
142	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