ecearth3/sources/tm5mp/base/phys_convec_ec2tm.F90

!###############################################################################
!
! Convert ecmwf convective fields to TM convective fields.
!
! Original code by Dirk Olivie, KNMI.
!
! Research notes by DO :
!
!  sign:
!   mfup : positive
!   dtup : positive
!   mfdo : negative
!   dtdo : positive
!                             
!  input mass fluxes are contaminated by profiles : 
!  the same for every grid point, different for each time step
!   mfup : not contaminated
!   dtup : between  0.      and -0.99e-19
!   mfdo : between -0.99e-6 and  0.99e-6
!   dtdo : betweeb  0.      and -0.99e-20
!
!  minumum values:
!   mfup : 1.e-6
!   dtup : 1.e-10
!   mfdo : 1.e-6
!   dtdo : 1.e-10
! 
!  remarks:
!  
!   1. start updraft is not always the lowest level : sometimes second or third
!  
!   2. end downdraft is not always lowest level : sometimes 6th level
!  
!   3. in the updraft cloud : delta_updraft = 1.e-4
!                             epsilon_updraft = 1.e-5
!      this automatically/naturally decreases the mass flux : GOOD GOOD GOOD
!  
!   4. there is a part incloud where the mass flus is constant; 
!      later over many layer (order of 10) the mass flux decreases
!   
!   5. detrainment updraft 2 or 3 levels higher (sometimes 15 levels : i = 17565) 
!      then the updraft mass flux
!  
!   6. only detrainment downdraft lowest level different from zero
!  
!   7. the cloudbase of updraft and downdraft correspond sometimes quite well
!
!###############################################################################
!
#define IF_NOTOK_RETURN(action) if (status/=0) then; write (gol,'("in ",a)') rname; call goErr; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; write (gol,'("in ",a)') rname; call goErr; action; return; end if
!
!###############################################################################


module phys_convec_ec2tm

  use GO, only : gol, goPr, goErr

  implicit none
  
  ! --- in/out -----------------------------
  
  private
  
  public  ::  ECconv_to_TMconv
  public  ::  convec_mfldet_to_entdet
  
  ! --- const ---------------------------------------
  
  character(len=*), parameter  ::  mname = 'phys_convec_ec2tm'

  
