source: flexpart.git/src/calcmatrix.f90 @ 02095e3

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 02095e3 was d8eed02, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Minor edits

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