C================================================================= C C PROGRAM netheat C C 正味加熱量の日変化図作成プログラム C C 履歴 2000/05/25 小高正嗣 C 2000/11/21 小高正嗣; D 論用に書き直し C================================================================= program netheat !----------------------------------------------------------------- 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 ( nxp = nx, nzp= 100) parameter ( NBND = 16 ) parameter ( NBNDSR1= 10 ) ! NIR band number 1(4.3μm) parameter ( NBNDSR2= 10 ) ! NIR band number 1(2.7μm) parameter ( NBNDSR3= 8 ) ! NIR band number 1(2.0μm) parameter ( DELNU = 25.0D0 ) ! band width parameter ( SIGMA = 5.67D-8 ) parameter ( PI = 3.14159263D0 ) parameter ( TSOL = 5760.0D0 ) ! surface temp. of Sun parameter ( DAY = 86400.0D0) parameter ( HOUR = 3600.0D0 ) parameter ( NT = 100 ) 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 \ TPOT ( -2:NXB, -2:NZB, NROT ), TEMP ( -2:NXB, -2:NZB ), \ QVAP ( -2:NXB, -2:NZB, NROT ) real*8 \ u ( -2:nxb, -2:nzb ), w ( -2:nxb, -2:nzb ) dimension \ FUMOMA(-2:NXB,2), FVMOMA(-2:NXB,2), FHEATA(-2:NXB,2), \ FEVAPA(-2:NXB,2), FWPREA(-2:NXB,2) dimension \ tempg ( -2:nxb,nzg ) DIMENSION \ FDZ ( -2:NZB ), FDZM ( -2:NZB ) DIMENSION \ DTIRAD ( -2:NXB, -2:NZB ), \ DTSRAD ( -2:NXB, -2:NZB ), \ DTQIRAD( -2:NXB, -2:NZB ), \ DTQSRAD( -2:NXB, -2:NZB ) DIMENSION \ OPL (-2:NZB), OPL_S(-2:NZB) dimension \ QTAU (-2:NXB,-2:NZB), QDTAU (-2:NXB,-2:NZB), \ QTAUIR1(-2:NXB,-2:NZB), QDTAUIR1(-2:NXB,-2:NZB), \ QTAUIR2(-2:NXB,-2:NZB), QDTAUIR2(-2:NXB,-2:NZB) DIMENSION \ TAUN (0:NZ+1,0:NZ+1,NBND), \ TAUN_SR1(0:NZ+1,0:NZ+1,NBNDSR1), \ TAUN_SR2(0:NZ+1,0:NZ+1,NBNDSR2), \ TAUN_SR3(0:NZ+1,0:NZ+1,NBNDSR3) dimension \ WNU(NBND), SNU(NBND), RNU(NBND) dimension \ WNU_SR1(NBNDSR1), SNU_SR1(NBNDSR1), RNU_SR1(NBNDSR1), \ WNU_SR2(NBNDSR2), SNU_SR2(NBNDSR2), RNU_SR2(NBNDSR2), \ WNU_SR3(NBNDSR3), SNU_SR3(NBNDSR3), RNU_SR3(NBNDSR3) dimension \ BBRAD ( -2:NXB, -2:NZB ), FXRADU( -2:NXB, -2:NZB ), \ FXRADD( -2:NXB, -2:NZB ), FXRADN( -2:NXB, -2:NZB ), \ BBRAD_S( -2:NXB ) dimension \ SXRADU (-2:NXB,-2:NZB), SXRADD (-2:NXB,-2:NZB), \ SXRADN (-2:NXB,-2:NZB), \ SXRADU1(-2:NXB,-2:NZB), SXRADD1(-2:NXB,-2:NZB), \ SXRADN1(-2:NXB,-2:NZB), \ SXRADU2(-2:NXB,-2:NZB), SXRADD2(-2:NXB,-2:NZB), \ SXRADN2(-2:NXB,-2:NZB), \ SXRADU3(-2:NXB,-2:NZB), SXRADD3(-2:NXB,-2:NZB), \ SXRADN3(-2:NXB,-2:NZB), \ SQRAD1 (-2:NXB,-2:NZB), \ SQRAD2 (-2:NXB,-2:NZB), \ SQRAD3 (-2:NXB,-2:NZB), \ F0(-2:NXB), FXRADG(-2:NXB,-2:NZB) DIMENSION \ FXRADU_NU(-2:NXB,-2:NZB,NBND), \ FXRADD_NU(-2:NXB,-2:NZB,NBND), \ SXRADU_NU1(-2:NXB,-2:NZB,NBNDSR1), \ SXRADD_NU1(-2:NXB,-2:NZB,NBNDSR1), \ SXRADU_NU3(-2:NXB,-2:NZB,NBNDSR3), \ SXRADD_NU3(-2:NXB,-2:NZB,NBNDSR3) dimension \ FXQRADU (-2:NXB,-2:NZB), FXQRADD (-2:NXB,-2:NZB), \ FXQRADN (-2:NXB,-2:NZB), \ FXQRADU1(-2:NXB,-2:NZB), FXQRADD1(-2:NXB,-2:NZB), \ FXQRADN1(-2:NXB,-2:NZB), \ FXQRADU2(-2:NXB,-2:NZB), FXQRADD2(-2:NXB,-2:NZB), \ FXQRADN2(-2:NXB,-2:NZB), \ QRAD1 (-2:NXB,-2:NZB), QRAD2 (-2:NXB,-2:NZB) dimension \ SXQRADU(-2:NXB,-2:NZB), SXQRADD(-2:NXB,-2:NZB), \ SXQRADN(-2:NXB,-2:NZB), SXSOL (-2:NXB,-2:NZB) dimension \ SCRADW1(-2:NXB,-2:NZB), SCRADW2(-2:NXB,-2:NZB) dimension \ ASR(0:NZ,0:NZ), BSR(0:NZ,0:NZ), CSR(0:NZ,0:NZ), \ DSR(0:NZ,0:NZ), RSR(0:NZ) , SSR(0:NZ) dimension TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB ) real*4 \ DTPRAD(nx+1,nz+1) real*4 \ qvapp ( nxp+1, nzp+1 ), TLEV1, TLEV2 real*4 \ YMIN1, YMAX1, YMIN2, YMAX2 character*2 mint character*2 lt character*2 date character*17 header character*100 dumy NAMELIST /CONST1/ GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT NAMELIST /CONST2/ CDRAG, VSFC0, PSFC !----------------------------------------------------------------- ! 物理定数の設定 !----------------------------------------------------------------- YMAX1 = 600.0 YMIN1 = -YMAX1 YMAX2 = 25 YMIN2 = -YMAX2 call SETCST O ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT ) CALL SETCS2 O ( CDRAG, VSFC0, TSFC(-2), QVSFC(-2), O FCOLI, O PSFC, C RMVAP, D NX, NXB ) call SETCSP C ( S00, EXCTL, OBLIC, OMGA0 ) call SETDUST O ( QEXT_SR, QEXT_IR1, QEXT_IR2, O QALB_SR, QALB_IR1, QALB_IR2, O QASM_SR, QASM_IR1, QASM_IR2, O QDENS , QREFF , QRMOD ) C C *** Delta-Eddington Approximation *** (2000/11/09) C QALB_SRE = QALB_DE( QALB_SR, QASM_SR ) QASM_SRE = QASM_DE( QASM_SR ) QALB_IR1 = QALB_DE( QALB_IR1, QASM_IR1 ) QASM_IR1 = QASM_DE( QASM_IR1 ) QALB_IR2 = QALB_DE( QALB_IR2, QASM_IR2 ) QASM_IR2 = QASM_DE( QASM_IR2 ) !----------------------------------------------------------------- open (11,file='exparam') read (11,*) irun read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, TIMES read (11,*) YLAT, ALS close(11) YLAT = YLAT * PI/180.0D0 ! transfer from degree to rad ALS = ALS * PI/180.0D0 ! transfer from degree to rad tinterv = nloop2 * ( dtime / 2.) 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 !----------------------------------------------------------------- ! 放射場の初期条件設定; バンドパラメターの設定 ! (Houghton, 1986;The physics of atmosphere, 2nd ed., ! Cambridge Univ. Press) DATA ( WNU(I), I=1,NBND ) 1 / 512.5D0, 537.5D0, 562.5D0, 587.5D0, 2 612.5D0, 637.5D0, 662.5D0, 687.5D0, 3 712.5D0, 737.5D0, 762.5D0, 787.5D0, 4 812.5D0, 837.5D0, 862.5D0, 887.5D0 / DATA ( SNU(I), I=1,NBND ) 1 / 1.952D-2, 2.785D-1, 5.495D-1, 5.331D0, 2 5.196D+2, 7.778D+3, 8.746D+4, 2.600D+4, 3 1.232D+3, 2.042D+2, 7.278D0 , 1.337D0, 4 3.974D-1, 1.280D-2, 2.501D-3, 3.937D-3 / DATA ( RNU(I), I=1,NBND ) 1 / 2.870D-1, 1.215D0 , 2.404D0 , 1.958D+1, 2 5.804D+1, 2.084D+2, 7.594D+2, 2.635D+2, 3 8.387D+1, 2.852D+1, 6.239D0 , 2.765D0 , 4 8.897D-1, 3.198D-1, 1.506D-1, 1.446D-1 / ! CO2 4.3 μm DATA ( WNU_SR1(I), I=1,NBNDSR1 ) 1 / 2212.5D0, 2237.5D0, 2262.5D0, 2287.5D0, 2 2312.5D0, 2337.5D0, 2362.5D0, 2387.5D0, 3 2412.5D0, 2437.5D0 / DATA ( SNU_SR1(I), I=1,NBNDSR1 ) 1 / 9.504D-1, 2.217D+2, 4.566D+3, 7.965D+3, 2 1.055D+5, 5.587D+5, 6.819D+5, 1.256D+4, 3 7.065D-2, 8.522D-2 / DATA ( RNU_SR1(I), I=1,NBNDSR1 ) 1 / 2.866D0 , 3.000D+1, 1.134D+2, 2.011D+2, 2 5.880D+2, 1.206D+3, 1.182D+3, 8.873D+1, 3 3.404D-1, 4.236D-1 / ! CO2 2.7 μm DATA ( WNU_SR2(I), I=1,NBNDSR2 ) 1 / 3150.0D0, 3250.0D0, 3350.0D0, 3450.0D0, 2 3550.0D0, 3650.0D0, 3750.0D0, 3850.0D0, 3 3950.0D0, 4050.0D0 / DATA ( SNU_SR2(I), I=1,NBNDSR2 ) 1 / 1.324D-1, 7.731D-2, 1.232D0 , 5.159D0, 2 4.299D+3, 1.543D+4, 1.649D+4, 1.180D-1, 3 1.464D-2, 1.251D-2 / DATA ( RNU_SR2(I), I=1,NBNDSR2 ) 1 / 9.836D-1, 4.900D-1, 2.952D0 , 7.639D0 , 2 1.914D+2, 3.245D+2, 2.722D+2, 9.535D-1, 3 2.601D-1, 2.021D-1 / ! CO2 2.0 μm DATA ( WNU_SR3(I), I=1,NBNDSR3 ) 1 / 4650.0D0, 4750.0D0, 4850.0D0, 4950.0D0, 2 5050.0D0, 5150.0D0, 5250.0D0, 5350.0D0 / DATA ( SNU_SR3(I), I=1,NBNDSR3 ) 1 / 2.185D-1, 2.040D0 , 1.197D+2, 4.829D+2, 2 8.778D+1, 8.346D+1, 8.518D-2, 4.951D-1 / DATA ( RNU_SR3(I), I=1,NBNDSR3 ) 1 / 1.916D0, 6.475D0 , 3.112D+1, 5.759D+1, 2 2.012D+1, 1.804D+1, 8.474D-1, 1.597D0 / !----------------------------------------------------------------- CALL SETGRD O ( ZGRID(-2), FDZ(-2), FDZM(-2), DX, D NZ, NZB ) call RBASIC I ( PRES0, PSFC , SNU , RNU , O OPL , TAUN , C GRAV , DELNU, D NZ , NZB , NBND ) !----------------------------------------------------------------- write(6,*) 'data number = ', nloop1 write(6,*) 'INPUT DATA BEGINING, END, & INTERVAL' write(6,*) 'itb ite interval ?' read (*,*) itb, ite, interv ilate = INT( TIMES / HOUR ) call sgpwsn read(*,*) iws call gropn(-iws) call sldiv( 'T', 1, 2 ) call swlset( 'LALT', .TRUE. ) call sglset( 'LFULL', .TRUE. ) call sglset( 'LCORNER', .FALSE. ) call udlset( 'LMSG', .FALSE. ) call uzfact(0.6) * do it = 1, ite read(34) TIME, ((TPOT(IX,IZ,1),ix=-2,nxb), IZ=-2,NZB ) read(35) TIME, ((qvap(ix,iz,1),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 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 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 ) header = 'day=' // date // ' ' // 'LT=' // lt // ':' // mint c************************************************************** !----------------------------------------------------------------- do IX = -2,NXB TSFC(IX) = TEMPG(ix,NZG) end do call CLFXSOL I ( S00 , ALS , YLAT , TIME, TIMES, O S0 , COSZ , C EXCTL, OBLIC , OMGA0, DAY , PI ) C C *** CO2 Infrared Heaating/Cooling *** C CALL CLRAD I ( TPOT , TAUN , PSFC , I 1, O DTIRAD, B TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ , C WNU , C TSFC , CPDRY , RDRY , DX , DELNU, GRAV, W TEMP , BBRAD , BBRAD_S, W FXRADU, FXRADD, FXRADN, FXRADU_NU, FXRADD_NU, D NX , NXB , NZ , NZB , NROT , NBND ) C C *** Dust Optical Depth Set up *** (2000/11/09) C C do ITQ = 1, NROT C do IZ = -2, NZB C do IX = -2, NXB C QVAP(IX,IZ,ITQ) = QVAP0(IZ) C end do C end do C end do call SETQTAU I ( DENS0, QVAP , I 1, O QTAU , QDTAU, O QTAUIR1, QDTAUIR1, QTAUIR2, QDTAUIR2, C QREFF , QDENS , QALB_SR, QASM_SR , C QEXT_SR, QEXT_IR1, QEXT_IR2, PI , G FDZ , DX, D NX , NXB , NZ , NZB , NROT ) C write(*,*) QTAU(0,0), QTAUIR1(0,0), QTAUIR2(0,0) C C *** Dust Infrared Radiative Heating/Cooling *** (2000/11/09) C For 5-11.6 micro m C DINU1 = 900.0D0 DINU2 = 2000.0D0 DELNUD = DELNU * 4.0 call CLDIRAD I ( TPOT , PSFC , 1, I QTAUIR1, QDTAUIR1, I DINU1 , DINU2 , DELNUD, O QRAD1 , B PRES0 , DENS0 , TEMP0, TPOT0, C QALB_IR1, QASM_IR1, C RDRY , CPDRY , TSFC , GRAV , PI , W TEMP , BBRAD , BBRAD_S, W FXQRADU1, FXQRADD1, FXQRADN1, W SCRADW1, SCRADW2, W ASR , BSR , CSR , DSR , W RSR , SSR , G FDZ , FDZM , DX , D NX , NXB , NZ , NZB , NROT ) C C For 20-200 micro m C DINU1 = 50.0D0 DINU2 = 500.0D0 DELNUD = DELNU * 2.0 call CLDIRAD I ( TPOT , PSFC , 1, I QTAUIR2, QDTAUIR2, I DINU1 , DINU2 , DELNUD, O QRAD2 , B PRES0 , DENS0 , TEMP0, TPOT0, C QALB_IR2, QASM_IR2, C RDRY , CPDRY , TSFC , GRAV , PI , W TEMP , BBRAD , BBRAD_S, W FXQRADU2, FXQRADD2, FXQRADN2, W SCRADW1, SCRADW2, W ASR , BSR , CSR , DSR , W RSR , SSR , G FDZ , FDZM , DX , D NX , NXB , NZ , NZB , NROT ) do iz = -2, NZB do ix = -2, NXB DTQIRAD(ix,iz) = QRAD1 (ix,iz) + QRAD2 (ix,iz) FXQRADU(ix,iz) = FXQRADU1(ix,iz) + FXQRADU2(ix,iz) FXQRADD(ix,iz) = FXQRADD1(ix,iz) + FXQRADD2(ix,iz) FXQRADN(ix,iz) = FXQRADN1(ix,iz) + FXQRADN2(ix,iz) end do end do C C *** Solar Heating/Cooling *** (2000/10/31) C if ( COSZ .gt. 0.0D0 ) then do IZ = -2, NZB OPL_S(IZ) = OPL(IZ) / COSZ end do call CLTAUN I ( OPL_S, SNU_SR1, RNU_SR1, PRES0, O TAUN_SR1, C DELNU, D NZ , NZB, NBNDSR1 ) call CLTAUN I ( OPL_S, SNU_SR2, RNU_SR2, PRES0, O TAUN_SR2, C DELNU*4.0, D NZ , NZB, NBNDSR2 ) call CLTAUN I ( OPL_S, SNU_SR3, RNU_SR3, PRES0, O TAUN_SR3, C DELNU*4.0, D NZ , NZB, NBNDSR3 ) C *** CO2 4.3μm *** call CLSRAD I ( TAUN_SR1, S0 , COSZ , O SXRADU1, SXRADD1, SXRADN1, SQRAD1, B TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ , C WNU_SR1, C CPDRY , TSOL , DX , SIGMA, C DELNU , GRAV , PI , W SXRADU_NU1, SXRADD_NU1, D NX , NXB, NZ , NZB , NROT, NBNDSR1 ) C *** CO2 2.7μm *** call CLSRAD I ( TAUN_SR2, S0 , COSZ , O SXRADU2, SXRADD2, SXRADN2, SQRAD2, B TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ , C WNU_SR2, C CPDRY , TSOL , DX , SIGMA, C DELNU*4.0 , GRAV , PI , W SXRADU_NU1, SXRADD_NU1, D NX , NXB, NZ , NZB , NBNDSR2 ) C *** CO2 2.0μm *** call CLSRAD I ( TAUN_SR3, S0 , COSZ , O SXRADU3, SXRADD3, SXRADN3, SQRAD3, B TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ , C WNU_SR3, C CPDRY , TSOL , DX , SIGMA, C DELNU*4.0 , GRAV , PI , W SXRADU_NU3, SXRADD_NU3, D NX , NXB, NZ , NZB , NBNDSR3 ) do iz = -2, NZB do ix = -2, NXB DTSRAD(ix,iz) = SQRAD1(ix,iz) \ + SQRAD2(ix,iz) + SQRAD3(ix,iz) SXRADD(ix,iz) = SXRADD1(ix,iz) \ + SXRADD2(ix,iz) + SXRADD3(ix,iz) end do end do C *** Dust *** (2000/11/09) call CLDSRAD I ( QTAU , QDTAU , S0 , COSZ , O DTQSRAD, SXSOL , B TPOT0 , TEMP0 , PRES0 , DENS0 , C ALBEDO, QALB_SRE, QASM_SRE, C DELNU , CPDRY , TSOL , SIGMA, C GRAV , PI , W SXQRADU, SXQRADD, SXQRADN, W SCRADW1, SCRADW2, W ASR , BSR , CSR , DSR , W RSR , SSR , G FDZ , FDZM , DX , D NX , NXB , NZ , NZB ) else do iz = -2, NZB do ix = -2, NXB SXRADU(ix,iz) = 0.0D0 SXRADD(ix,iz) = 0.0D0 SXRADN(ix,iz) = 0.0D0 SXQRADU(ix,iz)= 0.0D0 SXQRADD(ix,iz)= 0.0D0 SXQRADN(ix,iz)= 0.0D0 SXSOL (ix,iz) = 0.0D0 SQRAD1(ix,iz) = 0.0D0 SQRAD2(ix,iz) = 0.0D0 SQRAD3(ix,iz) = 0.0D0 DTSRAD(ix,iz) = 0.0D0 DTQSRAD(ix,iz)= 0.0D0 end do end do end if do iz = 0, NZ do ix = 0, NX DTPRAD(ix+1,iz+1) = ( C \ DTIRAD(ix,iz) + DTSRAD(ix,iz) + \ DTIRAD(ix,iz) + \ DTQIRAD(ix,iz) + DTQSRAD(ix,iz) ) * 86400.0 end do end do C==================================================================== call grfrm call grswnd ( 0.0, 51.2, 10.0, 20.0 ) C call grsvpt ( 0.15, 0.95, 0.15, 0.506 ) call grsvpt ( 0.15, 0.95, 0.18, 0.34 ) call grstrn(1) call uspfit call grstrf call udpset ( 'rsizet', 0.01 ) C call ussttl ( 'Radiative Heating [K/day]', C \ 'X[km]', '', 'Z[km]' ) call ussttl ( '', 'X[km]', '', 'Z[km]' ) c** 2-grid noize の除去 *************************************** do iz = 106,nz+1 C do iz = 1,nz+1 do ix = 1, nx+1 qvapp(ix,iz-105) = ( qvap(ix-2,iz-1,1) \ + 2.0 * qvap(ix-1,iz-1,1) \ + qvap(ix+1,iz-1,1) ) / 4.0 end do end do c************************************************************** do i = 1,4 TLEV2 = i + 1 TLEV1 = TLEV2 - 1 TLEV2 = TLEV2 * 1.0e-5 TLEV1 = TLEV1 * 1.0e-5 IPAT = 600 + i call uestlv(TLEV1, TLEV2, IPAT) end do TLEV2 = 10.0e-5 TLEV1 = 5.0e-5 IPAT = 605 call uestlv(TLEV1, TLEV2, IPAT) call uetone ( qvapp, nxp+1, nxp+1, nzp+1 ) call ueitlv call usdaxs call udpset ( 'ICYCLE', 5 ) call udgcla ( -100.0, 100.0, 10.0 ) call udcntz ( DTPRAD(1, 106), nx+1, nxp+1, 101, ibr, nbr3) c************************************************************** call uxsttl( 'T', header, 1. ) end if end if end do call grcls call FCLOSE stop end