| 1 | |
|---|
| 2 | ! $Header$ |
|---|
| 3 | |
|---|
| 4 | SUBROUTINE diagphy(airephy, tit, iprt, tops, topl, sols, soll, sens, evap, & |
|---|
| 5 | rain_fall, snow_fall, ts, d_etp_tot, d_qt_tot, d_ec_tot, fs_bound, & |
|---|
| 6 | fq_bound) |
|---|
| 7 | ! ====================================================================== |
|---|
| 8 | |
|---|
| 9 | ! Purpose: |
|---|
| 10 | ! Compute the thermal flux and the watter mass flux at the atmosphere |
|---|
| 11 | ! boundaries. Print them and also the atmospheric enthalpy change and |
|---|
| 12 | ! the atmospheric mass change. |
|---|
| 13 | |
|---|
| 14 | ! Arguments: |
|---|
| 15 | ! airephy-------input-R- grid area |
|---|
| 16 | ! tit---------input-A15- Comment to be added in PRINT (CHARACTER*15) |
|---|
| 17 | ! iprt--------input-I- PRINT level ( <=0 : no PRINT) |
|---|
| 18 | ! tops(klon)--input-R- SW rad. at TOA (W/m2), positive up. |
|---|
| 19 | ! topl(klon)--input-R- LW rad. at TOA (W/m2), positive down |
|---|
| 20 | ! sols(klon)--input-R- Net SW flux above surface (W/m2), positive up |
|---|
| 21 | ! (i.e. -1 * flux absorbed by the surface) |
|---|
| 22 | ! soll(klon)--input-R- Net LW flux above surface (W/m2), positive up |
|---|
| 23 | ! (i.e. flux emited - flux absorbed by the surface) |
|---|
| 24 | ! sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down |
|---|
| 25 | ! evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux |
|---|
| 26 | ! (kg/m2/s), positive up |
|---|
| 27 | ! rain_fall(klon) |
|---|
| 28 | ! --input-R- Liquid watter mass flux (kg/m2/s), positive down |
|---|
| 29 | ! snow_fall(klon) |
|---|
| 30 | ! --input-R- Solid watter mass flux (kg/m2/s), positive down |
|---|
| 31 | ! ts(klon)----input-R- Surface temperature (K) |
|---|
| 32 | ! d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy |
|---|
| 33 | ! change (W/m2) |
|---|
| 34 | ! d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass |
|---|
| 35 | ! change (kg/m2/s) |
|---|
| 36 | ! d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy |
|---|
| 37 | ! change (W/m2) |
|---|
| 38 | |
|---|
| 39 | ! fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2) |
|---|
| 40 | ! fq_bound---output-R- Watter mass flux at the atmosphere boundaries |
|---|
| 41 | ! (kg/m2/s) |
|---|
| 42 | |
|---|
| 43 | ! J.L. Dufresne, July 2002 |
|---|
| 44 | ! Version prise sur |
|---|
| 45 | ! ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd |
|---|
| 46 | ! le 25 Novembre 2002. |
|---|
| 47 | ! ====================================================================== |
|---|
| 48 | |
|---|
| 49 | USE dimphy |
|---|
| 50 | IMPLICIT NONE |
|---|
| 51 | |
|---|
| 52 | include "YOMCST.h" |
|---|
| 53 | include "YOETHF.h" |
|---|
| 54 | |
|---|
| 55 | ! Input variables |
|---|
| 56 | REAL airephy(klon) |
|---|
| 57 | CHARACTER *15 tit |
|---|
| 58 | INTEGER iprt |
|---|
| 59 | REAL tops(klon), topl(klon), sols(klon), soll(klon) |
|---|
| 60 | REAL sens(klon), evap(klon), rain_fall(klon), snow_fall(klon) |
|---|
| 61 | REAL ts(klon) |
|---|
| 62 | REAL d_etp_tot, d_qt_tot, d_ec_tot |
|---|
| 63 | ! Output variables |
|---|
| 64 | REAL fs_bound, fq_bound |
|---|
| 65 | |
|---|
| 66 | ! Local variables |
|---|
| 67 | REAL stops, stopl, ssols, ssoll |
|---|
| 68 | REAL ssens, sfront, slat |
|---|
| 69 | REAL airetot, zcpvap, zcwat, zcice |
|---|
| 70 | REAL rain_fall_tot, snow_fall_tot, evap_tot |
|---|
| 71 | |
|---|
| 72 | INTEGER i |
|---|
| 73 | |
|---|
| 74 | INTEGER pas |
|---|
| 75 | SAVE pas |
|---|
| 76 | DATA pas/0/ |
|---|
| 77 | !$OMP THREADPRIVATE(pas) |
|---|
| 78 | |
|---|
| 79 | pas = pas + 1 |
|---|
| 80 | stops = 0. |
|---|
| 81 | stopl = 0. |
|---|
| 82 | ssols = 0. |
|---|
| 83 | ssoll = 0. |
|---|
| 84 | ssens = 0. |
|---|
| 85 | sfront = 0. |
|---|
| 86 | evap_tot = 0. |
|---|
| 87 | rain_fall_tot = 0. |
|---|
| 88 | snow_fall_tot = 0. |
|---|
| 89 | airetot = 0. |
|---|
| 90 | |
|---|
| 91 | ! Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de |
|---|
| 92 | ! la glace, on travaille par difference a la chaleur specifique de l' |
|---|
| 93 | ! air sec. En effet, comme on travaille a niveau de pression donne, |
|---|
| 94 | ! toute variation de la masse d'un constituant est totalement |
|---|
| 95 | ! compense par une variation de masse d'air. |
|---|
| 96 | |
|---|
| 97 | zcpvap = rcpv - rcpd |
|---|
| 98 | zcwat = rcw - rcpd |
|---|
| 99 | zcice = rcs - rcpd |
|---|
| 100 | |
|---|
| 101 | DO i = 1, klon |
|---|
| 102 | stops = stops + tops(i)*airephy(i) |
|---|
| 103 | stopl = stopl + topl(i)*airephy(i) |
|---|
| 104 | ssols = ssols + sols(i)*airephy(i) |
|---|
| 105 | ssoll = ssoll + soll(i)*airephy(i) |
|---|
| 106 | ssens = ssens + sens(i)*airephy(i) |
|---|
| 107 | sfront = sfront + (evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice)* & |
|---|
| 108 | ts(i)*airephy(i) |
|---|
| 109 | evap_tot = evap_tot + evap(i)*airephy(i) |
|---|
| 110 | rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i) |
|---|
| 111 | snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i) |
|---|
| 112 | airetot = airetot + airephy(i) |
|---|
| 113 | END DO |
|---|
| 114 | stops = stops/airetot |
|---|
| 115 | stopl = stopl/airetot |
|---|
| 116 | ssols = ssols/airetot |
|---|
| 117 | ssoll = ssoll/airetot |
|---|
| 118 | ssens = ssens/airetot |
|---|
| 119 | sfront = sfront/airetot |
|---|
| 120 | evap_tot = evap_tot/airetot |
|---|
| 121 | rain_fall_tot = rain_fall_tot/airetot |
|---|
| 122 | snow_fall_tot = snow_fall_tot/airetot |
|---|
| 123 | |
|---|
| 124 | slat = rlvtt*rain_fall_tot + rlstt*snow_fall_tot |
|---|
| 125 | ! Heat flux at atm. boundaries |
|---|
| 126 | fs_bound = stops - stopl - (ssols+ssoll) + ssens + sfront + slat |
|---|
| 127 | ! Watter flux at atm. boundaries |
|---|
| 128 | fq_bound = evap_tot - rain_fall_tot - snow_fall_tot |
|---|
| 129 | |
|---|
| 130 | IF (iprt>=1) WRITE (6, 6666) tit, pas, fs_bound, d_etp_tot, fq_bound, & |
|---|
| 131 | d_qt_tot |
|---|
| 132 | |
|---|
| 133 | IF (iprt>=1) WRITE (6, 6668) tit, pas, d_etp_tot + d_ec_tot - fs_bound, & |
|---|
| 134 | d_qt_tot - fq_bound |
|---|
| 135 | |
|---|
| 136 | IF (iprt>=2) WRITE (6, 6667) tit, pas, stops, stopl, ssols, ssoll, ssens, & |
|---|
| 137 | slat, evap_tot, rain_fall_tot + snow_fall_tot |
|---|
| 138 | |
|---|
| 139 | RETURN |
|---|
| 140 | |
|---|
| 141 | 6666 FORMAT ('Phys. Flux Budget ', A15, 1I6, 2F8.2, 2(1PE13.5)) |
|---|
| 142 | 6667 FORMAT ('Phys. Boundary Flux ', A15, 1I6, 6F8.2, 2(1PE13.5)) |
|---|
| 143 | 6668 FORMAT ('Phys. Total Budget ', A15, 1I6, F8.2, 2(1PE13.5)) |
|---|
| 144 | |
|---|
| 145 | END SUBROUTINE diagphy |
|---|
| 146 | |
|---|
| 147 | ! ====================================================================== |
|---|
| 148 | SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, & |
|---|
| 149 | u, v, paprs, pplay, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec & |
|---|
| 150 | !#ifdef ISO |
|---|
| 151 | ! & , xt,xtl,xts,d_xtw, d_xtl, d_xts & |
|---|
| 152 | !#endif |
|---|
| 153 | & ) |
|---|
| 154 | ! ====================================================================== |
|---|
| 155 | |
|---|
| 156 | ! Purpose: |
|---|
| 157 | ! Calcul la difference d'enthalpie et de masse d'eau entre 2 appels, |
|---|
| 158 | ! et calcul le flux de chaleur et le flux d'eau necessaire a ces |
|---|
| 159 | ! changements. Ces valeurs sont moyennees sur la surface de tout |
|---|
| 160 | ! le globe et sont exprime en W/2 et kg/s/m2 |
|---|
| 161 | ! Outil pour diagnostiquer la conservation de l'energie |
|---|
| 162 | ! et de la masse dans la physique. Suppose que les niveau de |
|---|
| 163 | ! pression entre couche ne varie pas entre 2 appels. |
|---|
| 164 | |
|---|
| 165 | ! Plusieurs de ces diagnostics peuvent etre fait en parallele: les |
|---|
| 166 | ! bilans sont sauvegardes dans des tableaux indices. On parlera |
|---|
| 167 | ! "d'indice de diagnostic" |
|---|
| 168 | |
|---|
| 169 | |
|---|
| 170 | ! ====================================================================== |
|---|
| 171 | ! Arguments: |
|---|
| 172 | ! airephy-------input-R- grid area |
|---|
| 173 | ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15) |
|---|
| 174 | ! iprt----input-I- PRINT level ( <=1 : no PRINT) |
|---|
| 175 | ! idiag---input-I- indice dans lequel sera range les nouveaux |
|---|
| 176 | ! bilans d' entalpie et de masse |
|---|
| 177 | ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse |
|---|
| 178 | ! sont compare au bilan de d'enthalpie de masse de |
|---|
| 179 | ! l'indice numero idiag2 |
|---|
| 180 | ! Cas parriculier : si idiag2=0, pas de comparaison, on |
|---|
| 181 | ! sort directement les bilans d'enthalpie et de masse |
|---|
| 182 | ! dtime----input-R- time step (s) |
|---|
| 183 | ! t--------input-R- temperature (K) |
|---|
| 184 | ! q--------input-R- vapeur d'eau (kg/kg) |
|---|
| 185 | ! ql-------input-R- liquid watter (kg/kg) |
|---|
| 186 | ! qs-------input-R- solid watter (kg/kg) |
|---|
| 187 | ! u--------input-R- vitesse u |
|---|
| 188 | ! v--------input-R- vitesse v |
|---|
| 189 | ! paprs----input-R- pression a intercouche (Pa) |
|---|
| 190 | ! pplay----input-R- pression au milieu de couche (Pa) |
|---|
| 191 | |
|---|
| 192 | ! the following total value are computed by UNIT of earth surface |
|---|
| 193 | |
|---|
| 194 | ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy |
|---|
| 195 | ! change (J/m2) during one time step (dtime) for the whole |
|---|
| 196 | ! atmosphere (air, watter vapour, liquid and solid) |
|---|
| 197 | ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the |
|---|
| 198 | ! total watter (kg/m2) change during one time step (dtime), |
|---|
| 199 | ! d_qw------output-R- same, for the watter vapour only (kg/m2/s) |
|---|
| 200 | ! d_ql------output-R- same, for the liquid watter only (kg/m2/s) |
|---|
| 201 | ! d_qs------output-R- same, for the solid watter only (kg/m2/s) |
|---|
| 202 | ! d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column |
|---|
| 203 | |
|---|
| 204 | ! other (COMMON...) |
|---|
| 205 | ! RCPD, RCPV, .... |
|---|
| 206 | |
|---|
| 207 | ! J.L. Dufresne, July 2002 |
|---|
| 208 | ! ====================================================================== |
|---|
| 209 | |
|---|
| 210 | USE dimphy |
|---|
| 211 | !#ifdef ISO |
|---|
| 212 | ! USE infotrac_phy, ONLY: niso,ntraciso |
|---|
| 213 | !#endif |
|---|
| 214 | ! isotopes inutiles ici |
|---|
| 215 | IMPLICIT NONE |
|---|
| 216 | |
|---|
| 217 | include "YOMCST.h" |
|---|
| 218 | include "YOETHF.h" |
|---|
| 219 | |
|---|
| 220 | ! Input variables |
|---|
| 221 | REAL airephy(klon) |
|---|
| 222 | CHARACTER *15 tit |
|---|
| 223 | INTEGER iprt, idiag, idiag2 |
|---|
| 224 | REAL dtime |
|---|
| 225 | REAL t(klon, klev), q(klon, klev), ql(klon, klev), qs(klon, klev) |
|---|
| 226 | REAL u(klon, klev), v(klon, klev) |
|---|
| 227 | REAL paprs(klon, klev+1), pplay(klon, klev) |
|---|
| 228 | ! Output variables |
|---|
| 229 | REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec |
|---|
| 230 | |
|---|
| 231 | ! Local variables |
|---|
| 232 | |
|---|
| 233 | REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot, & |
|---|
| 234 | qs_tot, ec_tot |
|---|
| 235 | ! h_vcol_tot-- total enthalpy of vertical air column |
|---|
| 236 | ! (air with watter vapour, liquid and solid) (J/m2) |
|---|
| 237 | ! h_dair_tot-- total enthalpy of dry air (J/m2) |
|---|
| 238 | ! h_qw_tot---- total enthalpy of watter vapour (J/m2) |
|---|
| 239 | ! h_ql_tot---- total enthalpy of liquid watter (J/m2) |
|---|
| 240 | ! h_qs_tot---- total enthalpy of solid watter (J/m2) |
|---|
| 241 | ! qw_tot------ total mass of watter vapour (kg/m2) |
|---|
| 242 | ! ql_tot------ total mass of liquid watter (kg/m2) |
|---|
| 243 | ! qs_tot------ total mass of solid watter (kg/m2) |
|---|
| 244 | ! ec_tot------ total cinetic energy (kg/m2) |
|---|
| 245 | |
|---|
| 246 | REAL zairm(klon, klev) ! layer air mass (kg/m2) |
|---|
| 247 | REAL zqw_col(klon) |
|---|
| 248 | REAL zql_col(klon) |
|---|
| 249 | REAL zqs_col(klon) |
|---|
| 250 | REAL zec_col(klon) |
|---|
| 251 | REAL zh_dair_col(klon) |
|---|
| 252 | REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon) |
|---|
| 253 | |
|---|
| 254 | !#ifdef ISO |
|---|
| 255 | ! real xt(ntraciso,klon,klev), xtl(ntraciso,klon,klev), & |
|---|
| 256 | ! & xts(ntraciso,klon,klev) |
|---|
| 257 | ! real xtw_tot(ntraciso), xtl_tot(ntraciso), xts_tot(ntraciso) |
|---|
| 258 | ! real d_xtw(ntraciso), d_xtl(ntraciso), d_xts(ntraciso) |
|---|
| 259 | ! REAL zxtw_col(ntraciso,klon) |
|---|
| 260 | ! REAL zxtl_col(ntraciso,klon) |
|---|
| 261 | ! REAL zxts_col(ntraciso,klon) |
|---|
| 262 | ! integer ixt |
|---|
| 263 | !#endif |
|---|
| 264 | REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs |
|---|
| 265 | |
|---|
| 266 | REAL airetot, zcpvap, zcwat, zcice |
|---|
| 267 | |
|---|
| 268 | INTEGER i, k |
|---|
| 269 | |
|---|
| 270 | INTEGER ndiag ! max number of diagnostic in parallel |
|---|
| 271 | PARAMETER (ndiag=10) |
|---|
| 272 | INTEGER pas(ndiag) |
|---|
| 273 | SAVE pas |
|---|
| 274 | DATA pas/ndiag*0/ |
|---|
| 275 | !$OMP THREADPRIVATE(pas) |
|---|
| 276 | |
|---|
| 277 | REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag), & |
|---|
| 278 | h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag), & |
|---|
| 279 | qs_pre(ndiag), ec_pre(ndiag) |
|---|
| 280 | SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre, h_qs_pre, qw_pre, ql_pre, & |
|---|
| 281 | qs_pre, ec_pre |
|---|
| 282 | !$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre) |
|---|
| 283 | !$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre) |
|---|
| 284 | !#ifdef ISO |
|---|
| 285 | ! REAL xtw_pre(ntraciso,ndiag), xtl_pre(ntraciso,ndiag), & |
|---|
| 286 | ! & xts_pre(ntraciso,ndiag) |
|---|
| 287 | ! SAVE xtw_pre, xtl_pre, xts_pre |
|---|
| 288 | ! !$OMP THREADPRIVATE(xtw_pre, xtl_pre, xts_pre) |
|---|
| 289 | !#endif |
|---|
| 290 | ! ====================================================================== |
|---|
| 291 | |
|---|
| 292 | #ifdef ISOVERIF |
|---|
| 293 | write(*,*) 'diagetp tmp 293' |
|---|
| 294 | #endif |
|---|
| 295 | DO k = 1, klev |
|---|
| 296 | DO i = 1, klon |
|---|
| 297 | ! layer air mass |
|---|
| 298 | zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg |
|---|
| 299 | END DO |
|---|
| 300 | END DO |
|---|
| 301 | |
|---|
| 302 | ! Reset variables |
|---|
| 303 | DO i = 1, klon |
|---|
| 304 | zqw_col(i) = 0. |
|---|
| 305 | zql_col(i) = 0. |
|---|
| 306 | zqs_col(i) = 0. |
|---|
| 307 | zec_col(i) = 0. |
|---|
| 308 | zh_dair_col(i) = 0. |
|---|
| 309 | zh_qw_col(i) = 0. |
|---|
| 310 | zh_ql_col(i) = 0. |
|---|
| 311 | zh_qs_col(i) = 0. |
|---|
| 312 | !#ifdef ISO |
|---|
| 313 | ! do ixt=1,ntraciso |
|---|
| 314 | ! zxtw_col(ixt,i)=0. |
|---|
| 315 | ! zxtl_col(ixt,i)=0. |
|---|
| 316 | ! zxts_col(ixt,i)=0. |
|---|
| 317 | ! enddo |
|---|
| 318 | !#endif |
|---|
| 319 | END DO |
|---|
| 320 | |
|---|
| 321 | zcpvap = rcpv |
|---|
| 322 | zcwat = rcw |
|---|
| 323 | zcice = rcs |
|---|
| 324 | |
|---|
| 325 | ! Compute vertical sum for each atmospheric column |
|---|
| 326 | ! ================================================ |
|---|
| 327 | DO k = 1, klev |
|---|
| 328 | DO i = 1, klon |
|---|
| 329 | ! Watter mass |
|---|
| 330 | zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k) |
|---|
| 331 | zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k) |
|---|
| 332 | zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k) |
|---|
| 333 | !#ifdef ISO |
|---|
| 334 | ! do ixt=1,ntraciso |
|---|
| 335 | ! zxtw_col(ixt,i) = zxtw_col(ixt,i) + xt(ixt,i,k)*zairm(i,k) |
|---|
| 336 | ! zxtl_col(ixt,i) = zxtl_col(ixt,i) + xtl(ixt,i,k)*zairm(i,k) |
|---|
| 337 | ! zxts_col(ixt,i) = zxts_col(ixt,i) + xts(ixt,i,k)*zairm(i,k) |
|---|
| 338 | ! enddo |
|---|
| 339 | !#endif |
|---|
| 340 | ! Cinetic Energy |
|---|
| 341 | zec_col(i) = zec_col(i) + 0.5*(u(i,k)**2+v(i,k)**2)*zairm(i, k) |
|---|
| 342 | ! Air enthalpy |
|---|
| 343 | zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-q(i,k)-ql(i,k)-qs(i,k))* & |
|---|
| 344 | zairm(i, k)*t(i, k) |
|---|
| 345 | zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k) |
|---|
| 346 | zh_ql_col(i) = zh_ql_col(i) + zcwat*ql(i, k)*zairm(i, k)*t(i, k) - & |
|---|
| 347 | rlvtt*ql(i, k)*zairm(i, k) |
|---|
| 348 | zh_qs_col(i) = zh_qs_col(i) + zcice*qs(i, k)*zairm(i, k)*t(i, k) - & |
|---|
| 349 | rlstt*qs(i, k)*zairm(i, k) |
|---|
| 350 | |
|---|
| 351 | END DO |
|---|
| 352 | END DO |
|---|
| 353 | |
|---|
| 354 | ! Mean over the planete surface |
|---|
| 355 | ! ============================= |
|---|
| 356 | qw_tot = 0. |
|---|
| 357 | ql_tot = 0. |
|---|
| 358 | qs_tot = 0. |
|---|
| 359 | ec_tot = 0. |
|---|
| 360 | h_vcol_tot = 0. |
|---|
| 361 | h_dair_tot = 0. |
|---|
| 362 | h_qw_tot = 0. |
|---|
| 363 | h_ql_tot = 0. |
|---|
| 364 | h_qs_tot = 0. |
|---|
| 365 | airetot = 0. |
|---|
| 366 | !#ifdef ISO |
|---|
| 367 | ! do ixt=1,ntraciso |
|---|
| 368 | ! xtw_tot(ixt) = 0. |
|---|
| 369 | ! xtl_tot(ixt) = 0. |
|---|
| 370 | ! xts_tot(ixt) = 0 |
|---|
| 371 | ! enddo |
|---|
| 372 | !#endif |
|---|
| 373 | |
|---|
| 374 | DO i = 1, klon |
|---|
| 375 | qw_tot = qw_tot + zqw_col(i)*airephy(i) |
|---|
| 376 | ql_tot = ql_tot + zql_col(i)*airephy(i) |
|---|
| 377 | qs_tot = qs_tot + zqs_col(i)*airephy(i) |
|---|
| 378 | ec_tot = ec_tot + zec_col(i)*airephy(i) |
|---|
| 379 | h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i) |
|---|
| 380 | h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i) |
|---|
| 381 | h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i) |
|---|
| 382 | h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i) |
|---|
| 383 | airetot = airetot + airephy(i) |
|---|
| 384 | !#ifdef ISO |
|---|
| 385 | ! do ixt=1,ntraciso |
|---|
| 386 | ! xtw_tot(ixt) = xtw_tot(ixt) + zxtw_col(ixt,i)*airephy(i) |
|---|
| 387 | ! xtl_tot(ixt) = xtl_tot(ixt) + zxtl_col(ixt,i)*airephy(i) |
|---|
| 388 | ! xts_tot(ixt) = xts_tot(ixt) + zxts_col(ixt,i)*airephy(i) |
|---|
| 389 | ! enddo |
|---|
| 390 | !#endif |
|---|
| 391 | END DO |
|---|
| 392 | |
|---|
| 393 | qw_tot = qw_tot/airetot |
|---|
| 394 | ql_tot = ql_tot/airetot |
|---|
| 395 | qs_tot = qs_tot/airetot |
|---|
| 396 | ec_tot = ec_tot/airetot |
|---|
| 397 | h_dair_tot = h_dair_tot/airetot |
|---|
| 398 | h_qw_tot = h_qw_tot/airetot |
|---|
| 399 | h_ql_tot = h_ql_tot/airetot |
|---|
| 400 | h_qs_tot = h_qs_tot/airetot |
|---|
| 401 | |
|---|
| 402 | h_vcol_tot = h_dair_tot + h_qw_tot + h_ql_tot + h_qs_tot |
|---|
| 403 | !#ifdef ISO |
|---|
| 404 | ! do ixt=1,ntraciso |
|---|
| 405 | ! xtw_tot(ixt) = xtw_tot(ixt)/airetot |
|---|
| 406 | ! xtl_tot(ixt) = xtl_tot(ixt)/airetot |
|---|
| 407 | ! xts_tot(ixt) = xts_tot(ixt)/airetot |
|---|
| 408 | ! enddo |
|---|
| 409 | !#endif |
|---|
| 410 | |
|---|
| 411 | |
|---|
| 412 | ! Compute the change of the atmospheric state compare to the one |
|---|
| 413 | ! stored in "idiag2", and convert it in flux. THis computation |
|---|
| 414 | ! is performed IF idiag2 /= 0 and IF it is not the first CALL |
|---|
| 415 | ! for "idiag" |
|---|
| 416 | ! =================================== |
|---|
| 417 | |
|---|
| 418 | IF ((idiag2>0) .AND. (pas(idiag2)/=0)) THEN |
|---|
| 419 | d_h_vcol = (h_vcol_tot-h_vcol_pre(idiag2))/dtime |
|---|
| 420 | d_h_dair = (h_dair_tot-h_dair_pre(idiag2))/dtime |
|---|
| 421 | d_h_qw = (h_qw_tot-h_qw_pre(idiag2))/dtime |
|---|
| 422 | d_h_ql = (h_ql_tot-h_ql_pre(idiag2))/dtime |
|---|
| 423 | d_h_qs = (h_qs_tot-h_qs_pre(idiag2))/dtime |
|---|
| 424 | d_qw = (qw_tot-qw_pre(idiag2))/dtime |
|---|
| 425 | d_ql = (ql_tot-ql_pre(idiag2))/dtime |
|---|
| 426 | d_qs = (qs_tot-qs_pre(idiag2))/dtime |
|---|
| 427 | d_ec = (ec_tot-ec_pre(idiag2))/dtime |
|---|
| 428 | d_qt = d_qw + d_ql + d_qs |
|---|
| 429 | !#ifdef ISO |
|---|
| 430 | ! do ixt=1,ntraciso |
|---|
| 431 | ! d_xtw(ixt) = (xtw_tot(ixt) - xtw_pre(ixt,idiag2) )/dtime |
|---|
| 432 | ! d_xtl(ixt) = (xtl_tot(ixt) - xtl_pre(ixt,idiag2) )/dtime |
|---|
| 433 | ! d_xts(ixt) = (xts_tot(ixt) - xts_pre(ixt,idiag2) )/dtime |
|---|
| 434 | ! enddo |
|---|
| 435 | !#endif |
|---|
| 436 | ELSE |
|---|
| 437 | d_h_vcol = 0. |
|---|
| 438 | d_h_dair = 0. |
|---|
| 439 | d_h_qw = 0. |
|---|
| 440 | d_h_ql = 0. |
|---|
| 441 | d_h_qs = 0. |
|---|
| 442 | d_qw = 0. |
|---|
| 443 | d_ql = 0. |
|---|
| 444 | d_qs = 0. |
|---|
| 445 | d_ec = 0. |
|---|
| 446 | d_qt = 0. |
|---|
| 447 | !#ifdef ISO |
|---|
| 448 | ! do ixt=1,ntraciso |
|---|
| 449 | ! d_xtw(ixt) = 0. |
|---|
| 450 | ! d_xtl(ixt) = 0. |
|---|
| 451 | ! d_xts(ixt) = 0. |
|---|
| 452 | ! enddo |
|---|
| 453 | !#endif |
|---|
| 454 | END IF |
|---|
| 455 | |
|---|
| 456 | IF (iprt>=2) THEN |
|---|
| 457 | WRITE (6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs |
|---|
| 458 | 9000 FORMAT ('Phys. Watter Mass Budget (kg/m2/s)', A15, 1I6, 10(1PE14.6)) |
|---|
| 459 | WRITE (6, 9001) tit, pas(idiag), d_h_vcol |
|---|
| 460 | 9001 FORMAT ('Phys. Enthalpy Budget (W/m2) ', A15, 1I6, 10(F8.2)) |
|---|
| 461 | WRITE (6, 9002) tit, pas(idiag), d_ec |
|---|
| 462 | 9002 FORMAT ('Phys. Cinetic Energy Budget (W/m2) ', A15, 1I6, 10(F8.2)) |
|---|
| 463 | END IF |
|---|
| 464 | |
|---|
| 465 | ! Store the new atmospheric state in "idiag" |
|---|
| 466 | |
|---|
| 467 | pas(idiag) = pas(idiag) + 1 |
|---|
| 468 | h_vcol_pre(idiag) = h_vcol_tot |
|---|
| 469 | h_dair_pre(idiag) = h_dair_tot |
|---|
| 470 | h_qw_pre(idiag) = h_qw_tot |
|---|
| 471 | h_ql_pre(idiag) = h_ql_tot |
|---|
| 472 | h_qs_pre(idiag) = h_qs_tot |
|---|
| 473 | qw_pre(idiag) = qw_tot |
|---|
| 474 | ql_pre(idiag) = ql_tot |
|---|
| 475 | qs_pre(idiag) = qs_tot |
|---|
| 476 | ec_pre(idiag) = ec_tot |
|---|
| 477 | !#ifdef ISO |
|---|
| 478 | ! do ixt=1,ntraciso |
|---|
| 479 | ! xtw_pre(ixt,idiag) = xtw_tot(ixt) |
|---|
| 480 | ! xtl_pre(ixt,idiag) = xtl_tot(ixt) |
|---|
| 481 | ! xts_pre(ixt,idiag) = xts_tot(ixt) |
|---|
| 482 | ! enddo |
|---|
| 483 | !#endif |
|---|
| 484 | |
|---|
| 485 | RETURN |
|---|
| 486 | END SUBROUTINE diagetpq |
|---|
