!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module BasicSet
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: basicset.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!$B4pK\>l$NCM$r<hF@$9$k(B. 
!   * basicset_file_init: $B4pK\>l$NCM$r(B netCDF $B%U%!%$%k$+$i<hF@(B
!   * basicset_calc_init: $B4pK\>l$N>pJs$r(B Namelist $B$+$i<hF@$7$FCM$r7W;;(B
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
                                                                 !=end

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

  !== save $BB0@-(B
  save
                                                                 !=begin
  !== Public Interface
  real(8), allocatable   :: ss_CpBasicZ(:,:)       !$BHfG.(B
  real(8), allocatable   :: ss_DensBasicZ(:,:)     !$BL)EY(B
  real(8), allocatable   :: ss_PressBasicZ(:,:)    !$BL5<!8505NO(B
  real(8), allocatable   :: ss_ExnerBasicZ(:,:)    !$BL5<!8505NO(B
  real(8), allocatable   :: ss_TempBasicZ(:,:)     !$B29EY(B
  real(8), allocatable   :: ss_PotTempBasicZ(:,:)  !$B290L(B
  real(8), allocatable   :: ss_VelSoundBasicZ(:,:) !$B2;B.(B
  real(8), allocatable   :: ss_MixRtVBasicZ(:,:) !$B?e>x5$:.9gHf(B
  real(8), allocatable   :: ss_MixRtCBasicZ(:,:) !$B1@N3:.9gHf(B
  real(8), allocatable   :: ss_MixRtRBasicZ(:,:)  !$B1+?e:.9gHf(B
                                                                 !=end

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module
  !
  !$B4pK\>l$NCM$r(B netCDF $B%U%!%$%k$+$i<hF@$9$k(B. 
  !netCDF $B$+$iCM$r<hF@$9$k$?$a$K(B gt4f90io $B$rMxMQ$9$k(B. 
  !
  subroutine BasicSet_file_init(basicfile)

    !=== Dependency
    use dc_trace, only: BeginSub, EndSub
    use gt4_history
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax

    !=== Input
    character(*), intent(in) :: basicfile
                                                                 !=end
    
    character(20)  :: name            !$BJQ?tL>(B

    call BeginSub("basicset_file_init", &
&                 fmt="%c",       &
&                 c1="Get basic state variables from netCDF file.")
    
    !--- $BG[Ns$N3d$jEv$F(B
    allocate(  &
         & ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtVBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtCBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtRBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)  )
    
                                                                 !=begin
    !=== Get a Value from netCDF File
    name = "CpBasicZ"
    call HistoryGet( basicfile, name, ss_CpBasicZ )

    name = "PressBasicZ"
    call HistoryGet( basicfile, name, ss_PressBasicZ )

    name = "ExnerBasicZ"
    call HistoryGet( basicfile, name, ss_ExnerBasicZ )

    name = "TempBasicZ"
    call HistoryGet( basicfile, name, ss_TempBasicZ )

    name = "PotTempBasicZ"
    call HistoryGet( basicfile, name, ss_PotTempBasicZ )

    name = "DensBasicZ"
    call HistoryGet( basicfile, name, ss_DensBasicZ )

    name = "VelSoundBasicZ"
    call HistoryGet( basicfile, name, ss_VelSoundBasicZ ) 

    name = "MixRtVBasicZ"
    call HistoryGet( basicfile, name, ss_MixRtVBasicZ ) 

    name = "MixRtCBasicZ"
    call HistoryGet( basicfile, name, ss_MixRtCBasicZ ) 

    name = "MixRtRBasicZ"
    call HistoryGet( basicfile, name, ss_MixRtRBasicZ ) 
                                                                 !=end   

    call EndSub("basicset_file_init")

  end subroutine BasicSet_file_init


                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !$BE57?E*$J4pK\>l$r7W;;$9$k$?$a$N%W%m%0%i%`(B basic $B$N=i4|CM$r(B
  !$BFI$_9~$_(B, $B4pK\>l$NCM$r7W;;$9$k$?$a$N%b%8%e!<%k(B
  !
  !$BE57?E*$J4pK\>l$N$r7W;;$9$k$?$a$N%5%V%k!<%A%s(B. 
  !$B%?%$%W$H$7$F$O0J2<$,@_Dj2DG=(B. 
  !   * $B29EY0lMM(B or $B290L0lMM(B
  !   * $B05NO0lMM(B or $B@E?e05J?9U(B
  !
  subroutine basicset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: DimXMin, DimXMax, DimZMin, DimZMax, &
         &                s_Z, RegXMin, RegZMin
    use physset,    only: Grav,     &! $B=ENO(B     [m/s]
         &                TempSfc,  &! $BI8=`29EY(B [K]
         &                PressSfc, &! $BI8=`05NO(B [Pa]
         &                GasRDry,  &! $B4%Ag@.J,$N5$BNDj?t(B [J/K kg]
         &                CpDry,    &! $B4%Ag@.J,$NJ?6QDj05HfG.(B [J/K kg]
         &                CvDry,    &! $B4%Ag@.J,$NJ?6QDj@QHfG.(B [J/K kg]
         &                epsilon    ! $B4%Ag@.J,$H<>=a@.J,$NJ,;RNLHf(B
    
    !=== Input
    character(*), intent(in) :: cfgfile  !$B@_Dj%U%!%$%k(B(NAMELIST)
                                                                 !=end
    !=== Work
    real(8)       :: LapseRate
    character(80) :: VaperType
    character(80) :: TempType
    character(80) :: PressType
    real(8)  :: ss_CvBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_Z(DimXMin:DimXMax, DimZMin:DimZMax)          ! 2-D $B:BI8(B
    real(8)  :: MixRtVBase
    real(8)  :: HumVaper
    real(8)  :: MolFrIni
    integer  :: k
                                                                 !=begin
    !=== NAMELIST 
    NAMELIST /basicset/ TempType, LapseRate, PressType, VaperType, HumVaper, MolFrIni
                                                                 !=end

    call BeginSub("basicset_calc_init", &
&                 fmt="%c",       &
&                 c1="calculate the basic state variables.")

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

    !=== $BJQ?t$N=i4|2=(B
    allocate( ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
         &    ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    ss_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_MixRtVBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_MixRtCBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_MixRtRBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)  )

    ss_CpBasicZ = 0.0d0
    ss_CvBasicZ = 0.0d0
    ss_DensBasicZ = 0.0d0
    ss_PressBasicZ = 0.0d0
    ss_ExnerBasicZ = 0.0d0
    ss_TempBasicZ = 0.0d0
    ss_PotTempBasicZ = 0.0d0
    ss_VelSoundBasicZ = 0.0d0
    ss_MixRtVBasicZ = 0.0d0
    ss_MixRtCBasicZ = 0.0d0
    ss_MixRtRBasicZ = 0.0d0
    
    !--- $B:BI8$r=i4|2=(B
    do k = DimZMin, DimZMax
       ss_Z(:,k) = s_Z(k)
    end do
    
    !--------------------------------------------
    !--- $BHfG.$r7h$a$k(B
    !---   $BHfG.$O4%AgBg5$$rG0F,$K0lMM$H$9$k(B. 
    !--------------------------------------------
    ss_CpBasicZ = CpDry
    ss_CvBasicZ = CvDry
  
    !--------------------------------------------
    !--- $B29EY(B/$B290L(B, $BL5<!8505NO$N4pK\>l(B
    !--------------------------------------------
    if (TempType == "uniform" .AND. PressType == "uniform") then 
       
       !--- $B29EY>l$O(B TempSfc $B$G0lMM(B
       ss_TempBasicZ = TempSfc
       
       !--- $BCOI=LL05NO$G0lDj(B. $BDj5A$h$jCM$O(B 1 $B$H$J$k(B. 
       ss_ExnerBasicZ = 1.0d0
       
       !--- $B290L(B
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ 
       
    elseif (TempType == "uniform" .AND. PressType == "hydrostatic") then 
       
       !--- $B29EY>l$O(B TempSfc $B$G0lMM(B
       ss_TempBasicZ = TempSfc
       
       !--- $BL5<!8505NO$H290L$rO"N)$7$F2r$/(B. 
       ss_ExnerBasicZ = &
            & dexp( - Grav * ss_Z / ( ss_CpBasicZ * ss_TempBasicZ ) )

       !--- $BL5<!8505NO$rMxMQ$7$F290L$r7W;;(B
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ
       
    elseif (TempType == "adiabat" .AND. PressType == "uniform") then 
       
       !--- $BCOI=LL$G$N290L(B(= TempSfc)$B$G0lDj(B
       ss_PotTempBasicZ = TempSfc
       
       !--- $BCOI=LL05NO$G0lDj(B. $BDj5A$h$jCM$O(B 1 $B$H$J$k(B. 
       ss_ExnerBasicZ = 1.0d0
       
       !--- $B29EY$N7W;;(B
       ss_TempBasicZ = ss_PotTempBasicZ * ss_ExnerBasicZ

    elseif (TempType == "adiabat" .AND. PressType == "hydrostatic") then 
       
       !--- $BCOI=LL$G$N290L(B(= TempSfc)$B$G0lDj(B
       ss_PotTempBasicZ = TempSfc
       
       !--- $BL5<!8505NO$r7W;;(B. $B290L0lDj(B
       ss_ExnerBasicZ = &
            & - Grav * ss_Z / ( ss_CpBasicZ * ss_PotTempBasicZ ) + 1.0d0

       !--- $B29EY$N7W;;(B
       ss_TempBasicZ = ss_PotTempBasicZ * ss_ExnerBasicZ
       
    elseif (TempType == "decrease" .AND. PressType == "hydrostatic") then 
       
       !--- $B29EY>l$O(B TempSfc $B$G0lMM(B
       ss_TempBasicZ = TempSfc - LapseRate * ss_Z 
       
       !--- $BL5<!8505NO$H290L$rO"N)$7$F2r$/(B. 
       ss_ExnerBasicZ = &
            & dexp( - Grav * ss_Z / ( ss_CpBasicZ * ss_TempBasicZ ) )

       !--- $BL5<!8505NO$rMxMQ$7$F290L$r7W;;(B
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ
       
    else
       call MessageNotify("Message", "basicset_calc_init", &
&                         "unknown TempType, or PressType")
    end if
    
    !--------------------------------------------
    !---- $BL)EY(B
    !--------------------------------------------
    ss_DensBasicZ = &
         & PressSfc * ( ss_ExnerBasicZ ** ( ss_CvBasicZ / GasRDry ) ) &
         &   / ( GasRDry * ss_PotTempBasicZ )
    
    !--------------------------------------------
    !---- $B05NO(B
    !--------------------------------------------
    ss_PressBasicZ = PressSfc * ( ss_ExnerBasicZ ** ( ss_CpBasicZ / GasRDry ) )
    
    !--------------------------------------------    
    !--- $BJ?6Q>l$N2;B.(B
    !--------------------------------------------
    ss_VelSoundBasicZ = &
         & ss_CpBasicZ * GasRDry * ss_ExnerBasicZ * ss_PotTempBasicZ  &
         &   / ss_CvBasicZ
    ss_VelSoundBasicZ = sqrt( ss_VelSoundBasicZ )


    !-------------------------------------------
    !---- $B?e>x5$:.9gHf(B
    !------ $B7h$^$C$?29EY$KBP$9$kK0OB?e>x5$:.9gHf$N<>EY(B HumVaper $B$K$9$k(B
    !------ $BK\Mh$O4X?t$r8F$S$?$$$H$3$m$@$,(B, $B%3%s%Q%$%k$N=gHV$,LdBj$K$J$k(B 
    !-------------------------------------------
    if ( VaperType == "decrease" ) then 
       MixRtVBase = &
            &   epsilon * 610.78d0 &
            & * dexp( 17.269d0 * ( ss_TempBasicZ(RegXMin,RegZMin) - 273.16d0 )   &
            &                  / ( ss_TempBasicZ(RegXMin,RegZMin) - 35.86d0  ) ) &
            & / ss_PressBasicZ(RegXMin,RegZMin)
       
       ss_MixRtVBasicZ = MixRtVBase * ( 1.0d0 - 5.0d-5 * ss_Z )
       
    else
       ss_MixRtVBasicZ = &
            & HumVaper &
            & * epsilon * 610.78d0 &
            & * dexp( 17.269d0 * ( ss_TempBasicZ - 273.16d0 )   &
            &                  / ( ss_TempBasicZ - 35.86d0  ) ) &
            & / ss_PressBasicZ
    end if




!!!
!!! ECCM $B$+$i$D$J$0>l9g(B
!!!
    if (         TempType  == "eccm" &
         & .AND. PressType == "eccm" &
         & .AND. VaperType == "eccm"   ) then 

       !--------------------------------------------
       !---- $B29EY(B, $B05NO(B, $B?e>x5$:.9gHf(B
       !--------------------------------------------
       do k = DimXMin, DimXMax
          call eccm( MolFrIni, HumVaper, &
               & ss_TempBasicZ(k,:), ss_PressBasicZ(k,:), ss_MixRtVBasicZ(k,:) )
       end do
       

       !--------------------------------------------
       !---- $B%(%/%9%J!<4X?t(B
       !--------------------------------------------
       ss_ExnerBasicZ = ( ss_PressBasicZ / PressSfc )**( GasRDry / CpDry )


       !--------------------------------------------
       !---- $B290L(B
       !--------------------------------------------
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ


       !--------------------------------------------
       !---- $BL)EY(B
       !--------------------------------------------
       ss_DensBasicZ = &
            & PressSfc * ( ss_ExnerBasicZ ** ( ss_CvBasicZ / GasRDry ) ) &
            &   / ( GasRDry * ss_PotTempBasicZ )
       
       !--------------------------------------------    
       !--- $BJ?6Q>l$N2;B.(B
       !--------------------------------------------
       ss_VelSoundBasicZ = &
            & ss_CpBasicZ * GasRDry * ss_ExnerBasicZ * ss_PotTempBasicZ  &
            &   / ss_CvBasicZ
       ss_VelSoundBasicZ = sqrt( ss_VelSoundBasicZ )
       
    end if





    !==== $B3NG'(B
!    write(*,*) "ss_CpBasicZ ", maxval( ss_CpBasicZ )
!    write(*,*) "ss_DensBasicZ ", maxval( ss_DensBasicZ )
!    write(*,*) "ss_ExnerBasicZ ", maxval( ss_ExnerBasicZ )
!    write(*,*) "ss_PotTempBasicZ ", maxval( ss_PotTempBasicZ )
!    write(*,*) "ss_VelSoundBasicZ ", maxval( ss_VelSoundBasicZ )
    

    call EndSub("basicset_calc_init")

  end subroutine basicset_calc_init
  
end module BasicSet
