calcmatrix.f90 5.09 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

4
subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
Matthias Langer's avatar
 
Matthias Langer committed
5
6
7
8
9
10
11
12
13
  !                        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                                *
  !                                                                            *
14
15
  !*****************************************************************************
  ! Changes:                                                                   *
Matthias Langer's avatar
 
Matthias Langer committed
16
17
18
19
  !  changed by C. Forster, November 2003 - February 2004                      *
  !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
  !  the convective redistribution matrix for the particles                    *
  !                                                                            *
20
21
22
23
24
25
  !   Unified ECMWF and GFS builds                                             *
  !   Marian Harustak, 12.5.2017                                               *
  !     - Merged calcmatrix and calcmatrix_gfs into one routine using if-then  *
  !       for meteo-type dependent code                                        *
  !*****************************************************************************
  !                                                                            *
Matthias Langer's avatar
 
Matthias Langer committed
26
27
28
  !  lconv        indicates whether there is convection in this cell, or not   *
  !  delt         time step for convection [s]                                 *
  !  cbmf         cloud base mass flux                                         *
29
  !  metdata_format format of metdata (ecmwf/gfs)                              *
Matthias Langer's avatar
 
Matthias Langer committed
30
31
32
33
34
35
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod
  use conv_mod
36
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
37
38
39
40

  implicit none

  real :: rlevmass,summe
41
  integer :: metdata_format
Matthias Langer's avatar
 
Matthias Langer committed
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

  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
72
    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
Matthias Langer's avatar
 
Matthias Langer committed
73
74
    pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
    phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
75
76
77
    else
      phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
    endif
Matthias Langer's avatar
 
Matthias Langer committed
78
79
80
81
82
83
    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.
Espen Sollum's avatar
Espen Sollum committed
84
85
    end do
  end do
Matthias Langer's avatar
 
Matthias Langer committed
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


  !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