[16] | 1 | !*********************************************************************** |
---|
| 2 | !* Copyright 2012,2013 * |
---|
| 3 | !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * |
---|
| 4 | !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* |
---|
| 5 | !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * |
---|
| 6 | !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * |
---|
| 7 | !* * |
---|
| 8 | !* This file is part of FLEXPART WRF * |
---|
| 9 | !* * |
---|
| 10 | !* FLEXPART is free software: you can redistribute it and/or modify * |
---|
| 11 | !* it under the terms of the GNU General Public License as published by* |
---|
| 12 | !* the Free Software Foundation, either version 3 of the License, or * |
---|
| 13 | !* (at your option) any later version. * |
---|
| 14 | !* * |
---|
| 15 | !* FLEXPART is distributed in the hope that it will be useful, * |
---|
| 16 | !* but WITHOUT ANY WARRANTY; without even the implied warranty of * |
---|
| 17 | !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
---|
| 18 | !* GNU General Public License for more details. * |
---|
| 19 | !* * |
---|
| 20 | !* You should have received a copy of the GNU General Public License * |
---|
| 21 | !* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
---|
| 22 | !*********************************************************************** |
---|
| 23 | subroutine convmix_kfeta(itime) |
---|
| 24 | ! i |
---|
| 25 | !************************************************************** |
---|
| 26 | ! handles all the calculations related to convective mixing |
---|
| 27 | ! Petra Seibert, Bernd C. Krueger, Feb 2001 |
---|
| 28 | ! nested grids included, Bernd C. Krueger, May 2001 |
---|
| 29 | ! |
---|
| 30 | ! Changes by Caroline Forster, April 2004 - February 2005: |
---|
| 31 | ! convmix called every lsynctime seconds |
---|
| 32 | ! CHANGES by A. Stohl: |
---|
| 33 | ! various run-time optimizations - February 2005 |
---|
| 34 | |
---|
| 35 | ! CHANGED on 10.10.2007 save convective mass fluxes, update them every dt_conv |
---|
| 36 | |
---|
| 37 | ! CHANGED by Weiguo WANG 13 Aug, 2007, use KFeta CU convection scheme |
---|
| 38 | ! |
---|
| 39 | ! Changes by J. Brioude: particles sorting is much more efficient. |
---|
| 40 | ! |
---|
| 41 | ! input for kftea cumulus scheme |
---|
| 42 | ! u1d - 1-d u wind velocity profile |
---|
| 43 | ! v1d - 1-d v wind velocity profile |
---|
| 44 | ! t1d - 1-D temperture (K) |
---|
| 45 | ! qv1d- 1-D water vapor mixin gratio (kg/kg) |
---|
| 46 | ! p1d - 1-D pressure profile (pa) |
---|
| 47 | ! rho1d-1-D density profile (kg/m3) |
---|
| 48 | ! w0avg1d - 1-D vertical velocity (m/s) |
---|
| 49 | ! all above are defined at T-point or P-poit |
---|
| 50 | ! dz1d - dz between full levels |
---|
| 51 | ! delx - grid size of column (m) |
---|
| 52 | ! dt - integraiton time step (s) |
---|
| 53 | ! cudt - cumulus activation time interval (min) |
---|
| 54 | ! kts - starting point in z for convection calculation |
---|
| 55 | ! kte - ending point |
---|
| 56 | |
---|
| 57 | ! output |
---|
| 58 | ! umf - updraft mass flux |
---|
| 59 | ! uer - updraft entrainment flux |
---|
| 60 | ! udr - updraft detrainment flux |
---|
| 61 | ! dmf - downdraft mass flux |
---|
| 62 | ! der - downdraft entrainemnt flux |
---|
| 63 | ! ddr - downdraft detrainment flux |
---|
| 64 | ! cu_top1 -- top of cumulus cloud (index? for height) |
---|
| 65 | ! cu_bot1 -- bottom of cumulus cloud (index) |
---|
| 66 | |
---|
| 67 | !************************************************************** |
---|
| 68 | |
---|
| 69 | use flux_mod |
---|
| 70 | use par_mod |
---|
| 71 | use com_mod |
---|
| 72 | use conv_mod |
---|
| 73 | |
---|
| 74 | implicit none |
---|
| 75 | |
---|
| 76 | integer :: igr,igrold, ipart, itime, ix, j, inest |
---|
| 77 | integer :: ipconv |
---|
| 78 | integer :: jy, kpart, ktop, ngrid,kz,kzp,a |
---|
| 79 | integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests) |
---|
| 80 | ! itime [s] current time |
---|
| 81 | ! igrid(maxpart) horizontal grid position of each particle |
---|
| 82 | ! igridn(maxpart,maxnests) dto. for nested grids |
---|
| 83 | ! ipoint(maxpart) pointer to access particles according to grid position |
---|
| 84 | |
---|
| 85 | logical :: lconv,warm_rain |
---|
| 86 | real :: x, y, xtn,ytn, ztold, delt |
---|
| 87 | real :: dt1,dt2,dtt,dummy |
---|
| 88 | real :: duma, dumz(nuvzmax+1) |
---|
| 89 | integer :: mind1,mind2 |
---|
| 90 | ! dt1,dt2,dtt,mind1,mind2 variables used for time interpolation |
---|
| 91 | integer :: itage,nage,duminc |
---|
| 92 | real,parameter :: eps=nxmax/3.e5 |
---|
| 93 | |
---|
| 94 | integer :: i,k,numberp(maxpart) |
---|
| 95 | !-- for KFeta |
---|
| 96 | real, dimension(nuvzmax):: u1d,v1d,t1d,dz1d,qv1d,p1d, & |
---|
| 97 | rho1d,w0avg1d,umf,uer,udr, & |
---|
| 98 | dmf,der,ddr,zh |
---|
| 99 | real :: cudt,delx,dt,cu_bot1,cu_top1,zp,fmix |
---|
| 100 | real, dimension(nuvzmax+1):: umfzf,dmfzf,zf |
---|
| 101 | integer :: kts,kte,if_update |
---|
| 102 | ! |
---|
| 103 | |
---|
| 104 | ! monitoring variables |
---|
| 105 | real :: sumconv,sumall,sumpart |
---|
| 106 | integer :: sumpartgrid(1000000) |
---|
| 107 | |
---|
| 108 | ! print *, "IN convmix_kfeta" |
---|
| 109 | ! write(*,'(//a,a//)') |
---|
| 110 | ! & '*** Stopping in subr. convmix ***', |
---|
| 111 | ! & ' This is not implemented for FLEXPART_WRF' |
---|
| 112 | ! stop |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | ! Calculate auxiliary variables for time interpolation |
---|
| 116 | !***************************************************** |
---|
| 117 | |
---|
| 118 | delt=real(abs(lsynctime)) |
---|
| 119 | |
---|
| 120 | ! dt_conv is given from input namelist |
---|
| 121 | ! dt_conv=3600.0 |
---|
| 122 | if_update=0 |
---|
| 123 | if ( mod(real(itime),dt_conv) .eq. 0 ) if_update=1 |
---|
| 124 | ! print*,'conv itime',itime,dt_conv,if_update |
---|
| 125 | |
---|
| 126 | ! delt=dt_conv |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | dt1=real(itime-memtime(1)) |
---|
| 130 | dt2=real(memtime(2)-itime) |
---|
| 131 | dtt=1./(dt1+dt2) |
---|
| 132 | mind1=memind(1) |
---|
| 133 | mind2=memind(2) |
---|
| 134 | |
---|
| 135 | lconv = .false. |
---|
| 136 | |
---|
| 137 | ! for KFeta |
---|
| 138 | ! warm_rain=.true. ! depends on mp_physics in WRF, may add an option in the future |
---|
| 139 | cudt = 10.0 ! cumulus para is called every 10 min in a time step, if dt < cudt, call once |
---|
| 140 | |
---|
| 141 | kts=1 |
---|
| 142 | kte=nuvz-1 |
---|
| 143 | |
---|
| 144 | ! if no particles are present return after initialization |
---|
| 145 | !******************************************************** |
---|
| 146 | |
---|
| 147 | if (numpart.le.0) return |
---|
| 148 | |
---|
| 149 | ! Assign igrid and igridn, which are pseudo grid numbers indicating particles |
---|
| 150 | ! that are outside the part of the grid under consideration |
---|
| 151 | ! (e.g. particles near the poles or particles in other nests). |
---|
| 152 | ! Do this for all nests but use only the innermost nest; for all others |
---|
| 153 | ! igrid shall be -1 |
---|
| 154 | ! Also, initialize index vector ipoint |
---|
| 155 | !************************************************************************ |
---|
| 156 | ! print*,'step1' |
---|
| 157 | do ipart=1,numpart |
---|
| 158 | igrid(ipart)=-1 |
---|
| 159 | do j=numbnests,1,-1 |
---|
| 160 | igridn(ipart,j)=-1 |
---|
| 161 | enddo |
---|
| 162 | ipoint(ipart)=ipart |
---|
| 163 | ! do not consider particles that are (yet) not part of simulation |
---|
| 164 | if (itra1(ipart).ne.itime) goto 20 |
---|
| 165 | x = xtra1(ipart) |
---|
| 166 | y = ytra1(ipart) |
---|
| 167 | |
---|
| 168 | ! Determine which nesting level to be used |
---|
| 169 | !********************************************************** |
---|
| 170 | |
---|
| 171 | ngrid=0 |
---|
| 172 | do j=numbnests,1,-1 |
---|
| 173 | if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. & |
---|
| 174 | y.gt.yln(j) .and. y.lt.yrn(j) ) then |
---|
| 175 | ngrid=j |
---|
| 176 | goto 23 |
---|
| 177 | endif |
---|
| 178 | end do |
---|
| 179 | 23 continue |
---|
| 180 | |
---|
| 181 | ! Determine nested grid coordinates |
---|
| 182 | !********************************** |
---|
| 183 | |
---|
| 184 | if (ngrid.gt.0) then |
---|
| 185 | ! nested grids |
---|
| 186 | xtn=(x-xln(ngrid))*xresoln(ngrid) |
---|
| 187 | ytn=(y-yln(ngrid))*yresoln(ngrid) |
---|
| 188 | ix=nint(xtn) |
---|
| 189 | jy=nint(ytn) |
---|
| 190 | igridn(ipart,ngrid) = 1 + jy*nxn(ngrid) + ix |
---|
| 191 | else if(ngrid.eq.0) then |
---|
| 192 | ! mother grid |
---|
| 193 | ix=nint(x) |
---|
| 194 | jy=nint(y) |
---|
| 195 | igrid(ipart) = 1 + jy*nx + ix |
---|
| 196 | endif |
---|
| 197 | |
---|
| 198 | 20 continue |
---|
| 199 | end do |
---|
| 200 | |
---|
| 201 | sumpart = 0. |
---|
| 202 | sumconv = 0. |
---|
| 203 | |
---|
| 204 | ! print*,'step2' |
---|
| 205 | |
---|
| 206 | !************************************************************************************** |
---|
| 207 | ! 1. Now, do everything for the mother domain and, later, for all of the nested domains |
---|
| 208 | ! While all particles have to be considered for redistribution, the Emanuel convection |
---|
| 209 | ! scheme only needs to be called once for every grid column where particles are present. |
---|
| 210 | ! Therefore, particles are sorted according to their grid position. Whenever a new grid |
---|
| 211 | ! cell is encountered by looping through the sorted particles, the convection scheme is called. |
---|
| 212 | !************************************************************************************** |
---|
| 213 | delx = dx |
---|
| 214 | |
---|
| 215 | ! sort particles according to horizontal position and calculate index vector IPOINT |
---|
| 216 | |
---|
| 217 | call sort2(numpart,igrid,ipoint) |
---|
| 218 | ! print*,'step2 after sort',minval(igrid),maxval(igrid),numpart |
---|
| 219 | |
---|
| 220 | ! count particle # in each column |
---|
| 221 | ! do 40 i=abs(minval(igrid)),maxval(igrid) |
---|
| 222 | ! sumpart=0. |
---|
| 223 | ! do 41 k=1,numpart |
---|
| 224 | ! if(igrid(k) .eq. i) then |
---|
| 225 | ! sumpart=sumpart+1 |
---|
| 226 | ! endif |
---|
| 227 | !41 continue |
---|
| 228 | ! do 42 k=1,numparc Changes by J. Brioude: the sort of particles is much more efficient.t |
---|
| 229 | !42 if(igrid(k) .eq. i) numberp(k)=int(sumpart) |
---|
| 230 | !40 continue |
---|
| 231 | ! JB |
---|
| 232 | if (maxval(igrid).gt.1000000) then |
---|
| 233 | print*,'too much x and y grid. modify convmix_kfeta.f' |
---|
| 234 | stop |
---|
| 235 | endif |
---|
| 236 | do k=1,1000000 |
---|
| 237 | sumpartgrid(k)=0 |
---|
| 238 | enddo |
---|
| 239 | do k=1,numpart |
---|
| 240 | if (igrid(k).gt.0) sumpartgrid(igrid(k))=sumpartgrid(igrid(k))+1 |
---|
| 241 | enddo |
---|
| 242 | do k=1,numpart |
---|
| 243 | if (igrid(k).gt.0) then |
---|
| 244 | numberp(k)=sumpartgrid(igrid(k)) |
---|
| 245 | else |
---|
| 246 | numberp(k)=0 |
---|
| 247 | endif |
---|
| 248 | enddo |
---|
| 249 | |
---|
| 250 | ! print*,'step3',numpart |
---|
| 251 | |
---|
| 252 | ! Now visit all grid columns where particles are present |
---|
| 253 | ! by going through the sorted particles |
---|
| 254 | |
---|
| 255 | igrold = -1 |
---|
| 256 | a=0 |
---|
| 257 | do kpart=1,numpart |
---|
| 258 | igr = igrid(kpart) |
---|
| 259 | if (igr .eq. -1 .or. numberp(kpart).le.20 & |
---|
| 260 | ! if (igr .eq. -1 |
---|
| 261 | ) goto 50 |
---|
| 262 | ipart = ipoint(kpart) |
---|
| 263 | |
---|
| 264 | ! sumall = sumall + 1 |
---|
| 265 | !c For one column, we only need to compute 1D met once |
---|
| 266 | |
---|
| 267 | if (igr .ne. igrold) then |
---|
| 268 | sumconv=sumconv+1 |
---|
| 269 | ! we are in a new grid column |
---|
| 270 | jy = (igr-1)/nx |
---|
| 271 | ix = igr - jy*nx - 1 |
---|
| 272 | a=a+1 |
---|
| 273 | ! print*,'a',a |
---|
| 274 | ! Interpolate all meteorological data needed for the convection scheme |
---|
| 275 | |
---|
| 276 | do kz=1,nuvz-1 ! nconvlev+1 |
---|
| 277 | ! FLEXPART_WRF - used add_sfc_level for the shifting |
---|
| 278 | ! for W, it is not shifted, make sure w is 'true' vertical velocity! |
---|
| 279 | |
---|
| 280 | kzp = kz + add_sfc_level |
---|
| 281 | u1d(kz)=(u_wrf(ix,jy,kzp,mind1)*dt2+ & |
---|
| 282 | u_wrf(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 283 | v1d(kz)=(v_wrf(ix,jy,kzp,mind1)*dt2+ & |
---|
| 284 | v_wrf(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 285 | t1d(kz)=(tth(ix,jy,kzp,mind1)*dt2+ & |
---|
| 286 | tth(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 287 | qv1d(kz)=(qvh(ix,jy,kzp,mind1)*dt2+ & |
---|
| 288 | qvh(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 289 | p1d(kz)=(pph(ix,jy,kzp,mind1)*dt2+ & |
---|
| 290 | pph(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 291 | dz1d(kz)=(zzh(ix,jy,kzp+1,mind1)*dt2+ & |
---|
| 292 | zzh(ix,jy,kzp+1,mind2)*dt1)*dtt- & |
---|
| 293 | (zzh(ix,jy,kzp,mind1)*dt2+ & |
---|
| 294 | zzh(ix,jy,kzp,mind2)*dt1)*dtt |
---|
| 295 | w0avg1d(kz)=(w_wrf(ix,jy,kz,mind1)*dt2+ & |
---|
| 296 | w_wrf(ix,jy,kz,mind2)*dt1)*dtt+ & |
---|
| 297 | (w_wrf(ix,jy,kz+1,mind1)*dt2+ & |
---|
| 298 | w_wrf(ix,jy,kz+1,mind2)*dt1)*dtt |
---|
| 299 | w0avg1d(kz)=0.5*w0avg1d(kz) |
---|
| 300 | rho1d(kz)=p1d(kz)/ & |
---|
| 301 | (t1d(kz)*(1.0+0.608*qv1d(kz))) & |
---|
| 302 | /287.0 |
---|
| 303 | |
---|
| 304 | ! write(*,'(1x,I10,10F10.2)')kz,u1d(kz),v1d(kz),w0avg1d(kz), |
---|
| 305 | ! & t1d(kz),qv1d(kz),p1d(kz)/100,dz1d(kz) |
---|
| 306 | |
---|
| 307 | |
---|
| 308 | enddo |
---|
| 309 | |
---|
| 310 | ! -- old convective mass fluxes |
---|
| 311 | do k=kts,kte |
---|
| 312 | umf(k)=umf3(ix,jy,k) |
---|
| 313 | uer(k)=uer3(ix,jy,k) |
---|
| 314 | udr(k)=udr3(ix,jy,k) |
---|
| 315 | dmf(k)=dmf3(ix,jy,k) |
---|
| 316 | der(k)=der3(ix,jy,k) |
---|
| 317 | ddr(k)=ddr3(ix,jy,k) |
---|
| 318 | enddo |
---|
| 319 | cu_top1=cu_top(ix,jy) |
---|
| 320 | cu_bot1=cu_bot(ix,jy) |
---|
| 321 | |
---|
| 322 | |
---|
| 323 | ! write(*,*)'1-D wind' |
---|
| 324 | |
---|
| 325 | ! Calculate convection flux, updrought flux, entrainment, detrainment flux |
---|
| 326 | ! downdraft flux, entrainment ,detrainment flux |
---|
| 327 | warm_rain=.false. |
---|
| 328 | if (mp_physics .eq. 1) warm_rain = .true. |
---|
| 329 | ! -- Update fluxes |
---|
| 330 | if (if_update .eq. 1 ) then !! if_update |
---|
| 331 | ! write(*,*)'update convective fluxes, itime=',itime/3600. |
---|
| 332 | ! print*,u1d(4:8),v1d(4:8),t1d(4:8),dz1d(4:8),qv1d(4:8) |
---|
| 333 | ! print*,p1d(4:8),rho1d(4:8),w0avg1d(4:8) |
---|
| 334 | ! print*,cudt,delx,dt_conv,warm_rain |
---|
| 335 | ! print*,'attend' |
---|
| 336 | ! pause |
---|
| 337 | CALL KF_ETA(nuvzmax,u1d,v1d,t1d,dz1d,qv1d,p1d, & ! IN |
---|
| 338 | rho1d,w0avg1d,cudt,delx,dt_conv,warm_rain,kts,kte, & ! IN |
---|
| 339 | umf,uer,udr,dmf,der,ddr,cu_bot1,cu_top1) ! OUT |
---|
| 340 | dummy=0. |
---|
| 341 | do k=kts,kte |
---|
| 342 | umf3(ix,jy,k)=umf(k) |
---|
| 343 | uer3(ix,jy,k)=uer(k) |
---|
| 344 | udr3(ix,jy,k)=udr(k) |
---|
| 345 | dmf3(ix,jy,k)=dmf(k) |
---|
| 346 | der3(ix,jy,k)=der(k) |
---|
| 347 | ddr3(ix,jy,k)=ddr(k) |
---|
| 348 | dummy=dummy+umf(k)+uer(k)+udr(k) |
---|
| 349 | enddo |
---|
| 350 | if (dummy.gt.0.) then |
---|
| 351 | ! print*,'dummy',dummy |
---|
| 352 | duminc=1 |
---|
| 353 | else |
---|
| 354 | duminc=0 |
---|
| 355 | endif |
---|
| 356 | cu_top(ix,jy)=cu_top1 |
---|
| 357 | cu_bot(ix,jy)=cu_bot1 |
---|
| 358 | ! if (cu_top1.gt.1.) print*,'cu h',cu_bot1,cu_top1 |
---|
| 359 | endif !! if_update |
---|
| 360 | ! if (a.gt.2000) then |
---|
| 361 | ! print*,'after kf_eta',a |
---|
| 362 | ! a=0 |
---|
| 363 | ! endif |
---|
| 364 | |
---|
| 365 | ! write(*,*)'ix,jy=',ix,jy,itime |
---|
| 366 | ! write(*,*)'previous column part#=',sumpart |
---|
| 367 | ! write(*,*)'FLUX,k,umf,uer,udr,dmf,der,ddr' |
---|
| 368 | ! write(*,*)'cu_bot1,cu_top1=',cu_bot1,cu_top1 |
---|
| 369 | ! if (cu_top1 .lt. cu_bot1) write(*,*)'umf=', umf(1),umf(10) |
---|
| 370 | |
---|
| 371 | ! do kz=kts,kte |
---|
| 372 | ! write(*,'(1x,I10,10E10.2)')kz,umf(kz),uer(kz),udr(kz),dmf(kz) |
---|
| 373 | ! & ,der(kz),ddr(kz),p1d(kz)/100 |
---|
| 374 | |
---|
| 375 | ! enddo |
---|
| 376 | |
---|
| 377 | sumpart=0 |
---|
| 378 | IF (cu_top1 .gt. cu_bot1+1 ) then ! lconv |
---|
| 379 | ! write(250+lconvection,*)'-1',itime,ix,jy |
---|
| 380 | lconv = .true. |
---|
| 381 | |
---|
| 382 | ! Prepare data for redistributing particle |
---|
| 383 | |
---|
| 384 | CALL pre_redist_kf(nuvzmax,nuvz,umf,dmf,dz1d,p1d,delx,delt, & ! IN |
---|
| 385 | cu_bot1,cu_top1, & ! IN |
---|
| 386 | zf,zh, & ! OUT |
---|
| 387 | umfzf,dmfzf,fmix) ! OUT |
---|
| 388 | |
---|
| 389 | else |
---|
| 390 | lconv= .false. |
---|
| 391 | |
---|
| 392 | ENDIF ! lconv |
---|
| 393 | |
---|
| 394 | |
---|
| 395 | igrold = igr |
---|
| 396 | ktop = 0 |
---|
| 397 | endif |
---|
| 398 | |
---|
| 399 | sumpart=sumpart+1 |
---|
| 400 | ! treat particle only if column has convection |
---|
| 401 | if (lconv .eqv. .true.) then |
---|
| 402 | ! assign new vertical position to particle |
---|
| 403 | ! ztold=ztra1(ipart) |
---|
| 404 | zp=ztra1(ipart) |
---|
| 405 | !C write(*,*)'befrore convection zp= ',zp |
---|
| 406 | |
---|
| 407 | ! write(*,*)'part No =',sumpart |
---|
| 408 | CALL redist_kf(lconvection,ldirect,delt,delx, & ! IN |
---|
| 409 | dz1d,nzmax, nz,umf,uer,udr,dmf, & ! IN |
---|
| 410 | der,ddr,rho1d, & ! IN |
---|
| 411 | zf,zh, & ! IN |
---|
| 412 | umfzf,dmfzf,fmix, & ! IN |
---|
| 413 | zp) ! IN/OUT |
---|
| 414 | |
---|
| 415 | if (zp .lt. 0.0) zp=-1.0*zp |
---|
| 416 | if (zp .gt. height(nz)-0.5) & |
---|
| 417 | zp = height(nz)-0.5 |
---|
| 418 | ! if (abs(zp-ztra1(ipart)) .ge. 1e-5) |
---|
| 419 | ! &write(250+lconvection,*)ztra1(ipart),zp,zp-ztra1(ipart) |
---|
| 420 | ! if (duminc.eq.1) |
---|
| 421 | ! + print*,'true conv',dummy,zp-ztra1(ipart),cu_top1-cu_bot1 |
---|
| 422 | |
---|
| 423 | ztra1(ipart) = zp |
---|
| 424 | !C write(*,*)'after convection, zp=',zp |
---|
| 425 | |
---|
| 426 | !C OLD call redist(ipart,ktop,ipconv) |
---|
| 427 | ! if (ipconv.le.0) sumconv = sumconv+1 |
---|
| 428 | |
---|
| 429 | ! Calculate the gross fluxes across layer interfaces |
---|
| 430 | !*************************************************** |
---|
| 431 | |
---|
| 432 | if (iflux.eq.1) then |
---|
| 433 | itage=abs(itra1(ipart)-itramem(ipart)) |
---|
| 434 | do nage=1,nageclass |
---|
| 435 | if (itage.lt.lage(nage)) goto 37 |
---|
| 436 | enddo |
---|
| 437 | 37 continue |
---|
| 438 | ! print*,'step 4' |
---|
| 439 | |
---|
| 440 | if (nage.le.nageclass) & |
---|
| 441 | call calcfluxes(nage,ipart,real(xtra1(ipart)), & |
---|
| 442 | real(ytra1(ipart)),ztold) |
---|
| 443 | endif |
---|
| 444 | |
---|
| 445 | endif !(lconv .eqv. .true) |
---|
| 446 | enddo |
---|
| 447 | 50 continue |
---|
| 448 | |
---|
| 449 | ! write(*,*)'total convective columns=',sumconv, |
---|
| 450 | ! & 'time=', 1.0*itime/3600 |
---|
| 451 | |
---|
| 452 | !*********************************************************************************** |
---|
| 453 | ! 2. Nested domains |
---|
| 454 | !*********************************************************************************** |
---|
| 455 | |
---|
| 456 | ! sort particles according to horizontal position and calculate index vector IPOINT |
---|
| 457 | |
---|
| 458 | do inest=1,numbnests |
---|
| 459 | delx = dxn(inest) |
---|
| 460 | if (delx .le. 10000.0) goto 70 ! for small grid size, no need to do convection |
---|
| 461 | do ipart=1,numpart |
---|
| 462 | ipoint(ipart)=ipart |
---|
| 463 | igrid(ipart) = igridn(ipart,inest) |
---|
| 464 | enddo |
---|
| 465 | call sort2(numpart,igrid,ipoint) |
---|
| 466 | |
---|
| 467 | ! print*,'step in nest' |
---|
| 468 | ! Now visit all grid columns where particles are present |
---|
| 469 | ! by going through the sorted particles |
---|
| 470 | |
---|
| 471 | igrold = -1 |
---|
| 472 | do kpart=1,numpart |
---|
| 473 | igr = igrid(kpart) |
---|
| 474 | if (igr .eq. -1) goto 60 |
---|
| 475 | ipart = ipoint(kpart) |
---|
| 476 | ! sumall = sumall + 1 |
---|
| 477 | if (igr .ne. igrold) then |
---|
| 478 | ! we are in a new grid column |
---|
| 479 | jy = (igr-1)/nxn(inest) |
---|
| 480 | ix = igr - jy*nxn(inest) - 1 |
---|
| 481 | |
---|
| 482 | ! Interpolate all meteorological data needed for the convection scheme |
---|
| 483 | |
---|
| 484 | ! Interpolate all meteorological data needed for the convection scheme |
---|
| 485 | |
---|
| 486 | do kz=1,nuvz ! nconvlev+1 |
---|
| 487 | ! FLEXPART_WRF - used add_sfc_level for the shifting |
---|
| 488 | ! for W, it is not shifted, make sure w is 'true' vertical velocity! |
---|
| 489 | |
---|
| 490 | kzp = kz + add_sfc_level |
---|
| 491 | u1d(kz)=(un_wrf(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 492 | un_wrf(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 493 | v1d(kz)=(vn_wrf(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 494 | vn_wrf(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 495 | t1d(kz)=(tthn(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 496 | tthn(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 497 | qv1d(kz)=(qvhn(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 498 | qvhn(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 499 | p1d(kz)=(pphn(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 500 | pphn(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 501 | dz1d(kz)=(zzhn(ix,jy,kzp+1,mind1,inest)*dt2+ & |
---|
| 502 | zzhn(ix,jy,kzp+1,mind2,inest)*dt1)*dtt- & |
---|
| 503 | (zzhn(ix,jy,kzp,mind1,inest)*dt2+ & |
---|
| 504 | zzhn(ix,jy,kzp,mind2,inest)*dt1)*dtt |
---|
| 505 | w0avg1d(kz)=(wn_wrf(ix,jy,kz,mind1,inest)*dt2+ & |
---|
| 506 | wn_wrf(ix,jy,kz,mind2,inest)*dt1)*dtt+ & |
---|
| 507 | (wn_wrf(ix,jy,kz+1,mind1,inest)*dt2+ & |
---|
| 508 | wn_wrf(ix,jy,kz+1,mind2,inest)*dt1)*dtt |
---|
| 509 | w0avg1d(kz)=0.5*w0avg1d(kz) |
---|
| 510 | rho1d(kz)=p1d(kz)/ & |
---|
| 511 | (t1d(kz)*(1.0+0.608*qv1d(kz))) & |
---|
| 512 | /287.0 |
---|
| 513 | |
---|
| 514 | enddo |
---|
| 515 | |
---|
| 516 | !C Old convective mass fluxes |
---|
| 517 | do k=kts,kte |
---|
| 518 | umf(k)=umf3n(ix,jy,k,inest) |
---|
| 519 | uer(k)=uer3n(ix,jy,k,inest) |
---|
| 520 | udr(k)=udr3n(ix,jy,k,inest) |
---|
| 521 | dmf(k)=dmf3n(ix,jy,k,inest) |
---|
| 522 | der(k)=der3n(ix,jy,k,inest) |
---|
| 523 | ddr(k)=ddr3n(ix,jy,k,inest) |
---|
| 524 | enddo |
---|
| 525 | cu_top1=cu_topn(ix,jy,inest) |
---|
| 526 | cu_bot1=cu_botn(ix,jy,inest) |
---|
| 527 | |
---|
| 528 | |
---|
| 529 | |
---|
| 530 | !alculate convection flux, updrought flux, entrainment, detrainment flux |
---|
| 531 | ! downdraft flux, entrainment ,detrainment flux |
---|
| 532 | warm_rain = .false. |
---|
| 533 | if (mp_physicsn(inest) .eq. 1 ) warm_rain = .true. |
---|
| 534 | |
---|
| 535 | if (if_update .eq. 1 ) then !!! update |
---|
| 536 | CALL KF_ETA(nuvzmax,u1d,v1d,t1d,dz1d,qv1d,p1d, & ! IN |
---|
| 537 | rho1d,w0avg1d,cudt,delx,dt_conv,warm_rain,kts,kte, & ! IN |
---|
| 538 | umf,uer,udr,dmf,der,ddr,cu_bot1,cu_top1) ! OUT |
---|
| 539 | do k=kts,kte |
---|
| 540 | umf3n(ix,jy,k,inest)=umf(k) |
---|
| 541 | uer3n(ix,jy,k,inest)=uer(k) |
---|
| 542 | udr3n(ix,jy,k,inest)=udr(k) |
---|
| 543 | dmf3n(ix,jy,k,inest)=dmf(k) |
---|
| 544 | der3n(ix,jy,k,inest)=der(k) |
---|
| 545 | ddr3n(ix,jy,k,inest)=ddr(k) |
---|
| 546 | enddo |
---|
| 547 | cu_topn(ix,jy,inest)=cu_top1 |
---|
| 548 | cu_botn(ix,jy,inest)=cu_bot1 |
---|
| 549 | endif !!! update |
---|
| 550 | |
---|
| 551 | |
---|
| 552 | |
---|
| 553 | IF (cu_top1 .gt. cu_bot1) then ! lconv |
---|
| 554 | |
---|
| 555 | lconv = .true. |
---|
| 556 | |
---|
| 557 | ! Prepare data for redistributing particle |
---|
| 558 | |
---|
| 559 | CALL pre_redist_kf(nuvzmax,nuvz,umf,dmf,dz1d,p1d,delx,delt, & ! IN |
---|
| 560 | cu_bot1,cu_top1, & ! IN |
---|
| 561 | zf,zh, & ! OUT |
---|
| 562 | umfzf,dmfzf,fmix) ! OUT |
---|
| 563 | |
---|
| 564 | else |
---|
| 565 | lconv = .false. |
---|
| 566 | ENDIF ! lconv |
---|
| 567 | |
---|
| 568 | |
---|
| 569 | igrold = igr |
---|
| 570 | ktop = 0 |
---|
| 571 | endif |
---|
| 572 | |
---|
| 573 | ! treat particle only if column has convection |
---|
| 574 | if (lconv .eqv. .true.) then |
---|
| 575 | ! assign new vertical position to particle |
---|
| 576 | |
---|
| 577 | ! ztold=ztra1(ipart) |
---|
| 578 | zp=ztra1(ipart) |
---|
| 579 | |
---|
| 580 | CALL redist_kf(lconvection,ldirect,delt,delx, & ! IN |
---|
| 581 | dz1d,nzmax, nz,umf,uer,udr,dmf, & ! IN |
---|
| 582 | der,ddr,rho1d, & ! IN |
---|
| 583 | zf,zh, & ! IN |
---|
| 584 | umfzf,dmfzf, & ! IN |
---|
| 585 | zp) ! IN/OUT |
---|
| 586 | |
---|
| 587 | if (zp .lt. 0.0) zp=-1.0*zp |
---|
| 588 | if (zp .gt. height(nz)-0.5) & |
---|
| 589 | zp = height(nz)-0.5 |
---|
| 590 | ztra1(ipart) = zp |
---|
| 591 | |
---|
| 592 | ! Calculate the gross fluxes across layer interfaces |
---|
| 593 | !*************************************************** |
---|
| 594 | |
---|
| 595 | if (iflux.eq.1) then |
---|
| 596 | itage=abs(itra1(ipart)-itramem(ipart)) |
---|
| 597 | do nage=1,nageclass |
---|
| 598 | if (itage.lt.lage(nage)) goto 47 |
---|
| 599 | enddo |
---|
| 600 | 47 continue |
---|
| 601 | |
---|
| 602 | if (nage.le.nageclass) & |
---|
| 603 | call calcfluxes(nage,ipart,real(xtra1(ipart)), & |
---|
| 604 | real(ytra1(ipart)),ztold) |
---|
| 605 | endif |
---|
| 606 | |
---|
| 607 | endif !(lconv .eqv. .true.) |
---|
| 608 | |
---|
| 609 | enddo |
---|
| 610 | 60 continue |
---|
| 611 | enddo |
---|
| 612 | 70 continue !inest - loop |
---|
| 613 | ! print*,'end of convmix' |
---|
| 614 | !-------------------------------------------------------------------------- |
---|
| 615 | ! write(*,*)'############################################' |
---|
| 616 | ! write(*,*)'TIME=', |
---|
| 617 | ! & itime |
---|
| 618 | ! write(*,*)'fraction of particles under convection', |
---|
| 619 | ! & sumconv/(sumall+0.001) |
---|
| 620 | ! write(*,*)'total number of particles', |
---|
| 621 | ! & sumall |
---|
| 622 | ! write(*,*)'number of particles under convection', |
---|
| 623 | ! & sumconv |
---|
| 624 | ! write(*,*)'############################################' |
---|
| 625 | |
---|
| 626 | return |
---|
| 627 | end subroutine convmix_kfeta |
---|