C*********************************************************************
C
C     PROGRAM fxmom2
C
C     地表摩擦・顕熱フラックス日変化図作成プログラム
C
C     履歴 2000/11/21 小高正嗣; D 論用に書き直し
C          2000/12/02 小高正嗣; 高時間分解能用
C
C*********************************************************************
      program fxmom2
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = NZG )
      parameter ( NT  = 100 )
      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 )

      dimension
     \     u    ( -2:nxb, -2:nzb ), v    ( -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 )

      dimension
     \     FUMOMA(-2:NXB,2), FVMOMA(-2:NXB,2), FHEATA(-2:NXB,2), 
     \     FEVAPA(-2:NXB,2), FWPREA(-2:NXB,2)

      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 
     \     xplot ( nx+1 ), stress( nx+1 )

      real*4
     \     YMIN1,  YMAX1,  YMIN2,  YMAX2


      real*4 
     \     stress_MEAN(0:NT), stress_MAX(0:NT), stress_SD(0:NT),
     \     FHEAT(0:NT), TIMEL(0:NT)
  
      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

!-----------------------------------------------------------------      
      YMIN1 = 0.0
      YMAX1 = 0.04
      YMIN2 = -5.0
      YMAX2 = 20.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 SETCST( grav, psd, amdry, rdry, cpdry, rmvap, alatent )

      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

C      ilate = INT( TIMES / HOUR )

      do it = 0, nloop1-1
         
         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(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)
         read(40) time, ((tempg(ix,iz),ix=-2,nxb), iz=1,nzg)

         read(72) TIME, FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA

         do IX = -2, NXB
            TSFC(IX) = tempg(IX, NZG)
         end do

!- 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
                 
         call CLDRAG
     I      ( uu    , vv    , ttpot, 
     I        1, 
     I        TPOT0, 
     O        RICHD, CDRAGV, 
     G        ZGRID, 
     C        GRAV , TSFC ,
     D        NX, NXB, NZB, NROT )


         do ix = 1,nx+1
            vsfc = SQRT( u(ix-1,0)**2 + v(ix-1,0)**2 )
            stress(ix) = d0p(0) * cdragv(ix-1) * VSFC **2
         end do

         
         STRESS_MAX(it) = 0.0
         STRESS_MEAN(it) = 0.0

         do ix = 1, NX
            STRESS_MAX(it) = 
     \           MAX( STRESS_MAX(it), stress(ix))
            STRESS_MEAN(it) = 
     \           STRESS_MEAN(it) + stress(ix) / NX
         end do

         STRESS_SD(it) = 0.0
         do ix = 1, NX
            STRESS_SD(it) = STRESS_SD(it) +
     \           ( stress(ix) - STRESS_MEAN(it) )**2/NX 
         end do
         STRESS_SD(it)=SQRT(STRESS_SD(it))

         FHEAT(it) = 0.0
         do ix = 1, NX
            FHEAT(it) = FHEAT(it) + 
     \           FHEATA(IX, 1) / tinterv / NX 
         end do


         TIMEL(it) = (MOD ( TIME, DAY ) + TIMES ) / HOUR

      end do

      call sgpwsn
      read(*,*) iws

C      CALL SWISET('IWIDTH',  500)
C      CALL SWISET('IHEIGHT', 500)
C      call swlset( 'LDUMP', .TRUE. )

      call gropn(iws)
      call sglset( 'LCORNER', .FALSE. )

      call grfrm
      CALL sglset( 'LCHAR', .FALSE. )
      call grswnd( TIMEL(0), TIMEL(nloop1-1), YMIN1, YMAX1 )   
      call grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Local Time', '[Hour]', 
     \     'Surface Stress', '[N/m^2]' )
      call usdaxs
      call uulinz(nloop1, TIMEL, STRESS_MAX, 1, 3  )
      call uulinz(nloop1, TIMEL, STRESS_MEAN, 2, 3  )
      call uulinz(nloop1, TIMEL, STRESS_SD, 3, 3  )
               
      call uxsttl( 'T', dheader, 1. )      
c      call uxsttl( 'T', pheader, 1. )      
     
      call grcls                                                        

      call FCLOSE


      end

