phylmd/CV30_routines/cv30_closure.f90

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(:) ! (ncum) {2 <= icb <= nl - 3}
    integer, intent(in):: 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(:) ! (ncum) {2 <= icb <= nl - 3}
18	integer, intent(in):: inb(:) ! (ncum)
19	real pbase(klon)
20	real p(:, :) ! (klon, klev)
21	real, intent(in):: ph(:, :) ! (ncum, klev + 1)
22	real tv(klon, klev), buoy(klon, klev)
23
24	! input/output:
25	real sig(klon, klev), w0(klon, klev)
26
27	! output:
28	real cape(klon)
29	real m(klon, klev)
30
31	! Local:
32	integer ncum
33	integer i, j, k, icbmax
34	real deltap, fac, w, amu
35	real dtmin(klon, klev), sigold(klon, klev)
36
37	!-------------------------------------------------------
38
39	ncum = size(icb)
40
41	! Initialization
42
43	do k=1, nl
44	do i=1, ncum
45	m(i, k)=0.0
46	enddo
47	enddo
48
49	! Reset sig(i) and w0(i) for i>inb and i<icb
50
51	! update sig and w0 above LNB:
52
53	do k=1, nl-1
54	do i=1, ncum
55	if ((inb(i) < (nl-1)).and.(k >= (inb(i) + 1)))then
56	sig(i, k)=beta*sig(i, k) &
57	+ 2.alphabuoy(i, inb(i))*ABS(buoy(i, inb(i)))
58	sig(i, k)=AMAX1(sig(i, k), 0.0)
59	w0(i, k)=beta*w0(i, k)
60	endif
61	end do
62	end do
63
64	! compute icbmax:
65
66	icbmax=2
67	do i=1, ncum
68	icbmax=MAX(icbmax, icb(i))
69	end do
70
71	! update sig and w0 below cloud base:
72
73	do k=1, icbmax
74	do i=1, ncum
75	if (k <= icb(i))then
76	sig(i, k)=betasig(i, k)-2.alphabuoy(i, icb(i))buoy(i, icb(i))
77	sig(i, k)=amax1(sig(i, k), 0.0)
78	w0(i, k)=beta*w0(i, k)
79	endif
80	end do
81	end do
82
83	! Reset fractional areas of updrafts and w0 at initial time
84	! and after 10 time steps of no convection
85
86	do k=1, nl-1
87	do i=1, ncum
88	if (sig(i, klev) < 1.5.or.sig(i, klev) > 12.0)then
89	sig(i, k)=0.0
90	w0(i, k)=0.0
91	endif
92	end do
93	end do
94
95	! Calculate convective available potential energy (cape),
96	! vertical velocity (w), fractional area covered by
97	! undilute updraft (sig), and updraft mass flux (m)
98
99	do i=1, ncum
100	cape(i)=0.0
101	end do
102
103	! compute dtmin (minimum buoyancy between ICB and given level k):
104
105	do i=1, ncum
106	do k=1, nl
107	dtmin(i, k)=100.0
108	enddo
109	enddo
110
111	do i=1, ncum
112	do k=1, nl
113	do j=minorig, nl
114	if ((k >= (icb(i) + 1)).and.(k <= inb(i)).and. &
115	(j >= icb(i)).and.(j <= (k-1)))then
116	dtmin(i, k)=AMIN1(dtmin(i, k), buoy(i, j))
117	endif
118	end do
119	end do
120	end do
121
122	! The interval on which cape is computed starts at pbase:
123
124	do k=1, nl
125	do i=1, ncum
126	if ((k >= (icb(i) + 1)).and.(k <= inb(i))) then
127	deltap = MIN(pbase(i), ph(i, k-1))-MIN(pbase(i), ph(i, k))
128	cape(i)=cape(i) + rdbuoy(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)=betasig(i, k) + alphadtmin(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)facSQRT(cape(i)) + beta*w0(i, k)
137	amu=0.5(sig(i, k) + sigold(i, k))w
138	m(i, k)=amu0.007p(i, k)*(ph(i, k)-ph(i, k + 1))/tv(i, k)
139	w0(i, k)=w
140	endif
141
142	end do
143	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