calcmatrix.f90 4.98 KB
Newer Older
1
subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
Matthias Langer's avatar
 
Matthias Langer committed
2
3
4
5
6
7
8
9
10
  !                        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                                *
  !                                                                            *
11
12
  !*****************************************************************************
  ! Changes:                                                                   *
Matthias Langer's avatar
 
Matthias Langer committed
13
14
15
16
  !  changed by C. Forster, November 2003 - February 2004                      *
  !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
  !  the convective redistribution matrix for the particles                    *
  !                                                                            *
17
18
19
20
21
22
  !   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
23
24
25
  !  lconv        indicates whether there is convection in this cell, or not   *
  !  delt         time step for convection [s]                                 *
  !  cbmf         cloud base mass flux                                         *
26
  !  metdata_format format of metdata (ecmwf/gfs)                              *
Matthias Langer's avatar
 
Matthias Langer committed
27
28
29
30
31
32
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod
  use conv_mod
33
  use class_gribfile
Matthias Langer's avatar
 
Matthias Langer committed
34
35
36
37

  implicit none

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

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


  !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