C********************************************************************* C C PROGRAM xzcont_ptm_3 C C x-z平面等値線図作成プログラム(時刻, 条件指定) C C 履歴 1997/07/24 中島健介 C 1997/07/25 小高正嗣 C 1997/07/26 小高正嗣 C 1997/07/26 小高正嗣 C 1998/09/17 小高正嗣 C 1998/10/02 小高正嗣 C 1999/01/10 小高正嗣 C 1999/02/04 小高正嗣 C 1999/02/07 小高正嗣 C C********************************************************************* program xzcont_ptm_3 *----------------------------------------------------- implicit real*8 ( a-h, o-z ) *----------------------------------------------------- #include "grid_size_M.f" parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) ) *----------------------------------------------------- integer ibr(nbr3) *----------------------------------------------------- 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*8 \ tpot ( -2:nxb, -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 ) real*4 \ tpotp ( nx+1, nz+1 ), \ tpall ( nx+1, nz+1 ), \ tpave ( nz+1 ), \ tpdev ( nx+1, nz+1 ) character*6 stime c character*2 stime character*11 ttime character*100 dumy character*8 ltime *----------------------------------------------------- real*4 xvpt real*4 yvpt real*4 tpdev_p(nz+1, nz+1) *----------------------------------------------------- parameter ( nzprof = nz ) parameter ( nzp = nz+1, ntp = 120 ) CHARACTER*10 NMZPRF( NZPROF ) character*5 nmrun character*4 from character*1 chara character*1 pz real*4 tend real*4 dtpad0,dqvad0 real*4 zp(nzp) real*4 tptend(nzp),tpmean(nzp),dtpadv(nzp), & dtpdif(nzp),dtpdnl(nzp), & dtpdi0(nzp),dtprad(nzp),tperror(nzp) real*8 qvtend(nzp),qvmean(nzp),dqvadv(nzp), & dqvdif(nzp),dqvdnl(nzp), & dqvdi0(nzp),dqvnfl(nzp),qverror(nzp) real*4 qvmean_4(nzp) real*4 itime(100) real*4 heat(100) real*4 sflux_h real*4 timeh(264) real*4 heath(264) real*4 tinterv integer fnum real*4 para1 integer begin,end,interv integer nloop1 C===================================================================== xvpt = 0.05 yvpt = 0.4 open (11,file='exparam.1') read (11,*) irun read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu tinterv = nloop2 / 2.0 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, & INTERVAL' write(6,*) 'itb ite interval ?' read (*,*) itb, ite, interv iflag = 0 call sgpwsn read(*,*) iws call gropn(iws) call swlset( 'LALT', .TRUE. ) call sglset( 'LFULL', .TRUE. ) call sglset( 'LCORNER', .FALSE. ) call uzfact(0.5) call UDLSET( 'LMSG', .FALSE. ) * do it = 1, ite read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb) 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 へ************************************ c time = time - 86400.*2 c time = time / 3600. + 6 c************************************************************** write(stime,100) INT(time) 100 format(I6.6) c 100 format(I2.2) c ltime = 'LT=' // stime // ':00' ttime = 't=' // stime // 'sec' do iz = 1,nz+1 do ix = 1, nx+1 tpotp(ix,iz) = tpot(ix-1,iz-1) tpall(ix,iz) = tpot(ix-1,iz-1) + tpot0(iz-1) end do end do * do iz = 1, nz+1 tpave(iz) = 0. do ix = 1,nx+1 tpave(iz) = tpave(iz) + tpall(ix,iz) /(nx+1) end do end do do iz = 1,nz+1 do ix = 1, nx+1 tpdev(ix,iz) = tpall(ix,iz) - tpave(iz) end do end do do iz = 1,nz+1 do ix = 1, nz+1 tpdev_p(ix,iz) = tpdev(ix,iz) end do end do if (iflag .eq. 4) then iflag = 1 xvpt = 0.05 else iflag = iflag + 1 end if if ( iflag .eq. 1 ) then call grfrm else call grfig end if call grswnd ( 0.0, 10.0, 0.0, 10.0 ) call grsvpt ( xvpt, xvpt+0.2, 0.48, 0.68 ) call grstrn (1) call uspfit call grstrf call udpset ( 'rsizet', 0.01 ) c call ussttl ( 'anomaly THETA[K]', 'X[km]', '', 'Z[km]' ) call udpset ( 'ICYCLE', 2 ) call uestlv ( -10., 0., 603) call uetone ( tpdev_p, nz+1, nz+1, nz+1 ) call ueitlv if ( iflag .eq. 1 ) then call ussttl ( '', '', '', 'Z[km]' ) call usdaxs else if ( iflag .eq. 4 ) then c call ussttl ( '', 'X[km]', '', '' ) call usdaxs else call slpvpr( 3 ) call uxaxdv( 'B', 1., 2. ) call uxaxdv( 'T', 1., 2. ) call uyaxdv( 'L', 1., 2. ) call uyaxdv( 'R', 1., 2. ) end if call udgcla ( -4., 4., 0.5 ) call udcntz ( tpdev_p, nz+1, nz+1, nz+1, ibr, nbr3) call uxsttl ( 'T', ttime, 1. ) c call uxsttl ( 'T', ltime, 1. ) call uzinit * * xvpt = xvpt + 0.23 end if end if end do call FCLOSE iflag = 0 xvpt = 0.05 do it = 1, ite fnum = 50 + it open(fnum) read(fnum,'(a2,a15,g15.7,a6,g15.7,a4,9a15)') & chara, nmrun, t_begin, from, t_end, & pz, (nmzprf(j),j=40,48) do iz = 1, nzp read(fnum,'(10g15.7)') zp(iz), & tptend(iz), & tpmean(iz), & dtpadv(iz), & dtpdif(iz), & dtpdnl(iz), & dtpad0, & dtpdi0(iz), & dtprad(iz), & tperror(iz) end do close(fnum) do iz = 1, nzp zp(iz) = zp(iz) /1000. end do if ( it .ge. itb .and. it .le. ite ) then if ( mod( (it-itb), interv ) .eq. 0 ) then if (iflag .eq. 4) then iflag = 1 xvpt = 0.05 else iflag = iflag + 1 end if call grfig call grswnd( 230., 270., 0., 5. ) call grsvpt ( xvpt, xvpt+0.2, 0.25, 0.45 ) call grstrn(1) call grstrf if ( iflag .eq. 1 ) then c call ussttl ( '', '', '', 'Z[km]' ) call usdaxs else if ( iflag .eq. 4 ) then c call ussttl ( '', 'THETA[K]', '', '' ) call usdaxs else call usdaxs c call slpvpr( 3 ) c call uxaxdv( 'B', 2., 10. ) c call uxaxdv( 'T', 2., 10. ) c call uyaxdv( 'L', 1., 2. ) c call uyaxdv( 'R', 1., 2. ) end if call uulin(nzp/2, tpmean, zp) xvpt = xvpt + 0.23 end if end if end do * * iflag = 0 xvpt = 0.05 do it = 1, ite fnum = 50 + it open(fnum) read(fnum,'(a2,a15,g15.7,a6,g15.7,a4,9a15)') & chara, nmrun, t_begin, from, t_end, & pz, (nmzprf(j),j=40,48) do iz = 1, nzp read(fnum,'(10g15.7)') zp(iz), & tptend(iz), & tpmean(iz), & dtpadv(iz), & dtpdif(iz), & dtpdnl(iz), & dtpad0, & dtpdi0(iz), & dtprad(iz), & tperror(iz) end do close(fnum) do iz = 1, nzp zp(iz) = zp(iz) /1000. end do if ( it .ge. itb .and. it .le. ite ) then if ( mod( (it-itb), interv ) .eq. 0 ) then if (iflag .eq. 4) then iflag = 1 xvpt = 0.05 else iflag = iflag + 1 end if call grfig call grswnd( 230., 270., 0., 5. ) call grsvpt ( xvpt, xvpt+0.2, 0.25, 0.45 ) call grstrn(1) call grstrf if ( iflag .eq. 1 ) then c call ussttl ( '', '', '', 'Z[km]' ) call usdaxs else if ( iflag .eq. 4 ) then c call ussttl ( '', 'THETA[K]', '', '' ) call usdaxs else call usdaxs c call slpvpr( 3 ) c call uxaxdv( 'B', 2., 10. ) c call uxaxdv( 'T', 2., 10. ) c call uyaxdv( 'L', 1., 2. ) c call uyaxdv( 'R', 1., 2. ) end if call uulin(nzp/2, tpmean, zp) c--------------------------------------------------------- c read(13,*) (itime(j),j=1,nloop1),(heat(j),j=1,nloop1) c write(*,*) 't_begin,t_end' c read(*,*) begin,end c begin = begin / 1800 c end = end / 1800 cc begin = begin / tinterv cc end = end / tinterv c do i = 1, nloop1 cc heat(i) = heat(i)/tinterv c heat(i) = heat(i)/1800. c end do c interv = end - begin + 1 c sflux_h = 0. c do i = begin, end c sflux_h = sflux_h + heat(i) / interv c end do c sflux_h = sflux_h / ( cpdry * dens0(0) * z0p(0)*2 * 51200.) c dtpdif(1) = dtpdif(1) + sflux_h do iz = 1, nzp dtpdif(iz) = ( dtpdif(iz) + dtpdi0(iz) )* 86400. tptend(iz) = tptend(iz) * 86400. dtpadv(iz) = ( dtpadv(iz) + dtpdnl(iz) )* 86400. if ( dtpadv(iz) .gt. 100.0 ) then dtpadv(iz) = 100.0 else if ( dtpadv(iz) .lt. -100.0 ) then dtpadv(iz) = - 100.0 end if dtprad(iz) = dtprad(iz) * 86400. if ( dtprad(iz) .gt. 100.0 ) then dtprad(iz) = 100.0 end if end do call grfig call grswnd( -100., 100., 0., 5. ) call grsvpt ( xvpt, xvpt+0.2, 0.02, 0.22 ) call grstrn(1) call grstrf c call sglset( 'LCHAR', .TRUE. ) if ( iflag .eq. 1 ) then c call ussttl ( '', '', '', 'Z[km]' ) call usdaxs else if ( iflag .eq. 4 ) then call ussttl ( '', '[K/day]', '', '' ) call usdaxs else call usdaxs c call slpvpr( 3 ) c call uxaxdv( 'B', 2., 10. ) c call uxaxdv( 'T', 2., 10. ) c call uyaxdv( 'L', 1., 2. ) c call uyaxdv( 'R', 1., 2. ) end if c call sgsplc( 'tend' ) c call uulinz(nzp/2, tptend, zp, 1, 1) c call sgrset( 'FFCT', 0.8 ) c call sgsplc( 'adv' ) call uulinz(nzp/2, dtpadv, zp, 1, 1) c call sgsplc( 'dif' ) c call sgrset( 'FFCT', 0.65 ) call uulinz(nzp/2, dtpdif, zp, 2, 1) c call sgsplc( 'rad' ) c call sgrset( 'FFCT', 0.9 ) call uulinz(nzp/2, dtprad, zp, 3, 1) c call sglset( 'LCHAR', .FALSE. ) xvpt = xvpt + 0.23 end if end if end do call grcls end