C********************************************************************* C C PROGRAM tgrand C C 地中・地表温度日変化図作成プログラム(時刻, 条件指定) C C 履歴 2000/09/22 小高正嗣 C 2000/11/21 小高正嗣; D 論用に書き直し C C********************************************************************* program tgrand !----------------------------------------------------------------- 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 \ TPOT (-2:NXB,-2:NZB), \ TEMPG (-2:NXB,NZG), ! 地中温度 \ GGRID (-1:NZGB), ! 鉛直座標 \ GDZ (0:NZGB), ! 鉛直格子間隔(無次元) \ GDZM (0:NZG) ! 鉛直格子間隔(半整数格子上) real*8 \ 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 \ GGRID_p(NZP), \ TEMPG_p(NZP) real*4 \ ZMIN , ZMAX, TMIN, TMAX real*4 \ TSFC(0:NT), TIMEL(0:NT), TEMPA(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 !----------------------------------------------------------------- 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 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 DATA ( GGRID(I), I=-1,NZGB ) 1 / -6.0D0 , 2 -6.0D0 , -4.0D0 , -2.7D0 , -1.8D0, -1.2D0, 3 -0.79D0, -0.53D0, -0.35D0, -0.2D0, -0.1D0, 4 0.0D0 , 0.0D0 / do iz = 1, nzp GGRID_p(iz)= GGRID(iz) end do ZMIN = GGRID_p(1) ZMAX = 0.0 TMIN = 160.0 TMAX = 290.0 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', 840) CALL SWISET('IHEIGHT', 420) call sglset( 'LFULL', .TRUE. ) call swlset( 'LDUMP', .TRUE. ) call gropn(iws) call sglset( 'LCORNER', .FALSE. ) CALL sglset( 'LCHAR', .TRUE. ) call sldiv ('Y', 2, 1 ) call grfrm CALL sglset( 'LCHAR', .TRUE. ) call grswnd( TMIN, TMAX, ZMIN, ZMAX ) call grsvpt( 0.2, 0.8, 0.2, 0.8) call grstrn( 1 ) call grstrf call ussttl( 'Grand Temperature', 'K', 'Depth', 'Z' ) call usdaxs li = 23 do it = 1, ite 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 c*** 時刻変換: 秒から LT へ************************************ time = time + times time = ( time - INT( time / day ) * day ) / hour if ( time .gt. 24 ) then time = time - 24 end if write(lt, '(I2.2)') INT(time) c************************************************************** do IZ = 1, NZG TEMPG_p(IZ) = 0.0 end do do IZ = 1, NZG do IX = 1, NX TEMPG_p(IZ) = TEMPG_p(IZ) + TEMPG(IX,IZ) / NX end do end do CALL sgrset( 'FFCT', 0.05 ) CALL sgsplc( lt ) call SGSPLI( li ) call sgplu(nzp, tempg_p, ggrid_p) li = li + 10 end if end if end do Call uxsttl( 'T', dheader, 1. ) C CALL uxsttl( 'T', pheader, 1. ) call FCLOSE call FOPEN ( irun+1 ) ilate = INT(TIMES / HOUR ) do it = 1, nloop1 read(34) TIME, ((TPOT(IX,IZ), IX=-2,NXB),IZ=-2,NZB) read(40) TIME, ((TEMPG(IX,IZ),IX=-2,NXB),IZ= 1,NZG) TEMPA(it+ilate) = 0.0 TSFC (it+ilate) = 0.0 do IX = 1, NX TEMPA(it+ilate) = TEMPA(it+ilate) + \ ( TPOT(IX,0) + TPOT0(0) ) / NX TSFC (it+ilate) = TSFC (it+ilate) + \ TEMPG(IX,NZG) / NX end do TIMEL(it+ilate) = (MOD ( TIME, DAY ) + TIMES ) / HOUR end do TEMPA(ilate) = TEMPA(nloop1+ilate) TSFC(ilate) = TSFC(nloop1+ilate) TIMEL(ilate) = TIMEL(nloop1+ilate) il = 0 do it = 1, nloop1+ilate if ( TIMEL(it) .ge. 24.0 ) then TIMEL(il) = TIMEL(it) - 24.0 TSFC (il) = TSFC(it) TEMPA(il) = TEMPA(it) il = il + 1 end if end do call grfrm call sglset( 'LCHAR', .FALSE. ) call grswnd( 0.0, 24.0, TMIN, TMAX ) call grsvpt( 0.2, 0.8, 0.2, 0.8) call grstrn( 1 ) call grstrf call ussttl( 'Local Time', 'Hour', \ 'Grand Surface Temp.', 'K' ) call usdaxs call SGSPLI( 22 ) call sgplu(nloop1+1, TIMEL, TSFC) call SGSPLI( 42 ) call sgplu(nloop1+1, TIMEL, TEMPA) call uxsttl( 'T', dheader, 1. ) c call uxsttl( 'T', pheader, 1. ) call grcls call FCLOSE end