scale_atmos_phy_mp_convert.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_phy_mp_convert
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  use scale_atmos_phy_mp_suzuki10, only: &
      nbin, &
      nccn, &
      nspc
  use scale_process, only: &
      prc_mpistop, &
      prc_masterrank, &
      prc_myrank
  use scale_comm, only: &
      comm_datatype
  use scale_const, only: &
      pi    => const_pi,    &
      eps   => const_eps,   &
      undef => const_undef, &
      dwatr => const_dwatr, &
      dice  => const_dice
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_mp_bulk2bin
  public :: atmos_phy_mp_bin2bulk

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
  !--- Indeces
  integer, parameter :: il = 1
  integer, parameter :: ic = 2
  integer, parameter :: ip = 3
  integer, parameter :: id = 4
  integer, parameter :: iss= 5
  integer, parameter :: ig = 6
  integer, parameter :: ih = 7

contains
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_bulk2bin(  &
       xs, xe,              & ! [IN]
       ys, ye,              & ! [IN]
       dims,                & ! [IN]
       it,                  & ! [IN]
       rank,                & ! [IN]
       handle,              & ! [IN]
       basename_org,        & ! [IN]
       dens_org,            & ! [IN]
       qtrc_org             ) ! [INOUT]
    use scale_grid_nest, only: &
       parent_imax,     &
       parent_jmax
    use scale_specfunc, only: &
       sf_gamma
    use gtool_file, only: &
       fileread
    use scale_atmos_phy_mp, only: &
       atmos_phy_mp_name, &
       qa_mp, &
       qs_mp
    use scale_atmos_hydrometeor, only: &
       i_qv
    implicit none

    integer,          intent(in) :: xs, xe
    integer,          intent(in) :: ys, ye
    integer,          intent(in) :: dims(:)
    integer,          intent(in) :: it
    integer,          intent(in) :: rank
    integer,          intent(in) :: handle
    character(len=*), intent(in) :: basename_org
    real(RP),         intent(in) :: dens_org(:,:,:)
    real(RP),         intent(inout) :: qtrc_org(:,:,:,:)
