C********************************************************************* C C PROGRAM fxmom2 C C 地表摩擦・顕熱フラックス日変化図作成プログラム C C 履歴 2000/11/21 小高正嗣; D 論用に書き直し C 2000/12/02 小高正嗣; 高時間分解能用 C 2001/07/15 小高正嗣; ながれ用 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), stress_cr(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*3 name character*100 dumy !----------------------------------------------------------------- YMIN1 = 0.0 YMAX1 = 0.05 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 CALL SWISET('IWIDTH', 420) CALL SWISET('IHEIGHT', 420) 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', 'Pa' ) call usdaxs call SGSPLI( 43 ) call sgplu ( nloop1, TIMEL, STRESS_MAX ) call SGSPLI( 93 ) call sgplu ( nloop1, TIMEL, STRESS_MEAN ) C call uulinz( nloop1, TIMEL, STRESS_SD, 3, 3 ) call uxsttl( 'T', dheader, 1. ) c call uxsttl( 'T', pheader, 1. ) do ix = 1, nloop1 stress_cr(ix) = 0.03 end do call SGSPLI( 772 ) call sgplu ( nloop1, timel, stress_cr ) do ix = 1, nloop1 stress_cr(ix) = 0.04 end do call SGSPLI( 22 ) call sgplu ( nloop1, timel, stress_cr ) call grcls call FCLOSE end