1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcmatrix(lconv,delt,cbmf)
! o i o
!*****************************************************************************
! *
! This subroutine calculates the matrix describing convective *
! redistribution of mass in a grid column, using the subroutine *
! convect43c.f provided by Kerry Emanuel. *
! *
! Petra Seibert, Bernd C. Krueger, 2000-2001 *
! *
! changed by C. Forster, November 2003 - February 2004 *
! array fmassfrac(nconvlevmax,nconvlevmax) represents *
! the convective redistribution matrix for the particles *
! *
! lconv indicates whether there is convection in this cell, or not *
! delt time step for convection [s] *
! cbmf cloud base mass flux *
! *
!*****************************************************************************
use par_mod
use com_mod
use conv_mod
implicit none
real :: rlevmass,summe
integer :: iflag, k, kk, kuvz
!1-d variables for convection
!variables for redistribution matrix
real :: cbmfold, precip, qprime
real :: tprime, wd, f_qvsat
real :: delt,cbmf
logical :: lconv
lconv = .false.
! calculate pressure at eta levels for use in convect
! and assign temp & spec. hum. to 1D workspace
! -------------------------------------------------------
! pconv(1) is the pressure at the first level above ground
! phconv(k) is the pressure between levels k-1 and k
! dpr(k) is the pressure difference "around" tconv(k)
! phconv(kmax) must also be defined 1/2 level above pconv(kmax)
! Therefore, we define k = kuvz-1 and let kuvz start from 2
! top layer cannot be used for convection because p at top of this layer is
! not given
phconv(1) = psconv
! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
do kuvz = 2,nuvz
k = kuvz-1
pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
dpr(k) = phconv(k) - phconv(kuvz)
qsconv(k) = f_qvsat( pconv(k), tconv(k) )
! initialize mass fractions
do kk=1,nconvlev
fmassfrac(k,kk)=0.
enddo
enddo
!note that Emanuel says it is important
!a. to set this =0. every grid point
!b. to keep this value in the calling programme in the iteration
! CALL CONVECTION
!******************
cbmfold = cbmf
! Convert pressures to hPa, as required by Emanuel scheme
!********************************************************
!!$ do k=1,nconvlev !old
do k=1,nconvlev+1 !bugfix
pconv_hpa(k)=pconv(k)/100.
phconv_hpa(k)=phconv(k)/100.
end do
phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
call convect(nconvlevmax, nconvlev, delt, iflag, &
precip, wd, tprime, qprime, cbmf)
! do not update fmassfrac and cloudbase massflux
! if no convection takes place or
! if a CFL criterion is violated in convect43c.f
if (iflag .ne. 1 .and. iflag .ne. 4) then
cbmf=cbmfold
goto 200
endif
! do not update fmassfrac and cloudbase massflux
! if the old and the new cloud base mass
! fluxes are zero
if (cbmf.le.0..and.cbmfold.le.0.) then
cbmf=cbmfold
goto 200
endif
! Update fmassfrac
! account for mass displaced from level k to level k
lconv = .true.
do k=1,nconvtop
rlevmass = dpr(k)/ga
summe = 0.
do kk=1,nconvtop
fmassfrac(k,kk) = delt*fmass(k,kk)
summe = summe + fmassfrac(k,kk)
end do
fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
end do
200 continue
end subroutine calcmatrix
|