!    real(RP),         intent(in) :: dens_org(dims(1)+2,dims(2),dims(3))
!    real(RP),         intent(inout) :: qtrc_org(dims(1)+2,dims(2),dims(3),11)

    real(RP)         :: sigma_sdf(4)
    real(RP)         :: r0_sdf(4)
    real(RP)         :: n0_sdf(4)
    real(RP)         :: rho_sdf(4)
    character(len=11), parameter :: fname_micpara="micpara.dat"
    character(len=H_SHORT)       :: mp_type_outer = "NONE"
    character(len=H_SHORT)       :: mp_type_inner = "NONE"

    namelist / param_atmos_phy_mp_bin2bulk / &
       sigma_sdf, &
       r0_sdf,    &
       n0_sdf,    &
       rho_sdf,   &
       mp_type_outer, &
       mp_type_inner

    real(RP), allocatable :: read3di(:,:,:)
    real(RP), allocatable :: qc_tmp(:,:,:)
    real(RP), allocatable :: qr_tmp(:,:,:)
    real(RP), allocatable :: qi_tmp(:,:,:)
    real(RP), allocatable :: qs_tmp(:,:,:)
    real(RP), allocatable :: qg_tmp(:,:,:)

    integer  :: fid_micpara
    real(RP) :: coef0, coef1, coef2
    real(RP) :: tmp_hyd, n_hyd, lambda_hyd
    real(RP) :: dummy( nbin ), radc( nbin )
    integer  :: k, i, j, iq, iqa, ierr
    integer  :: nnbin, nnspc, nn
    !---------------------------------------------------------------------------

    !--- define coefficients
    coef0 = 4.0_rp/3.0_rp*pi
    coef1 = 4.0_rp/3.0_rp*sqrt(pi/2.0_rp)

    !--- determine the parameters for interpolating SDF from qxx, Nxx of parent domain
    sigma_sdf(1) = 0.2_rp
    sigma_sdf(2) = 0.35_rp
    sigma_sdf(3) = 0.35_rp
    sigma_sdf(4) = 0.35_rp
    r0_sdf(1)    = 5.e-6_rp
    r0_sdf(2)    = 2.61e-6_rp
    r0_sdf(3)    = 5.e-6_rp
    r0_sdf(4)    = 2.61e-6_rp
    n0_sdf(1)    = 8.0e+6_rp
    n0_sdf(2)    = 0.0_rp
    n0_sdf(3)    = 3.0e+6_rp
    n0_sdf(4)    = 4.0e+6_rp
    rho_sdf(1)   = dwatr
    rho_sdf(2)   = dice
    rho_sdf(3)   = 100.0_rp
    rho_sdf(4)   = 400.0_rp

    !---- initiate
    allocate( read3di( parent_imax(handle), parent_jmax(handle), dims(1) ) )
    allocate( qc_tmp(dims(1)+2,dims(2),dims(3)) )
    allocate( qr_tmp(dims(1)+2,dims(2),dims(3)) )
    allocate( qi_tmp(dims(1)+2,dims(2),dims(3)) )
    allocate( qs_tmp(dims(1)+2,dims(2),dims(3)) )
    allocate( qg_tmp(dims(1)+2,dims(2),dims(3)) )

    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_mp_bin2bulk,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       write(*,*) 'xxx Not found namelist PARAM_ATMOS_PHY_MP_BIN2BULK. Check!'
       call prc_mpistop
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP_BIN2BULK. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_mp_bin2bulk)

    if( trim(mp_type_inner) /= "SUZUKI10" ) then
       write(*,*) 'xxx MP_TYPE_INNER should be SUZUKI10 Check!'
       write(*,*) 'Now MP_TYPE_INNER set as ', mp_type_inner
       call prc_mpistop
    endif

    !--- read micpara.dat (microphysical parameter) and broad cast
    fid_micpara = io_get_available_fid()
    !--- open parameter of cloud microphysics
    open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )

    !--- micpara.dat does not exist
    if( ierr == 0 ) then

      read( fid_micpara,* ) nnspc, nnbin

      if( nnbin /= nbin ) then
         write(*,*) 'xxx nbin in inc_tracer and nbin in micpara.dat is different check!'
         call prc_mpistop
      end if

      ! grid parameter
      if( io_l ) write(io_fid_log,*)  '*** Radius of cloud ****'
      do iq = 1, nbin
        read( fid_micpara,* ) nn, dummy( iq ), radc( iq )
        if( io_l ) write(io_fid_log,'(a,1x,i3,1x,a,1x,e15.7,1x,a)') &
                  "Radius of ", iq, "th cloud bin (bin center)= ", radc( iq ) , "[m]"
      end do

      close ( fid_micpara )

    else

      write(*,*) 'xxx micpara.dat does not exist. check!'
      call prc_mpistop

    endif

    call fileread( read3di(:,:,:), basename_org, "QV", it, rank )
    do k = 1, dims(1)
       qtrc_org(k+2,xs:xe,ys:ye,i_qv) = read3di(:,:,k)
    end do

    if( trim(mp_type_outer) == "KESSLER" ) then

      if( io_l ) write(io_fid_log,*)
      if( io_l ) write(io_fid_log,*) '+++ SDF of Bin model is created from'
      if( io_l ) write(io_fid_log,*) '+++ Kessler type Bulk microphysical model'

       call fileread( read3di(:,:,:), basename_org, "QC", it, rank )
       do k = 1, dims(1)
          qc_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       call fileread( read3di(:,:,:), basename_org, "QR", it, rank )
       do k = 1, dims(1)
          qr_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       do k = 1, dims(1)
       do i = xs, xe
       do j = ys, ye

         do iq = 1, nbin
            dummy(iq) = coef1 / sigma_sdf(1) * dwatr * radc( iq )**3 &
                      * exp(  &
                            - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
                            / sigma_sdf(1) / sigma_sdf(1) &
                            )
         enddo

         tmp_hyd = 0.0_rp
         do iq = 1, nbin
            tmp_hyd = tmp_hyd + dummy(iq)
         enddo

         coef2 = ( qc_tmp(k+2,i,j)+qr_tmp(k+2,i,j) ) &
               / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
         do iq = 1, nbin
           qtrc_org(k+2,i,j,qs_mp+iq) = coef2 * dummy(iq)
         enddo

       enddo
       enddo
       enddo

    elseif( trim(mp_type_outer) == "TOMITA08" ) then

       if( io_l ) write(io_fid_log,*)
       if( io_l ) write(io_fid_log,*) '+++ SDF of Bin model is created from'
       if( io_l ) write(io_fid_log,*) '+++ TOMITA08 Bulk microphysical model'

       call fileread( read3di(:,:,:), basename_org, "QC", it, rank )
       do k = 1, dims(1)
          qc_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do
