C*********************************************************************
C
C     PROGRAM vsurf2
C
C     地表風速・地表気温日変化図作成プログラム(VL1 との比較用)
C
C     履歴 2001/01/08 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      program vsurf2
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = NZG )
      parameter ( NT  = 200 )
      parameter ( DAY = 86400.0D0 )
      parameter ( HOUR= 3600.0D0 )
!-----------------------------------------------------------------      

      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 )

      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 )

      dimension
     \     u    ( -2:nxb, -2:nzb ),
     \     tpot ( -2:nxb, -2:nzb )


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

      real*4
     \     YMIN1,  YMAX1, YMIN2, YMAX2

      real*4 
     \     TIMEL(NT), USURF(NT), TEMP(NT)

      real*4
     \     USURF_MAX, USURF_MIN, USURF_MEAN, USURF_DEV,
     \     TEMP_MAX , TEMP_MIN , TEMP_MEAN , TEMP_DEV

      character*14 dheader
      character*18 pheader

      character*5 CLAT
      character*5 CLS
      character*2 mint
      character*2 lt
      character*2 date
      character*1 CSGI
      character*100 dumy

      character*100 note1, note2
!-----------------------------------------------------------------      
      YMIN1 = 0.0
      YMAX1 = 15.0

      YMIN2 = 235.0
      YMAX2 = 250.0

      open(10,file='exparam')
      read(10,*) IRUN
      read(10,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read(10,*) YLAT, ALS

      write(CLS,  '(F5.1)') ALS
      write(CLAT, '(F5.1)') YLAT
      pheader = 'Ls=' // CLS // ' ' // CSGI(172) // '=' // CLAT

      tinterv = nloop2 * ( dtime / 2. )

c*** 時刻変換: 日付設定 ************************************
      time = time0 + times
      write(date, '(I2.2)') INT( time/day ) + 1
      dheader = 'day=' // date 
c**************************************************************

      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

      write(6,*) 'data number = ', nloop1
      write(6,*) 'INPUT DATA BEGINING, END, & IX' 
      write(6,*) 'itb ite ix?'
      read (*,*) itb, ite, itx

      il = 1
      do it = itb, ite
         
         read(31) time, ((u(ix,iz),ix=-2,nxb), iz=-2,nzb)
         read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)

!- rm 2-grid noize ----------------------------------------------------
         if ( it .ge. itb ) then

         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
!- end rm 2-grid noize ------------------------------------------------


         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 ----------------------------------------------------


         USURF(il) = REAL( abs(u(itx,0)) )
         TEMP (il) = pai0p(0) * ( tpot(itx,0) + tpot0p(0) )

         TIMEL(il) = REAL( (MOD ( TIME, DAY ) + TIMES ) / HOUR )
         write(*,*) TIMEL(il), USURF(il), TEMP(il)

         TIMEL(il) = REAL
     \        ( MOD( (MOD ( TIME, DAY ) + TIMES ), HOUR ) / 60.0 )
         if ( TIMEL(il) .eq. 0.0 ) then
            TIMEL(il) = 60.0
         end if



         il = il + 1 
         end if

      end do

      ntmax = ite-itb+1

!- min, max, mean, sigam ----------------------------------------------

      USURF_MAX = usurf(1)
      USURF_MIN = usurf(1)
      USURF_MEAN = 0.0
      USURF_DEV = 0.0

      TEMP_MAX = temp(1)
      TEMP_MIN = temp(1)
      TEMP_MEAN = 0.0
      TEMP_DEV  = 0.0

      do it = 2, ntmax
         usurf_max = max( usurf_max, usurf(it) )
         usurf_min = min( usurf_min, usurf(it) )
         temp_max = max( temp_max, temp(it) )
         temp_min = min( temp_min, temp(it) )
      end do

      do it = 1, ntmax
         usurf_mean = usurf_mean + usurf(it) / ntmax
         temp_mean = temp_mean + temp(it) / ntmax
      end do

      do it = 1, ntmax
         usurf_dev = usurf_dev + abs( usurf_mean - usurf(it) ) / ntmax
         temp_dev  = temp_dev + abs( temp_mean - temp(it) ) / ntmax
      end do

      write(*,*) usurf_max, usurf_min, usurf_dev
      write(*,*) temp_max, temp_min, temp_dev

      write(note1,'(a5,f5.2,a6,f5.2,a7,f5.2,a7,f5.2)') 
     \     'Umax=', usurf_max, ' Umin=', usurf_min, 
     \     ' Umean=', usurf_mean, ' sigma=', usurf_dev

      write(note2,'(a5,f5.1,a6,f5.1,a7,f5.1,a7,f5.2)') 
     \     'Tmax=', temp_max, ' Tmin=', temp_min, 
     \     ' Tmean=', temp_mean, ' sigma=', temp_dev

!- min, max, mean, sigam ----------------------------------------------

      call sgpwsn
      read(*,*) iws

      CALL SWISET('IWIDTH',  840)
      CALL SWISET('IHEIGHT', 400)
      call swlset( 'LDUMP', .TRUE. )

      call gropn(iws)
      call sglset( 'LCORNER', .FALSE. )
      call swlset( 'LALT',  .TRUE. )
      call sglset( 'LFULL', .TRUE. )
      call uzfact(0.5)

      call grfrm
      call grswnd( TIMEL(1), TIMEL(ite-itb+1), YMIN1, YMAX1 )   
      call grsvpt( 0.15, 0.95, 0.30, 0.43)          
      call grstrn( 1 )
      call grstrf
      call ussttl( ' ', '', 'U', 'm/sec' )
      call usdaxs
      call uulinz(ntmax, TIMEL, USURF, 1, 3  )

      call uxsttl( 'T', dheader, 1. )      
      call uxsttl( 'B', note1, -1. )

      call grfig
      call grswnd( TIMEL(1), TIMEL(ite-itb+1), YMIN2, YMAX2 )   
      call grsvpt( 0.15, 0.95, 0.1, 0.23)          
      call grstrn( 1 )
      call grstrf
      call ussttl( ' ', 'Minute', 'Temp.', 'K' )
      call usdaxs
      call uulinz(ntmax, TIMEL, temp, 1, 3  )

      call uxsttl( 'B', note2, -1. )                     

      call grcls                                                        

      call FCLOSE


      end

