calcmatrix.f90 6.38 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
!**********************************************************************
! 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/>.   *
!**********************************************************************

22
subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
Matthias Langer's avatar
 
Matthias Langer committed
23
24
25
26
27
28
29
30
31
  !                        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                                *
  !                                                                            *
32
33
  !*****************************************************************************
  ! Changes:                                                                   *
Matthias Langer's avatar
 
Matthias Langer committed
34
35
36
37
  !  changed by C. Forster, November 2003 - February 2004                      *
  !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
  !  the convective redistribution matrix for the particles                    *
  !                                                                            *
38
39
40
41
42
43
  !   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
44
45
46
  !  lconv        indicates whether there is convection in this cell, or not   *
  !  delt         time step for convection [s]                                 *
  !  cbmf         cloud base mass flux                                         *
47
  !  metdata_format format of metdata (ecmwf/gfs)                              *
Matthias Langer's avatar
 
Matthias Langer committed
48
49
50
51
52
53
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod
  use conv_mod
54
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
55
56
57
58

  implicit none

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

  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
90
    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
Matthias Langer's avatar
 
Matthias Langer committed
91
92
    pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
    phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
93
94
95
    else
      phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
    endif
Matthias Langer's avatar
 
Matthias Langer committed
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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    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