C********************************************************************* C C PROGRAM xstress C C 地表摩擦水平分布図(時刻指定) C C 履歴 1999/08/03 小高正嗣 C 2000/11/21 小高正嗣; D 論用に書き直し C C********************************************************************* program xstress *----------------------------------------------------- implicit real*8 ( a-h, o-z ) *----------------------------------------------------- #include "grid_size_M.f" parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) ) *----------------------------------------------------- integer ibr(nbr3) *----------------------------------------------------- parameter ( day = 86400.0 ) parameter ( hour = 3600.0 ) parameter ( PI = 3.141592D0 ) 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 ), \ w ( -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 ) 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 \ up ( nx+1,nz+1 ), \ vp ( nx+1,nz+1 ), \ wp ( nx+1,nz+1 ), \ tpotp ( nx+1, nz+1 ), \ tpall ( nx+1, nz+1 ), \ cdragvp( nx+1 ), \ richdp( nx+1 ), \ tsfcp ( nx+1 ), \ fheat ( nx+1 ), \ xplot ( nx+1 ), \ stress( nx+1 ), \ stress_cr( nx+1 ) real*4 \ stress_1( nx+1 ), \ stress_2( nx+1 ), \ cdragvp_1( nx+1 ), \ cdragvp_2( nx+1 ), \ uplot(nx+1) character*2 mint character*2 lt character*8 ltime character*2 date character*14 header character*100 dumy *----------------------------------------------------- open (11,file='exparam') read (11,*) irun read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times read (11,*) ylat, als close(11) 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 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', 720) CALL SWISET('IHEIGHT', 500) call swlset( 'LDUMP', .TRUE. ) call gropn(iws) call swlset( 'LALT', .TRUE. ) call sglset( 'LFULL', .TRUE. ) call sglset( 'LCORNER', .FALSE. ) call uzfact(0.6) call sldiv ('T', 1, 2 ) do it = 1, ite 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(33) time, ((w(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) 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 do IX = -2, NXB TSFC(IX) = tempg(IX, NZG) end do c*** 時刻変換: 秒から LT へ************************************ time = time + times write(date,'(I2.2)') INT( time / day ) + 1 time = ( time - INT( time / day ) * day ) / hour if ( INT(time) .ge. 24 ) then time = time - 24 end if write(lt,'(I2.2)') INT(time) write(mint, '(I2.2)') INT( MOD(time,1)*60.0 ) ltime = 'LT=' // lt // ':' // mint header = 'day' // date // ' ' // ltime c************************************************************** !- 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 iz = 1,nz+1 do ix = 1, nx+1 up(ix,iz) = u(ix-1,iz-1) vp(ix,iz) = v(ix-1,iz-1) wp(ix,iz) = w(ix-1,iz-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 ix = 1,nx+1 cdragvp( ix )= cdragv(ix-1) richdp( ix ) = MIN(richd(ix-1),0.25) tsfcp ( ix ) = tsfc(ix-1) end do * * do ix = 1,nx+1 vsfc = SQRT( up(ix,1)**2 + vp(ix,1)**2 ) fheat(ix) = - cdragvp(ix) * VSFC * CPDRY *dens0(0) \ * ( tpall(ix,1)*PPAI0(0) - tsfcp(ix) ) stress(ix) = d0p(0) * cdragvp(ix) * VSFC **2 end do do iz = -2,nzb do ix = -2, nxb uu (ix,iz,1) = u(ix,iz) vv (ix,iz,1) = 10.0 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 cdragvp_1( ix )= cdragv(ix-1) vsfc = SQRT( up(ix,1)**2 + 10.0**2 ) stress_1(ix) = d0p(0) * cdragvp_1(ix) * VSFC **2 end do do iz = -2,nzb do ix = -2, nxb uu (ix,iz,1) = u(ix,iz)+10.0 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 cdragvp_2( ix )= cdragv(ix-1) vsfc = SQRT( (up(ix,1)+10.0)**2 + vp(ix,1)**2 ) stress_2(ix) = d0p(0) * cdragvp_2(ix) * VSFC **2 end do call grfrm call grswnd ( 0.0, 51.2, 0., 0.08 ) call grsvpt ( 0.15, 0.95, 0.1, 0.3 ) call grstrn (1) call uspfit call grstrf call udpset ( 'rsizet', 0.01 ) call ussttl ( '', 'X[km]', 'Stress', 'Pa' ) call usdaxs call uxsttl ( 'T', header, 1. ) call SGSPLI( 42 ) call sgplu ( nx+1, xplot, stress ) call SGSPLI( 33 ) call sgplu ( nx+1, xplot, stress_2 ) do ix = 1, nx+1 stress_cr(ix) = 0.03 end do call SGSPLI( 772 ) call sgplu ( nx+1, xplot, stress_cr ) do ix = 1, nx+1 stress_cr(ix) = 0.04 end do call SGSPLI( 22 ) call sgplu ( nx+1, xplot, stress_cr ) call grfrm call grswnd ( 0.0, 51.2, 0., 0.08 ) call grsvpt ( 0.15, 0.95, 0.1, 0.3 ) call grstrn (1) call uspfit call grstrf call udpset ( 'rsizet', 0.01 ) call ussttl ( '', 'X[km]', 'Stress', 'Pa' ) call usdaxs call uxsttl ( 'T', header, 1. ) call SGSPLI( 42 ) call sgplu ( nx+1, xplot, stress ) call SGSPLI( 33 ) call sgplu ( nx+1, xplot, stress_1 ) do ix = 1, nx+1 stress_cr(ix) = 0.03 end do call SGSPLI( 772 ) call sgplu ( nx+1, xplot, stress_cr ) do ix = 1, nx+1 stress_cr(ix) = 0.04 end do call SGSPLI( 22 ) call sgplu ( nx+1, xplot, stress_cr ) end if end if end do call grcls call FCLOSE end