C*********************************************************************
C
C     PROGRAM xstress
C
C     地表摩擦水平分布図(時刻指定)
C
C     履歴 1999/08/03 小高正嗣
C          2000/11/21 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      program xstress
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------

#include "grid_size_M.f"

      parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) )
*-----------------------------------------------------
      integer ibr(nbr3)
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
      parameter ( PI = 3.141592D0 ) 

      dimension
     \   zgrid ( -2:nzb ),
     \   temp0 ( -2:nzb ),
     \   tpot0 ( -2:nzb ),
     \   dens0 ( -2:nzb ),
     \   ppai0 ( -2:nzb ),
     \   pres0 ( -2:nzb ),
     \   rvap0 ( -2:nzb ),
     \   qvap0 ( -2:nzb )

      dimension
     \     u    ( -2:nxb, -2:nzb ),
     \     v    ( -2:nxb, -2:nzb ),
     \     w    ( -2:nxb, -2:nzb ),
     \     tpot ( -2:nxb, -2:nzb ),
     \     tempg( -2:nxb, nzg    )

      dimension
     \     uu   ( -2:NXB, -2:NZB, NROT ),
     \     vv   ( -2:NXB, -2:NZB, NROT ),
     \     ttpot( -2:NXB, -2:NZB, NROT )

      DIMENSION
     \     CDRAGV( -2:NXB ),
     \     RICHD( -2:NXB ),
     \     TSFC ( -2:NXB )

      real*4
     \    z0p ( 0:nz ), 
     \    p0p ( 0:nz ),
     \    d0p ( 0:nz ),
     \ temp0p ( 0:nz ),
     \ tpot0p ( 0:nz ),
     \  pai0p ( 0:nz ),
     \ rvap0p ( 0:nz ),
     \ qvap0p ( 0:nz )

      real*4 
     \     us ( nx+1 ),
     \     up ( nx+1,nz+1 ),
     \     vp ( nx+1,nz+1 ),
     \     wp ( nx+1,nz+1 ),
     \     tpotp ( nx+1, nz+1 ),
     \     tpall ( nx+1, nz+1 ),
     \     cdragvp( nx+1 ),
     \     richdp( nx+1 ),
     \     tsfcp ( nx+1 ),
     \     fheat ( nx+1 ),
     \     xplot ( nx+1 ),
     \     stress( nx+1 ),
     \     stress_cr( nx+1 )

      real*4
     \     stress_1( nx+1 ),
     \     stress_2( nx+1 ),
     \     cdragvp_1( nx+1 ),
     \     cdragvp_2( nx+1 ),
     \     uplot(nx+1)

      character*2 mint
      character*2 lt
      character*8 ltime
      character*2 date
      character*14 header

      character*100 dumy
*-----------------------------------------------------
      
      open (11,file='exparam')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read (11,*) ylat, als
      close(11)

      call FOPEN ( irun+1 )

      do iz = 1,9
         read(16,'(a100)') dumy
      end do

      do iz = -2, nzb
         read(16,*) ii, zgrid(iz), pres0(iz), dens0(iz), temp0(iz),
     \                     tpot0(iz), ppai0(iz), qvap0(iz), rvap0(iz)
      end do

      do iz = 0, nz
         z0p(iz)   = zgrid(iz)
         p0p(iz)   = pres0(iz)
         d0p(iz)   = dens0(iz)
         temp0p(iz)= temp0(iz)
         tpot0p(iz)= tpot0(iz)
         pai0p(iz) = ppai0(iz)
         rvap0p(iz)= rvap0(iz)
         qvap0p(iz)= qvap0(iz)
      end do

      do ix = 1, nx+1
         xplot(ix) = 0.1 * (ix - 1)
      end do

      write(6,*) 'data number = ', nloop1
      write(6,*) 'INPUT DATA BEGINING, END, & INTERVAL' 
      write(6,*) 'itb ite interval ?'
      read (*,*) itb, ite, interv

      call sgpwsn
      read(*,*) iws

      CALL SWISET('IWIDTH',  720)
      CALL SWISET('IHEIGHT', 500)

      call gropn(iws)
c      call gropn(-iws)
c      call sldiv( 'T', 2, 6 )
      call swlset( 'LALT',  .TRUE. )
      call sglset( 'LFULL', .TRUE. )
      call sglset( 'LCORNER', .FALSE. )
      call uzfact(0.6)

      do it = 1, ite

      read(31) time, ((u(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(32) time, ((v(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(33) time, ((w(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(40) time, ((tempg(ix,iz),ix=-2,nxb), iz=1,nzg)

      write(0,*) 'reading time = ', time

      if ( it .ge. itb   .and.  it .le. ite ) then
         if ( mod( (it-itb), interv ) .eq. 0 ) then

            write(0,*) 'writing time = ', time
            
            do IX = -2, NXB
               TSFC(IX) = tempg(IX, NZG)
            end do

c*** 時刻変換: 秒から LT へ************************************
            time = time + times
            write(date,'(I2.2)') INT( time / day ) + 1

            time = ( time - INT( time / day ) * day ) / hour 

            if ( INT(time) .ge. 24 ) then
               time = time - 24
            end if

            write(lt,'(I2.2)') INT(time)
            write(mint, '(I2.2)') INT( MOD(time,1)*60.0 )

            ltime = 'LT=' // lt // ':' // mint
            header = 'day' // date  // ' ' // ltime

c**************************************************************

!- rm 2-grid noize ----------------------------------------------------

            do iz = -2,nzb
            do ix = 0,nx
               uu(ix,iz,1) = ( u(ix-1,iz) + 
     \                         2.0* u(ix,iz) + u(ix+1,iz) )/4.0
               ttpot(ix,iz,1) = ( tpot(ix-1,iz) + 2.0*tpot(ix,iz)
     \                        + tpot(ix+1,iz) ) / 4.0 
            end do
            end do


            do iz = -2,NZB
               uu(-2,iz,1)  = u(nx-2,iz) 
               uu(-1,iz,1)  = u(nx-1,iz) 
               uu(0,iz,1)   = u(nx,iz) 
               uu(nx+1,iz,1)= u(1,iz) 
               uu(nx+2,iz,1)= u(2,iz) 
               ttpot(-2,iz,1)  = tpot(nx-2,iz) 
               ttpot(-1,iz,1)  = tpot(nx-1,iz) 
               ttpot(0,iz,1)   = tpot(nx,iz) 
               ttpot(nx+1,iz,1)= tpot(1,iz) 
               ttpot(nx+2,iz,1)= tpot(2,iz) 
            end do

            do iz = -2,nzb
            do ix = -2,nxb
               u(ix,iz) = uu(ix,iz,1)
               tpot(ix,iz) = ttpot(ix,iz,1)
            end do
            end do

!- rm 2-grid noize ----------------------------------------------------
            
            do iz = -2,nzb            
               do ix = -2, nxb
                  uu   (ix,iz,1) = u(ix,iz)
                  vv   (ix,iz,1) = v(ix,iz)
                  ttpot(ix,iz,1) = tpot(ix,iz)
               end do
            end do
                 
            do ix = 1, NX+1
               us(ix) = REAL(u(ix-1,0))
            end do

            call grfrm
            call grswnd ( 0.0, 51.2, -25., 25. )
            call grsvpt ( 0.15, 0.95, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( '', 'X[km]', 'U', '[m/sec]' )
            call usdaxs
            call uxsttl ( 'T', header, 1. )
            call uulinz ( nx+1, xplot, us, 1, 3 )


         end if
      end if

      end do

      call grcls

      call FCLOSE

      end


