!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module DisturbSet_calc
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: disturbset.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!$B>qMp>l$N@_Dj$rFI$_9~$_(B, $B>qMp>l$NCM$r7W;;$9$k(B. 
!
!$BE57?E*$J>qMp>l$H$7$F0J2<$,@_Dj2DG=(B. 
!   * $BL5<!8505NO(B: $B>qMp$J$7(B or X $BJ}8~%,%&%7%"%sJ,I[(B or Z $BJ}8~%,%&%7%"%sJ,I[(B or XZ $BJ}8~%,%&%7%"%sJ,I[(B 
!   * $B290L(B: $B>qMp$J$7(B or X $BJ}8~%,%&%7%"%sJ,I[(B or Z $BJ}8~%,%&%7%"%sJ,I[(B or XZ $BJ}8~%,%&%7%"%sJ,I[(B 
!   * $BB.EY(B: $B>qMp$J$7(B or X $BJ}8~B.EY0lDj(B or Z $BJ}8~B.EY0lDj(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
!
                                                                 !=end

module disturbset
  !== $B0EL[$N7?@k8@6X;_(B
  implicit none

  !== save $BB0@-(B
  save
                                                                 !=begin
  !== Public Interface
  real(8), allocatable  :: ss_PotTempDisturb(:, :) !$B290L(B
  real(8), allocatable  :: ss_ExnerDisturb(:, :)   !$B%(%/%9%J!<4X?t(B
  real(8), allocatable  :: fs_VelXDisturb(:, :)    !$BB.EY(B u
  real(8), allocatable  :: sf_VelZDisturb(:, :)    !$BB.EY(B w
  real(8), allocatable  :: ss_MixRtVDisturb(:, :) !$B?e>x5$:.9gHf(B
  real(8), allocatable  :: ss_MixRtCDisturb(:, :) !$B1@?e:.9gHf(B
  real(8), allocatable  :: ss_MixRtRDisturb(:, :)  !$B1+?e:.9gHf(B
                                                                 !=end

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !NAMELIST $B$+$i>qMp>l$N=i4|CM@8@.%W%m%0%i%`(B disturb $B$N@_Dj$rFI$_<h$k(B. 
  !
  subroutine disturbset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: s_X, s_Z, f_Z, Xmin, Xmax, Zmin, Zmax, &
         &                DimXMin, DimXMax, DimZMin, DimZMax
    use physset,    only: Pi
    
    !=== Input
    character(*), intent(in) :: cfgfile
                                                                 !=end
    !=== Work
    character(80) :: ExnerDisturbType !$BL5<!8505NO$N>qMp(B (DisturbType.f90 $B;2>H(B)
    real(8)       :: ExnerSigmaX      !$B>qMp$NI}(B (X $BJ}8~(B)
    real(8)       :: ExnerSigmaZ      !$B>qMp$NI}(B (Z $BJ}8~(B)
    real(8)       :: ExnerMax         !$B>qMp$N:GBgCM(B
    real(8)       :: ExnerCenterX     !$B>qMp$NCf?4(B (X $BJ}8~(B)
    real(8)       :: ExnerCenterZ     !$B>qMp$NCf?4(B (X $BJ}8~(B)
    
    character(80) :: PotTempDisturbType  !$B290L$N>qMp(B (DisturbType.f90 $B;2>H(B)
    real(8)       :: PotTempSigmaX       !$B>qMp$NI}(B (X $BJ}8~(B)
    real(8)       :: PotTempSigmaZ       !$B>qMp$NI}(B (X $BJ}8~(B)
    real(8)       :: PotTempMax          !$B>qMp$N:GBgCM(B
    real(8)       :: PotTempCenterX      !$B>qMp$NCf?4(B (X $BJ}8~(B)
    real(8)       :: PotTempCenterZ      !$B>qMp$NCf?4(B (X $BJ}8~(B)
    
    character(80) :: VelDisturbType   !$BB.EY$N>qMp(B (DisturbType.f90 $B;2>H(B)
    real(8)       :: VelX             !$BB.EY(B u 
    real(8)       :: VelZ             !$BB.EY(B w

    character(80) :: MixRtDisturbType   !$B?e>x5$NL$N>qMp(B (DisturbType.f90 $B;2>H(B)
    real(8)       :: DeltaMixRtV

    integer       :: i, k
    real(8)       :: mu
    real(8)       :: sigma

    real(8)       :: beta
    real(8)       :: Xr, Xc, Zr, Zc
    real(8)       :: DeltaPotTemp


                                                                 !=end
                                                                 !=begin
    !=== NAMELIST
    NAMELIST /disturbset/ &
         & ExnerDisturbType, ExnerMax, ExnerSigmaX, ExnerSigmaZ, &
         &                  ExnerCenterX, ExnerCenterZ, &
         & PotTempDisturbType, PotTempMax, PotTempSigmaX, PotTempSigmaZ, &
         &                  PotTempCenterX, PotTempCenterZ, &
         & VelDisturbType, VelX, VelZ, &
         & MixRtDisturbType
                                                                 !=end    
    call BeginSub("disturbset_calc_init", &
&                 fmt="%c",       &
&                 c1="Initialize values of disturbance. ")

    open (10, FILE=cfgfile)
    read(10, NML=disturbset)
    close(10)

    !==== $BJQ?t$NCM$N%A%'%C%/(B

    if ( ExnerSigmaX > 1 )  then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaX is > 1")
       stop
    end if

    if ( ExnerSigmaZ > 1 )  then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaZ is > 1")
       stop
    end if

    if ( PotTempSigmaX > 1 )  then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaX is > 1")
       stop
    end if

    if ( PotTempSigmaZ > 1 )  then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaZ is > 1")
       stop
    end if

    if ( ExnerCenterX > 1 ) then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterX is > 1")
       stop
    end if

    if ( ExnerCenterZ > 1 ) then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterZ is > 1")
       stop
    end if

    if ( PotTempCenterX > 1 ) then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterX is > 1")
       stop
    end if

    if ( PotTempCenterZ > 1 ) then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterZ is > 1")
       stop
    end if


    !==== $BG[Ns$N3d$jEv$F(B
    allocate( ss_PotTempDisturb(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_ExnerDisturb(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    fs_VelXDisturb(DimXMin:DimXMax, DimZMin:DimZMax),    &
         &    sf_VelZDisturb(DimXMin:DimXMax, DimZMin:DimZMax),    &  
         &    ss_MixRtVDisturb(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_MixRtCDisturb(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_MixRtRDisturb(DimXMin:DimXMax, DimZMin:DimZMax)  )
    !==== $B=i4|2=(B
    ss_ExnerDisturb = 0.0d0
    ss_PotTempDisturb = 0.0d0
    fs_VelXDisturb = 0.0d0
    sf_VelZDisturb = 0.0d0
    ss_MixRtVDisturb = 0.0d0
    ss_MixRtCDisturb = 0.0d0
    ss_MixRtRDisturb  = 0.0d0
    
    
    !-----------------------------------------------------------
    !--- $B%(%/%9%J!<4X?t$N>qMp(B
    !-----------------------------------------------------------
    select case (ExnerDisturbType)
       
    case ("Exner_Gauss_X")
       !--- X $BJ}8~$N%,%&%7%"%sJ,I[(B
       mu = ( Xmax - Xmin ) * ExnerCenterX
       sigma = ( Xmax - Xmin ) * ExnerSigmaX
       
       do i = DimXMin, DimXMax
          ss_ExnerDisturb(i,:) = &
               & ExnerMax * dexp(-((s_X(i) - mu) / sigma )**2.0d0 / 2.0d0) 
       end do
       
    case ("Exner_Gauss_Z")    
       !--- Z $BJ}8~$N%,%&%7%"%sJ,I[(B
       mu = ( Zmax - Zmin ) * ExnerCenterZ
       sigma = ( Zmax - Zmin ) * ExnerSigmaZ
       
       do k = DimZMin, DimZMax
          ss_ExnerDisturb(:,k) = &
               & ExnerMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("Exner_Gauss_XZ")
       !--- $BF1?41_>u$N%,%&%7%"%sJ,I[(B
       mu = ( Xmax - Xmin ) * ExnerCenterX
      sigma = ( Xmax - Xmin ) * ExnerSigmaX
       
       do i = DimXMin, DimXMax
          ss_ExnerDisturb(i,:) = &
               & ExnerMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
       mu = ( Zmax - Zmin ) * ExnerCenterZ
       sigma = ( Zmax - Zmin ) * ExnerSigmaZ
       
       do k = DimZMin, DimZMax
          ss_ExnerDisturb(:,k) = ss_ExnerDisturb(:,k) &
               & * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do

    end select

    !-----------------------------------------------------------
    !--- $B2>290L$N>qMp(B
    !-----------------------------------------------------------
    select case (PotTempDisturbType)
       
    case ("PotTemp_Gauss_X")
       !--- X $BJ}8~$N%,%&%7%"%sJ,I[(B
       mu = ( Xmax - Xmin ) * PotTempCenterX
       sigma = ( Xmax - Xmin ) * PotTempSigmaX
       
       do i = DimXMin, DimXMax
          ss_PotTempDisturb(i,:) = &
               & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("PotTemp_Gauss_Z")
       !--- Z $BJ}8~$N%,%&%7%"%sJ,I[(B
       mu = ( Zmax - Zmin ) * PotTempCenterX
       sigma = ( Zmax - Zmin ) * PotTempSigmaZ
       
       do k = DimZMin, DimZMax
          ss_PotTempDisturb(:,k) = &
               & PotTempMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("PotTemp_Gauss_XZ")     
       !--- $BF1?41_>u$N%,%&%7%"%sJ,I[(B
       mu = ( Xmax - Xmin ) * PotTempCenterX
       sigma = ( Xmax - Xmin ) * PotTempSigmaX
       
       do i = DimXMin, DimXMax
          ss_PotTempDisturb(i,:) = &
               & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
       mu = ( Zmax - Zmin ) * PotTempCenterZ
       sigma = ( Zmax - Zmin ) * PotTempSigmaZ
       
       do k = DimZMin, DimZMax
          ss_PotTempDisturb(:,k) =  &
               & ss_PotTempDisturb(:,k) &
               &   * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
!       where (ss_PotTempDisturb < PotTempMax * 1.0d-4)
!          ss_PotTempDisturb = 0.0d0 
!       end where
       
    case default
       
       ss_PotTempDisturb = 0.0d0
       
    end select

     
    !----------------------------------------------------------
    !--- $BB.EY(B
    !----------------------------------------------------------
    select case (VelDisturbType)
       
    case ("Zero")
       fs_VelXDisturb = 0.0d0
       sf_VelZDisturb = 0.0d0
       
    case ("VelX_test")
       fs_VelXDisturb = VelX
       fs_VelXDisturb(-5:10, :) = 0.0d0
       fs_VelXDisturb(90:105, :) = 0.0d0
       sf_VelZDisturb = 0.0d0
       
    case ("VelX_Uniform")
       !--- $BB.EY(B u $B0lDj(B
       fs_VelXDisturb = VelX
       sf_VelZDisturb = 0.0d0
       
    case ("VelZ_Uniform")
       !--- $BB.EY(B w $B0lDj(B
       fs_VelXDisturb = 0.0d0
       sf_VelZDisturb = VelZ
       
    case ("Solid_Rotation")
       !--- $BNN0hCf?4$N9dBN2sE>(B
       
       do i = DimXMin, DimXMax   
          fs_VelXDisturb(i,:) = &
               & 2.0d0*pi / 1000.0d0 * ( f_Z(:) - (Zmax + Zmin) / 2.0D0 )
       end do
       
       do k = DimZMin, DimZMax
          sf_VelZDisturb(:,k) = - 2.0d0 * pi / 1000.0d0 &
               &                 * ( s_X(:) - (Xmax + Xmin) / 2.0D0 )
       end do
       
    end select
    

!!!
!!! Klemp and Wilhelmson (1978) $B$N@_Dj(B (2005-09-26 $BDI2C(B)
!!!
    select case (PotTempDisturbType)
       case ( "KW78" )
          Xr = 1.0d4     
          Xc = 1.2d4
          Zr = 2.0d3
          Zc = 1.5d3
          DeltaPotTemp = 1.5d0

          do i = DimXMin, DimXMax
             do k = DimZMin, DimZMax          
                beta =  &
                     & ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
                     &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1 
                
                if ( beta < 1.0d0 ) then 
                   ss_PotTempDisturb(i,k) = &
                        & DeltaPotTemp * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
                end if

             end do
          end do
          
       end select

!!!
!!!$B?e>x5$:.9gHf$KBP$9$k@_Dj(B (2005-10-05)
!!!  $B8=:_$O2?$b$7$J$$@_Dj$K$J$C$F$$$k(B
!!!
    select case (MixRtDisturbType)
       case ( "KW78" )
          Xr = 1.0d4     
          Xc = 1.2d4
          Zr = 2.0d3
          Zc = 1.5d3
          DeltaMixRtV = 1.0d-4
          
          do i = DimXMin, DimXMax
             do k = DimZMin, DimZMax          
                beta =  &
                     & ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
                     &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1 
                
                if ( beta < 1.0d0 ) then 
                   ss_MixRtVDisturb(i,k) = &
                        & DeltaMixRtV * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
                end if
                
             end do
          end do
          
       end select
       
       
    call EndSub("disturbset_calc_init")
    
  end subroutine disturbset_calc_init
  
end module disturbset
