!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DisturbEnv
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: disturbenv.f90,v 1.1.1.1 2006/04/25 03:43:58 deepconv Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!$B%G%U%)%k%H$N>qMp$rM?$($k$?$a$N%k!<%A%s(B. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * $B@_Dj2DG=$J>qMp$N%?%$%W$rA}$d$9(B
!  * $B%(%i!<=hM}$K(B gf4f90io $B$rMxMQ(B
!

subroutine DisturbEnv( cfgfile,         &
  & pz_VelX,     xr_VelZ,               &
  & xz_Exner,    xz_PotTemp,   xza_MixRt )
  !
  !$B%G%U%)%k%H$N>qMp$rM?$($k$?$a$N%k!<%A%s(B. 
  !
  
  !$B%b%8%e!<%kFI$_9~$_(B
  use gridset,    only:DimXMin,         &! $BG[Ns$N(B X $BJ}8~$N2<8B(B
    &                  DimXMax,         &! $BG[Ns$N(B X $BJ}8~$N>e8B(B
    &                  DimZMin,         &! $BG[Ns$N(B Z $BJ}8~$N2<8B(B
    &                  DimZMax,         &! $BG[Ns$N(B Z $BJ}8~$N>e8B(B
    &                  MarginZ,          &! $B7W;;NN0h$N%^!<%8%s(B
    &                  RegZMin,         &! $BJ*M}NN0h$N(B Z $BJ}8~$N2<8B(B
    &                  DelZ,            &! Z $BJ}8~$N9o$_I}(B
    &                  SpcNum,          &! $B6E=L@-2=3X<o$N?t(B
    &                  XMin, XMax,      &! 
    &                  ZMin, ZMax,      &! 
    &                  s_X,             &! X $B:BI8<4(B($B%9%+%i!<3J;RE@(B)
    &                  s_Z               ! Z $B:BI8<4(B($B%9%+%i!<3J;RE@(B)
  use fileset,    only:RandomFile        ! $BMp?t%U%!%$%k(B
  use basicset,   only:                  &
    &                  SpcWetMolFr,     &!$B6E=L@.J,$N=i4|%b%kHf(B
    &                  MolWtWet,        &!$B6E=L@.J,$NJ,;RNL(B
    &                  MolWtDry,        &!$B4%Ag@.J,$NJ,;RNL(B
    &                  xz_TempBasicZ,   &! $B4pK\>l$N29EY(B
    &                  xz_PressBasicZ,  &! $B4pK\>l$N05NO(B
    &                  xza_MixRtBasicZ    ! $B4pK\>l$N:.9gHf(B
  use Boundary, only:  xza_BoundaryXCyc_xza , &!
    &                  xza_BoundaryZSym_xza 
  use ECCM,     only:  ECCM_MolFr


  !$B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !$BJQ?tDj5A(B
  character(*), intent(in) :: cfgfile
  real(8), intent(out)  :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !$B?eJ?IwB.$N>qMp@.J,(B
  real(8), intent(out)  :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax) 
                                    !$B1tD>IwB.$N>qMp@.J,(B 
  real(8), intent(out)  :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !$B%(%/%9%J!<4X?t$N>qMp@.J,(B 
  real(8), intent(out)  :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !$B290L$N>qMp@.J,(B 
  real(8), intent(out)  :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum)  
                                    !$B6E=L@.J,$N:.9gHf(B($B>qMp@.J,(B)
  real(8), parameter         :: Pi = 3.1415926535897932385d0 
                                    !$B1_<~N((B
  integer       :: i, k, s, n
  real(8)       :: Humidity         !$BAjBP<>EY(B
  real(8)       :: XcRate           !$B>qMp$NCf?40LCV(B($B?eJ?J}8~(B)$B$NNN0h$KBP$9$k3d9g(B
  real(8)       :: XrRate           !$B>qMp$NH>7B(B($B?eJ?J}8~(B)$B$NNN0h$KBP$9$k3d9g(B
  real(8)       :: ZcRate           !$B>qMp$NCf?40LCV(B($B1tD>J}8~(B)$B$NNN0h$KBP$9$k3d9g(B
  real(8)       :: ZrRate           !$B>qMp$NH>7B(B($B1tD>J}8~(B)$B$NNN0h$KBP$9$k3d9g(B
  real(8)       :: Xc               !$B>qMp$NCf?40LCV(B($B?eJ?J}8~(B)
  real(8)       :: Xr               !$B>qMp$NH>7B(B($B?eJ?J}8~(B)
  real(8)       :: Zc               !$B>qMp$NCf?40LCV(B($B1tD>J}8~(B)
  real(8)       :: Zr               !$B>qMp$NH>7B(B($B1tD>J}8~(B)
  real(8)       :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
                                    !$B>qMp$N:GBgCM$KBP$9$k3d9g(B
  real(8)       :: DelMax           !$B290L>qMp$N:GBgCM(B
  real(8)       :: Random           !$B%U%!%$%k$+$i<hF@$7$?Mp?t(B
  character(20) :: Type             !$B290L>qMp$N%?%$%W(B
  real(8)       :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                    !$B%b%kHf(B

  !-------------------------------------------------------------
  ! $B=i4|2=(B
  !-------------------------------------------------------------
  !NAMELIST $B%U%!%$%k$NFI$_9~$_(B
  NAMELIST /disturbset/ &
    & Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, &
    & Humidity

  open (10, FILE=cfgfile)
  read(10, NML=disturbset)
  close(10)
  
  !$B=i4|2=(B
  pz_VelX = 0.0d0
  xr_VelZ = 0.0d0
  xz_Exner = 0.0d0
  xz_PotTemp = 0.0d0
  xza_MixRt   = 0.0d0

  !-------------------------------------------------------------
  ! $B290L$N>qMp(B
  !-------------------------------------------------------------
  select case(Type)

  case ("Thermal")
    ! $B;XDj$5$l$?>l=j$K290L$N(B 1 $B$D$NG.E*%5!<%^%k$rM?$($k(B. 
    
    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        beta(i,k) =                                 &
          & (                                       &
          &      ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
          &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0  &
          &  ) ** 5.0d-1
      end do
    end do
    
    where ( beta < 1.0d0 )
      xz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
    end where


  case ("Thermal-Gauss")
    ! $B;XDj$5$l$?>l=j$K(B, $B%,%&%7%"%s$J290L$N>qMp$rM?$($k(B. 

    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    do i = DimXMin, DimXMax
      xz_PotTemp(i,:) = &
        & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 ) 
    end do
    
    do k = DimZMin, DimZMax
      xz_PotTemp(:,k) =  &
        & xz_PotTemp(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
    end do
    
    where ( xz_PotTemp < DelMax * 1.0d-4)
      xz_PotTemp = 0.0d0 
    end where
    
  case ("Exner-Gauss")
    ! $B;XDj$5$l$?>l=j$K(B, $B%,%&%7%"%s$J290L$N>qMp$rM?$($k(B. 

    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    do i = DimXMin, DimXMax
      xz_Exner(i,:) = &
        & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 ) 
    end do
    
    do k = DimZMin, DimZMax
      xz_Exner(:,k) =  &
        & xz_Exner(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
    end do
    
    where ( xz_Exner < DelMax * 1.0d-4)
      xz_Exner = 0.0d0 
    end where
    
  case ("Thermal-Random")
    ! $B:G2<AX$K%i%s%@%`$J>qMp$rM?$($k(B

    open(20,file=RandomFile)
    do i = DimXMin, DimXMax
      read(20,*) random
      xz_Exner(i, RegZMin+1) = DelMax * random
    end do
    close(20)
    
  case ("HS2001")
    ! Hueso and Sanchez-Lavega $B$rLO$7$?@_Dj(B

    i = ( DimXMax - DimXMin - 10) * 5.0d-1
    k = minloc( s_Z, 1, s_Z > 2.5d4 ) - MarginZ
    n = 5.0d3 / DelZ 
    write(*,*) 'DisturbEnv: ', i, k, n

    xz_PotTemp(i-n:i,k-n:k) = DelMax
    
  end select
  !-------------------------------------------------------------
  ! $B>x5$05$N@_Dj(B. 
  !-------------------------------------------------------------
  if (Humidity /= 0.0d0) then 
    do i = DimXMin, DimXMax      
      call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, xz_TempBasicZ(i,:), &
        &              xz_PressBasicZ(i,:), xz_MolFr(i,:,:) )
    end do
    
    !$B5$Aj$N%b%kHf$r:.9gHf$KJQ49(B
    do s = 1, SpcNum
      xza_MixRt(:,:,s) =                                                    &
        & xz_MolFr(:,:,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(:,:,s)
    end do
  end if
  
  
  !$BCM$,>.$5$/$J$j$9$.$J$$$h$&$K:GDcCM$rM?$($k(B
!  where (xza_MixRt <= 1.0d-20 )
!    xza_MixRt = 1.0d-20
!  end where
  
  !$B6-3&>r7o(B
  xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
  xza_MixRt = xza_BoundaryZSym_xza( xza_MixRt )
  
end subroutine DisturbEnv