contains


  ! =============================================================

  !
  ! o Level order is always top->down (ecmwf order)
  !
  ! o ECMWF mass fluxes mflu_ec and mfld_ec are defined for half levels.
  !   According to the original code, the mass fluxes read from the grib
  !   files are valid for the upper half level of a cell.
  !   The bottom values mflu_ec(lm) and mfld_ec(lm) should have been
  !   set to zero.
  !
  

  subroutine ECconv_to_TMconv( lm, zh_ec, &
                               mflu_ec, detu_ec, mfld_ec, detd_ec, &
                               entu   , detu   , entd   , detd   , &
                               status )

    ! --- in/out -------------------------------
    
    integer, intent(in)   ::  lm
    real, intent(in)      ::    zh_ec(0:lm)   ! geopot. height       m        (half lev)
    real, intent(in)      ::  mflu_ec(0:lm)   !   updr. mass flux    kg/m2/s  (half lev)   
    real, intent(in)      ::  detu_ec(1:lm)   !   updr. detr. rate   kg/m3/s  (full lev)
    real, intent(in)      ::  mfld_ec(0:lm)   ! downdr. mass flux    kg/m2/s  (half lev) 
    real, intent(in)      ::  detd_ec(1:lm)   ! downdr. detr. rate   kg/m3/s  (full lev)

    real, intent(out)     ::  entu(1:lm)      !   updr. entr.        kg/m2/s  (full lev)
    real, intent(out)     ::  detu(1:lm)      ! downdr. detr.        kg/m2/s  (full lev)
    real, intent(out)     ::  entd(1:lm)      !   updr. entr.        kg/m2/s  (full lev)
    real, intent(out)     ::  detd(1:lm)      ! downdr. detr.        kg/m2/s  (full lev)
    
    integer, intent(out)  ::  status

    ! --- const -----------------------------------------
  
    character(len=*), parameter  ::  rname = mname//'/ECconv_to_TMconv'

    ! --- local ------------------------------
    
    integer               :: l           ! level number                        

    real                  ::  mflu(0:lm)   !   updr. mass flux    kg/m2/s (half lev)   
    real                  ::  mfld(0:lm)   ! downdr. mass flux    kg/m2/s (half lev) 

    ! layer thickness
    real                  :: dz          ! m

    ! updraft top
    integer               :: uptop       ! layer
    logical               :: uptop_found
    
    ! downdraft top
    integer               :: dotop       ! layer 
    logical               :: dotop_found

    ! --- begin --------------------------------
    
    ! copy ecmwf arrays into local arrays
    mflu = mflu_ec
    detu = detu_ec
    entu = 0.0
    mfld = mfld_ec
    detd = detd_ec
    entd = 0.0

    ! removing small values due to gribbing
    where ( mflu <  1.e-6  ) mflu = 0.0
    where ( detu <  1.e-10 ) detu = 0.0
    where ( mfld > -1.e-6  ) mfld = 0.0
    where ( detd <  1.e-10 ) detd = 0.0

    ! height integral of detrainment rates
    do l = 1, lm
      dz = zh_ec(l-1) - zh_ec(l) 
      detu(l) = detu(l)*dz    ! kg/m3/s  ->  kg/m2/s
      detd(l) = detd(l)*dz    ! kg/m3/s  ->  kg/m2/s
    enddo

    ! find updraftop
    uptop = 0
    uptop_found = .false.
    do l = 1, lm 
      if ( mflu(l) > 0.0 ) then
        uptop       = l
        uptop_found = .true.
        exit
      end if 
    end do

    ! find downdrafttop
    dotop = 0
    dotop_found = .false.
    do l = 1, lm 
      if ( mfld(l) < 0.0 ) then
        dotop       = l
        dotop_found = .true.
        exit
      end if 
    end do

    ! updraft
    if ( uptop_found ) then 
       ! set updr. entr/detr to zero above uptop
       do l = 1, uptop-1
          entu(l) = 0.0
          detu(l) = 0.0
       enddo
       ! loop from uptop to bot:
       do l = uptop, lm
          ! entr  = out through top - in through bot + out by detr
          entu(l) =    mflu(l-1)    -     mflu(l)    +    detu(l)      
       enddo
    else
       ! no updr. at all
       mflu(:) = 0.0
       detu(:) = 0.0
       entu(:) = 0.0
    endif

    ! downdraft
    if ( dotop_found ) then 
       ! set updr. entr/detr to zero above dotop
       !AJS: BUG ? do l = 1, dotop+1
       do l = 1, dotop-1
          detd(l) = 0.0
          entd(l) = 0.0
       enddo  
       ! loop from dotop to bot:
       do l = dotop, lm
          ! entr  = out through top - in through bot + out by detr
          entd(l) =    mfld(l-1)    -     mfld(l)    +   detd(l)      
       enddo  
    else
       ! no downdr. at all
       mfld(:) = 0.0
       detd(:) = 0.0
       entd(:) = 0.0
    endif

    ! check on negative values:
    do l = 1, lm
       if ( entu(l) < 0.0 ) then 
          detu(l) = detu(l) - entu(l)
          entu(l) = 0.0
       endif
       if ( detu(l) < 0.0 ) then 
          entu(l) = entu(l) - detu(l)
          detu(l) = 0.0
       endif
       if ( entd(l) < 0.0 ) then 
          detd(l) = detd(l) - entd(l)
          entd(l) = 0.0
       endif
       if ( detd(l) < 0.0 ) then 
          entd(l) = entd(l) - detd(l)
          detd(l) = 0.0
       endif
    enddo

    ! ok
    status = 0

  end subroutine ECconv_to_TMconv 
  

  ! =============================================================

  !
  ! Generalized version of ECconv_to_TMconv .
  !
  ! o Level order is either 'u' upwards (tm order) or 'd' downwards (ecmwf order)
  !
  ! o ECMWF mass fluxes mflu_ec and mfld_ec are defined for half levels.
  !   According to the original code, the mass fluxes read from the grib
  !   files are valid for the upper half level of a cell.
  !   The bottom values mflu_ec(lm) and mfld_ec(lm) should have been
  !   set to zero.
  !
  

  subroutine convec_mfldet_to_entdet( updo, lm, zh_ec, &
                                      mflu_ec, detu_ec, mfld_ec, detd_ec, &
                                      entu   , detu   , entd   , detd   , &
                                      status )

    ! --- in/out -------------------------------
    
    character(len=1)      ::  updo
    integer, intent(in)   ::  lm
    real, intent(in)      ::    zh_ec(0:lm)   ! geopot. height       m        (half lev)
    real, intent(in)      ::  mflu_ec(0:lm)   !   updr. mass flux    kg/m2/s  (half lev)   
    real, intent(in)      ::  detu_ec(1:lm)   !   updr. detr. rate   kg/m3/s  (full lev)
    real, intent(in)      ::  mfld_ec(0:lm)   ! downdr. mass flux    kg/m2/s  (half lev) 
    real, intent(in)      ::  detd_ec(1:lm)   ! downdr. detr. rate   kg/m3/s  (full lev)

    real, intent(out)     ::  entu(1:lm)      !   updr. entr.        kg/m2/s  (full lev)
    real, intent(out)     ::  detu(1:lm)      ! downdr. detr.        kg/m2/s  (full lev)
    real, intent(out)     ::  entd(1:lm)      !   updr. entr.        kg/m2/s  (full lev)
    real, intent(out)     ::  detd(1:lm)      ! downdr. detr.        kg/m2/s  (full lev)
    
    integer, intent(out)  ::  status

    ! --- const -----------------------------------------
  
    character(len=*), parameter  ::  rname = mname//'/convec_mfldet_to_entdet'

    ! --- local ------------------------------
    
    integer               :: l           ! level number                        

    real                  ::  mflu(0:lm)   !   updr. mass flux    kg/m2/s (half lev)   
    real                  ::  mfld(0:lm)   ! downdr. mass flux    kg/m2/s (half lev) 

    ! layer thickness
    real                  :: dz          ! m

    ! updraft top
    integer               :: uptop       ! layer
    logical               :: uptop_found
    
    ! downdraft top
    integer               :: dotop       ! layer 
    logical               :: dotop_found

    ! --- begin --------------------------------
    
    ! copy ecmwf arrays into local arrays
    mflu = mflu_ec
    detu = detu_ec
    entu = 0.0
    mfld = mfld_ec
    detd = detd_ec
    entd = 0.0

    ! removing small values due to gribbing
    where ( mflu <  1.e-6  ) mflu = 0.0
    where ( detu <  1.e-10 ) detu = 0.0
    where ( mfld > -1.e-6  ) mfld = 0.0
    where ( detd <  1.e-10 ) detd = 0.0

    ! height integral of detrainment rates
    do l = 1, lm
      dz = abs( zh_ec(l) - zh_ec(l-1) )  ! m 
      detu(l) = detu(l)*dz    ! kg/m3/s  ->  kg/m2/s
      detd(l) = detd(l)*dz    ! kg/m3/s  ->  kg/m2/s
    end do
    
    ! levels upwards (TM) or downwards (EC)
    select case ( updo )
    
      ! ***
      
      case ( 'u', 'U' )

        ! find updraftop
        uptop = 0
        uptop_found = .false.
        ! loop over layers from top to bot:
        do l = lm, 1, -1
          ! flux through bottom ? then this is the top level
          if ( mflu(l-1) > 0.0 ) then
            uptop       = l
            uptop_found = .true.
            exit
          end if 
        end do

        ! compute entrainments of updraught
        if ( uptop_found ) then 
          ! loop from bot to uptop:
          do l = 1, uptop
             ! entr  = out through top - in through bot + out by detr
             entu(l) =     mflu(l)     -     mflu(l-1)  +    detu(l)      
          enddo
          ! set updr. entr/detr to zero above uptop
          do l = uptop+1, lm
             entu(l) = 0.0
             detu(l) = 0.0
          enddo
        else
          ! no updr. at all
          detu(:) = 0.0
          entu(:) = 0.0
        endif

        ! find downdrafttop
        dotop = 0
        dotop_found = .false.
        ! loop over layers from top to bot:
        do l = lm, 1, -1
          ! flux through bottom ? then this is the top level
          if ( mfld(l-1) < 0.0 ) then
            dotop       = l
            dotop_found = .true.
            exit
          end if 
        end do

        ! compute entrainments of downdraught
        if ( dotop_found ) then 
          ! loop from bot to dotop:
          do l = 1, dotop
             ! entr  = out through top - in through bot + out by detr
             entd(l) =    mfld(l)      -     mfld(l-1)  +   detd(l)      
          enddo  
          ! set updr. entr/detr to zero above dotop
          do l = dotop+1, lm
             detd(l) = 0.0
             entd(l) = 0.0
          enddo  
        else
          ! no downdr. at all
          detd(:) = 0.0
          entd(:) = 0.0
        endif
        
      ! ***
      
      case ( 'd', 'D' )

        ! find updraftop
        uptop = 0
        uptop_found = .false.
        ! loop over layers from top to bot:
        do l = 1, lm 
          ! flux through bottom ? then this is the top level
          if ( mflu(l) > 0.0 ) then
            uptop       = l
            uptop_found = .true.
            exit
          end if 
        end do

        ! compute entrainments of updraught
        if ( uptop_found ) then 
           ! set updr. entr/detr to zero above uptop
           do l = 1, uptop-1
              entu(l) = 0.0
              detu(l) = 0.0
           enddo
           ! loop from uptop to bot:
           do l = uptop, lm
              ! entr  = out through top - in through bot + out by detr
              entu(l) =    mflu(l-1)    -     mflu(l)    +    detu(l)      
           enddo
        else
           ! no updr. at all
           detu(:) = 0.0
           entu(:) = 0.0
        endif

        ! find downdrafttop
        dotop = 0
        dotop_found = .false.
        ! loop over layers from top to bot:
        do l = 1, lm 
          ! flux through bottom ? then this is the top level
          if ( mfld(l) < 0.0 ) then
            dotop       = l
            dotop_found = .true.
            exit
          end if 
        end do

        ! compute entrainments of downdraught
        if ( dotop_found ) then 
           ! set updr. entr/detr to zero above dotop
           do l = 1, dotop-1
              detd(l) = 0.0
              entd(l) = 0.0
           enddo  
           ! loop from dotop to bot:
           do l = dotop, lm
              ! entr  = out through top - in through bot + out by detr
              entd(l) =    mfld(l-1)    -     mfld(l)    +   detd(l)      
           enddo  
        else
           ! no downdr. at all
           detd(:) = 0.0
           entd(:) = 0.0
        endif
        
      ! ***
      
      case default
      
        write (gol,'("unsupported updo : ",a)') updo; call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      
    end select


    ! check on negative values:
    do l = 1, lm
       if ( entu(l) < 0.0 ) then 
          detu(l) = detu(l) - entu(l)
          entu(l) = 0.0
       endif
       if ( detu(l) < 0.0 ) then 
          entu(l) = entu(l) - detu(l)
          detu(l) = 0.0
       endif
       if ( entd(l) < 0.0 ) then 
          detd(l) = detd(l) - entd(l)
          entd(l) = 0.0
       endif
       if ( detd(l) < 0.0 ) then 
          entd(l) = entd(l) - detd(l)
          detd(l) = 0.0
       endif
    enddo


    ! ok
    status = 0

  end subroutine convec_mfldet_to_entdet
  

end module phys_convec_ec2tm