!       do k = 1, dims(1)
!       do i = xs, xe
!       do j = ys, ye
!          qtrc_tmp(k+2,i,j,I_QC) = qtrc_tmp(k+2,i,j,I_QC)*dens_org(k+2,i,j)
!       enddo
!       enddo
!       enddo

       call fileread( read3di(:,:,:), basename_org, "QR", it, rank )
       do k = 1, dims(1)
          qr_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       call fileread( read3di(:,:,:), basename_org, "QI", it, rank )
       do k = 1, dims(1)
          qi_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       call fileread( read3di(:,:,:), basename_org, "QS", it, rank )
       do k = 1, dims(1)
          qs_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       call fileread( read3di(:,:,:), basename_org, "QG", it, rank )
       do k = 1, dims(1)
          qg_tmp(k+2,xs:xe,ys:ye) = read3di(:,:,k)
       end do

       if( nspc == 1 ) then !--- put all hydrometeors to liquid (warm bin)

         do k = 1, dims(1)
         do i = xs, xe
         do j = ys, ye

           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
                        * exp(  &
                              - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
                              / sigma_sdf(1) / sigma_sdf(1) &
                              )
           enddo

           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo

           coef2 = ( &
                     qc_tmp(k+2,i,j)+qr_tmp(k+2,i,j) &
                   + qs_tmp(k+2,i,j)+qi_tmp(k+2,i,j) &
                   + qg_tmp(k+2,i,j) &
                   ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
           do iq = 1, nbin
             qtrc_org(k+2,i,j,qs_mp+(il-1)*nbin+iq) = coef2 * dummy(iq)
           enddo

         enddo
         enddo
         enddo

       elseif( nspc == 7 ) then  !--- put each hydrometeor to each category (ice bin)

         do k = 1, dims(1)
         do i = xs, xe
         do j = ys, ye

           !--- Rain and Cloud put into liquid bin (log-normal)
           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
                        * exp(  &
                              - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
                              / sigma_sdf(1) / sigma_sdf(1) &
                              )
           enddo

           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo

           coef2 = ( qc_tmp(k+2,i,j)+qr_tmp(k+2,i,j) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
           do iq = 1, nbin
             qtrc_org(k+2,i,j,qs_mp+(il-1)*nbin+iq) = coef2 * dummy(iq)
           enddo

           !--- Ice put into plate bin (log-normal)
           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(2) * rho_sdf(2) * radc( iq )**3 &
                        * exp(  &
                              - ( log( radc(iq) )-log( r0_sdf(2) ) )**2*0.5_rp &
                              / sigma_sdf(2) / sigma_sdf(2) &
                              )
           enddo

           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo

           coef2 = qi_tmp(k+2,i,j) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
           do iq = 1, nbin
             qtrc_org(k+2,i,j,qs_mp+(ip-1)*nbin+iq) = coef2 * dummy(iq)
           enddo

           !--- Snow put into snow bin (gamma)
           n_hyd = coef0 * n0_sdf(3) * rho_sdf(3)
           lambda_hyd = ( pi * rho_sdf(3) / 6.0_rp *n0_sdf(3) * sf_gamma(4.0_rp) &
                      / ( qs_tmp(k+2,i,j) &
                        + (0.50_rp-sign(0.50_rp,qs_tmp(k+2,i,j)-eps)) &
                        ) )**(0.25_rp)
           do iq = 1, nbin
              dummy(iq) = n_hyd * radc( iq )**3 &
                        * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
           enddo
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo

           coef2 = qs_tmp(k+2,i,j) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
           do iq = 1, nbin
             qtrc_org(k+2,i,j,qs_mp+(iss-1)*nbin+iq) = coef2 * dummy(iq)
           enddo

           !--- Graupel put into Graupel bin (gamma)
           n_hyd = coef0 * n0_sdf(4) * rho_sdf(4)
           lambda_hyd = ( pi * rho_sdf(4) / 6.0_rp *n0_sdf(4) * sf_gamma(4.0_rp) &
                      / ( qg_tmp(k+2,i,j) &
                        + (0.50_rp-sign(0.50_rp,qg_tmp(k+2,i,j)-eps)) &
                        ) )**(0.25_rp)
           do iq = 1, nbin
              dummy(iq) = n_hyd * radc( iq )**3 &
                        * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
           enddo
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo

           coef2 = qg_tmp(k+2,i,j) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
           do iq = 1, nbin
             qtrc_org(k+2,i,j,qs_mp+(ig-1)*nbin+iq) = coef2 * dummy(iq)
           enddo

         enddo
         enddo
         enddo

       endif

    elseif( trim(mp_type_outer) == "SN14" ) then

       write(*,*) 'SN14 is not supported for MP_TYPE_OUTER now'
       write(*,*) 'Please wait'
       call prc_mpistop

    elseif( trim(mp_type_outer) == "SUZUKI10" ) then

      if( io_l ) write(io_fid_log,*)
      if( io_l ) write(io_fid_log,*) '+++ SDF of Bin model is created directory'
      if( io_l ) write(io_fid_log,*) '+++ from Bin microphysical model'
      do iq = 1, qa_mp
         iqa = qs_mp + iq - 1
         call fileread( read3di(:,:,:), basename_org, atmos_phy_mp_name(iq), it, rank )
         do k = 1, dims(1)
            qtrc_org(k+2,xs:xe,ys:ye,iqa) = read3di(:,:,k)
         end do
      end do

    else

       write(*,*) 'MP_TYPE_OUTER should be KESSLER, TOMITA08, or SUZUKI10'
       write(*,*) 'Please check! Now MP_TYPE_OUTER set as ', mp_type_outer
       call prc_mpistop

    endif

    deallocate( read3di )
    deallocate( qc_tmp )
    deallocate( qr_tmp )
    deallocate( qi_tmp )
    deallocate( qs_tmp )
    deallocate( qg_tmp )

    return
  end subroutine atmos_phy_mp_bulk2bin

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_bin2bulk
    !---------------------------------------------------------------------------

    return
  end subroutine atmos_phy_mp_bin2bulk

end module scale_atmos_phy_mp_convert