!=================================================
! 3D cumulus model - kaminari
!
! FileName : base.f90
! Note     : $B4pK\>l7W;;%W%m%0%i%`(B
!            $BC10L$O(B MKS $B7O(B
!            $BL5<!8505NO(B = $B%(%/%9%J!<4X?t(B
! Author   : TAKAHASHI Koko
! Date     : 2003/10/21
!            2003/09/01 $B?75,:n@.(B
! 
!=================================================
program base

  !-------------------------------------------------
  ! $BJQ?t$J$I$N@_Dj(B
  !-------------------------------------------------

  use gtool_history
  implicit none

  ! $B6u4V2rA|EY@_Dj(B
  integer, parameter :: im = 5         ! $B3J;RE@$N@_Dj(B (x)
  integer, parameter :: jm = 5         ! $B3J;RE@$N@_Dj(B (y)
  integer, parameter :: km = 5         ! $B3J;RE@$N@_Dj(B (z)

  ! $B:BI8JQ?t(B
  real, parameter :: xmin = .0
  real, parameter :: xmax = 1e2
  real, parameter :: dx = (xmax-xmin)/im ! $B3J;R4V3V(B (x)
  real, parameter :: ymin = .0
  real, parameter :: ymax = 1e2
  real, parameter :: dy = (ymax-ymin)/jm ! $B3J;R4V3V(B (y)
  real            :: zs(0:im,0:jm)       ! $BCO7A(B (z_s)
  real, parameter :: zmin = .0            ! $B%b%G%kNN0h$NDl(B
  real, parameter :: zmax = 1e2         ! $B%b%G%kNN0h$N9b$5(B (H)
  real, parameter :: dz = (zmax-zmin)/km ! $B3J;R4V3V(B (y)
!  real            :: dz(0:im,0:jm,0:km)  ! $B3J;R4V3V(B (z^{*})

  !--- $B4pK\>l$NJQ?t(B
  real :: bpi(0:im,0:jm,0:km)            ! $BL5<!8505NO(B
  real :: bpr(0:im,0:jm,0:km)            ! $B05NO(B
  real :: brh(0:im,0:jm,0:km)            ! $BL)EY(B
  real :: btm(0:im,0:jm,0:km)            ! $B29EY(B
  real :: bth(0:im,0:jm,0:km)            ! $B290L(B
  real :: bvt(0:im,0:jm,0:km)            ! $B2>290L(B
  real :: bqv(0:im,0:jm,0:km)            ! $B?e>x5$$N:.9gHf(B
  real :: bqc(0:im,0:jm,0:km)            ! $B1@?e$N:.9gHf(B
  real :: bqr(0:im,0:jm,0:km)            ! $B1+?e$N:.9gHf(B

  !--- $B$=$NB>$NJQ?t(B
  real :: g12(0:im,0:jm,0:km)           ! $B%d%3%S%"%s$N(B1/2$B>h(B

  real :: x(0:im,0:jm,0:km)                        ! x $B:BI8(B
  real :: y(0:im,0:jm,0:km)                        ! y $B:BI8(B
  real :: z(0:im,0:jm,0:km)                        ! z^{*} $B:BI8(B
  real :: zn(0:im,0:jm,0:km)                       ! z $B:BI8(B

  !--- $BJ*M}%Q%i%a!<%?(B
  real, parameter :: bpr0 = 1e3        ! $B4pK\>l$NCOI=LL$G$N05NO(B
  real, parameter :: rd = 287           ! $B4%Ag6u5$$N5$BNDj?t(B
  real, parameter :: cp = 1952          ! $BDj05HfG.(B
  real, parameter :: cv = 1463          ! $BDj@QHfG.(B
  real, parameter :: grv = 9.8          ! $B=ENO(B

!  real, parameter :: rx = 50
!  real, parameter :: ry = 50

  !--- $B$=$NB>$NJQ?t(B         
  integer            :: i, j, k                   ! do $BJQ?t(B

  !-------------------------------------------------
  ! x-y $B:BI8$N@_Dj(B
  !-------------------------------------------------

  do i=2,im-1
     do j=2,jm-1
        do k=2,km-1
           x(i,j,k)= dx*(i-2)
           y(i,j,k)= dy*(j-2)
           zs(i,j) = 0
!          zs(i,j) = exp(-(x(i,j,k) - xmax/5)**2/(rx**2)    &
!                  &     - (y(i,j,k) - ymax/2)**2/(ry**2))
           if(k == 2) zn(i,j,k) = zs(i,j)
           if(k > 2) zn(i,j,k)= (zmax-zs(i,j))/km*(k-2)
           z(i,j,k) = zmax*(zn(i,j,k) - zs(i,j))            &
                    & /(zmax - zs(i,j))
        end do
     end do
  end do

  call output_gtool4_init
  call output_gtool4

  !-------------------------------------------------
  ! $B4pK\>l$N=i4|CM@_Dj(B
  !-------------------------------------------------

  !--- $B4pK\>l$G$NCOI=LL$G$NCM(B

  !--- $B4%Ag6u5$$N>l9g$N?e>x5$$N:.9gHf(B
!  bqv(0:im,0:jm,0:km) = .0
!  bqc(0:im,0:jm,0:km) = .0
!  bqr(0:im,0:jm,0:km) = .0

  !-------------------------------------------------
  ! $B4pK\>l$N7W;;(B
  !-------------------------------------------------

  !-- $B4pK\>l$N%G!<%?=q$-9~$_MQ%U%!%$%k$r3+$/(B

  do i = 2,im-1
     do j = 2,jm-1
        do k = 2,km-1
           !--- $B29EY(B: $BNc$($P(B
           btm(i,j,k) = 300 - 6*z(i,j,k)/bpr0
           !--- $B%U%!%$%k$+$iFI$_9~$`$J$i(B

           !--- $B%d%3%S%"%s$N(B1/2$B>h(B
           g12(i,j,k) = (zmax - zs(i,j))/zmax

           !--- $B<>=a6u5$$N>l9g$N?e>x5$$N:.9gHf(B
           bqv(i,j,k) = .0
           bqc(i,j,k) = .0
           bqr(i,j,k) = .0

           !--- $B05NO(B
           if(k==2) bpr(i,j,k) = bpr0
           if (k>2) bpr(i,j,k) = bpr0*bpi(i,j,k)**(cp/rd)

           !--- $B290L(B
           bth(i,j,k) = btm(i,j,k)/bpi(i,j,k)

           !--- $B2>290L(B
           bvt(i,j,k) = bth(i,j,k)*(1 + 0.61*bqv(i,j,k))

           !--- $BL)EY(B
           brh(i,j,k) = (bpr0/rd)                           &
                      & *(bpi(i,j,k)**(cv/cp)/bvt(i,j,k))

           !--- $B@E?e05J?9U$N<0(B -> z+dz $B$G$NL5<!8505NO(B
           if (k==2) bpi(i,j,k) = 1.
           if (k>2) bpi(i,j,k) = bpi(i,j,k-1)               &
                               & - 2*dz*g12(i,j,k)          &
                               &   *grv/(cp*bvt(i,j,k))

           !--- $B=PNO(B
           write(*,*) brh(i,j,k)
           call output_gtool4

!           dz(i,j,k) = z(i,j,k) - z(i,j,k-1)

        end do
     end do
  end do

  call output_gtool4_close
  stop

contains

!-----------------------------------------------------------
!--- $B%5%V%k!<%A%s(B

  !--- gtool4: $B%R%9%H%j!<=i4|2=(B
  subroutine output_gtool4_init
    call HistoryCreate(                                     & ! $B%R%9%H%j!<:n@.(B
          & file='base.nc',                                 &
          & title='3D cumulus model  base',                 &
          & source='kaminari - 3D cumulus model base',      &
          & institution='GFD_Dennou Club deepconv project', &
          & dims=(/'x','y','z'/), dimsizes=(/im,jm,km/),    &
          & longnames=(/'X-coordinate','Y-coordinate', 'Z-coordinate'/), &
          & units=(/'m','m','m'/),                          &
          & origin=0.0, interval=real(dx)                         &
          & )

    call HistoryPut('x',x(2:im-1,2,2))                        ! $BJQ?t=PNO(B
    call HistoryPut('y',y(2,2:jm-1,2))                        ! $BJQ?t=PNO(B
    call HistoryPut('z',z(2,2,2:km-1))                        ! $BJQ?t=PNO(B

    call HistoryAddVariable(                                & ! $BJQ?tDj5A(B
         & varname='bpr', dims=(/'x','y','z'/),             & ! $B05NO(B
         & longname='pressure', units='hPa', xtype='double')
    call HistoryAddVariable(                                &
         & varname='btm', dims=(/'x','y','z'/),             & ! $B29EY(B
         & longname='temperature', units='k', xtype='double')
    call HistoryAddVariable(                                &
         & varname='brh', dims=(/'x','y','z'/),             & ! $BL)EY(B
         & longname='density', units='kg/m3',           &
         & xtype='double')
    call HistoryAddVariable(                                &
         & varname='bqv', dims=(/'x','y','z'/),             & ! $B?e>x5$:.9gHf(B
         & longname='mixing ratio of water vapor',          &
         &  units='g/g', xtype='double')
    call HistoryAddVariable(                                &
         & varname='bqc', dims=(/'x','y','z'/),             & ! $B1@?e:.9gHf(B
         & longname='mixing ratio of cloud drops',          &
         & units='g/g', xtype='double')
    call HistoryAddVariable(                                &
         & varname='bqr', dims=(/'x','y','z'/),             & ! $B1+?e:.9gHf(B
         & longname='mixing ratio of rainwater',            &
         & units='g/g', xtype='double')
  end subroutine output_gtool4_init

  !--- gtool4: $B%U%!%$%k=PNO(B
  subroutine output_gtool4
    write(6,*) 'i,j,k = ', i,j,k
    call HistoryPut('bpr',real(bpr(i,j,k)))
    call HistoryPut('btm',real(btm(i,j,k)))
    call HistoryPut('brh',real(brh(i,j,k)))
    call HistoryPut('bqv',real(bqv(i,j,k)))
    call HistoryPut('bqc',real(bqc(i,j,k)))
    call HistoryPut('bqr',real(bqr(i,j,k)))
  end subroutine output_gtool4

  !--- gtool4: $B=*N;=hM}(B
  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program base
