Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:35:04 2011 FILE NAME: surfaceflux_diff.f90 PROGRAM NAME: surfaceflux_diff DIAGNOSTIC LIST LINE LEVEL( NO.): DIAGNOSTIC MESSAGE 122 vec ( 3): Unvectorized loop. 167 vec ( 4): Vectorized array expression. 168 vec ( 4): Vectorized array expression. 170 vec ( 4): Vectorized array expression. 170 vec ( 4): Vectorized array expression. 171 vec ( 4): Vectorized array expression. 171 vec ( 4): Vectorized array expression. 178 vec ( 4): Vectorized array expression. 178 vec ( 4): Vectorized array expression. 182 vec ( 3): Unvectorized loop. 183 vec ( 4): Vectorized array expression. 183 vec ( 4): Vectorized array expression. 191 vec ( 4): Vectorized array expression. 191 vec ( 4): Vectorized array expression. 192 vec ( 4): Vectorized array expression. 192 vec ( 4): Vectorized array expression. 195 vec ( 3): Unvectorized loop. Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:35:04 2011 FILE NAME: surfaceflux_diff.f90 PROGRAM NAME: surfaceflux_diff TRANSFORMATION LIST LINE FORTRAN STATEMENT 1 != Module HeatFlux 2 ! 3 ! Authors:: ODAKA Masatsugu 4 ! Version:: $Id: surfaceflux_diff.f90,v 1.8 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 ! 基づいて計算するモジュール. これは中島 (1994) で用いられた方法である. 13 ! 14 ! 熱フラックスを Fh と凝結物質のフラックス Fq とすると, 温度と凝結物質 15 ! の変化率 H, Q は 16 ! 17 ! H = Fh/Δz_1 18 ! Q = Fq/Δz_1 19 ! 20 ! と表される. ここで Δz_1 は最下層の格子間隔である. 21 ! 22 ! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって 23 ! 計算する. 24 ! 25 ! Fh = - Cdρ|V| * (π_1θ_1 - T_sfc) 26 ! Fh = - Cdρ|V| * (Q_1 - Q*(T_sfc)) 27 ! 28 ! ここで θ_1, Q_1, π_1 は最下層の温度と凝結物質の混合比および無次元 29 ! 圧力関数, T_sfc は下部境界の温度, Q*(T_sfc) は T_sfc で決まる飽和混 30 ! 合比である. 31 ! 32 ! バルク係数 Cd は一定とする. 無次元圧力関数 π_1 は基本場の値を用いる. 33 ! 風速値 |V| は 34 ! 35 ! V = ( V^2 + V_0^2 )^(1/2) 36 ! 37 ! と計算する. 38 ! 39 !== Error Handling 40 ! 41 !== Bugs 42 ! 43 !== Note 44 ! 45 ! 46 !== Future Plans 47 ! 48 ! 49 50 module Surfaceflux_diff 51 ! 52 !下部境界でのフラックスの計算モジュール 53 ! 54 55 !モジュール読み込み 56 use dc_types, only: DP, STRING 57 use dc_iounit, only: FileOpen 58 use dc_message, only: MessageNotify 59 use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 60 61 use mpi_wrapper,only: myrank 62 use gridset, only: imin, & !x 方向の配列の下限 63 & imax, & !x 方向の配列の上限 64 & jmin, & !y 方向の配列の下限 65 & jmax, & !y 方向の配列の上限 66 & kmin, & !z 方向の配列の下限 67 & kmax, & !z 方向の配列の上限 68 & nx, ny, nz, ncmax 69 use axesset, only: z_dz !z 方向の格子点間隔 70 use basicset, only: xyz_ExnerBZ !エクスナー関数の基本場 71 use composition,only : GasNum, SpcWetSymbol 72 use namelist_util, only: namelist_filename 73 use timeset, only: TimeN 74 use setmargin,only: SetMargin_xyz, SetMargin_xyzf 75 76 !暗黙の型宣言禁止 77 implicit none 78 79 !属性の指定 80 private 81 82 !関数を public に設定 83 public surfaceflux_diff_init 84 public surfaceflux_diff_forcing 85 86 !変数定義 87 real(DP), save :: Kappa = 800.0d0 88 89 contains 90 !!!------------------------------------------------------------------------!!! 91 subroutine Surfaceflux_Diff_init 92 ! 93 !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 94 ! 95 96 !暗黙の型宣言禁止 97 implicit none 98 99 !内部変数 100 integer :: l, unit 101 102 !--------------------------------------------------------------- 103 ! NAMELIST から情報を取得 104 ! 105 NAMELIST /surfaceflux_diff_nml/ Kappa 106 107 call FileOpen(unit, file=namelist_filename, mode='r') 108 read(unit, NML=surfaceflux_diff_nml) 109 close(unit) 110 111 if (myrank == 0) then 112 call MessageNotify( "M", "SurfaceFlux", "Kappa = %f", d=(/Kappa/)) 113 end if 114 115 call HistoryAutoAddVariable( & 116 & varname='PTempFlux', & 117 & dims=(/'x','y','z','t'/), & 118 & longname='surface flux of potential temperature', & 119 & units='kg.kg-1.s-1', & 120 & xtype='float') 121 122 do l = 1, ncmax 123 call HistoryAutoAddVariable( & 124 & varname=trim(SpcWetSymbol(l))//'_Flux', & 125 & dims=(/'x','y','z','t'/), & 126 & longname='Surface Flux term of ' & 127 & //trim(SpcWetSymbol(l))//' mixing ratio', & 128 & units='kg.kg-1.s-1', & 129 & xtype='float') 130 end do 131 132 end subroutine Surfaceflux_Diff_init 133 134 135 !!!------------------------------------------------------------------------!!! 136 subroutine Surfaceflux_Diff_forcing( & 137 & xyz_PTemp, xyzf_QMix, & 138 & xyz_DPTempDt, xyzf_DQMixDt & 139 & ) 140 ! 141 ! 下部境界からのフラックスによる温度の変化率を, 142 ! バルク方法に基づいて計算する. 143 ! 144 145 !暗黙の型宣言禁止 146 implicit none 147 148 !変数定義 149 real(DP), intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax) 150 !温位の擾乱成分 151 real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) 152 !温位の擾乱成分 153 real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax) 154 real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) 155 real(DP) :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax) 156 real(DP) :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax) 157 real(DP) :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax) 158 !地表面熱フラックス 159 real(DP) :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax) 160 161 integer :: kz !配列添字 162 integer :: l !ループ変数 163 164 ! 初期化 165 ! 166 kz = 1 167 xyz_Heatflux = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t333 = 1, (xyz_heatflux.DSC.U3 + 1 - xyz_heatflux.DSC.L3)*( . 1 xyz_heatflux.DSC.U2 + 1 - xyz_heatflux.DSC.L2)*( . 2 xyz_heatflux.DSC.U1 + 1 - xyz_heatflux.DSC.L1) . xyz_heatflux(xyz_heatflux.DSC.L1+t333-1,xyz_heatflux.DSC.L2, . 1 xyz_heatflux.DSC.L3) = 0.0000000000000000e+000 . end do 168 xyzf_QMixflux = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t342 = 1, xyzf_qmixflux.DSC.U4*(xyzf_qmixflux.DSC.U3 + 1 - . 1 xyzf_qmixflux.DSC.L3)*(xyzf_qmixflux.DSC.U2 + 1 - . 2 xyzf_qmixflux.DSC.L2)*(xyzf_qmixflux.DSC.U1 + 1 - . 3 xyzf_qmixflux.DSC.L1) . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t342-1,xyzf_qmixflux.DSC.L2, . 1 xyzf_qmixflux.DSC.L3,1) = 0.0000000000000000e+000 . end do 169 170 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,3) . !CDIR NODEP . do t356 = 1, J1 . !CDIR NODEP . do t358 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t358-1,t356-1+ . 1 xyz_dptempdt0.DSC.L2,t354+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t89+t358-1,t356-1+t91,t354+t93) . end do . end do . !CDIR NODEP . do t356 = J1 + 1, xyz_dptempdt0.DSC.U2 + 1 - . 1 xyz_dptempdt0.DSC.L2, 4 . !CDIR NODEP . do t358 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t358-1,t356-1+ . 1 xyz_dptempdt0.DSC.L2,t354+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t89+t358-1,t356-1+t91,t354+t93) . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t358-1,t356+ . 1 xyz_dptempdt0.DSC.L2,t354+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t89+t358-1,t356+t91,t354+t93) . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t358-1,t356+1+ . 1 xyz_dptempdt0.DSC.L2,t354+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t89+t358-1,t356+1+t91,t354+t93) . xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t358-1,t356+2+ . 1 xyz_dptempdt0.DSC.L2,t354+xyz_dptempdt0.DSC.L3) = . 2 xyz_dptempdt(t89+t358-1,t356+2+t91,t354+t93) . end do . end do . endif 171 xyzf_DQMixDt0 = xyzf_DQMixDt . if (xyzf_dqmixdt0.DSC.U2 + 1 - xyzf_dqmixdt0.DSC.L2 .gt. 0) then . J2 = and(xyzf_dqmixdt0.DSC.U2 + 1 - xyzf_dqmixdt0.DSC.L2,3) . !CDIR NODEP . do t370 = 1, J2 . !CDIR NODEP . do t372 = 1, xyzf_dqmixdt0.DSC.U1 + 1 - xyzf_dqmixdt0.DSC.L1 . xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t372-1,t370-1+ . 1 xyzf_dqmixdt0.DSC.L2,t368+xyzf_dqmixdt0.DSC.L3,t366+1) . 2 = xyzf_dqmixdt(t99+t372-1,t370-1+t101,t368+t103,t366+ . 3 1) . end do . end do . !CDIR NODEP . do t370 = J2 + 1, xyzf_dqmixdt0.DSC.U2 + 1 - . 1 xyzf_dqmixdt0.DSC.L2, 4 . !CDIR NODEP . do t372 = 1, xyzf_dqmixdt0.DSC.U1 + 1 - xyzf_dqmixdt0.DSC.L1 . xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t372-1,t370-1+ . 1 xyzf_dqmixdt0.DSC.L2,t368+xyzf_dqmixdt0.DSC.L3,t366+1) . 2 = xyzf_dqmixdt(t99+t372-1,t370-1+t101,t368+t103,t366+ . 3 1) . xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t372-1,t370+ . 1 xyzf_dqmixdt0.DSC.L2,t368+xyzf_dqmixdt0.DSC.L3,t366+1) . 2 = xyzf_dqmixdt(t99+t372-1,t370+t101,t368+t103,t366+1) . xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t372-1,t370+1+ . 1 xyzf_dqmixdt0.DSC.L2,t368+xyzf_dqmixdt0.DSC.L3,t366+1) . 2 = xyzf_dqmixdt(t99+t372-1,t370+1+t101,t368+t103,t366+ . 3 1) . xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t372-1,t370+2+ . 1 xyzf_dqmixdt0.DSC.L2,t368+xyzf_dqmixdt0.DSC.L3,t366+1) . 2 = xyzf_dqmixdt(t99+t372-1,t370+2+t101,t368+t103,t366+ . 3 1) . end do . end do . endif 172 173 !地表面熱フラックスによる加熱率を計算 174 ! * 単位は K/s 175 ! * エクスナー関数は基本場の値で代表させる. 176 ! * 格子点 xz では, 物理領域の最下端の添え字は kz = 1 177 178 xyz_HeatFlux(:,:,kz) = & . if (jmax + 1 - jmin .gt. 0) then . J3 = and(jmax + 1 - jmin,3) . !CDIR NODEP . do t382 = 1, J3 . D2 = 1.D0/(z_dz(1)*5.00000000000000e-001)** . 1 2.00000000000000e+000 . !CDIR NODEP . do t384 = 1, imax + 1 - imin . xyz_heatflux(xyz_heatflux.DSC.L1+t384-1,t382-1+ . 1 xyz_heatflux.DSC.L2,1) = -kappa*xyz_ptemp(t67+t384-1, . 2 t382-1+t69,1)*xyz_exnerbz(xyz_exnerbz.DSC.L1+t384-1, . 3 t382-1+xyz_exnerbz.DSC.L2,1)*D2 . end do . end do . !CDIR NODEP . do t382 = J3 + 1, jmax + 1 - jmin, 4 . D1 = z_dz(1) . !CDIR NODEP . do t384 = 1, imax + 1 - imin . xyz_heatflux(xyz_heatflux.DSC.L1+t384-1,t382-1+ . 1 xyz_heatflux.DSC.L2,1) = -kappa*xyz_ptemp(t67+t384-1, . 2 t382-1+t69,1)*xyz_exnerbz(xyz_exnerbz.DSC.L1+t384-1, . 3 t382-1+xyz_exnerbz.DSC.L2,1)/(D1*5.00000000000000e-001 . 4 )**2.00000000000000e+000 . xyz_heatflux(xyz_heatflux.DSC.L1+t384-1,t382+ . 1 xyz_heatflux.DSC.L2,1) = -kappa*xyz_ptemp(t67+t384-1, . 2 t382+t69,1)*xyz_exnerbz(xyz_exnerbz.DSC.L1+t384-1,t382 . 3 +xyz_exnerbz.DSC.L2,1)/(D1*5.00000000000000e-001)** . 4 2.00000000000000e+000 . xyz_heatflux(xyz_heatflux.DSC.L1+t384-1,t382+1+ . 1 xyz_heatflux.DSC.L2,1) = -kappa*xyz_ptemp(t67+t384-1, . 2 t382+1+t69,1)*xyz_exnerbz(xyz_exnerbz.DSC.L1+t384-1, . 3 t382+1+xyz_exnerbz.DSC.L2,1)/(D1*5.00000000000000e-001 . 4 )**2.00000000000000e+000 . xyz_heatflux(xyz_heatflux.DSC.L1+t384-1,t382+2+ . 1 xyz_heatflux.DSC.L2,1) = -kappa*xyz_ptemp(t67+t384-1, . 2 t382+2+t69,1)*xyz_exnerbz(xyz_exnerbz.DSC.L1+t384-1, . 3 t382+2+xyz_exnerbz.DSC.L2,1)/(D1*5.00000000000000e-001 . 4 )**2.00000000000000e+000 . end do . end do . endif 179 & - Kappa * xyz_PTemp(:,:,kz) * xyz_ExnerBZ(:,:,kz) & 180 & / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) 181 182 do l = 1, GasNum 183 xyzf_QMixFlux(:,:,kz,l) = & . if (jmax + 1 - jmin .gt. 0) then . J4 = and(jmax + 1 - jmin,3) . !CDIR NODEP . do t392 = 1, J4 . D4 = kappa/(z_dz(1)*5.00000000000000e-001)** . 1 2.00000000000000e+000 . !CDIR NODEP . do t394 = 1, imax + 1 - imin . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t394-1,t392-1+ . 1 xyzf_qmixflux.DSC.L2,1,l) = max( . 2 0.0000000000000000e+000,(-xyzf_qmix(t77+t394-1,t392-1+ . 3 t79,1,l)*D4)) . end do . end do . !CDIR NODEP . do t392 = J4 + 1, jmax + 1 - jmin, 4 . D3 = z_dz(1) . !CDIR NODEP . do t394 = 1, imax + 1 - imin . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t394-1,t392-1+ . 1 xyzf_qmixflux.DSC.L2,1,l) = max( . 2 0.0000000000000000e+000,(-kappa*xyzf_qmix(t77+t394-1, . 3 t392-1+t79,1,l)/(D3*5.00000000000000e-001)** . 4 2.00000000000000e+000)) . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t394-1,t392+ . 1 xyzf_qmixflux.DSC.L2,1,l) = max( . 2 0.0000000000000000e+000,(-kappa*xyzf_qmix(t77+t394-1, . 3 t392+t79,1,l)/(D3*5.00000000000000e-001)** . 4 2.00000000000000e+000)) . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t394-1,t392+1+ . 1 xyzf_qmixflux.DSC.L2,1,l) = max( . 2 0.0000000000000000e+000,(-kappa*xyzf_qmix(t77+t394-1, . 3 t392+1+t79,1,l)/(D3*5.00000000000000e-001)** . 4 2.00000000000000e+000)) . xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t394-1,t392+2+ . 1 xyzf_qmixflux.DSC.L2,1,l) = max( . 2 0.0000000000000000e+000,(-kappa*xyzf_qmix(t77+t394-1, . 3 t392+2+t79,1,l)/(D3*5.00000000000000e-001)** . 4 2.00000000000000e+000)) . end do . end do . endif 184 & max( & 185 & 0.0d0, & 186 & - Kappa * xyzf_QMix(:,:,kz,l) & 187 & / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) & 188 & ) 189 end do 190 191 xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux . if (xyz_dptempdt0.DSC.U2 + 1 - xyz_dptempdt0.DSC.L2 .gt. 0) then . J5 = and(xyz_dptempdt0.DSC.U2 + 1 - xyz_dptempdt0.DSC.L2,3) . !CDIR NODEP . do t402 = 1, J5 . !CDIR NODEP . do t404 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt(t89+t404-1,t402-1+t91,t400+t93) = . 1 xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t404-1,t402-1+ . 2 xyz_dptempdt0.DSC.L2,t400+xyz_dptempdt0.DSC.L3) + . 3 xyz_heatflux(xyz_heatflux.DSC.L1+t404-1,t402-1+ . 4 xyz_heatflux.DSC.L2,t400+xyz_heatflux.DSC.L3) . end do . end do . !CDIR NODEP . do t402 = J5 + 1, xyz_dptempdt0.DSC.U2 + 1 - . 1 xyz_dptempdt0.DSC.L2, 4 . !CDIR NODEP . do t404 = 1, xyz_dptempdt0.DSC.U1 + 1 - xyz_dptempdt0.DSC.L1 . xyz_dptempdt(t89+t404-1,t402-1+t91,t400+t93) = . 1 xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t404-1,t402-1+ . 2 xyz_dptempdt0.DSC.L2,t400+xyz_dptempdt0.DSC.L3) + . 3 xyz_heatflux(xyz_heatflux.DSC.L1+t404-1,t402-1+ . 4 xyz_heatflux.DSC.L2,t400+xyz_heatflux.DSC.L3) . xyz_dptempdt(t89+t404-1,t402+t91,t400+t93) = . 1 xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t404-1,t402+ . 2 xyz_dptempdt0.DSC.L2,t400+xyz_dptempdt0.DSC.L3) + . 3 xyz_heatflux(xyz_heatflux.DSC.L1+t404-1,t402+ . 4 xyz_heatflux.DSC.L2,t400+xyz_heatflux.DSC.L3) . xyz_dptempdt(t89+t404-1,t402+1+t91,t400+t93) = . 1 xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t404-1,t402+1+ . 2 xyz_dptempdt0.DSC.L2,t400+xyz_dptempdt0.DSC.L3) + . 3 xyz_heatflux(xyz_heatflux.DSC.L1+t404-1,t402+1+ . 4 xyz_heatflux.DSC.L2,t400+xyz_heatflux.DSC.L3) . xyz_dptempdt(t89+t404-1,t402+2+t91,t400+t93) = . 1 xyz_dptempdt0(xyz_dptempdt0.DSC.L1+t404-1,t402+2+ . 2 xyz_dptempdt0.DSC.L2,t400+xyz_dptempdt0.DSC.L3) + . 3 xyz_heatflux(xyz_heatflux.DSC.L1+t404-1,t402+2+ . 4 xyz_heatflux.DSC.L2,t400+xyz_heatflux.DSC.L3) . end do . end do . endif 192 xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux . if (xyzf_dqmixdt0.DSC.U2 + 1 - xyzf_dqmixdt0.DSC.L2 .gt. 0) then . J6 = and(xyzf_dqmixdt0.DSC.U2 + 1 - xyzf_dqmixdt0.DSC.L2,3) . !CDIR NODEP . do t419 = 1, J6 . !CDIR NODEP . do t421 = 1, xyzf_dqmixdt0.DSC.U1 + 1 - xyzf_dqmixdt0.DSC.L1 . xyzf_dqmixdt(t99+t421-1,t419-1+t101,t417+t103,t415+1) = . 1 xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t421-1,t419-1+ . 2 xyzf_dqmixdt0.DSC.L2,t417+xyzf_dqmixdt0.DSC.L3,t415+1) . 3 + xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t421-1,t419-1+ . 4 xyzf_qmixflux.DSC.L2,t417+xyzf_qmixflux.DSC.L3,t415+1) . end do . end do . !CDIR NODEP . do t419 = J6 + 1, xyzf_dqmixdt0.DSC.U2 + 1 - . 1 xyzf_dqmixdt0.DSC.L2, 4 . !CDIR NODEP . do t421 = 1, xyzf_dqmixdt0.DSC.U1 + 1 - xyzf_dqmixdt0.DSC.L1 . xyzf_dqmixdt(t99+t421-1,t419-1+t101,t417+t103,t415+1) = . 1 xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t421-1,t419-1+ . 2 xyzf_dqmixdt0.DSC.L2,t417+xyzf_dqmixdt0.DSC.L3,t415+1) . 3 + xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t421-1,t419-1+ . 4 xyzf_qmixflux.DSC.L2,t417+xyzf_qmixflux.DSC.L3,t415+1) . xyzf_dqmixdt(t99+t421-1,t419+t101,t417+t103,t415+1) = . 1 xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t421-1,t419+ . 2 xyzf_dqmixdt0.DSC.L2,t417+xyzf_dqmixdt0.DSC.L3,t415+1) . 3 + xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t421-1,t419+ . 4 xyzf_qmixflux.DSC.L2,t417+xyzf_qmixflux.DSC.L3,t415+1) . xyzf_dqmixdt(t99+t421-1,t419+1+t101,t417+t103,t415+1) = . 1 xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t421-1,t419+1+ . 2 xyzf_dqmixdt0.DSC.L2,t417+xyzf_dqmixdt0.DSC.L3,t415+1) . 3 + xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t421-1,t419+1+ . 4 xyzf_qmixflux.DSC.L2,t417+xyzf_qmixflux.DSC.L3,t415+1) . xyzf_dqmixdt(t99+t421-1,t419+2+t101,t417+t103,t415+1) = . 1 xyzf_dqmixdt0(xyzf_dqmixdt0.DSC.L1+t421-1,t419+2+ . 2 xyzf_dqmixdt0.DSC.L2,t417+xyzf_dqmixdt0.DSC.L3,t415+1) . 3 + xyzf_qmixflux(xyzf_qmixflux.DSC.L1+t421-1,t419+2+ . 4 xyzf_qmixflux.DSC.L2,t417+xyzf_qmixflux.DSC.L3,t415+1) . end do . end do . endif 193 194 call HistoryAutoPut(TimeN, 'PTempFlux', xyz_HeatFlux(1:nx,1:ny,1:nz)) 195 do l = 1, ncmax 196 call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Flux', xyzf_Qmixflux(1:nx,1:ny,1:nz,l)) 197 end do 198 199 ! Set margin 200 ! 201 call SetMargin_xyz(xyz_DPTempDt) 202 call SetMargin_xyzf(xyzf_DQMixDt) 203 204 end subroutine Surfaceflux_Diff_forcing 205 206 end module Surfaceflux_diff Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:35:04 2011 FILE NAME: surfaceflux_diff.f90 PROGRAM NAME: surfaceflux_diff FORMAT LIST LINE LOOP FORTRAN STATEMENT 1: != Module HeatFlux 2: ! 3: ! Authors:: ODAKA Masatsugu 4: ! Version:: $Id: surfaceflux_diff.f90,v 1.8 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: ! 基づいて計算するモジュール. これは中島 (1994) で用いられた方法である. 13: ! 14: ! 熱フラックスを Fh と凝結物質のフラックス Fq とすると, 温度と凝結物質 15: ! の変化率 H, Q は 16: ! 17: ! H = Fh/Δz_1 18: ! Q = Fq/Δz_1 19: ! 20: ! と表される. ここで Δz_1 は最下層の格子間隔である. 21: ! 22: ! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって 23: ! 計算する. 24: ! 25: ! Fh = - Cdρ|V| * (π_1θ_1 - T_sfc) 26: ! Fh = - Cdρ|V| * (Q_1 - Q*(T_sfc)) 27: ! 28: ! ここで θ_1, Q_1, π_1 は最下層の温度と凝結物質の混合比および無次元 29: ! 圧力関数, T_sfc は下部境界の温度, Q*(T_sfc) は T_sfc で決まる飽和混 30: ! 合比である. 31: ! 32: ! バルク係数 Cd は一定とする. 無次元圧力関数 π_1 は基本場の値を用いる. 33: ! 風速値 |V| は 34: ! 35: ! V = ( V^2 + V_0^2 )^(1/2) 36: ! 37: ! と計算する. 38: ! 39: !== Error Handling 40: ! 41: !== Bugs 42: ! 43: !== Note 44: ! 45: ! 46: !== Future Plans 47: ! 48: ! 49: 50: module Surfaceflux_diff 51: ! 52: !下部境界でのフラックスの計算モジュール 53: ! 54: 55: !モジュール読み込み 56: use dc_types, only: DP, STRING 57: use dc_iounit, only: FileOpen 58: use dc_message, only: MessageNotify 59: use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut 60: 61: use mpi_wrapper,only: myrank 62: use gridset, only: imin, & !x 方向の配列の下限 63: & imax, & !x 方向の配列の上限 64: & jmin, & !y 方向の配列の下限 65: & jmax, & !y 方向の配列の上限 66: & kmin, & !z 方向の配列の下限 67: & kmax, & !z 方向の配列の上限 68: & nx, ny, nz, ncmax 69: use axesset, only: z_dz !z 方向の格子点間隔 70: use basicset, only: xyz_ExnerBZ !エクスナー関数の基本場 71: use composition,only : GasNum, SpcWetSymbol 72: use namelist_util, only: namelist_filename 73: use timeset, only: TimeN 74: use setmargin,only: SetMargin_xyz, SetMargin_xyzf 75: 76: !暗黙の型宣言禁止 77: implicit none 78: 79: !属性の指定 80: private 81: 82: !関数を public に設定 83: public surfaceflux_diff_init 84: public surfaceflux_diff_forcing 85: 86: !変数定義 87: real(DP), save :: Kappa = 800.0d0 88: 89: contains 90: !!!------------------------------------------------------------------------!!! 91: subroutine Surfaceflux_Diff_init 92: ! 93: !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 94: ! 95: 96: !暗黙の型宣言禁止 97: implicit none 98: 99: !内部変数 100: integer :: l, unit 101: 102: !--------------------------------------------------------------- 103: ! NAMELIST から情報を取得 104: ! 105: NAMELIST /surfaceflux_diff_nml/ Kappa 106: 107: call FileOpen(unit, file=namelist_filename, mode='r') 108: read(unit, NML=surfaceflux_diff_nml) 109: close(unit) 110: 111: if (myrank == 0) then 112: call MessageNotify( "M", "SurfaceFlux", "Kappa = %f", d=(/Kappa/)) 113: end if 114: 115: call HistoryAutoAddVariable( & 116: & varname='PTempFlux', & 117: & dims=(/'x','y','z','t'/), & 118: & longname='surface flux of potential temperature', & 119: & units='kg.kg-1.s-1', & 120: & xtype='float') 121: 122: +------> do l = 1, ncmax 123: | call HistoryAutoAddVariable( & 124: | & varname=trim(SpcWetSymbol(l))//'_Flux', & 125: | & dims=(/'x','y','z','t'/), & 126: | & longname='Surface Flux term of ' & 127: | & //trim(SpcWetSymbol(l))//' mixing ratio', & 128: | & units='kg.kg-1.s-1', & 129: | & xtype='float') 130: +------ end do 131: 132: end subroutine Surfaceflux_Diff_init 133: 134: 135: !!!------------------------------------------------------------------------!!! 136: subroutine Surfaceflux_Diff_forcing( & 137: & xyz_PTemp, xyzf_QMix, & 138: & xyz_DPTempDt, xyzf_DQMixDt & 139: & ) 140: ! 141: ! 下部境界からのフラックスによる温度の変化率を, 142: ! バルク方法に基づいて計算する. 143: ! 144: 145: !暗黙の型宣言禁止 146: implicit none 147: 148: !変数定義 149: real(DP), intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax) 150: !温位の擾乱成分 151: real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) 152: !温位の擾乱成分 153: real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax) 154: real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) 155: real(DP) :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax) 156: real(DP) :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax) 157: real(DP) :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax) 158: !地表面熱フラックス 159: real(DP) :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax) 160: 161: integer :: kz !配列添字 162: integer :: l !ループ変数 163: 164: ! 初期化 165: ! 166: kz = 1 167: WW+==== xyz_Heatflux = 0.0d0 168: ++++=== xyzf_QMixflux = 0.0d0 169: 170: ++V==== xyz_DPTempDt0 = xyz_DPTempDt 171: +++V=== xyzf_DQMixDt0 = xyzf_DQMixDt 172: 173: !地表面熱フラックスによる加熱率を計算 174: ! * 単位は K/s 175: ! * エクスナー関数は基本場の値で代表させる. 176: ! * 格子点 xz では, 物理領域の最下端の添え字は kz = 1 177: 178: +V===== xyz_HeatFlux(:,:,kz) = & 179: & - Kappa * xyz_PTemp(:,:,kz) * xyz_ExnerBZ(:,:,kz) & 180: & / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) 181: 182: +------> do l = 1, GasNum 183: |+V==== xyzf_QMixFlux(:,:,kz,l) = & 184: | & max( & 185: | & 0.0d0, & 186: | & - Kappa * xyzf_QMix(:,:,kz,l) & 187: | & / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) & 188: | & ) 189: +------ end do 190: 191: ++V==== xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux 192: +++V=== xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux 193: 194: call HistoryAutoPut(TimeN, 'PTempFlux', xyz_HeatFlux(1:nx,1:ny,1:nz)) 195: +------> do l = 1, ncmax 196: | call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Flux', xyzf_Qmixflux(1:nx,1:ny,1:nz,l)) 197: +------ end do 198: 199: ! Set margin 200: ! 201: call SetMargin_xyz(xyz_DPTempDt) 202: call SetMargin_xyzf(xyzf_DQMixDt) 203: 204: end subroutine Surfaceflux_Diff_forcing 205: 206: end module Surfaceflux_diff