calcmatrix_gfs.f90 5.52 KB
Newer Older
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
!**********************************************************************
! 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                    *
  !                                                                            *
  !  Changes by C. Forster, November 2005, NCEP GFS version                    *
  !                                                                            *
  !  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
  do kuvz = 2,nuvz
    k = kuvz-1
    phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
    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
  end do

  !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