Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:34:40 2011 FILE NAME: radiation_heatbalance.f90 PROGRAM NAME: radiation_heatbalance DIAGNOSTIC LIST LINE LEVEL( NO.): DIAGNOSTIC MESSAGE 107 vec ( 1): Vectorized loop. 122 vec ( 4): Vectorized array expression. 123 vec ( 4): Vectorized array expression. 124 vec ( 4): Vectorized array expression. 125 vec ( 4): Vectorized array expression. 191 opt ( 11): Fused array assignments. :line 191 - 195 191 vec ( 4): Vectorized array expression. 191 vec ( 4): Vectorized array expression. 202 vec ( 1): Vectorized loop. 202 vec ( 1): Vectorized loop. 211 vec ( 1): Vectorized loop. 211 vec ( 1): Vectorized loop. 221 opt ( 11): Fused array assignments. :line 221 - 232 221 vec ( 4): Vectorized array expression. Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:34:40 2011 FILE NAME: radiation_heatbalance.f90 PROGRAM NAME: radiation_heatbalance TRANSFORMATION LIST LINE FORTRAN STATEMENT 1 != Module Radiation 2 ! 3 ! Authors:: YAMASHITA Tatsuya, SUGIYAMA Ko-ichiro 4 ! Version:: $Id: radiation_heatbalance.f90,v 1.6 2011-10-04 05:16:27 sugiyama Exp $ 5 ! Tag Name:: $Name: arare5-20111010 $ 6 ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. 7 ! License:: See COPYRIGHT[link:../../COPYRIGHT] 8 ! 9 !== Overview 10 ! 11 !モデルの放射過程を計算するためのパッケージ型モジュール 12 !具体的には以下の項を計算するための関数を格納する. 13 ! * 一様冷却 14 ! 15 !== Error Handling 16 ! 17 !== Bugs 18 ! 19 !== Note 20 ! 21 ! 22 !== Future Plans 23 ! 24 ! 25 26 module Radiation_HeatBalance 27 ! 28 !モデルの放射過程を計算するためのパッケージ型モジュール 29 !具体的には以下の項を計算するための関数を格納する. 30 ! * 一様冷却 31 ! 32 33 !モジュール読み込み 34 use dc_types, only: DP, STRING 35 use dc_iounit, only: FileOpen 36 use dc_message, only: MessageNotify 37 use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 38 39 use mpi_wrapper,only: myrank 40 use namelist_util, only: namelist_filename 41 use timeset, only: TimeN 42 use gridSet, only: imin, & !x 方向の配列の下限 43 & imax, & !x 方向の配列の上限 44 & jmin, & !y 方向の配列の下限 45 & jmax, & !y 方向の配列の上限 46 & kmin, & !z 方向の配列の下限 47 & kmax, & !z 方向の配列の上限 48 & nx, ny, nz 49 use axesset, only: z_Z, z_dz !Z 座標軸(スカラー格子点) 50 use constants,only: DayTime, & ! 1 日の長さ [s] 51 & PressSfc, & !地表面圧力 52 & CpDry, & !定圧比熱 53 & CvDry, & !定積比熱 54 & GasRDry !気体定数 55 use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場 56 & xyz_VelSoundBZ, &! 音速 57 & xyz_VPTempBZ, &! 仮温位 58 & xyz_PTempBZ !温位の基本場 59 use setmargin,only: SetMargin_xyz 60 61 !暗黙の型宣言禁止 62 implicit none 63 64 !private 属性 65 private 66 67 !関数を public にする. 68 public Radiation_heatbalance_init 69 public Radiation_heatbalance_forcing 70 71 !変数定義 72 real(DP),save :: RadCoolRate = 0.0d0 !一様放射加熱率 [K/day] 73 integer, save :: IdxHeatUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 74 integer, save :: IdxHeatDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 75 integer, save :: IdxCoolUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 76 integer, save :: IdxCoolDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 77 real(DP),save, allocatable :: xyz_RadHeightHeat(:,:,:) !放射加熱が存在する領域 78 real(DP),save, allocatable :: xyz_RadHeightCool(:,:,:) !放射加熱が存在する領域 79 80 contains 81 82 !!!------------------------------------------------------------------------!!! 83 subroutine Radiation_heatbalance_init 84 ! 85 !NAMELIST から放射強制の設定を取得 86 ! 87 88 !暗黙の型宣言禁止 89 implicit none 90 91 !入力変数 92 real(8) :: HeightHeatUp = 0.0d0 !放射強制を与える鉛直領域の上限 93 real(8) :: HeightHeatDown = 0.0d0 !放射強制を与える鉛直領域の上限 94 real(8) :: HeightCoolUp = 0.0d0 !放射強制を与える鉛直領域の上限 95 real(8) :: HeightCoolDown = 0.0d0 !放射強制を与える鉛直領域の上限 96 integer :: k !ループ変数 97 integer :: unit 98 99 ! NAMELIST から情報を取得 100 NAMELIST /radiation_heatbalance_nml/ & 101 & RadCoolRate, HeightHeatUp, HeightHeatDown, HeightCoolUp, HeightCoolDown 102 103 call FileOpen(unit, file=namelist_filename, mode='r') 104 read(unit, NML=radiation_heatbalance_nml) 105 close(unit) 106 107 do k = kmin, kmax-1 108 if ( z_Z(k) <= HeightHeatUp .AND. HeightHeatUp < z_Z(k+1) ) then 109 IdxHeatUp = k 110 elseif ( z_Z(k) <= HeightHeatDown .AND. HeightHeatDown < z_Z(k+1) ) then 111 IdxHeatDown = k 112 elseif ( z_Z(k) <= HeightCoolUp .AND. HeightCoolUp < z_Z(k+1) ) then 113 IdxCoolUp = k 114 elseif ( z_Z(k) <= HeightCoolDown .AND. HeightCoolDown < z_Z(k+1) ) then 115 IdxCoolDown = k 116 end if 117 end do 118 119 allocate( xyz_RadHeightHeat(imin:imax, jmin:jmax, kmin:kmax) ) 120 allocate( xyz_RadHeightCool(imin:imax, jmin:jmax, kmin:kmax) ) 121 122 xyz_RadHeightHeat = 1.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t245 = 1, (xyz_radheightheat.DSC.U3 - xyz_radheightheat.DSC.L3 . 1 + 1)*(t13 - t14 + 1)*(t11 - t12 + 1) . xyz_radheightheat(t12+t245-1,t14,xyz_radheightheat.DSC.L3) = . 1 1.00000000000000e+000 . end do 123 xyz_RadHeightHeat(:,:,IdxHeatDown:IdxHeatUp) = 1.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t254=1,(idxheatup+1-idxheatdown)*(t13-t14+1)*(t11-t12+1) . xyz_radheightheat(t12+t254-1,t14,idxheatdown) = . 1 1.00000000000000e+000 . end do 124 xyz_RadHeightCool = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t263 = 1, (xyz_radheightcool.DSC.U3 - xyz_radheightcool.DSC.L3 . 1 + 1)*(t2 - t3 + 1)*(t0 - t1 + 1) . xyz_radheightcool(t1+t263-1,t3,xyz_radheightcool.DSC.L3) = . 1 0.0000000000000000e+000 . end do 125 xyz_RadHeightCool(:,:,IdxCoolDown:IdxCoolUp) = 1.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t272=1,(idxcoolup+1-idxcooldown)*(t2-t3+1)*(t0-t1+1) . xyz_radheightcool(t1+t272-1,t3,idxcooldown) = . 1 1.00000000000000e+000 . end do 126 127 ! Output 128 ! 129 if (myrank == 0) then 130 call MessageNotify( "M", & 131 & "Radiation", "RadCoolRate = %f", d=(/RadCoolRate/)) 132 call MessageNotify( "M", & 133 & "Radiation", "HeightHeatUp = %f", d=(/HeightHeatUp/)) 134 call MessageNotify( "M", & 135 & "Radiation", "HeightHeatDown= %f", d=(/HeightHeatDown/)) 136 call MessageNotify( "M", & 137 & "Radiation", "HeightCoolUp = %f", d=(/HeightCoolUp/)) 138 call MessageNotify( "M", & 139 & "Radiation", "HeightCoolDown= %f", d=(/HeightCoolDown/)) 140 end if 141 142 call HistoryAutoAddVariable( & 143 & varname='PTempRad', & 144 & dims=(/'x','y','z','t'/), & 145 & longname='Radiation term of potential temperature', & 146 & units='K.s-1"', & 147 & xtype='float') 148 149 call HistoryAutoAddVariable( & 150 & varname='ExnerRad', & 151 & dims=(/'x','y','z','t'/), & 152 & longname='Radiation term of exner function', & 153 & units='s-1"', & 154 & xtype='float') 155 156 end subroutine Radiation_heatbalance_init 157 158 !!!------------------------------------------------------------------------!!! 159 subroutine radiation_heatbalance_forcing(xyz_Exner, xyz_PTemp, xyz_DPTempDt, xyz_DExnerDt) 160 ! 161 ! 温位の放射強制項. 162 ! 地表面から RadHeight で指定された高度までの間で空間的に一様な加熱, 163 ! RadHeight から RadHeight2 までの間で空間的な一様な冷却を与える. 164 ! 放射強制全体として加熱と冷却がバランスするように時々刻々一様加熱率 165 ! および一様冷却率を変化させる. 166 ! 冷却の振幅を与え, 加熱の振幅を調節する. 167 168 !暗黙の型宣言を禁止 169 implicit none 170 171 !変数定義 172 real(DP), intent(in) :: xyz_Exner(imin:imax, jmin:jmax, kmin:kmax) 173 real(DP), intent(in) :: xyz_PTemp(imin:imax, jmin:jmax, kmin:kmax) 174 real(DP), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax) 175 real(DP), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax) 176 real(DP) :: xyz_DPTempDt0(imin:imax, jmin:jmax, kmin:kmax) 177 real(DP) :: xyz_DExnerDt0(imin:imax, jmin:jmax, kmin:kmax) 178 real(DP) :: xyz_Rad(imin:imax, jmin:jmax, kmin:kmax) 179 real(DP) :: xyz_RadPI(imin:imax, jmin:jmax, kmin:kmax) 180 real(DP) :: xyz_DensSum(imin:imax, jmin:jmax, kmin:kmax) 181 !密度(基本場成分と擾乱成分の和) 182 real(DP) :: HeatSum !ある時刻での単位時間当たりの加熱量 183 real(DP) :: CoolSum !ある時刻での単位時間当たりの冷却量 184 real(DP) :: RadHeatRate 185 integer :: i, j, k 186 187 ! 初期化 188 ! 189 HeatSum = 0.0d0 190 CoolSum = 0.0d0 191 xyz_DPTempDt0 = xyz_DPTempDt . if (xyz_dptempdt0.DSC.U2 + 1 - xyz_dptempdt0.DSC.L2 .gt. 0) then . J1 = and(xyz_dptempdt0.DSC.U2 + 1 - xyz_dptempdt0.DSC.L2,1) . !CDIR NODEP . do t438 = 1, J1 . !CDIR NODEP . do t440 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t440-1,t438-1+ . 1 xyz_dptempdt0.DSC.L2,t436+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t142+t440-1,t438-1+t144,t436+t146) . xyz_dexnerdt0(xyz_dexnerdt0.DSC.L1+t440-1,t438-1+ . 1 xyz_dexnerdt0.DSC.L2,t436+xyz_dexnerdt0.DSC.L3) = . 2 xyz_dexnerdt(t152+t440-1,t438-1+t154,t436+t156) . xyz_denssum(xyz_denssum.DSC.L1+t440-1,t438-1+ . 1 xyz_denssum.DSC.L2,t436+xyz_denssum.DSC.L3) = presssfc . 2 *(xyz_exnerbz(t56+t440-1,t438-1+t58,t436+t60)+ . 3 xyz_exner(t122+t440-1,t438-1+t124,t436+t126))**(cvdry/ . 4 gasrdry)/(gasrdry*(xyz_ptempbz(xyz_ptempbz.DSC.L1+t440 . 5 -1,t438-1+xyz_ptempbz.DSC.L2,t436+xyz_ptempbz.DSC.L3)+ . 6 xyz_ptemp(t132+t440-1,t438-1+t134,t436+t136))) . end do . end do . !CDIR NODEP . do t438 = J1 + 1, xyz_dptempdt0.DSC.U2 + 1 - . 1 xyz_dptempdt0.DSC.L2, 2 . !CDIR NODEP . do t440 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t440-1,t438-1+ . 1 xyz_dptempdt0.DSC.L2,t436+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t142+t440-1,t438-1+t144,t436+t146) . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t440-1,t438+ . 1 xyz_dptempdt0.DSC.L2,t436+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t142+t440-1,t438+t144,t436+t146) . xyz_dexnerdt0(xyz_dexnerdt0.DSC.L1+t440-1,t438-1+ . 1 xyz_dexnerdt0.DSC.L2,t436+xyz_dexnerdt0.DSC.L3) = . 2 xyz_dexnerdt(t152+t440-1,t438-1+t154,t436+t156) . xyz_dexnerdt0(xyz_dexnerdt0.DSC.L1+t440-1,t438+ . 1 xyz_dexnerdt0.DSC.L2,t436+xyz_dexnerdt0.DSC.L3) = . 2 xyz_dexnerdt(t152+t440-1,t438+t154,t436+t156) . xyz_denssum(xyz_denssum.DSC.L1+t440-1,t438-1+ . 1 xyz_denssum.DSC.L2,t436+xyz_denssum.DSC.L3) = presssfc . 2 *(xyz_exnerbz(t56+t440-1,t438-1+t58,t436+t60)+ . 3 xyz_exner(t122+t440-1,t438-1+t124,t436+t126))**(cvdry/ . 4 gasrdry)/(gasrdry*(xyz_ptempbz(xyz_ptempbz.DSC.L1+t440 . 5 -1,t438-1+xyz_ptempbz.DSC.L2,t436+xyz_ptempbz.DSC.L3)+ . 6 xyz_ptemp(t132+t440-1,t438-1+t134,t436+t136))) . xyz_denssum(xyz_denssum.DSC.L1+t440-1,t438+ . 1 xyz_denssum.DSC.L2,t436+xyz_denssum.DSC.L3) = presssfc . 2 *(xyz_exnerbz(t56+t440-1,t438+t58,t436+t60)+xyz_exner( . 3 t122+t440-1,t438+t124,t436+t126))**(cvdry/gasrdry)/( . 4 gasrdry*(xyz_ptempbz(xyz_ptempbz.DSC.L1+t440-1,t438+ . 5 xyz_ptempbz.DSC.L2,t436+xyz_ptempbz.DSC.L3)+xyz_ptemp( . 6 t132+t440-1,t438+t134,t436+t136))) . end do . end do . endif 192 xyz_DExnerDt0 = xyz_DExnerDt 193 194 ! 密度の算出 195 xyz_DensSum = PressSfc & 196 & * (xyz_ExnerBZ + xyz_Exner)**( CvDry /GasRDry ) & 197 & / (GasRDry * (xyz_PTempBZ + xyz_PTemp ) ) 198 199 ! ある時刻での (加熱量/加熱率) を計算 200 do k = IdxHeatDown, IdxHeatUp 201 do j = 1, ny 202 do i = 1, nx 203 HeatSum = HeatSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) 204 end do 205 end do . if (ny .gt. 0) then . J2 = and(ny,3) . do j = 1, J2 . !CDIR NODEP . do i = 1, nx . heatsum = heatsum + z_dz(k)*cpdry*xyz_denssum(i,j,k) . end do . end do . do j = J2 + 1, ny, 4 . D1 = z_dz(k) . !CDIR NODEP . do i = 1, nx . heatsum = heatsum + D1*cpdry*xyz_denssum(i,j,k) + D1* . 1 cpdry*xyz_denssum(i,j+1,k) + D1*cpdry*xyz_denssum(i,j+ . 2 2,k) + D1*cpdry*xyz_denssum(i,j+3,k) . end do . end do . endif 206 end do 207 208 ! ある時刻での (冷却量/冷却率) を計算 209 do k = IdxCoolDown, IdxCoolUp 210 do j = 1, ny 211 do i = 1, nx 212 CoolSum = CoolSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) 213 end do 214 end do . if (ny .gt. 0) then . J3 = and(ny,3) . do j = 1, J3 . !CDIR NODEP . do i = 1, nx . coolsum = coolsum + z_dz(k)*cpdry*xyz_denssum(i,j,k) . end do . end do . do j = J3 + 1, ny, 4 . D2 = z_dz(k) . !CDIR NODEP . do i = 1, nx . coolsum = coolsum + D2*cpdry*xyz_denssum(i,j,k) + D2* . 1 cpdry*xyz_denssum(i,j+1,k) + D2*cpdry*xyz_denssum(i,j+ . 2 2,k) + D2*cpdry*xyz_denssum(i,j+3,k) . end do . end do . endif 215 end do 216 217 ! 加熱率の算出 218 ! RadCoolRate が負値であることに注意. 219 RadHeatRate = - RadCoolRate * CoolSum / HeatSum 220 221 xyz_Rad = & . !CDIR NODEP . do t473 = 1, t370 . xyz_rad(xyz_rad.DSC.L1+t473-1,t471+xyz_rad.DSC.L2,t469+ . 1 xyz_rad.DSC.L3) = xyz_radheightheat(t12+t473-1,t471+t14,t469 . 2 +t16)*radheatrate/(xyz_exnerbz(t56+t473-1,t471+t58,t469+t60) . 3 *daytime) + xyz_radheightcool(xyz_radheightcool.DSC.L1+t473- . 4 1,t471+xyz_radheightcool.DSC.L2,t469+ . 5 xyz_radheightcool.DSC.L3)*radcoolrate/(xyz_exnerbz(t56+t473- . 6 1,t471+t58,t469+t60)*daytime) . xyz_dptempdt(t142+t473-1,t471+t144,t469+t146) = xyz_dptempdt0( . 1 xyz_dptempdt0.DSC.L1+t473-1,t471+xyz_dptempdt0.DSC.L2,t469+ . 2 xyz_dptempdt0.DSC.L3) + xyz_rad(xyz_rad.DSC.L1+t473-1,t471+ . 3 xyz_rad.DSC.L2,t469+xyz_rad.DSC.L3) . xyz_radpi(xyz_radpi.DSC.L1+t473-1,t471+xyz_radpi.DSC.L2,t469+ . 1 xyz_radpi.DSC.L3) = xyz_velsoundbz(xyz_velsoundbz.DSC.L1+ . 2 t473-1,t471+xyz_velsoundbz.DSC.L2,t469+xyz_velsoundbz.DSC.L3 . 3 )**2.00000000000000e+000*xyz_rad(xyz_rad.DSC.L1+t473-1,t471+ . 4 xyz_rad.DSC.L2,t469+xyz_rad.DSC.L3)/(cpdry*xyz_vptempbz( . 5 xyz_vptempbz.DSC.L1+t473-1,t471+xyz_vptempbz.DSC.L2,t469+ . 6 xyz_vptempbz.DSC.L3)**2.00000000000000e+000) . xyz_dexnerdt(t152+t473-1,t471+t154,t469+t156) = xyz_dexnerdt0( . 1 xyz_dexnerdt0.DSC.L1+t473-1,t471+xyz_dexnerdt0.DSC.L2,t469+ . 2 xyz_dexnerdt0.DSC.L3) + xyz_radpi(xyz_radpi.DSC.L1+t473-1, . 3 t471+xyz_radpi.DSC.L2,t469+xyz_radpi.DSC.L3) . end do 222 & xyz_RadHeightHeat * RadHeatRate / ( xyz_ExnerBZ * DayTime ) & 223 & + xyz_RadHeightCool * RadCoolRate / ( xyz_ExnerBZ * DayTime ) 224 225 xyz_DPTempDt = xyz_DPTempDt0 + xyz_Rad 226 227 ! 圧力変化 228 ! 229 xyz_RadPI = (xyz_VelSoundBZ ** 2.0d0) * xyz_Rad & 230 & / (CpDry * xyz_VPTempBZ ** 2.0d0) 231 232 xyz_DExnerDt = xyz_DExnerDt0 + xyz_RadPI 233 234 ! output 235 ! 236 call HistoryAutoPut(TimeN, 'PTempRad', xyz_Rad(1:nx,1:ny,1:nz)) 237 call HistoryAutoPut(TimeN, 'ExnerRad', xyz_RadPI(1:nx,1:ny,1:nz)) 238 239 240 ! Set Margin 241 ! 242 call SetMargin_xyz(xyz_DPTempDt) 243 call SetMargin_xyz(xyz_DExnerDt) 244 245 end subroutine radiation_heatbalance_forcing 246 247 end module Radiation_HeatBalance Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:34:40 2011 FILE NAME: radiation_heatbalance.f90 PROGRAM NAME: radiation_heatbalance FORMAT LIST LINE LOOP FORTRAN STATEMENT 1: != Module Radiation 2: ! 3: ! Authors:: YAMASHITA Tatsuya, SUGIYAMA Ko-ichiro 4: ! Version:: $Id: radiation_heatbalance.f90,v 1.6 2011-10-04 05:16:27 sugiyama Exp $ 5: ! Tag Name:: $Name: arare5-20111010 $ 6: ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. 7: ! License:: See COPYRIGHT[link:../../COPYRIGHT] 8: ! 9: !== Overview 10: ! 11: !モデルの放射過程を計算するためのパッケージ型モジュール 12: !具体的には以下の項を計算するための関数を格納する. 13: ! * 一様冷却 14: ! 15: !== Error Handling 16: ! 17: !== Bugs 18: ! 19: !== Note 20: ! 21: ! 22: !== Future Plans 23: ! 24: ! 25: 26: module Radiation_HeatBalance 27: ! 28: !モデルの放射過程を計算するためのパッケージ型モジュール 29: !具体的には以下の項を計算するための関数を格納する. 30: ! * 一様冷却 31: ! 32: 33: !モジュール読み込み 34: use dc_types, only: DP, STRING 35: use dc_iounit, only: FileOpen 36: use dc_message, only: MessageNotify 37: use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 38: 39: use mpi_wrapper,only: myrank 40: use namelist_util, only: namelist_filename 41: use timeset, only: TimeN 42: use gridSet, only: imin, & !x 方向の配列の下限 43: & imax, & !x 方向の配列の上限 44: & jmin, & !y 方向の配列の下限 45: & jmax, & !y 方向の配列の上限 46: & kmin, & !z 方向の配列の下限 47: & kmax, & !z 方向の配列の上限 48: & nx, ny, nz 49: use axesset, only: z_Z, z_dz !Z 座標軸(スカラー格子点) 50: use constants,only: DayTime, & ! 1 日の長さ [s] 51: & PressSfc, & !地表面圧力 52: & CpDry, & !定圧比熱 53: & CvDry, & !定積比熱 54: & GasRDry !気体定数 55: use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場 56: & xyz_VelSoundBZ, &! 音速 57: & xyz_VPTempBZ, &! 仮温位 58: & xyz_PTempBZ !温位の基本場 59: use setmargin,only: SetMargin_xyz 60: 61: !暗黙の型宣言禁止 62: implicit none 63: 64: !private 属性 65: private 66: 67: !関数を public にする. 68: public Radiation_heatbalance_init 69: public Radiation_heatbalance_forcing 70: 71: !変数定義 72: real(DP),save :: RadCoolRate = 0.0d0 !一様放射加熱率 [K/day] 73: integer, save :: IdxHeatUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 74: integer, save :: IdxHeatDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 75: integer, save :: IdxCoolUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 76: integer, save :: IdxCoolDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 77: real(DP),save, allocatable :: xyz_RadHeightHeat(:,:,:) !放射加熱が存在する領域 78: real(DP),save, allocatable :: xyz_RadHeightCool(:,:,:) !放射加熱が存在する領域 79: 80: contains 81: 82: !!!------------------------------------------------------------------------!!! 83: subroutine Radiation_heatbalance_init 84: ! 85: !NAMELIST から放射強制の設定を取得 86: ! 87: 88: !暗黙の型宣言禁止 89: implicit none 90: 91: !入力変数 92: real(8) :: HeightHeatUp = 0.0d0 !放射強制を与える鉛直領域の上限 93: real(8) :: HeightHeatDown = 0.0d0 !放射強制を与える鉛直領域の上限 94: real(8) :: HeightCoolUp = 0.0d0 !放射強制を与える鉛直領域の上限 95: real(8) :: HeightCoolDown = 0.0d0 !放射強制を与える鉛直領域の上限 96: integer :: k !ループ変数 97: integer :: unit 98: 99: ! NAMELIST から情報を取得 100: NAMELIST /radiation_heatbalance_nml/ & 101: & RadCoolRate, HeightHeatUp, HeightHeatDown, HeightCoolUp, HeightCoolDown 102: 103: call FileOpen(unit, file=namelist_filename, mode='r') 104: read(unit, NML=radiation_heatbalance_nml) 105: close(unit) 106: 107: V------> do k = kmin, kmax-1 108: | if ( z_Z(k) <= HeightHeatUp .AND. HeightHeatUp < z_Z(k+1) ) then 109: | IdxHeatUp = k 110: | elseif ( z_Z(k) <= HeightHeatDown .AND. HeightHeatDown < z_Z(k+1) ) then 111: | IdxHeatDown = k 112: | elseif ( z_Z(k) <= HeightCoolUp .AND. HeightCoolUp < z_Z(k+1) ) then 113: | IdxCoolUp = k 114: | elseif ( z_Z(k) <= HeightCoolDown .AND. HeightCoolDown < z_Z(k+1) ) then 115: | IdxCoolDown = k 116: | end if 117: V------ end do 118: 119: allocate( xyz_RadHeightHeat(imin:imax, jmin:jmax, kmin:kmax) ) 120: allocate( xyz_RadHeightCool(imin:imax, jmin:jmax, kmin:kmax) ) 121: 122: WWW==== xyz_RadHeightHeat = 1.0d0 123: W**==== xyz_RadHeightHeat(:,:,IdxHeatDown:IdxHeatUp) = 1.0d0 124: ***==== xyz_RadHeightCool = 0.0d0 125: ***==== xyz_RadHeightCool(:,:,IdxCoolDown:IdxCoolUp) = 1.0d0 126: 127: ! Output 128: ! 129: if (myrank == 0) then 130: call MessageNotify( "M", & 131: & "Radiation", "RadCoolRate = %f", d=(/RadCoolRate/)) 132: call MessageNotify( "M", & 133: & "Radiation", "HeightHeatUp = %f", d=(/HeightHeatUp/)) 134: call MessageNotify( "M", & 135: & "Radiation", "HeightHeatDown= %f", d=(/HeightHeatDown/)) 136: call MessageNotify( "M", & 137: & "Radiation", "HeightCoolUp = %f", d=(/HeightCoolUp/)) 138: call MessageNotify( "M", & 139: & "Radiation", "HeightCoolDown= %f", d=(/HeightCoolDown/)) 140: end if 141: 142: call HistoryAutoAddVariable( & 143: & varname='PTempRad', & 144: & dims=(/'x','y','z','t'/), & 145: & longname='Radiation term of potential temperature', & 146: & units='K.s-1"', & 147: & xtype='float') 148: 149: call HistoryAutoAddVariable( & 150: & varname='ExnerRad', & 151: & dims=(/'x','y','z','t'/), & 152: & longname='Radiation term of exner function', & 153: & units='s-1"', & 154: & xtype='float') 155: 156: end subroutine Radiation_heatbalance_init 157: 158: !!!------------------------------------------------------------------------!!! 159: subroutine radiation_heatbalance_forcing(xyz_Exner, xyz_PTemp, xyz_DPTempDt, xyz_DExnerDt) 160: ! 161: ! 温位の放射強制項. 162: ! 地表面から RadHeight で指定された高度までの間で空間的に一様な加熱, 163: ! RadHeight から RadHeight2 までの間で空間的な一様な冷却を与える. 164: ! 放射強制全体として加熱と冷却がバランスするように時々刻々一様加熱率 165: ! および一様冷却率を変化させる. 166: ! 冷却の振幅を与え, 加熱の振幅を調節する. 167: 168: !暗黙の型宣言を禁止 169: implicit none 170: 171: !変数定義 172: real(DP), intent(in) :: xyz_Exner(imin:imax, jmin:jmax, kmin:kmax) 173: real(DP), intent(in) :: xyz_PTemp(imin:imax, jmin:jmax, kmin:kmax) 174: real(DP), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax) 175: real(DP), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax) 176: real(DP) :: xyz_DPTempDt0(imin:imax, jmin:jmax, kmin:kmax) 177: real(DP) :: xyz_DExnerDt0(imin:imax, jmin:jmax, kmin:kmax) 178: real(DP) :: xyz_Rad(imin:imax, jmin:jmax, kmin:kmax) 179: real(DP) :: xyz_RadPI(imin:imax, jmin:jmax, kmin:kmax) 180: real(DP) :: xyz_DensSum(imin:imax, jmin:jmax, kmin:kmax) 181: !密度(基本場成分と擾乱成分の和) 182: real(DP) :: HeatSum !ある時刻での単位時間当たりの加熱量 183: real(DP) :: CoolSum !ある時刻での単位時間当たりの冷却量 184: real(DP) :: RadHeatRate 185: integer :: i, j, k 186: 187: ! 初期化 188: ! 189: HeatSum = 0.0d0 190: CoolSum = 0.0d0 191: **V----> xyz_DPTempDt0 = xyz_DPTempDt 192: ||| xyz_DExnerDt0 = xyz_DExnerDt 193: ||| 194: ||| ! 密度の算出 195: **V---- xyz_DensSum = PressSfc & 196: & * (xyz_ExnerBZ + xyz_Exner)**( CvDry /GasRDry ) & 197: & / (GasRDry * (xyz_PTempBZ + xyz_PTemp ) ) 198: 199: ! ある時刻での (加熱量/加熱率) を計算 200: +------> do k = IdxHeatDown, IdxHeatUp 201: |+-----> do j = 1, ny 202: ||V----> do i = 1, nx 203: ||| HeatSum = HeatSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) 204: ||V---- end do 205: |+----- end do 206: +------ end do 207: 208: ! ある時刻での (冷却量/冷却率) を計算 209: +------> do k = IdxCoolDown, IdxCoolUp 210: |+-----> do j = 1, ny 211: ||V----> do i = 1, nx 212: ||| CoolSum = CoolSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) 213: ||V---- end do 214: |+----- end do 215: +------ end do 216: 217: ! 加熱率の算出 218: ! RadCoolRate が負値であることに注意. 219: RadHeatRate = - RadCoolRate * CoolSum / HeatSum 220: 221: **V----> xyz_Rad = & 222: ||| & xyz_RadHeightHeat * RadHeatRate / ( xyz_ExnerBZ * DayTime ) & 223: ||| & + xyz_RadHeightCool * RadCoolRate / ( xyz_ExnerBZ * DayTime ) 224: ||| 225: ||| xyz_DPTempDt = xyz_DPTempDt0 + xyz_Rad 226: ||| 227: ||| ! 圧力変化 228: ||| ! 229: ||| xyz_RadPI = (xyz_VelSoundBZ ** 2.0d0) * xyz_Rad & 230: ||| & / (CpDry * xyz_VPTempBZ ** 2.0d0) 231: ||| 232: **V---- xyz_DExnerDt = xyz_DExnerDt0 + xyz_RadPI 233: 234: ! output 235: ! 236: call HistoryAutoPut(TimeN, 'PTempRad', xyz_Rad(1:nx,1:ny,1:nz)) 237: call HistoryAutoPut(TimeN, 'ExnerRad', xyz_RadPI(1:nx,1:ny,1:nz)) 238: 239: 240: ! Set Margin 241: ! 242: call SetMargin_xyz(xyz_DPTempDt) 243: call SetMargin_xyz(xyz_DExnerDt) 244: 245: end subroutine radiation_heatbalance_forcing 246: 247: end module Radiation_HeatBalance