| 1 | !WRF:MODEL_LAYER:PHYSICS |
|---|
| 2 | ! |
|---|
| 3 | |
|---|
| 4 | MODULE module_mp_gsfcgce |
|---|
| 5 | |
|---|
| 6 | USE module_wrf_error |
|---|
| 7 | |
|---|
| 8 | !JJS 1/3/2008 vvvvv |
|---|
| 9 | |
|---|
| 10 | ! common /bt/ |
|---|
| 11 | REAL, PRIVATE :: rd1, rd2, al, cp |
|---|
| 12 | |
|---|
| 13 | ! common /cont/ |
|---|
| 14 | REAL, PRIVATE :: c38, c358, c610, c149, & |
|---|
| 15 | c879, c172, c409, c76, & |
|---|
| 16 | c218, c580, c141 |
|---|
| 17 | ! common /b3cs/ |
|---|
| 18 | REAL, PRIVATE :: ag, bg, as, bs, & |
|---|
| 19 | aw, bw, bgh, bgq, & |
|---|
| 20 | bsh, bsq, bwh, bwq |
|---|
| 21 | |
|---|
| 22 | ! common /size/ |
|---|
| 23 | REAL, PRIVATE :: tnw, tns, tng, & |
|---|
| 24 | roqs, roqg, roqr |
|---|
| 25 | |
|---|
| 26 | ! common /bterv/ |
|---|
| 27 | REAL, PRIVATE :: zrc, zgc, zsc, & |
|---|
| 28 | vrc, vgc, vsc |
|---|
| 29 | |
|---|
| 30 | ! common /bsnw/ |
|---|
| 31 | REAL, PRIVATE :: alv, alf, als, t0, t00, & |
|---|
| 32 | avc, afc, asc, rn1, bnd1, & |
|---|
| 33 | rn2, bnd2, rn3, rn4, rn5, & |
|---|
| 34 | rn6, rn7, rn8, rn9, rn10, & |
|---|
| 35 | rn101,rn10a, rn11,rn11a, rn12 |
|---|
| 36 | |
|---|
| 37 | REAL, PRIVATE :: rn14, rn15,rn15a, rn16, rn17, & |
|---|
| 38 | rn17a,rn17b,rn17c, rn18, rn18a, & |
|---|
| 39 | rn19,rn19a,rn19b, rn20, rn20a, & |
|---|
| 40 | rn20b, bnd3, rn21, rn22, rn23, & |
|---|
| 41 | rn23a,rn23b, rn25,rn30a, rn30b, & |
|---|
| 42 | rn30c, rn31, beta, rn32 |
|---|
| 43 | |
|---|
| 44 | REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a |
|---|
| 45 | |
|---|
| 46 | ! common /rsnw1/ |
|---|
| 47 | REAL, PRIVATE :: rn10b, rn10c, rnn191, rnn192, rn30, & |
|---|
| 48 | rnn30a, rn33, rn331, rn332 |
|---|
| 49 | |
|---|
| 50 | ! |
|---|
| 51 | REAL, PRIVATE, DIMENSION( 31 ) :: aa1, aa2 |
|---|
| 52 | DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, & |
|---|
| 53 | .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, & |
|---|
| 54 | .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, & |
|---|
| 55 | .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, & |
|---|
| 56 | .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, & |
|---|
| 57 | .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, & |
|---|
| 58 | .4834e-6/ |
|---|
| 59 | DATA aa2/.4006, .4831, .5320, .5307, .5319, & |
|---|
| 60 | .5249, .4888, .3894, .4047, .4318, & |
|---|
| 61 | .4771, .5183, .5463, .5651, .5813, & |
|---|
| 62 | .5655, .5478, .5203, .4906, .4447, & |
|---|
| 63 | .4126, .3960, .4149, .4320, .4506, & |
|---|
| 64 | .4483, .4460, .4433, .4413, .4382, & |
|---|
| 65 | .4361/ |
|---|
| 66 | |
|---|
| 67 | !JJS 1/3/2008 ^^^^^ |
|---|
| 68 | |
|---|
| 69 | CONTAINS |
|---|
| 70 | |
|---|
| 71 | !------------------------------------------------------------------- |
|---|
| 72 | ! NASA/GSFC GCE |
|---|
| 73 | ! Tao et al, 2001, Meteo. & Atmos. Phy., 97-137 |
|---|
| 74 | !------------------------------------------------------------------- |
|---|
| 75 | ! SUBROUTINE gsfcgce( th, th_old, & |
|---|
| 76 | SUBROUTINE gsfcgce( th, & |
|---|
| 77 | qv, ql, & |
|---|
| 78 | qr, qi, & |
|---|
| 79 | qs, & |
|---|
| 80 | ! qvold, qlold, & |
|---|
| 81 | ! qrold, qiold, & |
|---|
| 82 | ! qsold, & |
|---|
| 83 | rho, pii, p, dt_in, z, & |
|---|
| 84 | ht, dz8w, grav, & |
|---|
| 85 | rhowater, rhosnow, & |
|---|
| 86 | itimestep, & |
|---|
| 87 | ids,ide, jds,jde, kds,kde, & ! domain dims |
|---|
| 88 | ims,ime, jms,jme, kms,kme, & ! memory dims |
|---|
| 89 | its,ite, jts,jte, kts,kte, & ! tile dims |
|---|
| 90 | rainnc, rainncv, & |
|---|
| 91 | snownc, snowncv, sr, & |
|---|
| 92 | graupelnc, graupelncv, & |
|---|
| 93 | ! f_qg, qg, pgold & |
|---|
| 94 | f_qg, qg, & |
|---|
| 95 | ihail, ice2 & |
|---|
| 96 | ) |
|---|
| 97 | |
|---|
| 98 | !------------------------------------------------------------------- |
|---|
| 99 | IMPLICIT NONE |
|---|
| 100 | !------------------------------------------------------------------- |
|---|
| 101 | ! |
|---|
| 102 | ! JJS 2/15/2005 |
|---|
| 103 | ! |
|---|
| 104 | INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , & |
|---|
| 105 | ims,ime, jms,jme, kms,kme , & |
|---|
| 106 | its,ite, jts,jte, kts,kte |
|---|
| 107 | INTEGER, INTENT(IN ) :: itimestep, ihail, ice2 |
|---|
| 108 | |
|---|
| 109 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 110 | INTENT(INOUT) :: & |
|---|
| 111 | th, & |
|---|
| 112 | qv, & |
|---|
| 113 | ql, & |
|---|
| 114 | qr, & |
|---|
| 115 | qi, & |
|---|
| 116 | qs, & |
|---|
| 117 | qg |
|---|
| 118 | ! |
|---|
| 119 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 120 | INTENT(IN ) :: & |
|---|
| 121 | ! th_old, & |
|---|
| 122 | ! qvold, & |
|---|
| 123 | ! qlold, & |
|---|
| 124 | ! qrold, & |
|---|
| 125 | ! qiold, & |
|---|
| 126 | ! qsold, & |
|---|
| 127 | ! qgold, & |
|---|
| 128 | rho, & |
|---|
| 129 | pii, & |
|---|
| 130 | p, & |
|---|
| 131 | dz8w, & |
|---|
| 132 | z |
|---|
| 133 | |
|---|
| 134 | REAL, DIMENSION( ims:ime , jms:jme ), & |
|---|
| 135 | INTENT(INOUT) :: rainnc, & |
|---|
| 136 | rainncv, & |
|---|
| 137 | snownc, & |
|---|
| 138 | snowncv, & |
|---|
| 139 | sr, & |
|---|
| 140 | graupelnc, & |
|---|
| 141 | graupelncv |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht |
|---|
| 147 | |
|---|
| 148 | REAL, INTENT(IN ) :: dt_in, & |
|---|
| 149 | grav, & |
|---|
| 150 | rhowater, & |
|---|
| 151 | rhosnow |
|---|
| 152 | |
|---|
| 153 | LOGICAL, INTENT(IN), OPTIONAL :: F_QG |
|---|
| 154 | |
|---|
| 155 | ! LOCAL VAR |
|---|
| 156 | |
|---|
| 157 | |
|---|
| 158 | !jjs INTEGER :: min_q, max_q |
|---|
| 159 | |
|---|
| 160 | !jjs REAL, DIMENSION( its:ite , jts:jte ) & |
|---|
| 161 | !jjs :: rain, snow, graupel,ice |
|---|
| 162 | |
|---|
| 163 | ! |
|---|
| 164 | ! INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id |
|---|
| 165 | INTEGER :: itaobraun, istatmin, new_ice_sat, id |
|---|
| 166 | |
|---|
| 167 | INTEGER :: i, j, k |
|---|
| 168 | INTEGER :: iskip, ih, icount, ibud, i12h |
|---|
| 169 | REAL :: hour |
|---|
| 170 | REAL , PARAMETER :: cmin=1.e-20 |
|---|
| 171 | REAL :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2 |
|---|
| 172 | REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: & |
|---|
| 173 | th1, qv1, ql1, qr1, qi1, qs1, qg1 |
|---|
| 174 | |
|---|
| 175 | LOGICAL :: flag_qg |
|---|
| 176 | |
|---|
| 177 | ! |
|---|
| 178 | !c ihail = 0 for graupel, for tropical region |
|---|
| 179 | !c ihail = 1 for hail, for mid-lat region |
|---|
| 180 | |
|---|
| 181 | ! itaobraun: 0 for Tao's constantis, 1 for Braun's constants |
|---|
| 182 | !c if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6 |
|---|
| 183 | !c if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8 |
|---|
| 184 | itaobraun = 1 |
|---|
| 185 | |
|---|
| 186 | !c ice2 = 0 for 3 ice --- ice, snow and graupel/hail |
|---|
| 187 | !c ice2 = 1 for 2 ice --- ice and snow only |
|---|
| 188 | !c ice2 = 2 for 2 ice --- ice and graupel only, use ihail = 0 only |
|---|
| 189 | |
|---|
| 190 | i12h=int(86400./dt_in) |
|---|
| 191 | ! if (itimestep.eq.1.or.itimestep.eq.1441) then |
|---|
| 192 | ! if (mod(itimestep,1440).eq.1) then |
|---|
| 193 | if (mod(itimestep,i12h).eq.1) then |
|---|
| 194 | write(6,*) 'ihail=',ihail,' ice2=',ice2 |
|---|
| 195 | if (ice2.eq.0) then |
|---|
| 196 | write(6,*) 'Running 3-ice scheme in GSFCGCE with' |
|---|
| 197 | if (ihail.eq.0) then |
|---|
| 198 | write(6,*) ' ice, snow and graupel' |
|---|
| 199 | else if (ihail.eq.1) then |
|---|
| 200 | write(6,*) ' ice, snow and hail' |
|---|
| 201 | else |
|---|
| 202 | write(6,*) 'ihail has to be either 1 or 0' |
|---|
| 203 | stop |
|---|
| 204 | endif !ihail |
|---|
| 205 | else if (ice2.eq.1) then |
|---|
| 206 | write(6,*) 'Running 2-ice scheme in GSFCGCE with' |
|---|
| 207 | write(6,*) ' ice and snow' |
|---|
| 208 | else if (ice2.eq.2) then |
|---|
| 209 | write(6,*) 'Running 2-ice scheme in GSFCGCE with' |
|---|
| 210 | write(6,*) ' ice and graupel' |
|---|
| 211 | else |
|---|
| 212 | write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, or 2' |
|---|
| 213 | stop |
|---|
| 214 | endif !ice2 |
|---|
| 215 | endif !itimestep |
|---|
| 216 | |
|---|
| 217 | !c new_ice_sat = 0, 1 or 2 |
|---|
| 218 | new_ice_sat = 2 |
|---|
| 219 | |
|---|
| 220 | !c istatmin |
|---|
| 221 | istatmin = 180 |
|---|
| 222 | |
|---|
| 223 | !c id = 0 without in-line staticstics |
|---|
| 224 | !c id = 1 with in-line staticstics |
|---|
| 225 | id = 0 |
|---|
| 226 | |
|---|
| 227 | !c ibud = 0 no calculation of dth, dqv, dqrest and dqall |
|---|
| 228 | !c ibud = 1 yes |
|---|
| 229 | ibud = 0 |
|---|
| 230 | |
|---|
| 231 | !jjs dt=dt_in |
|---|
| 232 | !jjs rhoe_s=1.29 |
|---|
| 233 | ! |
|---|
| 234 | ! IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN |
|---|
| 235 | ! CALL wrf_error_fatal3 ( "module_mp_lin.b" , 130 , 'module_mp_lin: Improper use of Lin et al scheme; no ice phase. Please chose another one.') |
|---|
| 236 | ! ENDIF |
|---|
| 237 | |
|---|
| 238 | ! calculte fallflux and precipiation in MKS system |
|---|
| 239 | |
|---|
| 240 | call fall_flux(dt_in, qr, qi, qs, qg, p, & |
|---|
| 241 | rho, z, dz8w, ht, rainnc, & |
|---|
| 242 | rainncv, grav,itimestep, & |
|---|
| 243 | rhowater, rhosnow, & |
|---|
| 244 | snownc, snowncv, sr, & |
|---|
| 245 | graupelnc, graupelncv, & |
|---|
| 246 | ihail, ice2, & |
|---|
| 247 | ims,ime, jms,jme, kms,kme, & ! memory dims |
|---|
| 248 | its,ite, jts,jte, kts,kte ) ! tile dims |
|---|
| 249 | !----------------------------------------------------------------------- |
|---|
| 250 | |
|---|
| 251 | !c set up constants used internally in GCE |
|---|
| 252 | |
|---|
| 253 | call consat_s (ihail, itaobraun) |
|---|
| 254 | |
|---|
| 255 | |
|---|
| 256 | !c Negative values correction |
|---|
| 257 | |
|---|
| 258 | iskip = 1 |
|---|
| 259 | ! i12h=int(86400./dt_in) |
|---|
| 260 | |
|---|
| 261 | if (iskip.eq.0) then |
|---|
| 262 | call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 263 | itimestep,1, & |
|---|
| 264 | its,ite,jts,jte,kts,kte) |
|---|
| 265 | call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 266 | itimestep,2, & |
|---|
| 267 | its,ite,jts,jte,kts,kte) |
|---|
| 268 | call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 269 | itimestep,3, & |
|---|
| 270 | its,ite,jts,jte,kts,kte) |
|---|
| 271 | call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 272 | itimestep,4, & |
|---|
| 273 | its,ite,jts,jte,kts,kte) |
|---|
| 274 | call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 275 | itimestep,5, & |
|---|
| 276 | its,ite,jts,jte,kts,kte) |
|---|
| 277 | call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, & |
|---|
| 278 | itimestep,6, & |
|---|
| 279 | its,ite,jts,jte,kts,kte) |
|---|
| 280 | if (mod(itimestep,i12h).eq.1) then |
|---|
| 281 | print *,'negative correction used in gsfcgce mp ' |
|---|
| 282 | endif |
|---|
| 283 | endif ! iskip |
|---|
| 284 | |
|---|
| 285 | !c microphysics in GCE |
|---|
| 286 | |
|---|
| 287 | call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin, & |
|---|
| 288 | new_ice_sat, id, & |
|---|
| 289 | ! th, th_old, qv, ql, qr, & |
|---|
| 290 | th, qv, ql, qr, & |
|---|
| 291 | qi, qs, qg, & |
|---|
| 292 | ! qvold, qlold, qrold, & |
|---|
| 293 | ! qiold, qsold, qgold, & |
|---|
| 294 | rho, pii, p, itimestep, & |
|---|
| 295 | ids,ide, jds,jde, kds,kde, & ! domain dims |
|---|
| 296 | ims,ime, jms,jme, kms,kme, & ! memory dims |
|---|
| 297 | its,ite, jts,jte, kts,kte & ! tile dims |
|---|
| 298 | ) |
|---|
| 299 | |
|---|
| 300 | END SUBROUTINE gsfcgce |
|---|
| 301 | |
|---|
| 302 | !----------------------------------------------------------------------- |
|---|
| 303 | SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p, & |
|---|
| 304 | rho, z, dz8w, topo, rainnc, & |
|---|
| 305 | rainncv, grav, itimestep, & |
|---|
| 306 | rhowater, rhosnow, & |
|---|
| 307 | snownc, snowncv, sr, & |
|---|
| 308 | graupelnc, graupelncv, & |
|---|
| 309 | ihail, ice2, & |
|---|
| 310 | ims,ime, jms,jme, kms,kme, & ! memory dims |
|---|
| 311 | its,ite, jts,jte, kts,kte ) ! tile dims |
|---|
| 312 | !----------------------------------------------------------------------- |
|---|
| 313 | ! adopted from Jiun-Dar Chern's codes for Purdue Regional Model |
|---|
| 314 | ! adopted by Jainn J. Shi, 6/10/2005 |
|---|
| 315 | !----------------------------------------------------------------------- |
|---|
| 316 | |
|---|
| 317 | IMPLICIT NONE |
|---|
| 318 | INTEGER, INTENT(IN ) :: ihail, ice2, & |
|---|
| 319 | ims,ime, jms,jme, kms,kme, & |
|---|
| 320 | its,ite, jts,jte, kts,kte |
|---|
| 321 | INTEGER, INTENT(IN ) :: itimestep |
|---|
| 322 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 323 | INTENT(INOUT) :: qr, qi, qs, qg |
|---|
| 324 | REAL, DIMENSION( ims:ime , jms:jme ), & |
|---|
| 325 | INTENT(INOUT) :: rainnc, rainncv, & |
|---|
| 326 | snownc, snowncv, sr, & |
|---|
| 327 | graupelnc, graupelncv |
|---|
| 328 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 329 | INTENT(IN ) :: rho, z, dz8w, p |
|---|
| 330 | |
|---|
| 331 | REAL, INTENT(IN ) :: dt, grav, rhowater, rhosnow |
|---|
| 332 | |
|---|
| 333 | |
|---|
| 334 | REAL, DIMENSION( ims:ime , jms:jme ), & |
|---|
| 335 | INTENT(IN ) :: topo |
|---|
| 336 | |
|---|
| 337 | ! temperary vars |
|---|
| 338 | |
|---|
| 339 | REAL, DIMENSION( kts:kte ) :: sqrhoz |
|---|
| 340 | REAL :: tmp1, term0 |
|---|
| 341 | REAL :: pptrain, pptsnow, & |
|---|
| 342 | pptgraul, pptice |
|---|
| 343 | REAL, DIMENSION( kts:kte ) :: qrz, qiz, qsz, qgz, & |
|---|
| 344 | zz, dzw, prez, rhoz, & |
|---|
| 345 | orhoz |
|---|
| 346 | |
|---|
| 347 | |
|---|
| 348 | INTEGER :: k, i, j |
|---|
| 349 | ! |
|---|
| 350 | |
|---|
| 351 | REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vti |
|---|
| 352 | |
|---|
| 353 | REAL :: dtb, pi, consta, constc, gambp4, & |
|---|
| 354 | gamdp4, gam4pt5, gam4bbar |
|---|
| 355 | |
|---|
| 356 | ! Lin |
|---|
| 357 | REAL , PARAMETER :: xnor = 8.0e6 |
|---|
| 358 | ! REAL , PARAMETER :: xnos = 3.0e6 |
|---|
| 359 | REAL , PARAMETER :: xnos = 1.6e7 ! Tao's value |
|---|
| 360 | REAL , PARAMETER :: & |
|---|
| 361 | constb = 0.8, constd = 0.25, o6 = 1./6., & |
|---|
| 362 | cdrag = 0.6 |
|---|
| 363 | ! Lin |
|---|
| 364 | ! REAL , PARAMETER :: xnoh = 4.0e4 |
|---|
| 365 | REAL , PARAMETER :: xnoh = 2.0e5 ! Tao's value |
|---|
| 366 | REAL , PARAMETER :: rhohail = 917. |
|---|
| 367 | |
|---|
| 368 | ! Hobbs |
|---|
| 369 | REAL , PARAMETER :: xnog = 4.0e6 |
|---|
| 370 | REAL , PARAMETER :: rhograul = 400. |
|---|
| 371 | REAL , PARAMETER :: abar = 19.3, bbar = 0.37, & |
|---|
| 372 | p0 = 1.0e5 |
|---|
| 373 | |
|---|
| 374 | REAL , PARAMETER :: rhoe_s = 1.29 |
|---|
| 375 | |
|---|
| 376 | ! for terminal velocity flux |
|---|
| 377 | INTEGER :: min_q, max_q |
|---|
| 378 | REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz |
|---|
| 379 | LOGICAL :: notlast |
|---|
| 380 | |
|---|
| 381 | ! if (itimestep.eq.1) then |
|---|
| 382 | ! write(6, *) 'in fall_flux' |
|---|
| 383 | ! write(6, *) 'ims=', ims, ' ime=', ime |
|---|
| 384 | ! write(6, *) 'jms=', jms, ' jme=', jme |
|---|
| 385 | ! write(6, *) 'kms=', kms, ' kme=', kme |
|---|
| 386 | ! write(6, *) 'its=', its, ' ite=', ite |
|---|
| 387 | ! write(6, *) 'jts=', jts, ' jte=', jte |
|---|
| 388 | ! write(6, *) 'kts=', kts, ' kte=', kte |
|---|
| 389 | ! write(6, *) 'dt=', dt |
|---|
| 390 | ! write(6, *) 'ihail=', ihail |
|---|
| 391 | ! write(6, *) 'ICE2=', ICE2 |
|---|
| 392 | ! write(6, *) 'dt=', dt |
|---|
| 393 | ! endif |
|---|
| 394 | |
|---|
| 395 | !----------------------------------------------------------------------- |
|---|
| 396 | ! This program calculates precipitation fluxes due to terminal velocities. |
|---|
| 397 | !----------------------------------------------------------------------- |
|---|
| 398 | |
|---|
| 399 | dtb=dt |
|---|
| 400 | pi=acos(-1.) |
|---|
| 401 | consta=2115.0*0.01**(1-constb) |
|---|
| 402 | ! constc=152.93*0.01**(1-constd) |
|---|
| 403 | constc=78.63*0.01**(1-constd) |
|---|
| 404 | |
|---|
| 405 | ! Gamma function |
|---|
| 406 | gambp4=ggamma(constb+4.) |
|---|
| 407 | gamdp4=ggamma(constd+4.) |
|---|
| 408 | gam4pt5=ggamma(4.5) |
|---|
| 409 | gam4bbar=ggamma(4.+bbar) |
|---|
| 410 | |
|---|
| 411 | !*********************************************************************** |
|---|
| 412 | ! Calculate precipitation fluxes due to terminal velocities. |
|---|
| 413 | !*********************************************************************** |
|---|
| 414 | ! |
|---|
| 415 | !- Calculate termianl velocity (vt?) of precipitation q?z |
|---|
| 416 | !- Find maximum vt? to determine the small delta t |
|---|
| 417 | |
|---|
| 418 | j_loop: do j = jts, jte |
|---|
| 419 | i_loop: do i = its, ite |
|---|
| 420 | |
|---|
| 421 | pptrain = 0. |
|---|
| 422 | pptsnow = 0. |
|---|
| 423 | pptgraul = 0. |
|---|
| 424 | pptice = 0. |
|---|
| 425 | |
|---|
| 426 | do k = kts, kte |
|---|
| 427 | qrz(k)=qr(i,k,j) |
|---|
| 428 | rhoz(k)=rho(i,k,j) |
|---|
| 429 | orhoz(k)=1./rhoz(k) |
|---|
| 430 | prez(k)=p(i,k,j) |
|---|
| 431 | sqrhoz(k)=sqrt(rhoe_s/rhoz(k)) |
|---|
| 432 | zz(k)=z(i,k,j) |
|---|
| 433 | dzw(k)=dz8w(i,k,j) |
|---|
| 434 | enddo !k |
|---|
| 435 | |
|---|
| 436 | DO k = kts, kte |
|---|
| 437 | qiz(k)=qi(i,k,j) |
|---|
| 438 | ENDDO |
|---|
| 439 | |
|---|
| 440 | DO k = kts, kte |
|---|
| 441 | qsz(k)=qs(i,k,j) |
|---|
| 442 | ENDDO |
|---|
| 443 | |
|---|
| 444 | IF (ice2 .eq. 0) THEN |
|---|
| 445 | DO k = kts, kte |
|---|
| 446 | qgz(k)=qg(i,k,j) |
|---|
| 447 | ENDDO |
|---|
| 448 | ELSE |
|---|
| 449 | DO k = kts, kte |
|---|
| 450 | qgz(k)=0. |
|---|
| 451 | ENDDO |
|---|
| 452 | ENDIF |
|---|
| 453 | |
|---|
| 454 | |
|---|
| 455 | ! |
|---|
| 456 | !-- rain |
|---|
| 457 | ! |
|---|
| 458 | t_del_tv=0. |
|---|
| 459 | del_tv=dtb |
|---|
| 460 | notlast=.true. |
|---|
| 461 | DO while (notlast) |
|---|
| 462 | ! |
|---|
| 463 | min_q=kte |
|---|
| 464 | max_q=kts-1 |
|---|
| 465 | ! |
|---|
| 466 | do k=kts,kte-1 |
|---|
| 467 | if (qrz(k) .gt. 1.0e-8) then |
|---|
| 468 | min_q=min0(min_q,k) |
|---|
| 469 | max_q=max0(max_q,k) |
|---|
| 470 | tmp1=sqrt(pi*rhowater*xnor/rhoz(k)/qrz(k)) |
|---|
| 471 | tmp1=sqrt(tmp1) |
|---|
| 472 | vtr(k)=consta*gambp4*sqrhoz(k)/tmp1**constb |
|---|
| 473 | vtr(k)=vtr(k)/6. |
|---|
| 474 | if (k .eq. 1) then |
|---|
| 475 | del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k)) |
|---|
| 476 | else |
|---|
| 477 | del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k)) |
|---|
| 478 | endif |
|---|
| 479 | else |
|---|
| 480 | vtr(k)=0. |
|---|
| 481 | endif |
|---|
| 482 | enddo |
|---|
| 483 | |
|---|
| 484 | if (max_q .ge. min_q) then |
|---|
| 485 | ! |
|---|
| 486 | !- Check if the summation of the small delta t >= big delta t |
|---|
| 487 | ! (t_del_tv) (del_tv) (dtb) |
|---|
| 488 | |
|---|
| 489 | t_del_tv=t_del_tv+del_tv |
|---|
| 490 | ! |
|---|
| 491 | if ( t_del_tv .ge. dtb ) then |
|---|
| 492 | notlast=.false. |
|---|
| 493 | del_tv=dtb+del_tv-t_del_tv |
|---|
| 494 | endif |
|---|
| 495 | |
|---|
| 496 | ! use small delta t to calculate the qrz flux |
|---|
| 497 | ! termi is the qrz flux pass in the grid box through the upper boundary |
|---|
| 498 | ! termo is the qrz flux pass out the grid box through the lower boundary |
|---|
| 499 | ! |
|---|
| 500 | fluxin=0. |
|---|
| 501 | do k=max_q,min_q,-1 |
|---|
| 502 | fluxout=rhoz(k)*vtr(k)*qrz(k) |
|---|
| 503 | flux=(fluxin-fluxout)/rhoz(k)/dzw(k) |
|---|
| 504 | ! tmpqrz=qrz(k) |
|---|
| 505 | qrz(k)=qrz(k)+del_tv*flux |
|---|
| 506 | qrz(k)=amax1(0.,qrz(k)) |
|---|
| 507 | qr(i,k,j)=qrz(k) |
|---|
| 508 | fluxin=fluxout |
|---|
| 509 | enddo |
|---|
| 510 | if (min_q .eq. 1) then |
|---|
| 511 | pptrain=pptrain+fluxin*del_tv |
|---|
| 512 | else |
|---|
| 513 | qrz(min_q-1)=qrz(min_q-1)+del_tv* & |
|---|
| 514 | fluxin/rhoz(min_q-1)/dzw(min_q-1) |
|---|
| 515 | qr(i,min_q-1,j)=qrz(min_q-1) |
|---|
| 516 | endif |
|---|
| 517 | ! |
|---|
| 518 | else |
|---|
| 519 | notlast=.false. |
|---|
| 520 | endif |
|---|
| 521 | ENDDO |
|---|
| 522 | |
|---|
| 523 | ! |
|---|
| 524 | !-- snow |
|---|
| 525 | ! |
|---|
| 526 | t_del_tv=0. |
|---|
| 527 | del_tv=dtb |
|---|
| 528 | notlast=.true. |
|---|
| 529 | |
|---|
| 530 | DO while (notlast) |
|---|
| 531 | ! |
|---|
| 532 | min_q=kte |
|---|
| 533 | max_q=kts-1 |
|---|
| 534 | ! |
|---|
| 535 | do k=kts,kte-1 |
|---|
| 536 | if (qsz(k) .gt. 1.0e-8) then |
|---|
| 537 | min_q=min0(min_q,k) |
|---|
| 538 | max_q=max0(max_q,k) |
|---|
| 539 | tmp1=sqrt(pi*rhosnow*xnos/rhoz(k)/qsz(k)) |
|---|
| 540 | tmp1=sqrt(tmp1) |
|---|
| 541 | vts(k)=constc*gamdp4*sqrhoz(k)/tmp1**constd |
|---|
| 542 | vts(k)=vts(k)/6. |
|---|
| 543 | if (k .eq. 1) then |
|---|
| 544 | del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k)) |
|---|
| 545 | else |
|---|
| 546 | del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k)) |
|---|
| 547 | endif |
|---|
| 548 | else |
|---|
| 549 | vts(k)=0. |
|---|
| 550 | endif |
|---|
| 551 | enddo |
|---|
| 552 | |
|---|
| 553 | if (max_q .ge. min_q) then |
|---|
| 554 | ! |
|---|
| 555 | ! |
|---|
| 556 | !- Check if the summation of the small delta t >= big delta t |
|---|
| 557 | ! (t_del_tv) (del_tv) (dtb) |
|---|
| 558 | |
|---|
| 559 | t_del_tv=t_del_tv+del_tv |
|---|
| 560 | |
|---|
| 561 | if ( t_del_tv .ge. dtb ) then |
|---|
| 562 | notlast=.false. |
|---|
| 563 | del_tv=dtb+del_tv-t_del_tv |
|---|
| 564 | endif |
|---|
| 565 | |
|---|
| 566 | ! use small delta t to calculate the qsz flux |
|---|
| 567 | ! termi is the qsz flux pass in the grid box through the upper boundary |
|---|
| 568 | ! termo is the qsz flux pass out the grid box through the lower boundary |
|---|
| 569 | ! |
|---|
| 570 | fluxin=0. |
|---|
| 571 | do k=max_q,min_q,-1 |
|---|
| 572 | fluxout=rhoz(k)*vts(k)*qsz(k) |
|---|
| 573 | flux=(fluxin-fluxout)/rhoz(k)/dzw(k) |
|---|
| 574 | qsz(k)=qsz(k)+del_tv*flux |
|---|
| 575 | qsz(k)=amax1(0.,qsz(k)) |
|---|
| 576 | qs(i,k,j)=qsz(k) |
|---|
| 577 | fluxin=fluxout |
|---|
| 578 | enddo |
|---|
| 579 | if (min_q .eq. 1) then |
|---|
| 580 | pptsnow=pptsnow+fluxin*del_tv |
|---|
| 581 | else |
|---|
| 582 | qsz(min_q-1)=qsz(min_q-1)+del_tv* & |
|---|
| 583 | fluxin/rhoz(min_q-1)/dzw(min_q-1) |
|---|
| 584 | qs(i,min_q-1,j)=qsz(min_q-1) |
|---|
| 585 | endif |
|---|
| 586 | ! |
|---|
| 587 | else |
|---|
| 588 | notlast=.false. |
|---|
| 589 | endif |
|---|
| 590 | |
|---|
| 591 | ENDDO |
|---|
| 592 | |
|---|
| 593 | ! |
|---|
| 594 | ! ice2=0 --- with hail/graupel |
|---|
| 595 | ! ice2=1 --- without hail/graupel |
|---|
| 596 | ! |
|---|
| 597 | if (ice2.eq.0) then |
|---|
| 598 | ! |
|---|
| 599 | !-- If IHAIL=1, use hail. |
|---|
| 600 | !-- If IHAIL=0, use graupel. |
|---|
| 601 | ! |
|---|
| 602 | ! if (ihail .eq. 1) then |
|---|
| 603 | ! xnog = xnoh |
|---|
| 604 | ! rhograul = rhohail |
|---|
| 605 | ! endif |
|---|
| 606 | |
|---|
| 607 | t_del_tv=0. |
|---|
| 608 | del_tv=dtb |
|---|
| 609 | notlast=.true. |
|---|
| 610 | ! |
|---|
| 611 | DO while (notlast) |
|---|
| 612 | ! |
|---|
| 613 | min_q=kte |
|---|
| 614 | max_q=kts-1 |
|---|
| 615 | ! |
|---|
| 616 | do k=kts,kte-1 |
|---|
| 617 | if (qgz(k) .gt. 1.0e-8) then |
|---|
| 618 | if (ihail .eq. 1) then |
|---|
| 619 | ! for hail, based on Lin et al (1983) |
|---|
| 620 | min_q=min0(min_q,k) |
|---|
| 621 | max_q=max0(max_q,k) |
|---|
| 622 | tmp1=sqrt(pi*rhohail*xnoh/rhoz(k)/qgz(k)) |
|---|
| 623 | tmp1=sqrt(tmp1) |
|---|
| 624 | term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag) |
|---|
| 625 | vtg(k)=gam4pt5*term0*sqrt(1./tmp1) |
|---|
| 626 | vtg(k)=vtg(k)/6. |
|---|
| 627 | if (k .eq. 1) then |
|---|
| 628 | del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k)) |
|---|
| 629 | else |
|---|
| 630 | del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k)) |
|---|
| 631 | endif !k |
|---|
| 632 | else |
|---|
| 633 | ! added by JJS |
|---|
| 634 | ! for graupel, based on RH (1984) |
|---|
| 635 | min_q=min0(min_q,k) |
|---|
| 636 | max_q=max0(max_q,k) |
|---|
| 637 | tmp1=sqrt(pi*rhograul*xnog/rhoz(k)/qgz(k)) |
|---|
| 638 | tmp1=sqrt(tmp1) |
|---|
| 639 | tmp1=tmp1**bbar |
|---|
| 640 | tmp1=1./tmp1 |
|---|
| 641 | term0=abar*gam4bbar/6. |
|---|
| 642 | vtg(k)=term0*tmp1*(p0/prez(k))**0.4 |
|---|
| 643 | if (k .eq. 1) then |
|---|
| 644 | del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k)) |
|---|
| 645 | else |
|---|
| 646 | del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k)) |
|---|
| 647 | endif !k |
|---|
| 648 | endif !ihail |
|---|
| 649 | else |
|---|
| 650 | vtg(k)=0. |
|---|
| 651 | endif !qgz |
|---|
| 652 | enddo !k |
|---|
| 653 | |
|---|
| 654 | if (max_q .ge. min_q) then |
|---|
| 655 | ! |
|---|
| 656 | ! |
|---|
| 657 | !- Check if the summation of the small delta t >= big delta t |
|---|
| 658 | ! (t_del_tv) (del_tv) (dtb) |
|---|
| 659 | |
|---|
| 660 | t_del_tv=t_del_tv+del_tv |
|---|
| 661 | |
|---|
| 662 | if ( t_del_tv .ge. dtb ) then |
|---|
| 663 | notlast=.false. |
|---|
| 664 | del_tv=dtb+del_tv-t_del_tv |
|---|
| 665 | endif |
|---|
| 666 | |
|---|
| 667 | ! use small delta t to calculate the qgz flux |
|---|
| 668 | ! termi is the qgz flux pass in the grid box through the upper boundary |
|---|
| 669 | ! termo is the qgz flux pass out the grid box through the lower boundary |
|---|
| 670 | ! |
|---|
| 671 | fluxin=0. |
|---|
| 672 | do k=max_q,min_q,-1 |
|---|
| 673 | fluxout=rhoz(k)*vtg(k)*qgz(k) |
|---|
| 674 | flux=(fluxin-fluxout)/rhoz(k)/dzw(k) |
|---|
| 675 | qgz(k)=qgz(k)+del_tv*flux |
|---|
| 676 | qgz(k)=amax1(0.,qgz(k)) |
|---|
| 677 | qg(i,k,j)=qgz(k) |
|---|
| 678 | fluxin=fluxout |
|---|
| 679 | enddo |
|---|
| 680 | if (min_q .eq. 1) then |
|---|
| 681 | pptgraul=pptgraul+fluxin*del_tv |
|---|
| 682 | else |
|---|
| 683 | qgz(min_q-1)=qgz(min_q-1)+del_tv* & |
|---|
| 684 | fluxin/rhoz(min_q-1)/dzw(min_q-1) |
|---|
| 685 | qg(i,min_q-1,j)=qgz(min_q-1) |
|---|
| 686 | endif |
|---|
| 687 | ! |
|---|
| 688 | else |
|---|
| 689 | notlast=.false. |
|---|
| 690 | endif |
|---|
| 691 | ! |
|---|
| 692 | ENDDO |
|---|
| 693 | ENDIF !ice2 |
|---|
| 694 | ! |
|---|
| 695 | !-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL |
|---|
| 696 | ! |
|---|
| 697 | |
|---|
| 698 | t_del_tv=0. |
|---|
| 699 | del_tv=dtb |
|---|
| 700 | notlast=.true. |
|---|
| 701 | ! |
|---|
| 702 | DO while (notlast) |
|---|
| 703 | ! |
|---|
| 704 | min_q=kte |
|---|
| 705 | max_q=kts-1 |
|---|
| 706 | ! |
|---|
| 707 | do k=kts,kte-1 |
|---|
| 708 | if (qiz(k) .gt. 1.0e-8) then |
|---|
| 709 | min_q=min0(min_q,k) |
|---|
| 710 | max_q=max0(max_q,k) |
|---|
| 711 | vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16 ! Heymsfield and Donner |
|---|
| 712 | if (k .eq. 1) then |
|---|
| 713 | del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k)) |
|---|
| 714 | else |
|---|
| 715 | del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k)) |
|---|
| 716 | endif |
|---|
| 717 | else |
|---|
| 718 | vti(k)=0. |
|---|
| 719 | endif |
|---|
| 720 | enddo |
|---|
| 721 | |
|---|
| 722 | if (max_q .ge. min_q) then |
|---|
| 723 | ! |
|---|
| 724 | ! |
|---|
| 725 | !- Check if the summation of the small delta t >= big delta t |
|---|
| 726 | ! (t_del_tv) (del_tv) (dtb) |
|---|
| 727 | |
|---|
| 728 | t_del_tv=t_del_tv+del_tv |
|---|
| 729 | |
|---|
| 730 | if ( t_del_tv .ge. dtb ) then |
|---|
| 731 | notlast=.false. |
|---|
| 732 | del_tv=dtb+del_tv-t_del_tv |
|---|
| 733 | endif |
|---|
| 734 | |
|---|
| 735 | ! use small delta t to calculate the qiz flux |
|---|
| 736 | ! termi is the qiz flux pass in the grid box through the upper boundary |
|---|
| 737 | ! termo is the qiz flux pass out the grid box through the lower boundary |
|---|
| 738 | ! |
|---|
| 739 | |
|---|
| 740 | fluxin=0. |
|---|
| 741 | do k=max_q,min_q,-1 |
|---|
| 742 | fluxout=rhoz(k)*vti(k)*qiz(k) |
|---|
| 743 | flux=(fluxin-fluxout)/rhoz(k)/dzw(k) |
|---|
| 744 | qiz(k)=qiz(k)+del_tv*flux |
|---|
| 745 | qiz(k)=amax1(0.,qiz(k)) |
|---|
| 746 | qi(i,k,j)=qiz(k) |
|---|
| 747 | fluxin=fluxout |
|---|
| 748 | enddo |
|---|
| 749 | if (min_q .eq. 1) then |
|---|
| 750 | pptice=pptice+fluxin*del_tv |
|---|
| 751 | else |
|---|
| 752 | qiz(min_q-1)=qiz(min_q-1)+del_tv* & |
|---|
| 753 | fluxin/rhoz(min_q-1)/dzw(min_q-1) |
|---|
| 754 | qi(i,min_q-1,j)=qiz(min_q-1) |
|---|
| 755 | endif |
|---|
| 756 | ! |
|---|
| 757 | else |
|---|
| 758 | notlast=.false. |
|---|
| 759 | endif |
|---|
| 760 | ! |
|---|
| 761 | ENDDO !notlast |
|---|
| 762 | |
|---|
| 763 | ! prnc(i,j)=prnc(i,j)+pptrain |
|---|
| 764 | ! psnowc(i,j)=psnowc(i,j)+pptsnow |
|---|
| 765 | ! pgrauc(i,j)=pgrauc(i,j)+pptgraul |
|---|
| 766 | ! picec(i,j)=picec(i,j)+pptice |
|---|
| 767 | ! |
|---|
| 768 | |
|---|
| 769 | ! write(6,*) 'i=',i,' j=',j,' ', pptrain, pptsnow, pptgraul, pptice |
|---|
| 770 | ! call flush(6) |
|---|
| 771 | |
|---|
| 772 | snowncv(i,j) = pptsnow |
|---|
| 773 | snownc(i,j) = snownc(i,j) + pptsnow |
|---|
| 774 | graupelncv(i,j) = pptgraul |
|---|
| 775 | graupelnc(i,j) = graupelnc(i,j) + pptgraul |
|---|
| 776 | RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice |
|---|
| 777 | RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice |
|---|
| 778 | sr(i,j) = 0. |
|---|
| 779 | if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j) |
|---|
| 780 | |
|---|
| 781 | ENDDO i_loop |
|---|
| 782 | ENDDO j_loop |
|---|
| 783 | |
|---|
| 784 | ! if (itimestep.eq.6480) then |
|---|
| 785 | ! write(51,*) 'in the end of fallflux, itimestep=',itimestep |
|---|
| 786 | ! do j=jts,jte |
|---|
| 787 | ! do i=its,ite |
|---|
| 788 | ! if (rainnc(i,j).gt.400.) then |
|---|
| 789 | ! write(50,*) 'i=',i,' j=',j,' rainnc=',rainnc |
|---|
| 790 | ! endif |
|---|
| 791 | ! enddo |
|---|
| 792 | ! enddo |
|---|
| 793 | ! endif |
|---|
| 794 | |
|---|
| 795 | END SUBROUTINE fall_flux |
|---|
| 796 | |
|---|
| 797 | !---------------------------------------------------------------- |
|---|
| 798 | REAL FUNCTION ggamma(X) |
|---|
| 799 | !---------------------------------------------------------------- |
|---|
| 800 | IMPLICIT NONE |
|---|
| 801 | !---------------------------------------------------------------- |
|---|
| 802 | REAL, INTENT(IN ) :: x |
|---|
| 803 | REAL, DIMENSION(8) :: B |
|---|
| 804 | INTEGER ::j, K1 |
|---|
| 805 | REAL ::PF, G1TO2 ,TEMP |
|---|
| 806 | |
|---|
| 807 | DATA B/-.577191652,.988205891,-.897056937,.918206857, & |
|---|
| 808 | -.756704078,.482199394,-.193527818,.035868343/ |
|---|
| 809 | |
|---|
| 810 | PF=1. |
|---|
| 811 | TEMP=X |
|---|
| 812 | DO 10 J=1,200 |
|---|
| 813 | IF (TEMP .LE. 2) GO TO 20 |
|---|
| 814 | TEMP=TEMP-1. |
|---|
| 815 | 10 PF=PF*TEMP |
|---|
| 816 | 100 FORMAT(//,5X,'module_gsfcgce: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5) |
|---|
| 817 | WRITE(wrf_err_message,100)X |
|---|
| 818 | CALL wrf_error_fatal(wrf_err_message) |
|---|
| 819 | 20 G1TO2=1. |
|---|
| 820 | TEMP=TEMP - 1. |
|---|
| 821 | DO 30 K1=1,8 |
|---|
| 822 | 30 G1TO2=G1TO2 + B(K1)*TEMP**K1 |
|---|
| 823 | ggamma=PF*G1TO2 |
|---|
| 824 | |
|---|
| 825 | END FUNCTION ggamma |
|---|
| 826 | |
|---|
| 827 | !----------------------------------------------------------------------- |
|---|
| 828 | !c Correction of negative values |
|---|
| 829 | SUBROUTINE negcor ( X, rho, dz8w, & |
|---|
| 830 | ims,ime, jms,jme, kms,kme, & ! memory dims |
|---|
| 831 | itimestep, ics, & |
|---|
| 832 | its,ite, jts,jte, kts,kte ) ! tile dims |
|---|
| 833 | !----------------------------------------------------------------------- |
|---|
| 834 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 835 | INTENT(INOUT) :: X |
|---|
| 836 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 837 | INTENT(IN ) :: rho, dz8w |
|---|
| 838 | integer, INTENT(IN ) :: itimestep, ics |
|---|
| 839 | |
|---|
| 840 | !c Local variables |
|---|
| 841 | ! REAL, DIMENSION( kts:kte ) :: Y1, Y2 |
|---|
| 842 | REAL :: A0, A1, A2 |
|---|
| 843 | |
|---|
| 844 | A1=0. |
|---|
| 845 | A2=0. |
|---|
| 846 | do k=kts,kte |
|---|
| 847 | do j=jts,jte |
|---|
| 848 | do i=its,ite |
|---|
| 849 | A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j) |
|---|
| 850 | A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j) |
|---|
| 851 | enddo |
|---|
| 852 | enddo |
|---|
| 853 | enddo |
|---|
| 854 | |
|---|
| 855 | ! A1=0.0 |
|---|
| 856 | ! A2=0.0 |
|---|
| 857 | ! do k=kts,kte |
|---|
| 858 | ! A1=A1+Y1(k) |
|---|
| 859 | ! A2=A2+Y2(k) |
|---|
| 860 | ! enddo |
|---|
| 861 | |
|---|
| 862 | A0=0.0 |
|---|
| 863 | |
|---|
| 864 | if (A1.NE.0.0.and.A1.GT.A2) then |
|---|
| 865 | A0=(A1-A2)/A1 |
|---|
| 866 | |
|---|
| 867 | if (mod(itimestep,540).eq.0) then |
|---|
| 868 | if (ics.eq.1) then |
|---|
| 869 | write(61,*) 'kms=',kms,' kme=',kme,' kts=',kts,' kte=',kte |
|---|
| 870 | write(61,*) 'jms=',jms,' jme=',jme,' jts=',jts,' jte=',jte |
|---|
| 871 | write(61,*) 'ims=',ims,' ime=',ime,' its=',its,' ite=',ite |
|---|
| 872 | endif |
|---|
| 873 | if (ics.eq.1) then |
|---|
| 874 | write(61,*) 'qv timestep=',itimestep |
|---|
| 875 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 876 | else if (ics.eq.2) then |
|---|
| 877 | write(61,*) 'ql timestep=',itimestep |
|---|
| 878 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 879 | else if (ics.eq.3) then |
|---|
| 880 | write(61,*) 'qr timestep=',itimestep |
|---|
| 881 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 882 | else if (ics.eq.4) then |
|---|
| 883 | write(61,*) 'qi timestep=',itimestep |
|---|
| 884 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 885 | else if (ics.eq.5) then |
|---|
| 886 | write(61,*) 'qs timestep=',itimestep |
|---|
| 887 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 888 | else if (ics.eq.6) then |
|---|
| 889 | write(61,*) 'qg timestep=',itimestep |
|---|
| 890 | write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0 |
|---|
| 891 | else |
|---|
| 892 | write(61,*) 'wrong cloud specieis number' |
|---|
| 893 | endif |
|---|
| 894 | endif |
|---|
| 895 | |
|---|
| 896 | do k=kts,kte |
|---|
| 897 | do j=jts,jte |
|---|
| 898 | do i=its,ite |
|---|
| 899 | X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0) |
|---|
| 900 | enddo |
|---|
| 901 | enddo |
|---|
| 902 | enddo |
|---|
| 903 | endif |
|---|
| 904 | |
|---|
| 905 | END SUBROUTINE negcor |
|---|
| 906 | |
|---|
| 907 | SUBROUTINE consat_s (ihail,itaobraun) |
|---|
| 908 | |
|---|
| 909 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 910 | ! c |
|---|
| 911 | ! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c |
|---|
| 912 | ! squall-type convective line. J. Atmos. Sci., 46, 177-202. c |
|---|
| 913 | ! c |
|---|
| 914 | ! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c |
|---|
| 915 | ! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c |
|---|
| 916 | ! c |
|---|
| 917 | ! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c |
|---|
| 918 | ! Model. Part I: Model description. Terrestrial, Atmospheric and c |
|---|
| 919 | ! Oceanic Sciences, 4, 35-72. c |
|---|
| 920 | ! c |
|---|
| 921 | ! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c |
|---|
| 922 | ! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c |
|---|
| 923 | ! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c |
|---|
| 924 | ! radiation and surface processes in the Goddard Cumulus Ensemble c |
|---|
| 925 | ! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c |
|---|
| 926 | ! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c |
|---|
| 927 | ! c |
|---|
| 928 | ! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c |
|---|
| 929 | ! Rutledge, and J. Simpson, 2007: Improving simulations of c |
|---|
| 930 | ! convective system from TRMM LBA: Easterly and Westerly regimes. c |
|---|
| 931 | ! J. Atmos. Sci., 64, 1141-1164. c |
|---|
| 932 | ! c |
|---|
| 933 | ! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c |
|---|
| 934 | ! c |
|---|
| 935 | ! Implemented into WRF by Roger Shi 2006/2007 c |
|---|
| 936 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 937 | |
|---|
| 938 | ! itaobraun=0 ! see Tao and Simpson (1993) |
|---|
| 939 | ! itaobraun=1 ! see Tao et al. (2003) |
|---|
| 940 | |
|---|
| 941 | integer :: itaobraun |
|---|
| 942 | real :: cn0 |
|---|
| 943 | |
|---|
| 944 | !JJS 1/3/2008 vvvvv |
|---|
| 945 | !JJS the following common blocks have been moved to the top of |
|---|
| 946 | !JJS module_mp_gsfcgce_driver_instat.F |
|---|
| 947 | ! |
|---|
| 948 | ! real, dimension (1:31) :: a1, a2 |
|---|
| 949 | ! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, & |
|---|
| 950 | ! .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, & |
|---|
| 951 | ! .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, & |
|---|
| 952 | ! .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, & |
|---|
| 953 | ! .6513e-6,.5956e-6,.5333e-6,.4834e-6/ |
|---|
| 954 | ! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, & |
|---|
| 955 | ! .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, & |
|---|
| 956 | ! .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, & |
|---|
| 957 | ! .4382,.4361/ |
|---|
| 958 | !JJS 1/3/2008 ^^^^^ |
|---|
| 959 | |
|---|
| 960 | |
|---|
| 961 | ! ****************************************************************** |
|---|
| 962 | !JJS |
|---|
| 963 | al = 2.5e10 |
|---|
| 964 | cp = 1.004e7 |
|---|
| 965 | rd1 = 1.e-3 |
|---|
| 966 | rd2 = 2.2 |
|---|
| 967 | !JJS |
|---|
| 968 | cpi=4.*atan(1.) |
|---|
| 969 | cpi2=cpi*cpi |
|---|
| 970 | grvt=980. |
|---|
| 971 | cd1=6.e-1 |
|---|
| 972 | cd2=4.*grvt/(3.*cd1) |
|---|
| 973 | tca=2.43e3 |
|---|
| 974 | dwv=.226 |
|---|
| 975 | dva=1.718e-4 |
|---|
| 976 | amw=18.016 |
|---|
| 977 | ars=8.314e7 |
|---|
| 978 | scv=2.2904487 |
|---|
| 979 | t0=273.16 |
|---|
| 980 | t00=238.16 |
|---|
| 981 | alv=2.5e10 |
|---|
| 982 | alf=3.336e9 |
|---|
| 983 | als=2.8336e10 |
|---|
| 984 | avc=alv/cp |
|---|
| 985 | afc=alf/cp |
|---|
| 986 | asc=als/cp |
|---|
| 987 | rw=4.615e6 |
|---|
| 988 | cw=4.187e7 |
|---|
| 989 | ci=2.093e7 |
|---|
| 990 | c76=7.66 |
|---|
| 991 | c358=35.86 |
|---|
| 992 | c172=17.26939 |
|---|
| 993 | c409=4098.026 |
|---|
| 994 | c218=21.87456 |
|---|
| 995 | c580=5807.695 |
|---|
| 996 | c610=6.1078e3 |
|---|
| 997 | c149=1.496286e-5 |
|---|
| 998 | c879=8.794142 |
|---|
| 999 | c141=1.4144354e7 |
|---|
| 1000 | !*** DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY |
|---|
| 1001 | !*** DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION |
|---|
| 1002 | !********** HAIL OR GRAUPEL PARAMETERS ********** |
|---|
| 1003 | if(ihail .eq. 1) then |
|---|
| 1004 | roqg=.9 |
|---|
| 1005 | ag=sqrt(cd2*roqg) |
|---|
| 1006 | bg=.5 |
|---|
| 1007 | tng=.002 |
|---|
| 1008 | else |
|---|
| 1009 | roqg=.4 |
|---|
| 1010 | ag=351.2 |
|---|
| 1011 | ! AG=372.3 ! if ice913=1 6/15/02 tao's |
|---|
| 1012 | bg=.37 |
|---|
| 1013 | tng=.04 |
|---|
| 1014 | endif |
|---|
| 1015 | !********** SNOW PARAMETERS ********** |
|---|
| 1016 | !ccshie 6/15/02 tao's |
|---|
| 1017 | ! TNS=1. |
|---|
| 1018 | ! TNS=.08 ! if ice913=1, tao's |
|---|
| 1019 | tns=.16 ! if ice913=0, tao's |
|---|
| 1020 | roqs=.1 |
|---|
| 1021 | ! AS=152.93 |
|---|
| 1022 | as=78.63 |
|---|
| 1023 | ! BS=.25 |
|---|
| 1024 | bs=.11 |
|---|
| 1025 | !********** RAIN PARAMETERS ********** |
|---|
| 1026 | aw=2115. |
|---|
| 1027 | bw=.8 |
|---|
| 1028 | roqr=1. |
|---|
| 1029 | tnw=.08 |
|---|
| 1030 | !***************************************************************** |
|---|
| 1031 | bgh=.5*bg |
|---|
| 1032 | bsh=.5*bs |
|---|
| 1033 | bwh=.5*bw |
|---|
| 1034 | bgq=.25*bg |
|---|
| 1035 | bsq=.25*bs |
|---|
| 1036 | bwq=.25*bw |
|---|
| 1037 | !**********GAMMA FUNCTION CALCULATIONS************* |
|---|
| 1038 | ga3b = gammagce(3.+bw) |
|---|
| 1039 | ga4b = gammagce(4.+bw) |
|---|
| 1040 | ga6b = gammagce(6.+bw) |
|---|
| 1041 | ga5bh = gammagce((5.+bw)/2.) |
|---|
| 1042 | ga3g = gammagce(3.+bg) |
|---|
| 1043 | ga4g = gammagce(4.+bg) |
|---|
| 1044 | ga5gh = gammagce((5.+bg)/2.) |
|---|
| 1045 | ga3d = gammagce(3.+bs) |
|---|
| 1046 | ga4d = gammagce(4.+bs) |
|---|
| 1047 | ga5dh = gammagce((5.+bs)/2.) |
|---|
| 1048 | !CCCCC LIN ET AL., 1983 OR LORD ET AL., 1984 CCCCCCCCCCCCCCCCC |
|---|
| 1049 | ac1=aw |
|---|
| 1050 | !JJS |
|---|
| 1051 | ac2=ag |
|---|
| 1052 | ac3=as |
|---|
| 1053 | !JJS |
|---|
| 1054 | bc1=bw |
|---|
| 1055 | cc1=as |
|---|
| 1056 | dc1=bs |
|---|
| 1057 | zrc=(cpi*roqr*tnw)**0.25 |
|---|
| 1058 | zsc=(cpi*roqs*tns)**0.25 |
|---|
| 1059 | zgc=(cpi*roqg*tng)**0.25 |
|---|
| 1060 | vrc=aw*ga4b/(6.*zrc**bw) |
|---|
| 1061 | vsc=as*ga4d/(6.*zsc**bs) |
|---|
| 1062 | vgc=ag*ga4g/(6.*zgc**bg) |
|---|
| 1063 | ! **************************** |
|---|
| 1064 | ! RN1=1.E-3 |
|---|
| 1065 | rn1=9.4e-15 ! 6/15/02 tao's |
|---|
| 1066 | bnd1=6.e-4 |
|---|
| 1067 | rn2=1.e-3 |
|---|
| 1068 | ! BND2=1.25E-3 |
|---|
| 1069 | ! BND2=1.5E-3 ! if ice913=1 6/15/02 tao's |
|---|
| 1070 | bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's |
|---|
| 1071 | rn3=.25*cpi*tns*cc1*ga3d |
|---|
| 1072 | esw=1. |
|---|
| 1073 | rn4=.25*cpi*esw*tns*cc1*ga3d |
|---|
| 1074 | ! ERI=1. |
|---|
| 1075 | eri=.1 ! 6/17/02 tao's ice913=0 (not 1) |
|---|
| 1076 | rn5=.25*cpi*eri*tnw*ac1*ga3b |
|---|
| 1077 | ! AMI=1./(24.*4.19E-10) |
|---|
| 1078 | ami=1./(24.*6.e-9) ! 6/15/02 tao's |
|---|
| 1079 | rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami |
|---|
| 1080 | ! ESR=1. ! also if ice913=1 for tao's |
|---|
| 1081 | esr=.5 ! 6/15/02 for ice913=0 tao's |
|---|
| 1082 | rn7=cpi2*esr*tnw*tns*roqs |
|---|
| 1083 | esr=1. |
|---|
| 1084 | rn8=cpi2*esr*tnw*tns*roqr |
|---|
| 1085 | rn9=cpi2*tns*tng*roqs |
|---|
| 1086 | rn10=2.*cpi*tns |
|---|
| 1087 | rn101=.31*ga5dh*sqrt(cc1) |
|---|
| 1088 | rn10a=als*als/rw |
|---|
| 1089 | !JJS |
|---|
| 1090 | rn10b=alv/tca |
|---|
| 1091 | rn10c=ars/(dwv*amw) |
|---|
| 1092 | !JJS |
|---|
| 1093 | rn11=2.*cpi*tns/alf |
|---|
| 1094 | rn11a=cw/alf |
|---|
| 1095 | ! AMI50=1.51e-7 |
|---|
| 1096 | ami50=3.84e-6 ! 6/15/02 tao's |
|---|
| 1097 | ! AMI40=2.41e-8 |
|---|
| 1098 | ami40=3.08e-8 ! 6/15/02 tao's |
|---|
| 1099 | eiw=1. |
|---|
| 1100 | ! UI50=20. |
|---|
| 1101 | ui50=100. ! 6/15/02 tao's |
|---|
| 1102 | ri50=2.*5.e-3 |
|---|
| 1103 | cmn=1.05e-15 |
|---|
| 1104 | rn12=cpi*eiw*ui50*ri50**2 |
|---|
| 1105 | |
|---|
| 1106 | do 10 k=1,31 |
|---|
| 1107 | y1=1.-aa2(k) |
|---|
| 1108 | rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1) |
|---|
| 1109 | rn12a(k)=rn13(k)/ami50 |
|---|
| 1110 | rn12b(k)=aa1(k)*ami50**aa2(k) |
|---|
| 1111 | rn25a(k)=aa1(k)*cmn**aa2(k) |
|---|
| 1112 | 10 continue |
|---|
| 1113 | |
|---|
| 1114 | egw=1. |
|---|
| 1115 | rn14=.25*cpi*egw*tng*ga3g*ag |
|---|
| 1116 | egi=.1 |
|---|
| 1117 | rn15=.25*cpi*egi*tng*ga3g*ag |
|---|
| 1118 | egi=1. |
|---|
| 1119 | rn15a=.25*cpi*egi*tng*ga3g*ag |
|---|
| 1120 | egr=1. |
|---|
| 1121 | rn16=cpi2*egr*tng*tnw*roqr |
|---|
| 1122 | rn17=2.*cpi*tng |
|---|
| 1123 | rn17a=.31*ga5gh*sqrt(ag) |
|---|
| 1124 | rn17b=cw-ci |
|---|
| 1125 | rn17c=cw |
|---|
| 1126 | apri=.66 |
|---|
| 1127 | bpri=1.e-4 |
|---|
| 1128 | bpri=0.5*bpri ! 6/17/02 tao's |
|---|
| 1129 | rn18=20.*cpi2*bpri*tnw*roqr |
|---|
| 1130 | rn18a=apri |
|---|
| 1131 | rn19=2.*cpi*tng/alf |
|---|
| 1132 | rn19a=.31*ga5gh*sqrt(ag) |
|---|
| 1133 | rn19b=cw/alf |
|---|
| 1134 | ! |
|---|
| 1135 | rnn191=.78 |
|---|
| 1136 | rnn192=.31*ga5gh*sqrt(ac2/dva) |
|---|
| 1137 | ! |
|---|
| 1138 | rn20=2.*cpi*tng |
|---|
| 1139 | rn20a=als*als/rw |
|---|
| 1140 | rn20b=.31*ga5gh*sqrt(ag) |
|---|
| 1141 | bnd3=2.e-3 |
|---|
| 1142 | rn21=1.e3*1.569e-12/0.15 |
|---|
| 1143 | erw=1. |
|---|
| 1144 | rn22=.25*cpi*erw*ac1*tnw*ga3b |
|---|
| 1145 | rn23=2.*cpi*tnw |
|---|
| 1146 | rn23a=.31*ga5bh*sqrt(ac1) |
|---|
| 1147 | rn23b=alv*alv/rw |
|---|
| 1148 | |
|---|
| 1149 | |
|---|
| 1150 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1151 | !cc |
|---|
| 1152 | !cc "c0" in routine "consat" (2d), "consatrh" (3d) |
|---|
| 1153 | !cc if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6 |
|---|
| 1154 | !cc if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8 |
|---|
| 1155 | |
|---|
| 1156 | if (itaobraun .eq. 0) then |
|---|
| 1157 | cn0=1.e-8 |
|---|
| 1158 | beta=-.6 |
|---|
| 1159 | elseif (itaobraun .eq. 1) then |
|---|
| 1160 | cn0=1.e-6 |
|---|
| 1161 | beta=-.46 |
|---|
| 1162 | endif |
|---|
| 1163 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1164 | ! CN0=1.E-6 |
|---|
| 1165 | ! CN0=1.E-8 ! 6/15/02 tao's |
|---|
| 1166 | ! BETA=-.46 |
|---|
| 1167 | ! BETA=-.6 ! 6/15/02 tao's |
|---|
| 1168 | |
|---|
| 1169 | rn25=cn0 |
|---|
| 1170 | rn30a=alv*als*amw/(tca*ars) |
|---|
| 1171 | rn30b=alv/tca |
|---|
| 1172 | rn30c=ars/(dwv*amw) |
|---|
| 1173 | rn31=1.e-17 |
|---|
| 1174 | |
|---|
| 1175 | rn32=4.*51.545e-4 |
|---|
| 1176 | ! |
|---|
| 1177 | rn30=2.*cpi*tng |
|---|
| 1178 | rnn30a=alv*alv*amw/(tca*ars) |
|---|
| 1179 | ! |
|---|
| 1180 | rn33=4.*tns |
|---|
| 1181 | rn331=.65 |
|---|
| 1182 | rn332=.44*sqrt(ac3/dva)*ga5dh |
|---|
| 1183 | ! |
|---|
| 1184 | |
|---|
| 1185 | return |
|---|
| 1186 | END SUBROUTINE consat_s |
|---|
| 1187 | |
|---|
| 1188 | SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, & |
|---|
| 1189 | new_ice_sat, id, & |
|---|
| 1190 | ptwrf, qvwrf, qlwrf, qrwrf, & |
|---|
| 1191 | qiwrf, qswrf, qgwrf, & |
|---|
| 1192 | rho_mks, pi_mks, p0_mks,itimestep, & |
|---|
| 1193 | ids,ide, jds,jde, kds,kde, & |
|---|
| 1194 | ims,ime, jms,jme, kms,kme, & |
|---|
| 1195 | its,ite, jts,jte, kts,kte & |
|---|
| 1196 | ) |
|---|
| 1197 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1198 | ! c |
|---|
| 1199 | ! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c |
|---|
| 1200 | ! squall-type convective line. J. Atmos. Sci., 46, 177-202. c |
|---|
| 1201 | ! c |
|---|
| 1202 | ! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c |
|---|
| 1203 | ! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c |
|---|
| 1204 | ! c |
|---|
| 1205 | ! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c |
|---|
| 1206 | ! Model. Part I: Model description. Terrestrial, Atmospheric and c |
|---|
| 1207 | ! Oceanic Sciences, 4, 35-72. c |
|---|
| 1208 | ! c |
|---|
| 1209 | ! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c |
|---|
| 1210 | ! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c |
|---|
| 1211 | ! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c |
|---|
| 1212 | ! radiation and surface processes in the Goddard Cumulus Ensemble c |
|---|
| 1213 | ! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c |
|---|
| 1214 | ! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c |
|---|
| 1215 | ! c |
|---|
| 1216 | ! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c |
|---|
| 1217 | ! Rutledge, and J. Simpson, 2007: Improving simulations of c |
|---|
| 1218 | ! convective system from TRMM LBA: Easterly and Westerly regimes. c |
|---|
| 1219 | ! J. Atmos. Sci., 64, 1141-1164. c |
|---|
| 1220 | ! c |
|---|
| 1221 | ! Tao, W.-K., J. J. Shi, S. Lang, C. Peters-Lidard, A. Hou, S. c |
|---|
| 1222 | ! Braun, and J. Simpson, 2007: New, improved bulk-microphysical c |
|---|
| 1223 | ! schemes for studying precipitation processes in WRF. Part I: c |
|---|
| 1224 | ! Comparisons with other schemes. to appear on Mon. Wea. Rev. C |
|---|
| 1225 | ! c |
|---|
| 1226 | ! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c |
|---|
| 1227 | ! c |
|---|
| 1228 | ! Implemented into WRF by Roger Shi 2006/2007 c |
|---|
| 1229 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1230 | ! |
|---|
| 1231 | ! COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES |
|---|
| 1232 | ! |
|---|
| 1233 | integer, parameter :: nt=2880, nt2=2*nt |
|---|
| 1234 | |
|---|
| 1235 | !cc using scott braun's way for pint, pidep computations |
|---|
| 1236 | integer :: itaobraun,ice2,ihail,new_ice_sat,id,istatmin |
|---|
| 1237 | integer :: itimestep |
|---|
| 1238 | real :: tairccri, cn0, dt |
|---|
| 1239 | !cc |
|---|
| 1240 | |
|---|
| 1241 | !JJS common/bxyz/ n,isec,nran,kt1,kt2 |
|---|
| 1242 | !JJS common/option/ lipps,ijkadv,istatmin,iwater,itoga,imlifting,lin, |
|---|
| 1243 | !JJS 1 irf,iadvh,irfg,ismg,id |
|---|
| 1244 | |
|---|
| 1245 | !JJS common/timestat/ ndt_stat,idq |
|---|
| 1246 | !JJS common/iice/ new_ice_sat |
|---|
| 1247 | !JJS common/bt/ dt,d2t,rijl2,dts,f5,rd1,rd2,bound,al,cp,ra,ck,ce,eps, |
|---|
| 1248 | !JJS 1 psfc,fcor,sec,aminut,rdt |
|---|
| 1249 | |
|---|
| 1250 | !JJS the following common blocks have been moved to the top of |
|---|
| 1251 | !JJS module_mp_gsfcgce_driver_instat.F |
|---|
| 1252 | |
|---|
| 1253 | ! common/bt/ rd1,rd2,al,cp |
|---|
| 1254 | ! |
|---|
| 1255 | ! |
|---|
| 1256 | ! common/bterv/ zrc,zgc,zsc,vrc,vgc,vsc |
|---|
| 1257 | ! common/size/ tnw,tns,tng,roqs,roqg,roqr |
|---|
| 1258 | ! common/cont/ c38,c358,c610,c149,c879,c172,c409,c76,c218,c580,c141 |
|---|
| 1259 | ! common/b3cs/ ag,bg,as,bs,aw,bw,bgh,bgq,bsh,bsq,bwh,bwq |
|---|
| 1260 | ! common/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2, & |
|---|
| 1261 | ! rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a, & |
|---|
| 1262 | ! rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17, & |
|---|
| 1263 | ! rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b, & |
|---|
| 1264 | ! bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b, & |
|---|
| 1265 | ! rn30c,rn31,beta,rn32 |
|---|
| 1266 | ! common/rsnw1/ rn10b,rn10c,rnn191,rnn192,rn30,rnn30a,rn33,rn331, & |
|---|
| 1267 | ! rn332 |
|---|
| 1268 | !JJS |
|---|
| 1269 | |
|---|
| 1270 | real, dimension (its:ite,jts:jte,kts:kte) :: fv |
|---|
| 1271 | real, dimension (its:ite,jts:jte,kts:kte) :: dpt, dqv |
|---|
| 1272 | real, dimension (its:ite,jts:jte,kts:kte) :: qcl, qrn, & |
|---|
| 1273 | qci, qcs, qcg |
|---|
| 1274 | !JJS 10/16/06 vvvv |
|---|
| 1275 | ! real dpt1(ims:ime,jms:jme,kms:kme) |
|---|
| 1276 | ! real dqv1(ims:ime,jms:jme,kms:kme), |
|---|
| 1277 | ! 1 qcl1(ims:ime,jms:jme,kms:kme) |
|---|
| 1278 | ! real qrn1(ims:ime,jms:jme,kms:kme), |
|---|
| 1279 | ! 1 qci1(ims:ime,jms:jme,kms:kme) |
|---|
| 1280 | ! real qcs1(ims:ime,jms:jme,kms:kme), |
|---|
| 1281 | ! 1 qcg1(ims:ime,jms:jme,kms:kme) |
|---|
| 1282 | !JJS 10/16/06 ^^^^ |
|---|
| 1283 | |
|---|
| 1284 | !JJS |
|---|
| 1285 | |
|---|
| 1286 | real, dimension (ims:ime, kms:kme, jms:jme) :: ptwrf, qvwrf |
|---|
| 1287 | real, dimension (ims:ime, kms:kme, jms:jme) :: qlwrf, qrwrf, & |
|---|
| 1288 | qiwrf, qswrf, qgwrf |
|---|
| 1289 | !JJS 10/16/06 vvvv |
|---|
| 1290 | ! real ptwrfold(ims:ime, kms:kme, jms:jme) |
|---|
| 1291 | ! real qvwrfold(ims:ime, kms:kme, jms:jme), |
|---|
| 1292 | ! 1 qlwrfold(ims:ime, kms:kme, jms:jme) |
|---|
| 1293 | ! real qrwrfold(ims:ime, kms:kme, jms:jme), |
|---|
| 1294 | ! 1 qiwrfold(ims:ime, kms:kme, jms:jme) |
|---|
| 1295 | ! real qswrfold(ims:ime, kms:kme, jms:jme), |
|---|
| 1296 | ! 1 qgwrfold(ims:ime, kms:kme, jms:jme) |
|---|
| 1297 | !JJS 10/16/06 ^^^^ |
|---|
| 1298 | |
|---|
| 1299 | !JJS in MKS |
|---|
| 1300 | real, dimension (ims:ime, kms:kme, jms:jme) :: rho_mks |
|---|
| 1301 | real, dimension (ims:ime, kms:kme, jms:jme) :: pi_mks |
|---|
| 1302 | real, dimension (ims:ime, kms:kme, jms:jme) :: p0_mks |
|---|
| 1303 | !JJS |
|---|
| 1304 | ! real, dimension (its:ite,jts:jte,kts:kte) :: ww1 |
|---|
| 1305 | ! real, dimension (its:ite,jts:jte,kts:kte) :: rsw |
|---|
| 1306 | ! real, dimension (its:ite,jts:jte,kts:kte) :: rlw |
|---|
| 1307 | |
|---|
| 1308 | !JJS COMMON /BADV/ |
|---|
| 1309 | real, dimension (its:ite,jts:jte) :: & |
|---|
| 1310 | vg, zg, & |
|---|
| 1311 | ps, pg, & |
|---|
| 1312 | prn, psn, & |
|---|
| 1313 | pwacs, wgacr, & |
|---|
| 1314 | pidep, pint, & |
|---|
| 1315 | qsi, ssi, & |
|---|
| 1316 | esi, esw, & |
|---|
| 1317 | qsw, pr, & |
|---|
| 1318 | ssw, pihom, & |
|---|
| 1319 | pidw, pimlt, & |
|---|
| 1320 | psaut, qracs, & |
|---|
| 1321 | psaci, psacw, & |
|---|
| 1322 | qsacw, praci, & |
|---|
| 1323 | pmlts, pmltg, & |
|---|
| 1324 | asss, y1, y2 |
|---|
| 1325 | !JJS Y2(its:ite,jts:jte), DDE(NB) |
|---|
| 1326 | |
|---|
| 1327 | !JJS COMMON/BSAT/ |
|---|
| 1328 | real, dimension (its:ite,jts:jte) :: & |
|---|
| 1329 | praut, pracw, & |
|---|
| 1330 | psfw, psfi, & |
|---|
| 1331 | dgacs, dgacw, & |
|---|
| 1332 | dgaci, dgacr, & |
|---|
| 1333 | pgacs, wgacs, & |
|---|
| 1334 | qgacw, wgaci, & |
|---|
| 1335 | qgacr, pgwet, & |
|---|
| 1336 | pgaut, pracs, & |
|---|
| 1337 | psacr, qsacr, & |
|---|
| 1338 | pgfr, psmlt, & |
|---|
| 1339 | pgmlt, psdep, & |
|---|
| 1340 | pgdep, piacr, & |
|---|
| 1341 | y5, scv, & |
|---|
| 1342 | tca, dwv, & |
|---|
| 1343 | egs, y3, & |
|---|
| 1344 | y4, ddb |
|---|
| 1345 | |
|---|
| 1346 | !JJS COMMON/BSAT1/ |
|---|
| 1347 | real, dimension (its:ite,jts:jte) :: & |
|---|
| 1348 | pt, qv, & |
|---|
| 1349 | qc, qr, & |
|---|
| 1350 | qi, qs, & |
|---|
| 1351 | qg, tair, & |
|---|
| 1352 | tairc, rtair, & |
|---|
| 1353 | dep, dd, & |
|---|
| 1354 | dd1, qvs, & |
|---|
| 1355 | dm, rq, & |
|---|
| 1356 | rsub1, col, & |
|---|
| 1357 | cnd, ern, & |
|---|
| 1358 | dlt1, dlt2, & |
|---|
| 1359 | dlt3, dlt4, & |
|---|
| 1360 | zr, vr, & |
|---|
| 1361 | zs, vs, & |
|---|
| 1362 | dbz, pssub, & |
|---|
| 1363 | pgsub, dda |
|---|
| 1364 | |
|---|
| 1365 | !JJS COMMON/B5/ |
|---|
| 1366 | real, dimension (its:ite,jts:jte,kts:kte) :: rho |
|---|
| 1367 | real, dimension (kts:kte) :: & |
|---|
| 1368 | tb, qb, rho1, & |
|---|
| 1369 | ta, qa, ta1, qa1, & |
|---|
| 1370 | coef, z1, z2, z3, & |
|---|
| 1371 | am, am1, ub, vb, & |
|---|
| 1372 | wb, ub1, vb1, rrho, & |
|---|
| 1373 | rrho1, wbx |
|---|
| 1374 | |
|---|
| 1375 | !JJS COMMON/B6/ |
|---|
| 1376 | real, dimension (its:ite,jts:jte,kts:kte) :: p0, pi, f0 |
|---|
| 1377 | real, dimension (kts:kte) :: & |
|---|
| 1378 | fd, fe, & |
|---|
| 1379 | st, sv, & |
|---|
| 1380 | sq, sc, & |
|---|
| 1381 | se, sqa |
|---|
| 1382 | |
|---|
| 1383 | !JJS COMMON/BRH1/ |
|---|
| 1384 | real, dimension (kts:kte) :: & |
|---|
| 1385 | srro, qrro, sqc, sqr, & |
|---|
| 1386 | sqi, sqs, sqg, stqc, & |
|---|
| 1387 | stqr, stqi, stqs, stqg |
|---|
| 1388 | real, dimension (nt) :: & |
|---|
| 1389 | tqc, tqr, tqi, tqs, tqg |
|---|
| 1390 | |
|---|
| 1391 | !JJS common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny) |
|---|
| 1392 | |
|---|
| 1393 | !JJS COMMON/BLS/ |
|---|
| 1394 | real, dimension (its:ite,jts:jte) :: & |
|---|
| 1395 | y0, ts0, qss0 |
|---|
| 1396 | |
|---|
| 1397 | !JJS COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4) |
|---|
| 1398 | integer, dimension (its:ite,jts:jte) :: it |
|---|
| 1399 | integer, dimension (its:ite,jts:jte, 4) :: ics |
|---|
| 1400 | |
|---|
| 1401 | !JJS COMMON/MICRO/ |
|---|
| 1402 | real, dimension (its:ite,kts:kte,jts:jte) :: & |
|---|
| 1403 | physc, physe, physd, & |
|---|
| 1404 | physs, physm, physf |
|---|
| 1405 | |
|---|
| 1406 | !23456789012345678901234567890123456789012345678901234567890123456789012 |
|---|
| 1407 | |
|---|
| 1408 | ! |
|---|
| 1409 | !JJS 1/3/2008 vvvvv |
|---|
| 1410 | !JJS the following common blocks have been moved to the top of |
|---|
| 1411 | !JJS module_mp_gsfcgce_driver.F |
|---|
| 1412 | |
|---|
| 1413 | ! real, dimension (31) :: aa1, aa2 |
|---|
| 1414 | ! data aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, & |
|---|
| 1415 | ! .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, & |
|---|
| 1416 | ! .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, & |
|---|
| 1417 | ! .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, & |
|---|
| 1418 | ! .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, & |
|---|
| 1419 | ! .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, & |
|---|
| 1420 | ! .4834e-6/ |
|---|
| 1421 | ! data aa2/.4006, .4831, .5320, .5307, .5319, & |
|---|
| 1422 | ! .5249, .4888, .3894, .4047, .4318, & |
|---|
| 1423 | ! .4771, .5183, .5463, .5651, .5813, & |
|---|
| 1424 | ! .5655, .5478, .5203, .4906, .4447, & |
|---|
| 1425 | ! .4126, .3960, .4149, .4320, .4506, & |
|---|
| 1426 | ! .4483, .4460, .4433, .4413, .4382, & |
|---|
| 1427 | ! .4361/ |
|---|
| 1428 | |
|---|
| 1429 | !JJS 1/3/2008 ^^^^^ |
|---|
| 1430 | |
|---|
| 1431 | ! |
|---|
| 1432 | save |
|---|
| 1433 | |
|---|
| 1434 | ! write(6, *) 'in gsfcgce_open_ice2 at itimestep=',itimestep |
|---|
| 1435 | ! write(6, *) 'ims=', ims, ' ime=', ime |
|---|
| 1436 | ! write(6, *) 'jms=', jms, ' jme=', jme |
|---|
| 1437 | ! write(6, *) 'kms=', kms, ' kme=', kme |
|---|
| 1438 | ! write(6, *) 'its=', its, ' ite=', ite |
|---|
| 1439 | ! write(6, *) 'jts=', jts, ' jte=', jte |
|---|
| 1440 | ! write(6, *) 'kts=', kts, ' kte=', kte |
|---|
| 1441 | ! write(6, *) 'ihail=', ihail |
|---|
| 1442 | ! write(6, *) 'itaobraun=',itaobraun |
|---|
| 1443 | ! write(6, *) 'ICE2=', ICE2 |
|---|
| 1444 | ! write(6, *) 'istatmin=',istatmin |
|---|
| 1445 | ! write(6, *) 'new_ice_sat=', new_ice_sat |
|---|
| 1446 | ! write(6, *) 'id=', id |
|---|
| 1447 | ! write(6, *) 'dt=', dt |
|---|
| 1448 | ! endif |
|---|
| 1449 | |
|---|
| 1450 | !JJS convert from mks to cgs, and move from WRF grid to GCE grid |
|---|
| 1451 | do k=kts,kte |
|---|
| 1452 | do j=jts,jte |
|---|
| 1453 | do i=its,ite |
|---|
| 1454 | rho(i,j,k)=rho_mks(i,k,j)*0.001 |
|---|
| 1455 | p0(i,j,k)=p0_mks(i,k,j)*10.0 |
|---|
| 1456 | pi(i,j,k)=pi_mks(i,k,j) |
|---|
| 1457 | dpt(i,j,k)=ptwrf(i,k,j) |
|---|
| 1458 | dqv(i,j,k)=qvwrf(i,k,j) |
|---|
| 1459 | qcl(i,j,k)=qlwrf(i,k,j) |
|---|
| 1460 | qrn(i,j,k)=qrwrf(i,k,j) |
|---|
| 1461 | qci(i,j,k)=qiwrf(i,k,j) |
|---|
| 1462 | qcs(i,j,k)=qswrf(i,k,j) |
|---|
| 1463 | qcg(i,j,k)=qgwrf(i,k,j) |
|---|
| 1464 | !JJS 10/16/06 vvvv |
|---|
| 1465 | ! dpt1(i,j,k)=ptwrfold(i,k,j) |
|---|
| 1466 | ! dqv1(i,j,k)=qvwrfold(i,k,j) |
|---|
| 1467 | ! qcl1(i,j,k)=qlwrfold(i,k,j) |
|---|
| 1468 | ! qrn1(i,j,k)=qrwrfold(i,k,j) |
|---|
| 1469 | ! qci1(i,j,k)=qiwrfold(i,k,j) |
|---|
| 1470 | ! qcs1(i,j,k)=qswrfold(i,k,j) |
|---|
| 1471 | ! qcg1(i,j,k)=qgwrfold(i,k,j) |
|---|
| 1472 | !JJS 10/16/06 ^^^^ |
|---|
| 1473 | enddo !i |
|---|
| 1474 | enddo !j |
|---|
| 1475 | enddo !k |
|---|
| 1476 | |
|---|
| 1477 | if (ice2 .eq. 2) ihail = 0 |
|---|
| 1478 | |
|---|
| 1479 | do k=kts,kte |
|---|
| 1480 | do j=jts,jte |
|---|
| 1481 | do i=its,ite |
|---|
| 1482 | fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k)) |
|---|
| 1483 | enddo !i |
|---|
| 1484 | enddo !j |
|---|
| 1485 | enddo !k |
|---|
| 1486 | !JJS |
|---|
| 1487 | |
|---|
| 1488 | ! |
|---|
| 1489 | ! ****** THREE CLASSES OF ICE-PHASE (LIN ET AL, 1983) ********* |
|---|
| 1490 | |
|---|
| 1491 | !JJS D22T=D2T |
|---|
| 1492 | !JJS IF(IJKADV .EQ. 0) THEN |
|---|
| 1493 | !JJS D2T=D2T |
|---|
| 1494 | !JJS ELSE |
|---|
| 1495 | d2t=dt |
|---|
| 1496 | !JJS ENDIF |
|---|
| 1497 | ! |
|---|
| 1498 | ! itaobraun=0 ! original pint and pidep & see Tao and Simpson 1993 |
|---|
| 1499 | itaobraun=1 ! see Tao et al. (2003) |
|---|
| 1500 | ! |
|---|
| 1501 | if ( itaobraun.eq.0 ) then |
|---|
| 1502 | cn0=1.e-8 |
|---|
| 1503 | !c beta=-.6 |
|---|
| 1504 | elseif ( itaobraun.eq.1 ) then |
|---|
| 1505 | cn0=1.e-6 |
|---|
| 1506 | ! cn0=1.e-8 ! special |
|---|
| 1507 | !c beta=-.46 |
|---|
| 1508 | endif |
|---|
| 1509 | !C TAO 2007 START |
|---|
| 1510 | ! ICE2=2 ! 2ICE with cloud ice and graupel (no snow) - r2ice=1, r2iceg=0. |
|---|
| 1511 | ! ICE2=1 ! 2ice with cloud ice and snow (no graupel) - r2iceg=1, r2ice=0. |
|---|
| 1512 | !c |
|---|
| 1513 | r2ice=1. |
|---|
| 1514 | r2iceg=1. |
|---|
| 1515 | if (ice2 .eq. 1) then |
|---|
| 1516 | r2ice=0. |
|---|
| 1517 | r2iceg=1. |
|---|
| 1518 | endif |
|---|
| 1519 | if (ice2 .eq. 2) then |
|---|
| 1520 | r2ice=1. |
|---|
| 1521 | r2iceg=0. |
|---|
| 1522 | endif |
|---|
| 1523 | !C TAO 2007 END |
|---|
| 1524 | cmin=1.e-19 |
|---|
| 1525 | cmin1=1.e-20 |
|---|
| 1526 | cmin2=1.e-12 |
|---|
| 1527 | ucor=3071.29/tnw**.75 |
|---|
| 1528 | ucos=687.97*roqs**.25/tns**.75 |
|---|
| 1529 | ucog=687.97*roqg**.25/tng**.75 |
|---|
| 1530 | uwet=4.464**.95 |
|---|
| 1531 | |
|---|
| 1532 | rijl2 = 1. / (ide-ids) / (jde-jds) |
|---|
| 1533 | |
|---|
| 1534 | !JJScap $doacross local(j,i) |
|---|
| 1535 | |
|---|
| 1536 | !JJS DO 1 J=1,JMAX |
|---|
| 1537 | !JJS DO 1 I=1,IMAX |
|---|
| 1538 | do j=jts,jte |
|---|
| 1539 | do i=its,ite |
|---|
| 1540 | it(i,j)=1 |
|---|
| 1541 | enddo |
|---|
| 1542 | enddo |
|---|
| 1543 | |
|---|
| 1544 | f2=rd1*d2t |
|---|
| 1545 | f3=rd2*d2t |
|---|
| 1546 | |
|---|
| 1547 | ft=dt/d2t |
|---|
| 1548 | rft=rijl2*ft |
|---|
| 1549 | a0=.5*istatmin*rijl2 |
|---|
| 1550 | rt0=1./(t0-t00) |
|---|
| 1551 | bw3=bw+3. |
|---|
| 1552 | bs3=bs+3. |
|---|
| 1553 | bg3=bg+3. |
|---|
| 1554 | bsh5=2.5+bsh |
|---|
| 1555 | bgh5=2.5+bgh |
|---|
| 1556 | bwh5=2.5+bwh |
|---|
| 1557 | bw6=bw+6. |
|---|
| 1558 | bs6=bs+6. |
|---|
| 1559 | betah=.5*beta |
|---|
| 1560 | r10t=rn10*d2t |
|---|
| 1561 | r11at=rn11a*d2t |
|---|
| 1562 | r19bt=rn19b*d2t |
|---|
| 1563 | r20t=-rn20*d2t |
|---|
| 1564 | r23t=-rn23*d2t |
|---|
| 1565 | r25a=rn25 |
|---|
| 1566 | |
|---|
| 1567 | ! ami50 for use in PINT |
|---|
| 1568 | ami50=3.76e-8 |
|---|
| 1569 | ami100=1.51e-7 |
|---|
| 1570 | ami40=2.41e-8 |
|---|
| 1571 | |
|---|
| 1572 | !C ****************************************************************** |
|---|
| 1573 | |
|---|
| 1574 | !JJS DO 1000 K=2,kles |
|---|
| 1575 | do 1000 k=kts,kte |
|---|
| 1576 | kp=k+1 |
|---|
| 1577 | !JJS tb0=ta1(k) |
|---|
| 1578 | !JJS qb0=qa1(k) |
|---|
| 1579 | tb0=0. |
|---|
| 1580 | qb0=0. |
|---|
| 1581 | |
|---|
| 1582 | do 2000 j=jts,jte |
|---|
| 1583 | do 2000 i=its,ite |
|---|
| 1584 | |
|---|
| 1585 | rp0=3.799052e3/p0(i,j,k) |
|---|
| 1586 | pi0=pi(i,j,k) |
|---|
| 1587 | pir=1./(pi(i,j,k)) |
|---|
| 1588 | pr0=1./p0(i,j,k) |
|---|
| 1589 | r00=rho(i,j,k) |
|---|
| 1590 | r0s=sqrt(rho(i,j,k)) |
|---|
| 1591 | !JJS RR0=RRHO(K) |
|---|
| 1592 | rr0=1./rho(i,j,k) |
|---|
| 1593 | !JJS RRS=SRRO(K) |
|---|
| 1594 | rrs=sqrt(rr0) |
|---|
| 1595 | !JJS RRQ=QRRO(K) |
|---|
| 1596 | rrq=sqrt(rrs) |
|---|
| 1597 | f0(i,j,k)=al/cp/pi(i,j,k) |
|---|
| 1598 | f00=f0(i,j,k) |
|---|
| 1599 | fv0=fv(i,j,k) |
|---|
| 1600 | fvs=sqrt(fv(i,j,k)) |
|---|
| 1601 | zrr=1.e5*zrc*rrq |
|---|
| 1602 | zsr=1.e5*zsc*rrq |
|---|
| 1603 | zgr=1.e5*zgc*rrq |
|---|
| 1604 | cp409=c409*pi0 |
|---|
| 1605 | cv409=c409*avc |
|---|
| 1606 | cp580=c580*pi0 |
|---|
| 1607 | cs580=c580*asc |
|---|
| 1608 | alvr=r00*alv |
|---|
| 1609 | afcp=afc*pir |
|---|
| 1610 | avcp=avc*pir |
|---|
| 1611 | ascp=asc*pir |
|---|
| 1612 | vrcf=vrc*fv0 |
|---|
| 1613 | vscf=vsc*fv0 |
|---|
| 1614 | vgcf=vgc*fv0 |
|---|
| 1615 | vgcr=vgc*rrs |
|---|
| 1616 | dwvp=c879*pr0 |
|---|
| 1617 | r3f=rn3*fv0 |
|---|
| 1618 | r4f=rn4*fv0 |
|---|
| 1619 | r5f=rn5*fv0 |
|---|
| 1620 | r6f=rn6*fv0 |
|---|
| 1621 | r7r=rn7*rr0 |
|---|
| 1622 | r8r=rn8*rr0 |
|---|
| 1623 | r9r=rn9*rr0 |
|---|
| 1624 | r101f=rn101*fvs |
|---|
| 1625 | r10ar=rn10a*r00 |
|---|
| 1626 | r11rt=rn11*rr0*d2t |
|---|
| 1627 | r12r=rn12*r00 |
|---|
| 1628 | r14r=rn14*rrs |
|---|
| 1629 | r14f=rn14*fv0 |
|---|
| 1630 | r15r=rn15*rrs |
|---|
| 1631 | r15ar=rn15a*rrs |
|---|
| 1632 | r15f=rn15*fv0 |
|---|
| 1633 | r15af=rn15a*fv0 |
|---|
| 1634 | r16r=rn16*rr0 |
|---|
| 1635 | r17r=rn17*rr0 |
|---|
| 1636 | r17aq=rn17a*rrq |
|---|
| 1637 | r17as=rn17a*fvs |
|---|
| 1638 | r18r=rn18*rr0 |
|---|
| 1639 | r19rt=rn19*rr0*d2t |
|---|
| 1640 | r19aq=rn19a*rrq |
|---|
| 1641 | r19as=rn19a*fvs |
|---|
| 1642 | r20bq=rn20b*rrq |
|---|
| 1643 | r20bs=rn20b*fvs |
|---|
| 1644 | r22f=rn22*fv0 |
|---|
| 1645 | r23af=rn23a*fvs |
|---|
| 1646 | r23br=rn23b*r00 |
|---|
| 1647 | r25rt=rn25*rr0*d2t |
|---|
| 1648 | r31r=rn31*rr0 |
|---|
| 1649 | r32rt=rn32*d2t*rrs |
|---|
| 1650 | |
|---|
| 1651 | !JJS DO 100 J=3,JLES |
|---|
| 1652 | !JJS DO 100 I=3,ILES |
|---|
| 1653 | pt(i,j)=dpt(i,j,k) |
|---|
| 1654 | qv(i,j)=dqv(i,j,k) |
|---|
| 1655 | qc(i,j)=qcl(i,j,k) |
|---|
| 1656 | qr(i,j)=qrn(i,j,k) |
|---|
| 1657 | qi(i,j)=qci(i,j,k) |
|---|
| 1658 | qs(i,j)=qcs(i,j,k) |
|---|
| 1659 | qg(i,j)=qcg(i,j,k) |
|---|
| 1660 | ! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0 |
|---|
| 1661 | if (qc(i,j) .le. cmin1) qc(i,j)=0.0 |
|---|
| 1662 | if (qr(i,j) .le. cmin1) qr(i,j)=0.0 |
|---|
| 1663 | if (qi(i,j) .le. cmin1) qi(i,j)=0.0 |
|---|
| 1664 | if (qs(i,j) .le. cmin1) qs(i,j)=0.0 |
|---|
| 1665 | if (qg(i,j) .le. cmin1) qg(i,j)=0.0 |
|---|
| 1666 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 1667 | tairc(i,j)=tair(i,j)-t0 |
|---|
| 1668 | zr(i,j)=zrr |
|---|
| 1669 | zs(i,j)=zsr |
|---|
| 1670 | zg(i,j)=zgr |
|---|
| 1671 | vr(i,j)=0.0 |
|---|
| 1672 | vs(i,j)=0.0 |
|---|
| 1673 | vg(i,j)=0.0 |
|---|
| 1674 | |
|---|
| 1675 | ! *** COMPUTE ZR,ZS,ZG,VR,VS,VG ***************************** |
|---|
| 1676 | |
|---|
| 1677 | if (qr(i,j) .gt. cmin1) then |
|---|
| 1678 | dd(i,j)=r00*qr(i,j) |
|---|
| 1679 | y1(i,j)=dd(i,j)**.25 |
|---|
| 1680 | zr(i,j)=zrc/y1(i,j) |
|---|
| 1681 | vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.) |
|---|
| 1682 | endif |
|---|
| 1683 | |
|---|
| 1684 | if (qs(i,j) .gt. cmin1) then |
|---|
| 1685 | dd(i,j)=r00*qs(i,j) |
|---|
| 1686 | y1(i,j)=dd(i,j)**.25 |
|---|
| 1687 | zs(i,j)=zsc/y1(i,j) |
|---|
| 1688 | vs(i,j)=max(vscf*dd(i,j)**bsq, 0.) |
|---|
| 1689 | endif |
|---|
| 1690 | |
|---|
| 1691 | if (qg(i,j) .gt. cmin1) then |
|---|
| 1692 | dd(i,j)=r00*qg(i,j) |
|---|
| 1693 | y1(i,j)=dd(i,j)**.25 |
|---|
| 1694 | zg(i,j)=zgc/y1(i,j) |
|---|
| 1695 | if(ihail .eq. 1) then |
|---|
| 1696 | vg(i,j)=max(vgcr*dd(i,j)**bgq, 0.) |
|---|
| 1697 | else |
|---|
| 1698 | vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.) |
|---|
| 1699 | endif |
|---|
| 1700 | endif |
|---|
| 1701 | |
|---|
| 1702 | if (qr(i,j) .le. cmin2) vr(i,j)=0.0 |
|---|
| 1703 | if (qs(i,j) .le. cmin2) vs(i,j)=0.0 |
|---|
| 1704 | if (qg(i,j) .le. cmin2) vg(i,j)=0.0 |
|---|
| 1705 | |
|---|
| 1706 | ! ****************************************************************** |
|---|
| 1707 | ! *** Y1 : DYNAMIC VISCOSITY OF AIR (U) |
|---|
| 1708 | ! *** DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI) |
|---|
| 1709 | ! *** TCA : THERMAL CONDUCTIVITY OF AIR (KA) |
|---|
| 1710 | ! *** Y2 : KINETIC VISCOSITY (V) |
|---|
| 1711 | |
|---|
| 1712 | y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.) |
|---|
| 1713 | dwv(i,j)=dwvp*tair(i,j)**1.81 |
|---|
| 1714 | tca(i,j)=c141*y1(i,j) |
|---|
| 1715 | scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333) |
|---|
| 1716 | !JJS 100 CONTINUE |
|---|
| 1717 | |
|---|
| 1718 | !* 1 * PSAUT : AUTOCONVERSION OF QI TO QS ***1** |
|---|
| 1719 | !* 3 * PSACI : ACCRETION OF QI TO QS ***3** |
|---|
| 1720 | !* 4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT) ***4** |
|---|
| 1721 | !* 5 * PRACI : ACCRETION OF QI BY QR ***5** |
|---|
| 1722 | !* 6 * PIACR : ACCRETION OF QR OR QG BY QI ***6** |
|---|
| 1723 | |
|---|
| 1724 | !JJS DO 125 J=3,JLES |
|---|
| 1725 | !JJS DO 125 I=3,ILES |
|---|
| 1726 | psaut(i,j)=0.0 |
|---|
| 1727 | psaci(i,j)=0.0 |
|---|
| 1728 | praci(i,j)=0.0 |
|---|
| 1729 | piacr(i,j)=0.0 |
|---|
| 1730 | psacw(i,j)=0.0 |
|---|
| 1731 | qsacw(i,j)=0.0 |
|---|
| 1732 | dd(i,j)=1./zs(i,j)**bs3 |
|---|
| 1733 | |
|---|
| 1734 | if (tair(i,j).lt.t0) then |
|---|
| 1735 | esi(i,j)=exp(.025*tairc(i,j)) |
|---|
| 1736 | psaut(i,j)=r2iceg*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0) |
|---|
| 1737 | psaci(i,j)=r2iceg*r3f*esi(i,j)*qi(i,j)*dd(i,j) |
|---|
| 1738 | !JJS 3/30/06 |
|---|
| 1739 | ! to cut water to snow accretion by half |
|---|
| 1740 | ! PSACW(I,J)=R4F*QC(I,J)*DD(I,J) |
|---|
| 1741 | psacw(i,j)=r2iceg*0.5*r4f*qc(i,j)*dd(i,j) |
|---|
| 1742 | !JJS 3/30/06 |
|---|
| 1743 | praci(i,j)=r2iceg*r5f*qi(i,j)/zr(i,j)**bw3 |
|---|
| 1744 | piacr(i,j)=r2iceg*r6f*qi(i,j)*(zr(i,j)**(-bw6)) |
|---|
| 1745 | !JJS PIACR(I,J)=R6F*QI(I,J)/ZR(I,J)**BW6 |
|---|
| 1746 | else |
|---|
| 1747 | qsacw(i,j)=r2iceg*r4f*qc(i,j)*dd(i,j) |
|---|
| 1748 | endif |
|---|
| 1749 | |
|---|
| 1750 | !* 21 * PRAUT AUTOCONVERSION OF QC TO QR **21** |
|---|
| 1751 | !* 22 * PRACW : ACCRETION OF QC BY QR **22** |
|---|
| 1752 | |
|---|
| 1753 | pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3 |
|---|
| 1754 | praut(i,j)=0.0 |
|---|
| 1755 | y1(i,j)=qc(i,j)-bnd3 |
|---|
| 1756 | if (y1(i,j).gt.0.0) then |
|---|
| 1757 | praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j)) |
|---|
| 1758 | endif |
|---|
| 1759 | |
|---|
| 1760 | !* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971) **12** |
|---|
| 1761 | !* 13 * PSFI : BERGERON PROCESSES FOR QS **13** |
|---|
| 1762 | |
|---|
| 1763 | psfw(i,j)=0.0 |
|---|
| 1764 | psfi(i,j)=0.0 |
|---|
| 1765 | pidep(i,j)=0.0 |
|---|
| 1766 | |
|---|
| 1767 | if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then |
|---|
| 1768 | y1(i,j)=max( min(tairc(i,j), -1.), -31.) |
|---|
| 1769 | it(i,j)=int(abs(y1(i,j))) |
|---|
| 1770 | y1(i,j)=rn12a(it(i,j)) |
|---|
| 1771 | y2(i,j)=rn12b(it(i,j)) |
|---|
| 1772 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1773 | psfw(i,j)=r2iceg* & |
|---|
| 1774 | max(d2t*y1(i,j)*(y2(i,j)+r12r*qc(i,j))*qi(i,j),0.0) |
|---|
| 1775 | rtair(i,j)=1./(tair(i,j)-c76) |
|---|
| 1776 | y2(i,j)=exp(c218-c580*rtair(i,j)) |
|---|
| 1777 | qsi(i,j)=rp0*y2(i,j) |
|---|
| 1778 | esi(i,j)=c610*y2(i,j) |
|---|
| 1779 | ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. |
|---|
| 1780 | r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.) |
|---|
| 1781 | ! R_NCI=min(1.e-8*EXP(-.6*TAIRC(I,J)),1.) ! use Tao's |
|---|
| 1782 | dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.) |
|---|
| 1783 | rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) |
|---|
| 1784 | y3(i,j)=1./tair(i,j) |
|---|
| 1785 | dd(i,j)=y3(i,j)*(rn30a*y3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j) |
|---|
| 1786 | y1(i,j)=206.18*ssi(i,j)/dd(i,j) |
|---|
| 1787 | pidep(i,j)=y1(i,j)*sqrt(r_nci*qi(i,j)/r00) |
|---|
| 1788 | dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t |
|---|
| 1789 | if(dm(i,j).gt.cmin2) then |
|---|
| 1790 | a2=1. |
|---|
| 1791 | if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then |
|---|
| 1792 | a2=dep(i,j)/pidep(i,j) |
|---|
| 1793 | pidep(i,j)=dep(i,j) |
|---|
| 1794 | endif |
|---|
| 1795 | psfi(i,j)=r2iceg*a2*.5*qi(i,j)*y1(i,j)/(sqrt(ami100) & |
|---|
| 1796 | -sqrt(ami40)) |
|---|
| 1797 | elseif(dm(i,j).lt.-cmin2) then |
|---|
| 1798 | ! |
|---|
| 1799 | ! SUBLIMATION TERMS USED ONLY WHEN SATURATION ADJUSTMENT FOR ICE |
|---|
| 1800 | ! IS TURNED OFF |
|---|
| 1801 | ! |
|---|
| 1802 | pidep(i,j)=0. |
|---|
| 1803 | psfi(i,j)=0. |
|---|
| 1804 | else |
|---|
| 1805 | pidep(i,j)=0. |
|---|
| 1806 | psfi(i,j)=0. |
|---|
| 1807 | endif |
|---|
| 1808 | endif |
|---|
| 1809 | |
|---|
| 1810 | !TTT***** QG=QG+MIN(PGDRY,PGWET) |
|---|
| 1811 | !* 9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET) ***9** |
|---|
| 1812 | !* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT) **14** |
|---|
| 1813 | !* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT) **16** |
|---|
| 1814 | |
|---|
| 1815 | if(qc(i,j)+qr(i,j).lt.1.e-4) then |
|---|
| 1816 | ee1=.01 |
|---|
| 1817 | else |
|---|
| 1818 | ee1=1. |
|---|
| 1819 | endif |
|---|
| 1820 | ee2=0.09 |
|---|
| 1821 | egs(i,j)=ee1*exp(ee2*tairc(i,j)) |
|---|
| 1822 | ! EGS(I,J)=0.1 ! 6/15/02 tao's |
|---|
| 1823 | if (tair(i,j).ge.t0) egs(i,j)=1.0 |
|---|
| 1824 | y1(i,j)=abs(vg(i,j)-vs(i,j)) |
|---|
| 1825 | y2(i,j)=zs(i,j)*zg(i,j) |
|---|
| 1826 | y3(i,j)=5./y2(i,j) |
|---|
| 1827 | y4(i,j)=.08*y3(i,j)*y3(i,j) |
|---|
| 1828 | y5(i,j)=.05*y3(i,j)*y4(i,j) |
|---|
| 1829 | dd(i,j)=y1(i,j)*(y3(i,j)/zs(i,j)**5+y4(i,j)/zs(i,j)**3 & |
|---|
| 1830 | +y5(i,j)/zs(i,j)) |
|---|
| 1831 | pgacs(i,j)=r2ice*r2iceg*r9r*egs(i,j)*dd(i,j) |
|---|
| 1832 | !JJS 1/3/06 from Steve and Chunglin |
|---|
| 1833 | if (ihail.eq.1) then |
|---|
| 1834 | dgacs(i,j)=pgacs(i,j) |
|---|
| 1835 | else |
|---|
| 1836 | dgacs(i,j)=0. |
|---|
| 1837 | endif |
|---|
| 1838 | !JJS 1/3/06 from Steve and Chunglin |
|---|
| 1839 | wgacs(i,j)=r2ice*r2iceg*r9r*dd(i,j) |
|---|
| 1840 | ! WGACS(I,J)=0. ! 6/15/02 tao's |
|---|
| 1841 | y1(i,j)=1./zg(i,j)**bg3 |
|---|
| 1842 | |
|---|
| 1843 | if(ihail .eq. 1) then |
|---|
| 1844 | dgacw(i,j)=r2ice*max(r14r*qc(i,j)*y1(i,j), 0.0) |
|---|
| 1845 | else |
|---|
| 1846 | dgacw(i,j)=r2ice*max(r14f*qc(i,j)*y1(i,j), 0.0) |
|---|
| 1847 | endif |
|---|
| 1848 | |
|---|
| 1849 | qgacw(i,j)=dgacw(i,j) |
|---|
| 1850 | y1(i,j)=abs(vg(i,j)-vr(i,j)) |
|---|
| 1851 | y2(i,j)=zr(i,j)*zg(i,j) |
|---|
| 1852 | y3(i,j)=5./y2(i,j) |
|---|
| 1853 | y4(i,j)=.08*y3(i,j)*y3(i,j) |
|---|
| 1854 | y5(i,j)=.05*y3(i,j)*y4(i,j) |
|---|
| 1855 | dd(i,j)=r16r*y1(i,j)*(y3(i,j)/zr(i,j)**5+y4(i,j)/zr(i,j)**3 & |
|---|
| 1856 | +y5(i,j)/zr(i,j)) |
|---|
| 1857 | dgacr(i,j)=r2ice*max(dd(i,j), 0.0) |
|---|
| 1858 | qgacr(i,j)=dgacr(i,j) |
|---|
| 1859 | |
|---|
| 1860 | if (tair(i,j).ge.t0) then |
|---|
| 1861 | dgacs(i,j)=0.0 |
|---|
| 1862 | wgacs(i,j)=0.0 |
|---|
| 1863 | dgacw(i,j)=0.0 |
|---|
| 1864 | dgacr(i,j)=0.0 |
|---|
| 1865 | else |
|---|
| 1866 | pgacs(i,j)=0.0 |
|---|
| 1867 | qgacw(i,j)=0.0 |
|---|
| 1868 | qgacr(i,j)=0.0 |
|---|
| 1869 | endif |
|---|
| 1870 | |
|---|
| 1871 | !*******PGDRY : DGACW+DGACI+DGACR+DGACS ****** |
|---|
| 1872 | !* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH) **15** |
|---|
| 1873 | !* 17 * PGWET : WET GROWTH OF QG **17** |
|---|
| 1874 | |
|---|
| 1875 | dgaci(i,j)=0.0 |
|---|
| 1876 | wgaci(i,j)=0.0 |
|---|
| 1877 | pgwet(i,j)=0.0 |
|---|
| 1878 | |
|---|
| 1879 | if (tair(i,j).lt.t0) then |
|---|
| 1880 | y1(i,j)=qi(i,j)/zg(i,j)**bg3 |
|---|
| 1881 | if (ihail.eq.1) then |
|---|
| 1882 | dgaci(i,j)=r2ice*r15r*y1(i,j) |
|---|
| 1883 | wgaci(i,j)=r2ice*r15ar*y1(i,j) |
|---|
| 1884 | ! WGACI(I,J)=0. ! 6/15/02 tao's |
|---|
| 1885 | else |
|---|
| 1886 | |
|---|
| 1887 | !JJS DGACI(I,J)=r2ice*R15F*Y1(I,J) |
|---|
| 1888 | dgaci(i,j)=0. |
|---|
| 1889 | wgaci(i,j)=r2ice*r15af*y1(i,j) |
|---|
| 1890 | ! WGACI(I,J)=0. ! 6/15/02 tao's |
|---|
| 1891 | endif |
|---|
| 1892 | ! |
|---|
| 1893 | if (tairc(i,j).ge.-50.) then |
|---|
| 1894 | if (alf+rn17c*tairc(i,j) .eq. 0.) then |
|---|
| 1895 | write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j) |
|---|
| 1896 | endif |
|---|
| 1897 | y1(i,j)=1./(alf+rn17c*tairc(i,j)) |
|---|
| 1898 | if (ihail.eq.1) then |
|---|
| 1899 | y3(i,j)=.78/zg(i,j)**2+r17aq*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 1900 | else |
|---|
| 1901 | y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 1902 | endif |
|---|
| 1903 | y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) & |
|---|
| 1904 | -tca(i,j)*tairc(i,j) |
|---|
| 1905 | dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) & |
|---|
| 1906 | +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j))) |
|---|
| 1907 | pgwet(i,j)=r2ice*max(dd(i,j), 0.0) |
|---|
| 1908 | endif |
|---|
| 1909 | endif |
|---|
| 1910 | !JJS 125 CONTINUE |
|---|
| 1911 | |
|---|
| 1912 | !******** HANDLING THE NEGATIVE CLOUD WATER (QC) ****************** |
|---|
| 1913 | !******** HANDLING THE NEGATIVE CLOUD ICE (QI) ****************** |
|---|
| 1914 | |
|---|
| 1915 | !JJS DO 150 J=3,JLES |
|---|
| 1916 | !JJS DO 150 I=3,ILES |
|---|
| 1917 | |
|---|
| 1918 | y1(i,j)=qc(i,j)/d2t |
|---|
| 1919 | psacw(i,j)=min(y1(i,j), psacw(i,j)) |
|---|
| 1920 | praut(i,j)=min(y1(i,j), praut(i,j)) |
|---|
| 1921 | pracw(i,j)=min(y1(i,j), pracw(i,j)) |
|---|
| 1922 | psfw(i,j)= min(y1(i,j), psfw(i,j)) |
|---|
| 1923 | dgacw(i,j)=min(y1(i,j), dgacw(i,j)) |
|---|
| 1924 | qsacw(i,j)=min(y1(i,j), qsacw(i,j)) |
|---|
| 1925 | qgacw(i,j)=min(y1(i,j), qgacw(i,j)) |
|---|
| 1926 | |
|---|
| 1927 | y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) & |
|---|
| 1928 | +dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t |
|---|
| 1929 | qc(i,j)=qc(i,j)-y1(i,j) |
|---|
| 1930 | |
|---|
| 1931 | if (qc(i,j) .lt. 0.0) then |
|---|
| 1932 | a1=1. |
|---|
| 1933 | if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1. |
|---|
| 1934 | psacw(i,j)=psacw(i,j)*a1 |
|---|
| 1935 | praut(i,j)=praut(i,j)*a1 |
|---|
| 1936 | pracw(i,j)=pracw(i,j)*a1 |
|---|
| 1937 | psfw(i,j)=psfw(i,j)*a1 |
|---|
| 1938 | dgacw(i,j)=dgacw(i,j)*a1 |
|---|
| 1939 | qsacw(i,j)=qsacw(i,j)*a1 |
|---|
| 1940 | qgacw(i,j)=qgacw(i,j)*a1 |
|---|
| 1941 | qc(i,j)=0.0 |
|---|
| 1942 | endif |
|---|
| 1943 | !c |
|---|
| 1944 | ! |
|---|
| 1945 | !******** SHED PROCESS (WGACR=PGWET-DGACW-WGACI-WGACS) |
|---|
| 1946 | !c |
|---|
| 1947 | wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j) |
|---|
| 1948 | y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j) |
|---|
| 1949 | if (pgwet(i,j).ge.y2(i,j)) then |
|---|
| 1950 | wgacr(i,j)=0.0 |
|---|
| 1951 | wgaci(i,j)=0.0 |
|---|
| 1952 | wgacs(i,j)=0.0 |
|---|
| 1953 | else |
|---|
| 1954 | dgacr(i,j)=0.0 |
|---|
| 1955 | dgaci(i,j)=0.0 |
|---|
| 1956 | dgacs(i,j)=0.0 |
|---|
| 1957 | endif |
|---|
| 1958 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1959 | !c |
|---|
| 1960 | y1(i,j)=qi(i,j)/d2t |
|---|
| 1961 | psaut(i,j)=min(y1(i,j), psaut(i,j)) |
|---|
| 1962 | psaci(i,j)=min(y1(i,j), psaci(i,j)) |
|---|
| 1963 | praci(i,j)=min(y1(i,j), praci(i,j)) |
|---|
| 1964 | psfi(i,j)= min(y1(i,j), psfi(i,j)) |
|---|
| 1965 | dgaci(i,j)=min(y1(i,j), dgaci(i,j)) |
|---|
| 1966 | wgaci(i,j)=min(y1(i,j), wgaci(i,j)) |
|---|
| 1967 | ! |
|---|
| 1968 | y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) & |
|---|
| 1969 | +dgaci(i,j)+wgaci(i,j))*d2t |
|---|
| 1970 | qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t |
|---|
| 1971 | |
|---|
| 1972 | if (qi(i,j).lt.0.0) then |
|---|
| 1973 | a2=1. |
|---|
| 1974 | if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1. |
|---|
| 1975 | psaut(i,j)=psaut(i,j)*a2 |
|---|
| 1976 | psaci(i,j)=psaci(i,j)*a2 |
|---|
| 1977 | praci(i,j)=praci(i,j)*a2 |
|---|
| 1978 | psfi(i,j)=psfi(i,j)*a2 |
|---|
| 1979 | dgaci(i,j)=dgaci(i,j)*a2 |
|---|
| 1980 | wgaci(i,j)=wgaci(i,j)*a2 |
|---|
| 1981 | qi(i,j)=0.0 |
|---|
| 1982 | endif |
|---|
| 1983 | ! |
|---|
| 1984 | dlt3(i,j)=0.0 |
|---|
| 1985 | dlt2(i,j)=0.0 |
|---|
| 1986 | ! |
|---|
| 1987 | |
|---|
| 1988 | ! DLT4(I,J)=1.0 |
|---|
| 1989 | ! if(qc(i,j) .gt. 5.e-4) dlt4(i,j)=0.0 |
|---|
| 1990 | ! if(qs(i,j) .le. 1.e-4) dlt4(i,j)=1.0 |
|---|
| 1991 | ! |
|---|
| 1992 | ! IF (TAIR(I,J).ge.T0) THEN |
|---|
| 1993 | ! DLT4(I,J)=0.0 |
|---|
| 1994 | ! ENDIF |
|---|
| 1995 | |
|---|
| 1996 | if (tair(i,j).lt.t0) then |
|---|
| 1997 | if (qr(i,j).lt.1.e-4) then |
|---|
| 1998 | dlt3(i,j)=1.0 |
|---|
| 1999 | dlt2(i,j)=1.0 |
|---|
| 2000 | endif |
|---|
| 2001 | if (qs(i,j).ge.1.e-4) then |
|---|
| 2002 | dlt2(i,j)=0.0 |
|---|
| 2003 | endif |
|---|
| 2004 | endif |
|---|
| 2005 | |
|---|
| 2006 | if (ice2 .eq. 1) then |
|---|
| 2007 | dlt3(i,j)=1.0 |
|---|
| 2008 | dlt2(i,j)=1.0 |
|---|
| 2009 | endif |
|---|
| 2010 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2011 | pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t |
|---|
| 2012 | ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) & |
|---|
| 2013 | +psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t |
|---|
| 2014 | ! PS(I,J)=(PSAUT(I,J)+PSACI(I,J)+dlt4(i,j)*PSACW(I,J) |
|---|
| 2015 | ! 1 +PSFW(I,J)+PSFI(I,J)+DLT3(I,J)*PRACI(I,J))*D2T |
|---|
| 2016 | pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) & |
|---|
| 2017 | +dgacw(i,j))*d2t |
|---|
| 2018 | ! PG(I,J)=((1.-DLT3(I,J))*PRACI(I,J)+DGACI(I,J)+WGACI(I,J) |
|---|
| 2019 | ! 1 +DGACW(I,J)+(1.-dlt4(i,j))*PSACW(I,J))*D2T |
|---|
| 2020 | |
|---|
| 2021 | !JJS 150 CONTINUE |
|---|
| 2022 | |
|---|
| 2023 | !* 7 * PRACS : ACCRETION OF QS BY QR ***7** |
|---|
| 2024 | !* 8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT) ***8** |
|---|
| 2025 | |
|---|
| 2026 | !JJS DO 175 J=3,JLES |
|---|
| 2027 | !JJS DO 175 I=3,ILES |
|---|
| 2028 | |
|---|
| 2029 | y1(i,j)=abs(vr(i,j)-vs(i,j)) |
|---|
| 2030 | y2(i,j)=zr(i,j)*zs(i,j) |
|---|
| 2031 | y3(i,j)=5./y2(i,j) |
|---|
| 2032 | y4(i,j)=.08*y3(i,j)*y3(i,j) |
|---|
| 2033 | y5(i,j)=.05*y3(i,j)*y4(i,j) |
|---|
| 2034 | pracs(i,j)=r2ice*r2iceg*r7r*y1(i,j)*(y3(i,j)/zs(i,j)**5 & |
|---|
| 2035 | +y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j)) |
|---|
| 2036 | psacr(i,j)=r2iceg*r8r*y1(i,j)*(y3(i,j)/zr(i,j)**5 & |
|---|
| 2037 | +y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j)) |
|---|
| 2038 | qsacr(i,j)=psacr(i,j) |
|---|
| 2039 | |
|---|
| 2040 | if (tair(i,j).ge.t0) then |
|---|
| 2041 | pracs(i,j)=0.0 |
|---|
| 2042 | psacr(i,j)=0.0 |
|---|
| 2043 | else |
|---|
| 2044 | qsacr(i,j)=0.0 |
|---|
| 2045 | endif |
|---|
| 2046 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2047 | !* 2 * PGAUT : AUTOCONVERSION OF QS TO QG ***2** |
|---|
| 2048 | !* 18 * PGFR : FREEZING OF QR TO QG **18** |
|---|
| 2049 | |
|---|
| 2050 | pgaut(i,j)=0.0 |
|---|
| 2051 | pgfr(i,j)=0.0 |
|---|
| 2052 | |
|---|
| 2053 | if (tair(i,j) .lt. t0) then |
|---|
| 2054 | ! Y1(I,J)=EXP(.09*TAIRC(I,J)) |
|---|
| 2055 | ! PGAUT(I,J)=r2iceg*max(RN2*Y1(I,J)*(QS(I,J)-BND2), 0.0) |
|---|
| 2056 | ! IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2*Y1(I,J)*(QS(I,J)-BND2),0.0) |
|---|
| 2057 | y2(i,j)=exp(rn18a*(t0-tair(i,j))) |
|---|
| 2058 | !JJS PGFR(I,J)=r2ice*max(R18R*(Y2(I,J)-1.)/ZR(I,J)**7., 0.0) |
|---|
| 2059 | ! pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* & |
|---|
| 2060 | ! (zr(i,j)**(-7.)), 0.0) |
|---|
| 2061 | ! modify to prevent underflow on some computers (JD) |
|---|
| 2062 | temp = 1./zr(i,j) |
|---|
| 2063 | temp = temp*temp*temp*temp*temp*temp*temp |
|---|
| 2064 | pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* & |
|---|
| 2065 | temp, 0.0) |
|---|
| 2066 | endif |
|---|
| 2067 | |
|---|
| 2068 | !JJS 175 CONTINUE |
|---|
| 2069 | |
|---|
| 2070 | !******** HANDLING THE NEGATIVE RAIN WATER (QR) ******************* |
|---|
| 2071 | !******** HANDLING THE NEGATIVE SNOW (QS) ******************* |
|---|
| 2072 | |
|---|
| 2073 | !JJS DO 200 J=3,JLES |
|---|
| 2074 | !JJS DO 200 I=3,ILES |
|---|
| 2075 | |
|---|
| 2076 | y1(i,j)=qr(i,j)/d2t |
|---|
| 2077 | y2(i,j)=-qg(i,j)/d2t |
|---|
| 2078 | piacr(i,j)=min(y1(i,j), piacr(i,j)) |
|---|
| 2079 | dgacr(i,j)=min(y1(i,j), dgacr(i,j)) |
|---|
| 2080 | wgacr(i,j)=min(y1(i,j), wgacr(i,j)) |
|---|
| 2081 | wgacr(i,j)=max(y2(i,j), wgacr(i,j)) |
|---|
| 2082 | psacr(i,j)=min(y1(i,j), psacr(i,j)) |
|---|
| 2083 | pgfr(i,j)= min(y1(i,j), pgfr(i,j)) |
|---|
| 2084 | del=0. |
|---|
| 2085 | if(wgacr(i,j) .lt. 0.) del=1. |
|---|
| 2086 | y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) & |
|---|
| 2087 | +psacr(i,j)+pgfr(i,j))*d2t |
|---|
| 2088 | qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-del*wgacr(i,j)*d2t |
|---|
| 2089 | if (qr(i,j) .lt. 0.0) then |
|---|
| 2090 | a1=1. |
|---|
| 2091 | if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1. |
|---|
| 2092 | piacr(i,j)=piacr(i,j)*a1 |
|---|
| 2093 | dgacr(i,j)=dgacr(i,j)*a1 |
|---|
| 2094 | if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1 |
|---|
| 2095 | pgfr(i,j)=pgfr(i,j)*a1 |
|---|
| 2096 | psacr(i,j)=psacr(i,j)*a1 |
|---|
| 2097 | qr(i,j)=0.0 |
|---|
| 2098 | endif |
|---|
| 2099 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2100 | prn(i,j)=d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j) & |
|---|
| 2101 | +wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j)) |
|---|
| 2102 | ps(i,j)=ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j) & |
|---|
| 2103 | +dlt2(i,j)*psacr(i,j)) |
|---|
| 2104 | pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j) |
|---|
| 2105 | y1(i,j)=qs(i,j)/d2t |
|---|
| 2106 | pgacs(i,j)=min(y1(i,j), pgacs(i,j)) |
|---|
| 2107 | dgacs(i,j)=min(y1(i,j), dgacs(i,j)) |
|---|
| 2108 | wgacs(i,j)=min(y1(i,j), wgacs(i,j)) |
|---|
| 2109 | pgaut(i,j)=min(y1(i,j), pgaut(i,j)) |
|---|
| 2110 | pracs(i,j)=min(y1(i,j), pracs(i,j)) |
|---|
| 2111 | psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) & |
|---|
| 2112 | +pgaut(i,j)+pracs(i,j)) |
|---|
| 2113 | qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j) |
|---|
| 2114 | |
|---|
| 2115 | if (qs(i,j).lt.0.0) then |
|---|
| 2116 | a2=1. |
|---|
| 2117 | if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1. |
|---|
| 2118 | pgacs(i,j)=pgacs(i,j)*a2 |
|---|
| 2119 | dgacs(i,j)=dgacs(i,j)*a2 |
|---|
| 2120 | wgacs(i,j)=wgacs(i,j)*a2 |
|---|
| 2121 | pgaut(i,j)=pgaut(i,j)*a2 |
|---|
| 2122 | pracs(i,j)=pracs(i,j)*a2 |
|---|
| 2123 | psn(i,j)=psn(i,j)*a2 |
|---|
| 2124 | qs(i,j)=0.0 |
|---|
| 2125 | endif |
|---|
| 2126 | ! |
|---|
| 2127 | !C PSN(I,J)=D2T*(PGACS(I,J)+DGACS(I,J)+WGACS(I,J) |
|---|
| 2128 | !c +PGAUT(I,J)+PRACS(I,J)) |
|---|
| 2129 | y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) & |
|---|
| 2130 | +dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j)) |
|---|
| 2131 | pt(i,j)=pt(i,j)+afcp*y2(i,j) |
|---|
| 2132 | qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j) |
|---|
| 2133 | |
|---|
| 2134 | !JJS 200 CONTINUE |
|---|
| 2135 | |
|---|
| 2136 | !* 11 * PSMLT : MELTING OF QS **11** |
|---|
| 2137 | !* 19 * PGMLT : MELTING OF QG TO QR **19** |
|---|
| 2138 | |
|---|
| 2139 | !JJS DO 225 J=3,JLES |
|---|
| 2140 | !JJS DO 225 I=3,ILES |
|---|
| 2141 | |
|---|
| 2142 | psmlt(i,j)=0.0 |
|---|
| 2143 | pgmlt(i,j)=0.0 |
|---|
| 2144 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2145 | |
|---|
| 2146 | if (tair(i,j).ge.t0) then |
|---|
| 2147 | tairc(i,j)=tair(i,j)-t0 |
|---|
| 2148 | y1(i,j)=tca(i,j)*tairc(i,j)-alvr*dwv(i,j) & |
|---|
| 2149 | *(rp0-(qv(i,j)+qb0)) |
|---|
| 2150 | y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5 |
|---|
| 2151 | dd(i,j)=r11rt*y1(i,j)*y2(i,j)+r11at*tairc(i,j) & |
|---|
| 2152 | *(qsacw(i,j)+qsacr(i,j)) |
|---|
| 2153 | psmlt(i,j)=r2iceg*max(0.0, min(dd(i,j), qs(i,j))) |
|---|
| 2154 | |
|---|
| 2155 | if(ihail.eq.1) then |
|---|
| 2156 | y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 2157 | else |
|---|
| 2158 | y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 2159 | endif |
|---|
| 2160 | |
|---|
| 2161 | dd1(i,j)=r19rt*y1(i,j)*y3(i,j)+r19bt*tairc(i,j) & |
|---|
| 2162 | *(qgacw(i,j)+qgacr(i,j)) |
|---|
| 2163 | pgmlt(i,j)=r2ice*max(0.0, min(dd1(i,j), qg(i,j))) |
|---|
| 2164 | pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j)) |
|---|
| 2165 | qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j) |
|---|
| 2166 | qs(i,j)=qs(i,j)-psmlt(i,j) |
|---|
| 2167 | qg(i,j)=qg(i,j)-pgmlt(i,j) |
|---|
| 2168 | endif |
|---|
| 2169 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2170 | !* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00) **24** |
|---|
| 2171 | !* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00) **25** |
|---|
| 2172 | !* 26 * PIMLT : MELTING OF QI TO QC (T >= T0) **26** |
|---|
| 2173 | |
|---|
| 2174 | if (qc(i,j).le.cmin1) qc(i,j)=0.0 |
|---|
| 2175 | if (qi(i,j).le.cmin1) qi(i,j)=0.0 |
|---|
| 2176 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2177 | |
|---|
| 2178 | if(tair(i,j).le.t00) then |
|---|
| 2179 | pihom(i,j)=qc(i,j) |
|---|
| 2180 | else |
|---|
| 2181 | pihom(i,j)=0.0 |
|---|
| 2182 | endif |
|---|
| 2183 | if(tair(i,j).ge.t0) then |
|---|
| 2184 | pimlt(i,j)=qi(i,j) |
|---|
| 2185 | else |
|---|
| 2186 | pimlt(i,j)=0.0 |
|---|
| 2187 | endif |
|---|
| 2188 | pidw(i,j)=0.0 |
|---|
| 2189 | |
|---|
| 2190 | if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then |
|---|
| 2191 | tairc(i,j)=tair(i,j)-t0 |
|---|
| 2192 | y1(i,j)=max( min(tairc(i,j), -1.), -31.) |
|---|
| 2193 | it(i,j)=int(abs(y1(i,j))) |
|---|
| 2194 | y2(i,j)=aa1(it(i,j)) |
|---|
| 2195 | y3(i,j)=aa2(it(i,j)) |
|---|
| 2196 | y4(i,j)=exp(abs(beta*tairc(i,j))) |
|---|
| 2197 | y5(i,j)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j) |
|---|
| 2198 | pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j)) |
|---|
| 2199 | endif |
|---|
| 2200 | |
|---|
| 2201 | y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j) |
|---|
| 2202 | pt(i,j)=pt(i,j)+afcp*y1(i,j)+ascp*(pidep(i,j))*d2t |
|---|
| 2203 | qv(i,j)=qv(i,j)-(pidep(i,j))*d2t |
|---|
| 2204 | qc(i,j)=qc(i,j)-y1(i,j) |
|---|
| 2205 | qi(i,j)=qi(i,j)+y1(i,j) |
|---|
| 2206 | |
|---|
| 2207 | !* 31 * PINT : INITIATION OF QI **31** |
|---|
| 2208 | !* 32 * PIDEP : DEPOSITION OF QI **32** |
|---|
| 2209 | ! |
|---|
| 2210 | ! CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR |
|---|
| 2211 | ! THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER |
|---|
| 2212 | ! CONCENTRATION EXCEEDS THAT ALREADY PRESENT. |
|---|
| 2213 | !* 31 * pint : initiation of qi **31** |
|---|
| 2214 | !* 32 * pidep : deposition of qi **32** |
|---|
| 2215 | pint(i,j)=0.0 |
|---|
| 2216 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2217 | if ( itaobraun.eq.1 ) then |
|---|
| 2218 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2219 | if (tair(i,j) .lt. t0) then |
|---|
| 2220 | ! if (qi(i,j) .le. cmin) qi(i,j)=0. |
|---|
| 2221 | if (qi(i,j) .le. cmin2) qi(i,j)=0. |
|---|
| 2222 | tairc(i,j)=tair(i,j)-t0 |
|---|
| 2223 | rtair(i,j)=1./(tair(i,j)-c76) |
|---|
| 2224 | y2(i,j)=exp(c218-c580*rtair(i,j)) |
|---|
| 2225 | qsi(i,j)=rp0*y2(i,j) |
|---|
| 2226 | esi(i,j)=c610*y2(i,j) |
|---|
| 2227 | ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. |
|---|
| 2228 | ami50=3.76e-8 |
|---|
| 2229 | |
|---|
| 2230 | !ccif ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6 |
|---|
| 2231 | !ccif ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8 |
|---|
| 2232 | |
|---|
| 2233 | y1(i,j)=1./tair(i,j) |
|---|
| 2234 | |
|---|
| 2235 | !cc insert a restriction on ice collection that ice collection |
|---|
| 2236 | !cc should be stopped at -30 c (with cn0=1.e-6, beta=-.46) |
|---|
| 2237 | |
|---|
| 2238 | tairccri=tairc(i,j) ! in degree c |
|---|
| 2239 | if(tairccri.le.-30.) tairccri=-30. |
|---|
| 2240 | |
|---|
| 2241 | ! y2(i,j)=exp(betah*tairc(i,j)) |
|---|
| 2242 | y2(i,j)=exp(betah*tairccri) |
|---|
| 2243 | y3(i,j)=sqrt(qi(i,j)) |
|---|
| 2244 | dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) & |
|---|
| 2245 | /esi(i,j) |
|---|
| 2246 | pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0) |
|---|
| 2247 | |
|---|
| 2248 | r_nci=min(cn0*exp(beta*tairc(i,j)),1.) |
|---|
| 2249 | ! r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.) |
|---|
| 2250 | |
|---|
| 2251 | dd(i,j)=max(1.e-9*r_nci/r00-qi(i,j)*1.e-9/ami50, 0.) |
|---|
| 2252 | dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0) |
|---|
| 2253 | rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) |
|---|
| 2254 | dep(i,j)=dm(i,j)/(1.+rsub1(i,j)) |
|---|
| 2255 | pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.) |
|---|
| 2256 | |
|---|
| 2257 | ! pint(i,j)=min(pint(i,j), dep(i,j)) |
|---|
| 2258 | pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j)) |
|---|
| 2259 | |
|---|
| 2260 | ! if (pint(i,j) .le. cmin) pint(i,j)=0. |
|---|
| 2261 | if (pint(i,j) .le. cmin2) pint(i,j)=0. |
|---|
| 2262 | pt(i,j)=pt(i,j)+ascp*pint(i,j) |
|---|
| 2263 | qv(i,j)=qv(i,j)-pint(i,j) |
|---|
| 2264 | qi(i,j)=qi(i,j)+pint(i,j) |
|---|
| 2265 | endif |
|---|
| 2266 | endif ! if ( itaobraun.eq.1 ) |
|---|
| 2267 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2268 | ! |
|---|
| 2269 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2270 | if ( itaobraun.eq.0 ) then |
|---|
| 2271 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2272 | if (tair(i,j) .lt. t0) then |
|---|
| 2273 | if (qi(i,j) .le. cmin1) qi(i,j)=0. |
|---|
| 2274 | tairc(i,j)=tair(i,j)-t0 |
|---|
| 2275 | dd(i,j)=r31r*exp(beta*tairc(i,j)) |
|---|
| 2276 | rtair(i,j)=1./(tair(i,j)-c76) |
|---|
| 2277 | y2(i,j)=exp(c218-c580*rtair(i,j)) |
|---|
| 2278 | qsi(i,j)=rp0*y2(i,j) |
|---|
| 2279 | esi(i,j)=c610*y2(i,j) |
|---|
| 2280 | ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. |
|---|
| 2281 | dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.) |
|---|
| 2282 | rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) |
|---|
| 2283 | dep(i,j)=dm(i,j)/(1.+rsub1(i,j)) |
|---|
| 2284 | pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.) |
|---|
| 2285 | y1(i,j)=1./tair(i,j) |
|---|
| 2286 | y2(i,j)=exp(betah*tairc(i,j)) |
|---|
| 2287 | y3(i,j)=sqrt(qi(i,j)) |
|---|
| 2288 | dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b) & |
|---|
| 2289 | +rn10c*tair(i,j)/esi(i,j) |
|---|
| 2290 | pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.) |
|---|
| 2291 | pint(i,j)=pint(i,j)+pidep(i,j) |
|---|
| 2292 | pint(i,j)=min(pint(i,j),dep(i,j)) |
|---|
| 2293 | !c if (pint(i,j) .le. cmin2) pint(i,j)=0. |
|---|
| 2294 | pt(i,j)=pt(i,j)+ascp*pint(i,j) |
|---|
| 2295 | qv(i,j)=qv(i,j)-pint(i,j) |
|---|
| 2296 | qi(i,j)=qi(i,j)+pint(i,j) |
|---|
| 2297 | endif |
|---|
| 2298 | endif ! if ( itaobraun.eq.0 ) |
|---|
| 2299 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2300 | |
|---|
| 2301 | !JJS 225 CONTINUE |
|---|
| 2302 | |
|---|
| 2303 | !***** TAO ET AL (1989) SATURATION TECHNIQUE *********************** |
|---|
| 2304 | |
|---|
| 2305 | if (new_ice_sat .eq. 0) then |
|---|
| 2306 | |
|---|
| 2307 | !JJS DO 250 J=3,JLES |
|---|
| 2308 | !JJS DO 250 I=3,ILES |
|---|
| 2309 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2310 | cnd(i,j)=rt0*(tair(i,j)-t00) |
|---|
| 2311 | dep(i,j)=rt0*(t0-tair(i,j)) |
|---|
| 2312 | y1(i,j)=1./(tair(i,j)-c358) |
|---|
| 2313 | y2(i,j)=1./(tair(i,j)-c76) |
|---|
| 2314 | qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) |
|---|
| 2315 | qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) |
|---|
| 2316 | dd(i,j)=cp409*y1(i,j)*y1(i,j) |
|---|
| 2317 | dd1(i,j)=cp580*y2(i,j)*y2(i,j) |
|---|
| 2318 | if (qc(i,j).le.cmin) qc(i,j)=cmin |
|---|
| 2319 | if (qi(i,j).le.cmin) qi(i,j)=cmin |
|---|
| 2320 | if (tair(i,j).ge.t0) then |
|---|
| 2321 | dep(i,j)=0.0 |
|---|
| 2322 | cnd(i,j)=1. |
|---|
| 2323 | qi(i,j)=0.0 |
|---|
| 2324 | endif |
|---|
| 2325 | |
|---|
| 2326 | if (tair(i,j).lt.t00) then |
|---|
| 2327 | cnd(i,j)=0.0 |
|---|
| 2328 | dep(i,j)=1. |
|---|
| 2329 | qc(i,j)=0.0 |
|---|
| 2330 | endif |
|---|
| 2331 | |
|---|
| 2332 | y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j) |
|---|
| 2333 | ! if (qc(i,j) .ge. cmin .or. qi(i,j) .ge. cmin) then |
|---|
| 2334 | y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j)) |
|---|
| 2335 | y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j)) |
|---|
| 2336 | y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j) |
|---|
| 2337 | qvs(i,j)=y1(i,j)+y2(i,j) |
|---|
| 2338 | rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j)) |
|---|
| 2339 | cnd(i,j)=cnd(i,j)*rsub1(i,j) |
|---|
| 2340 | dep(i,j)=dep(i,j)*rsub1(i,j) |
|---|
| 2341 | if (qc(i,j).le.cmin) qc(i,j)=0. |
|---|
| 2342 | if (qi(i,j).le.cmin) qi(i,j)=0. |
|---|
| 2343 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2344 | !c ****** condensation or evaporation of qc ****** |
|---|
| 2345 | |
|---|
| 2346 | cnd(i,j)=max(-qc(i,j),cnd(i,j)) |
|---|
| 2347 | |
|---|
| 2348 | !c ****** deposition or sublimation of qi ****** |
|---|
| 2349 | |
|---|
| 2350 | dep(i,j)=max(-qi(i,j),dep(i,j)) |
|---|
| 2351 | |
|---|
| 2352 | pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j) |
|---|
| 2353 | qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j) |
|---|
| 2354 | qc(i,j)=qc(i,j)+cnd(i,j) |
|---|
| 2355 | qi(i,j)=qi(i,j)+dep(i,j) |
|---|
| 2356 | !JJS 250 continue |
|---|
| 2357 | endif |
|---|
| 2358 | |
|---|
| 2359 | if (new_ice_sat .eq. 1) then |
|---|
| 2360 | |
|---|
| 2361 | !JJS DO J=3,JLES |
|---|
| 2362 | !JJS DO I=3,ILES |
|---|
| 2363 | |
|---|
| 2364 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2365 | cnd(i,j)=rt0*(tair(i,j)-t00) |
|---|
| 2366 | dep(i,j)=rt0*(t0-tair(i,j)) |
|---|
| 2367 | y1(i,j)=1./(tair(i,j)-c358) |
|---|
| 2368 | y2(i,j)=1./(tair(i,j)-c76) |
|---|
| 2369 | qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) |
|---|
| 2370 | qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) |
|---|
| 2371 | dd(i,j)=cp409*y1(i,j)*y1(i,j) |
|---|
| 2372 | dd1(i,j)=cp580*y2(i,j)*y2(i,j) |
|---|
| 2373 | y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j) |
|---|
| 2374 | y1(i,j)=rt0*(tair(i,j)-t00)*qsw(i,j) |
|---|
| 2375 | y2(i,j)=rt0*(t0-tair(i,j))*qsi(i,j) |
|---|
| 2376 | ! IF (QC(I,J).LE.CMIN) QC(I,J)=CMIN |
|---|
| 2377 | ! IF (QI(I,J).LE.CMIN) QI(I,J)=CMIN |
|---|
| 2378 | |
|---|
| 2379 | if (tair(i,j).ge.t0) then |
|---|
| 2380 | ! QI(I,J)=0.0 |
|---|
| 2381 | dep(i,j)=0.0 |
|---|
| 2382 | cnd(i,j)=1. |
|---|
| 2383 | y2(i,j)=0. |
|---|
| 2384 | y1(i,j)=qsw(i,j) |
|---|
| 2385 | endif |
|---|
| 2386 | if (tair(i,j).lt.t00) then |
|---|
| 2387 | cnd(i,j)=0.0 |
|---|
| 2388 | dep(i,j)=1. |
|---|
| 2389 | y2(i,j)=qsi(i,j) |
|---|
| 2390 | y1(i,j)=0. |
|---|
| 2391 | ! QC(I,J)=0.0 |
|---|
| 2392 | endif |
|---|
| 2393 | |
|---|
| 2394 | ! Y1(I,J)=QC(I,J)*QSW(I,J)/(QC(I,J)+QI(I,J)) |
|---|
| 2395 | ! Y2(I,J)=QI(I,J)*QSI(I,J)/(QC(I,J)+QI(I,J)) |
|---|
| 2396 | |
|---|
| 2397 | y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j) |
|---|
| 2398 | qvs(i,j)=y1(i,j)+y2(i,j) |
|---|
| 2399 | rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j)) |
|---|
| 2400 | cnd(i,j)=cnd(i,j)*rsub1(i,j) |
|---|
| 2401 | dep(i,j)=dep(i,j)*rsub1(i,j) |
|---|
| 2402 | ! IF (QC(I,J).LE.CMIN) QC(I,J)=0. |
|---|
| 2403 | ! IF (QI(I,J).LE.CMIN) QI(I,J)=0. |
|---|
| 2404 | |
|---|
| 2405 | !C ****** CONDENSATION OR EVAPORATION OF QC ****** |
|---|
| 2406 | |
|---|
| 2407 | cnd(i,j)=max(-qc(i,j),cnd(i,j)) |
|---|
| 2408 | |
|---|
| 2409 | !C ****** DEPOSITION OR SUBLIMATION OF QI ****** |
|---|
| 2410 | |
|---|
| 2411 | dep(i,j)=max(-qi(i,j),dep(i,j)) |
|---|
| 2412 | |
|---|
| 2413 | pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j) |
|---|
| 2414 | qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j) |
|---|
| 2415 | qc(i,j)=qc(i,j)+cnd(i,j) |
|---|
| 2416 | qi(i,j)=qi(i,j)+dep(i,j) |
|---|
| 2417 | !JJS ENDDO |
|---|
| 2418 | !JJS ENDDO |
|---|
| 2419 | endif |
|---|
| 2420 | |
|---|
| 2421 | !c |
|---|
| 2422 | ! |
|---|
| 2423 | if (new_ice_sat .eq. 2) then |
|---|
| 2424 | !JJS do j=3,jles |
|---|
| 2425 | !JJS do i=3,iles |
|---|
| 2426 | dep(i,j)=0.0 |
|---|
| 2427 | cnd(i,j)=0.0 |
|---|
| 2428 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2429 | if (tair(i,j) .ge. 253.16) then |
|---|
| 2430 | y1(i,j)=1./(tair(i,j)-c358) |
|---|
| 2431 | qsw(i,j)=rp0*exp(c172-c409*y1(i,j)) |
|---|
| 2432 | dd(i,j)=cp409*y1(i,j)*y1(i,j) |
|---|
| 2433 | dm(i,j)=qv(i,j)+qb0-qsw(i,j) |
|---|
| 2434 | cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j)) |
|---|
| 2435 | !c ****** condensation or evaporation of qc ****** |
|---|
| 2436 | cnd(i,j)=max(-qc(i,j), cnd(i,j)) |
|---|
| 2437 | pt(i,j)=pt(i,j)+avcp*cnd(i,j) |
|---|
| 2438 | qv(i,j)=qv(i,j)-cnd(i,j) |
|---|
| 2439 | qc(i,j)=qc(i,j)+cnd(i,j) |
|---|
| 2440 | endif |
|---|
| 2441 | if (tair(i,j) .le. 258.16) then |
|---|
| 2442 | !c cnd(i,j)=0.0 |
|---|
| 2443 | y2(i,j)=1./(tair(i,j)-c76) |
|---|
| 2444 | qsi(i,j)=rp0*exp(c218-c580*y2(i,j)) |
|---|
| 2445 | dd1(i,j)=cp580*y2(i,j)*y2(i,j) |
|---|
| 2446 | dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(i,j)*qsi(i,j)) |
|---|
| 2447 | !c ****** deposition or sublimation of qi ****** |
|---|
| 2448 | dep(i,j)=max(-qi(i,j),dep(i,j)) |
|---|
| 2449 | pt(i,j)=pt(i,j)+ascp*dep(i,j) |
|---|
| 2450 | qv(i,j)=qv(i,j)-dep(i,j) |
|---|
| 2451 | qi(i,j)=qi(i,j)+dep(i,j) |
|---|
| 2452 | endif |
|---|
| 2453 | !JJS enddo |
|---|
| 2454 | !JJS enddo |
|---|
| 2455 | endif |
|---|
| 2456 | |
|---|
| 2457 | !c |
|---|
| 2458 | ! |
|---|
| 2459 | !* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS **10** |
|---|
| 2460 | !* 20 * PGSUB : SUBLIMATION OF QG **20** |
|---|
| 2461 | |
|---|
| 2462 | !JJS DO 275 J=3,JLES |
|---|
| 2463 | !JJS DO 275 I=3,ILES |
|---|
| 2464 | |
|---|
| 2465 | psdep(i,j)=0.0 |
|---|
| 2466 | pssub(i,j)=0.0 |
|---|
| 2467 | pgsub(i,j)=0.0 |
|---|
| 2468 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2469 | |
|---|
| 2470 | if(tair(i,j).lt.t0) then |
|---|
| 2471 | if(qs(i,j).lt.cmin1) qs(i,j)=0.0 |
|---|
| 2472 | if(qg(i,j).lt.cmin1) qg(i,j)=0.0 |
|---|
| 2473 | rtair(i,j)=1./(tair(i,j)-c76) |
|---|
| 2474 | qsi(i,j)=rp0*exp(c218-c580*rtair(i,j)) |
|---|
| 2475 | ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1. |
|---|
| 2476 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2477 | y1(i,j)=r10ar/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) & |
|---|
| 2478 | *qsi(i,j)) |
|---|
| 2479 | y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5 |
|---|
| 2480 | psdep(i,j)=r10t*ssi(i,j)*y2(i,j)/y1(i,j) |
|---|
| 2481 | pssub(i,j)=psdep(i,j) |
|---|
| 2482 | psdep(i,j)=r2iceg*max(psdep(i,j), 0.) |
|---|
| 2483 | pssub(i,j)=r2iceg*max(-qs(i,j), min(pssub(i,j), 0.)) |
|---|
| 2484 | |
|---|
| 2485 | if(ihail.eq.1) then |
|---|
| 2486 | y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 2487 | else |
|---|
| 2488 | y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5 |
|---|
| 2489 | endif |
|---|
| 2490 | |
|---|
| 2491 | pgsub(i,j)=r2ice*r20t*ssi(i,j)*y2(i,j)/y1(i,j) |
|---|
| 2492 | dm(i,j)=qv(i,j)+qb0-qsi(i,j) |
|---|
| 2493 | rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j) |
|---|
| 2494 | |
|---|
| 2495 | ! ******** DEPOSITION OR SUBLIMATION OF QS ********************** |
|---|
| 2496 | |
|---|
| 2497 | y1(i,j)=dm(i,j)/(1.+rsub1(i,j)) |
|---|
| 2498 | psdep(i,j)=r2iceg*min(psdep(i,j),max(y1(i,j),0.)) |
|---|
| 2499 | y2(i,j)=min(y1(i,j),0.) |
|---|
| 2500 | pssub(i,j)=r2iceg*max(pssub(i,j),y2(i,j)) |
|---|
| 2501 | |
|---|
| 2502 | ! ******** SUBLIMATION OF QG *********************************** |
|---|
| 2503 | |
|---|
| 2504 | dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.) |
|---|
| 2505 | pgsub(i,j)=r2ice*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.)) |
|---|
| 2506 | |
|---|
| 2507 | if(qc(i,j)+qi(i,j).gt.1.e-5) then |
|---|
| 2508 | dlt1(i,j)=1. |
|---|
| 2509 | else |
|---|
| 2510 | dlt1(i,j)=0. |
|---|
| 2511 | endif |
|---|
| 2512 | |
|---|
| 2513 | psdep(i,j)=dlt1(i,j)*psdep(i,j) |
|---|
| 2514 | pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j) |
|---|
| 2515 | pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j) |
|---|
| 2516 | |
|---|
| 2517 | pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j)) |
|---|
| 2518 | qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j) |
|---|
| 2519 | qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j) |
|---|
| 2520 | qg(i,j)=qg(i,j)-pgsub(i,j) |
|---|
| 2521 | endif |
|---|
| 2522 | |
|---|
| 2523 | !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) **23** |
|---|
| 2524 | |
|---|
| 2525 | ern(i,j)=0.0 |
|---|
| 2526 | |
|---|
| 2527 | if(qr(i,j).gt.0.0) then |
|---|
| 2528 | tair(i,j)=(pt(i,j)+tb0)*pi0 |
|---|
| 2529 | rtair(i,j)=1./(tair(i,j)-c358) |
|---|
| 2530 | qsw(i,j)=rp0*exp(c172-c409*rtair(i,j)) |
|---|
| 2531 | ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0 |
|---|
| 2532 | dm(i,j)=qv(i,j)+qb0-qsw(i,j) |
|---|
| 2533 | rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j) |
|---|
| 2534 | dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0) |
|---|
| 2535 | y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5 |
|---|
| 2536 | y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) & |
|---|
| 2537 | *qsw(i,j)) |
|---|
| 2538 | !cccc |
|---|
| 2539 | ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j) |
|---|
| 2540 | ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.)) |
|---|
| 2541 | pt(i,j)=pt(i,j)-avcp*ern(i,j) |
|---|
| 2542 | qv(i,j)=qv(i,j)+ern(i,j) |
|---|
| 2543 | qr(i,j)=qr(i,j)-ern(i,j) |
|---|
| 2544 | endif |
|---|
| 2545 | |
|---|
| 2546 | ! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0 |
|---|
| 2547 | if (qc(i,j) .le. cmin1) qc(i,j)=0. |
|---|
| 2548 | if (qr(i,j) .le. cmin1) qr(i,j)=0. |
|---|
| 2549 | if (qi(i,j) .le. cmin1) qi(i,j)=0. |
|---|
| 2550 | if (qs(i,j) .le. cmin1) qs(i,j)=0. |
|---|
| 2551 | if (qg(i,j) .le. cmin1) qg(i,j)=0. |
|---|
| 2552 | dpt(i,j,k)=pt(i,j) |
|---|
| 2553 | dqv(i,j,k)=qv(i,j) |
|---|
| 2554 | qcl(i,j,k)=qc(i,j) |
|---|
| 2555 | qrn(i,j,k)=qr(i,j) |
|---|
| 2556 | qci(i,j,k)=qi(i,j) |
|---|
| 2557 | qcs(i,j,k)=qs(i,j) |
|---|
| 2558 | qcg(i,j,k)=qg(i,j) |
|---|
| 2559 | |
|---|
| 2560 | !JJS 275 CONTINUE |
|---|
| 2561 | |
|---|
| 2562 | scc=0. |
|---|
| 2563 | see=0. |
|---|
| 2564 | |
|---|
| 2565 | ! DO 110 J=3,JLES |
|---|
| 2566 | ! DO 110 I=3,ILES |
|---|
| 2567 | ! DD(I,J)=MAX(-CND(I,J), 0.) |
|---|
| 2568 | ! CND(I,J)=MAX(CND(I,J), 0.) |
|---|
| 2569 | ! DD1(I,J)=MAX(-DEP(I,J), 0.) |
|---|
| 2570 | |
|---|
| 2571 | !ccshie 2/21/02 shie follow tao |
|---|
| 2572 | !CC for reference QI(I,J)=QI(I,J)-Y2(I,J)+PIDEP(I,J)*D2T |
|---|
| 2573 | !CC for reference QV(I,J)=QV(I,J)-(PIDEP(I,J))*D2T |
|---|
| 2574 | |
|---|
| 2575 | !c DEP(I,J)=MAX(DEP(I,J), 0.) |
|---|
| 2576 | ! DEP(I,J)=MAX(DEP(I,J), 0.)+PIDEP(I,J)*D2T |
|---|
| 2577 | ! SCC=SCC+CND(I,J) |
|---|
| 2578 | ! SEE=SEE+DD(I,J)+ERN(I,J) |
|---|
| 2579 | |
|---|
| 2580 | ! 110 CONTINUE |
|---|
| 2581 | |
|---|
| 2582 | ! SC(K)=SCC+SC(K) |
|---|
| 2583 | ! SE(K)=SEE+SE(K) |
|---|
| 2584 | |
|---|
| 2585 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2586 | !c henry: please take a look (start) |
|---|
| 2587 | !JJS modified by JJS on 5/1/2007 vvvvv |
|---|
| 2588 | |
|---|
| 2589 | !JJS do 305 j=3,jles |
|---|
| 2590 | !JJS do 305 i=3,iles |
|---|
| 2591 | dd(i,j)=max(-cnd(i,j), 0.) |
|---|
| 2592 | cnd(i,j)=max(cnd(i,j), 0.) |
|---|
| 2593 | dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t |
|---|
| 2594 | dep(i,j)=max(dep(i,j), 0.) |
|---|
| 2595 | !JJS 305 continue |
|---|
| 2596 | |
|---|
| 2597 | !JJS do 306 j=3,jles |
|---|
| 2598 | !JJS do 306 i=3,iles |
|---|
| 2599 | !JJS scc=scc+cnd(i,j) |
|---|
| 2600 | !JJS see=see+(dd(i,j)+ern(i,j)) |
|---|
| 2601 | ! |
|---|
| 2602 | !JJS sddd=sddd+(dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j)) |
|---|
| 2603 | !JJS ssss=ssss+dd1(i,j) |
|---|
| 2604 | !JJS |
|---|
| 2605 | ! shhh=shhh+rsw(i,j,k)*d2t |
|---|
| 2606 | ! sccc=sccc+rlw(i,j,k)*d2t |
|---|
| 2607 | !jjs |
|---|
| 2608 | !JJS smmm=smmm+(psmlt(i,j)+pgmlt(i,j)+pimlt(i,j)) |
|---|
| 2609 | !JJS sfff=sfff+d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) |
|---|
| 2610 | !JJS 1 +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j)) |
|---|
| 2611 | !JJS 2 -qracs(i,j)+pihom(i,j)+pidw(i,j) |
|---|
| 2612 | |
|---|
| 2613 | |
|---|
| 2614 | sccc=cnd(i,j) |
|---|
| 2615 | seee=dd(i,j)+ern(i,j) |
|---|
| 2616 | sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j) |
|---|
| 2617 | ssss=dd1(i,j) + pgsub(i,j) |
|---|
| 2618 | smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j) |
|---|
| 2619 | sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) & |
|---|
| 2620 | +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) & |
|---|
| 2621 | +wgacr(i,j))+pihom(i,j)+pidw(i,j) |
|---|
| 2622 | |
|---|
| 2623 | physc(i,k,j) = avcp * sccc / d2t |
|---|
| 2624 | physe(i,k,j) = avcp * seee / d2t |
|---|
| 2625 | physd(i,k,j) = ascp * sddd / d2t |
|---|
| 2626 | physs(i,k,j) = ascp * ssss / d2t |
|---|
| 2627 | physm(i,k,j) = afcp * smmm / d2t |
|---|
| 2628 | physf(i,k,j) = afcp * sfff / d2t |
|---|
| 2629 | |
|---|
| 2630 | !JJS modified by JJS on 5/1/2007 ^^^^^ |
|---|
| 2631 | |
|---|
| 2632 | 2000 continue |
|---|
| 2633 | |
|---|
| 2634 | 1000 continue |
|---|
| 2635 | |
|---|
| 2636 | !JJS **************************************************************** |
|---|
| 2637 | !JJS convert from GCE grid back to WRF grid |
|---|
| 2638 | do k=kts,kte |
|---|
| 2639 | do j=jts,jte |
|---|
| 2640 | do i=its,ite |
|---|
| 2641 | ptwrf(i,k,j) = dpt(i,j,k) |
|---|
| 2642 | qvwrf(i,k,j) = dqv(i,j,k) |
|---|
| 2643 | qlwrf(i,k,j) = qcl(i,j,k) |
|---|
| 2644 | qrwrf(i,k,j) = qrn(i,j,k) |
|---|
| 2645 | qiwrf(i,k,j) = qci(i,j,k) |
|---|
| 2646 | qswrf(i,k,j) = qcs(i,j,k) |
|---|
| 2647 | qgwrf(i,k,j) = qcg(i,j,k) |
|---|
| 2648 | enddo !i |
|---|
| 2649 | enddo !j |
|---|
| 2650 | enddo !k |
|---|
| 2651 | |
|---|
| 2652 | ! **************************************************************** |
|---|
| 2653 | |
|---|
| 2654 | return |
|---|
| 2655 | END SUBROUTINE saticel_s |
|---|
| 2656 | |
|---|
| 2657 | !JJS |
|---|
| 2658 | !JJS REAL FUNCTION GAMMA(X) |
|---|
| 2659 | !JJS Y=GAMMLN(X) |
|---|
| 2660 | !JJS GAMMA=EXP(Y) |
|---|
| 2661 | !JJS RETURN |
|---|
| 2662 | !JJS END |
|---|
| 2663 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2664 | !JJS real function GAMMLN (xx) |
|---|
| 2665 | real function gammagce (xx) |
|---|
| 2666 | !********************************************************************** |
|---|
| 2667 | real*8 cof(6),stp,half,one,fpf,x,tmp,ser |
|---|
| 2668 | data cof,stp / 76.18009173,-86.50532033,24.01409822, & |
|---|
| 2669 | -1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 / |
|---|
| 2670 | data half,one,fpf / .5, 1., 5.5 / |
|---|
| 2671 | ! |
|---|
| 2672 | x=xx-one |
|---|
| 2673 | tmp=x+fpf |
|---|
| 2674 | tmp=(x+half)*log(tmp)-tmp |
|---|
| 2675 | ser=one |
|---|
| 2676 | do j=1,6 |
|---|
| 2677 | x=x+one |
|---|
| 2678 | ser=ser+cof(j)/x |
|---|
| 2679 | enddo !j |
|---|
| 2680 | gammln=tmp+log(stp*ser) |
|---|
| 2681 | !JJS |
|---|
| 2682 | gammagce=exp(gammln) |
|---|
| 2683 | !JJS |
|---|
| 2684 | return |
|---|
| 2685 | END FUNCTION gammagce |
|---|
| 2686 | |
|---|
| 2687 | END MODULE module_mp_gsfcgce |
|---|
