| 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 | |
|---|
| 24 | ! THIS CODE IS TO REDISTRIBUTE PARTICLES INVOLVED IN UPDRAFT OR/AND DOWDRAFT |
|---|
| 25 | ! TWO OPTIONS: 2- simply well-mixed inside updraft |
|---|
| 26 | ! mixing is based on convective mass flux, downdraft is not considered |
|---|
| 27 | ! 3- particle positions based on entrainment/detrainment rates |
|---|
| 28 | ! if they are available |
|---|
| 29 | ! 0- no convection |
|---|
| 30 | ! 1- has been used by old code |
|---|
| 31 | |
|---|
| 32 | ! SUBROUTINE NAME: redist_kf |
|---|
| 33 | ! CALLED by convmix_kf.f |
|---|
| 34 | ! meteoroloy data are from pre_redist_kf.f that is called inside convmix_kf.f |
|---|
| 35 | ! |
|---|
| 36 | ! INPUT: |
|---|
| 37 | ! dx -- horizontal grid size |
|---|
| 38 | ! nuvzmax -- # of layer for mass flux |
|---|
| 39 | ! umf -- updraft mass flux |
|---|
| 40 | ! uer -- updraft entrament flux |
|---|
| 41 | ! udr -- updraft detrainment flux |
|---|
| 42 | ! dmf -- downdraft mass flux |
|---|
| 43 | ! der -- downdraft entrainment flux |
|---|
| 44 | ! ddr -- downdraft detrainment flux |
|---|
| 45 | ! rho1d -- air density |
|---|
| 46 | ! dz1d-- delt z between full levels |
|---|
| 47 | ! ldirect -- flag for forward (+1) or backward (-1) run |
|---|
| 48 | ! delt -- time step (s) |
|---|
| 49 | ! umfzf -- normalized up-mass-flux weighted distance, |
|---|
| 50 | ! dmfzf -- normalized down-mass-flux weighted distance,(not used) |
|---|
| 51 | ! fmix -- estimated fraction of particles that will be mixed (for LCONVECTION=2) |
|---|
| 52 | ! zf -- distance above ground at full levels |
|---|
| 53 | ! zh -- """"""""""""""""""""""""half levels |
|---|
| 54 | |
|---|
| 55 | ! IN/OUTPUT: |
|---|
| 56 | ! zp -- particle z position (IN/OUT) |
|---|
| 57 | |
|---|
| 58 | Subroutine redist_kf(mix_option,ldirect,delt,dx, & ! IN |
|---|
| 59 | dz1d,nuvzmax, nuvz,umf,uer,udr,dmf, & ! IN |
|---|
| 60 | der,ddr,rho1d, & ! IN |
|---|
| 61 | zf,zh, & ! IN |
|---|
| 62 | umfzf,dmfzf,fmix, & ! IN |
|---|
| 63 | zp) ! IN/OUT |
|---|
| 64 | |
|---|
| 65 | |
|---|
| 66 | IMPLICIT NONE |
|---|
| 67 | integer :: nuvzmax,nuvz,ldirect,i,j,k,kk,mix_option |
|---|
| 68 | integer,parameter :: well_mix=2 |
|---|
| 69 | integer,parameter :: prob_mix=3 |
|---|
| 70 | real,dimension(nuvzmax):: umf,uer,udr,zh,dz1d |
|---|
| 71 | real,dimension(nuvzmax):: dmf,der,ddr,rho1d |
|---|
| 72 | real,dimension(nuvzmax+1):: zf,umfzf,dmfzf |
|---|
| 73 | real :: zp,delt,dx,fup,fdown,ftotal,totalmass |
|---|
| 74 | real :: rn,rn1,rn2,ddz,fde,w_sub,ran3,fmix |
|---|
| 75 | real :: uptop,downtop ! top of updraft and downdraft |
|---|
| 76 | ! data well_mix,prob_mix/2,3/ |
|---|
| 77 | |
|---|
| 78 | ! Check if particle position is below convection cloud top |
|---|
| 79 | ! find top height of updraft |
|---|
| 80 | uptop=0.0 |
|---|
| 81 | do i=nuvz,1,-1 |
|---|
| 82 | if (umf(i) .gt. 1e-20) then |
|---|
| 83 | if (i .eq. nuvz) uptop = zh(nuvz) |
|---|
| 84 | if (i .lt. nuvz) uptop = zh(i+1) |
|---|
| 85 | goto 81 |
|---|
| 86 | endif |
|---|
| 87 | enddo |
|---|
| 88 | 81 continue ! updraft top |
|---|
| 89 | |
|---|
| 90 | downtop=0.0 |
|---|
| 91 | do i=nuvz,1,-1 |
|---|
| 92 | if (abs(dmf(i)) .gt. 1e-20) then |
|---|
| 93 | if (i .eq. nuvz) downtop=zh(nuvz) |
|---|
| 94 | if (i .lt. nuvz) downtop=zh(i+1) |
|---|
| 95 | goto 82 |
|---|
| 96 | endif |
|---|
| 97 | enddo |
|---|
| 98 | 82 continue ! downdraft top |
|---|
| 99 | |
|---|
| 100 | |
|---|
| 101 | if (zp .ge. uptop) goto 89 ! no convective adjustment |
|---|
| 102 | |
|---|
| 103 | DO i=1,nuvz |
|---|
| 104 | if (zp .ge. zf(i) .and. zp .lt. zf(i+1)) goto 90 |
|---|
| 105 | ENDDO |
|---|
| 106 | 90 kk=i ! kk is grid # of the particle before reposition |
|---|
| 107 | totalmass=rho1d(kk)*dz1d(kk)*dx*dx ! air mass kg |
|---|
| 108 | |
|---|
| 109 | |
|---|
| 110 | !+++++++++ |
|---|
| 111 | ! NOw start to reposition particle |
|---|
| 112 | ! |
|---|
| 113 | ! Simply mix all particles by assigning a random position within cloud |
|---|
| 114 | ! In this case, backward run is not wroking since this reposition process is not |
|---|
| 115 | ! reversible. |
|---|
| 116 | |
|---|
| 117 | IF (mix_option .eq. well_mix) then ! Well-mixed |
|---|
| 118 | |
|---|
| 119 | !! only consider updraft (usually >> downdraft flux) |
|---|
| 120 | !! Choose a random # evenly distributed in [0,1] |
|---|
| 121 | rn=ran3(88) |
|---|
| 122 | |
|---|
| 123 | if (rn .le. fmix) then ! inside cloud |
|---|
| 124 | write(*,*)'well mixed, fmix=,rn=',fmix,rn |
|---|
| 125 | ! --|-- umfzf(j), zf(j) |
|---|
| 126 | ! --|X- rn2 (particle position) |
|---|
| 127 | ! --|-- umfzf(j-1),zf(j-1) |
|---|
| 128 | |
|---|
| 129 | rn2=ran3(881) |
|---|
| 130 | |
|---|
| 131 | do j=2,nuvz+1 |
|---|
| 132 | if(umfzf(j) .ge. rn2) then |
|---|
| 133 | zp = zf(j)-(umfzf(j)-rn)/ & |
|---|
| 134 | (umfzf(j)-umfzf(j-1))*dz1d(j-1) |
|---|
| 135 | goto 92 |
|---|
| 136 | endif |
|---|
| 137 | enddo |
|---|
| 138 | 92 continue |
|---|
| 139 | write(*,*)'old position=',kk, & |
|---|
| 140 | 'repositioned at k=,umfzf(j),rn',j,umfzf(j),rn |
|---|
| 141 | endif ! inside cloud |
|---|
| 142 | ENDIF ! Well-mixed |
|---|
| 143 | |
|---|
| 144 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 145 | !!!!!!!!!!!!!!!! reposition based on probabilities of entrainment and detrainment |
|---|
| 146 | !!!!!!!!!!!!!!!! FORWARD treatment |
|---|
| 147 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 148 | IF (mix_option .eq. prob_mix .and. ldirect .gt. 0)then !! prob_mix+ forward |
|---|
| 149 | fup = abs(uer(kk))*delt/totalmass |
|---|
| 150 | fdown=abs(der(kk))*delt/totalmass |
|---|
| 151 | ftotal=fup+fdown ! ignore downward flux if 0*fdown |
|---|
| 152 | if (ftotal .le. 1e-10) goto 89 ! return if no cloud |
|---|
| 153 | ! pick a random number to see if particle is inside cloud |
|---|
| 154 | rn = ran3(88) |
|---|
| 155 | |
|---|
| 156 | |
|---|
| 157 | if (rn .le. ftotal) then ! in |
|---|
| 158 | rn1=ran3(881) |
|---|
| 159 | write(*,*)'inside cloud, kk,ftotal,rn=',kk,ftotal,rn |
|---|
| 160 | if (rn1 .le. fup/ftotal) then ! updraft |
|---|
| 161 | ! parcel will move upward till all particles inside are detrained |
|---|
| 162 | ! umf(k+1)=umf(k)+uer(k+1)-udr(k+1) |
|---|
| 163 | ! |
|---|
| 164 | ! ---|--- |
|---|
| 165 | ! ---|+++ |
|---|
| 166 | ! ---|--- zf(j+2) |
|---|
| 167 | ! ---|+++ udr(j+1), half level, partile is detrained here |
|---|
| 168 | ! ---|--- zf(j+1) |
|---|
| 169 | ! ........ |
|---|
| 170 | ! ---|+++ uer(kk), particle is entrained here |
|---|
| 171 | write(*,*)'updraft','fup=',fup,uer(kk)*delt,totalmass |
|---|
| 172 | write(*,*)'fdown=',fdown,fdown/ftotal |
|---|
| 173 | do j=kk,nuvz-1 |
|---|
| 174 | rn=ran3(882) |
|---|
| 175 | if (umf(j) .eq. 0.0) goto 98 |
|---|
| 176 | fde = abs(udr(j+1)/umf(j)) ! detrainment probability at level j+1 |
|---|
| 177 | if (rn .lt. fde .or. zh(j+1) .ge. uptop) goto 98 ! being detrtained from air parcel |
|---|
| 178 | ! or reach cloud top |
|---|
| 179 | enddo |
|---|
| 180 | 98 zp=zf(j+1) + rn*(zf(j+2)-zf(j+1)) ! particle is detrained between zf(j+1) and zf(j+2) |
|---|
| 181 | write(*,*)'detrain at',j+1,rn,zp,'old=',zf(kk) |
|---|
| 182 | else !up/down |
|---|
| 183 | ! move downward till detrained or get to ground |
|---|
| 184 | do j=kk,2,-1 |
|---|
| 185 | rn=ran3(883) |
|---|
| 186 | if(dmf(j) .eq. 0.0) goto 102 |
|---|
| 187 | fde = abs(ddr(j-1)/dmf(j)) ! detrainment at level j-1 |
|---|
| 188 | if(rn .lt. fde .or. j .eq. 2) goto 102 |
|---|
| 189 | enddo |
|---|
| 190 | 102 zp=zf(j-1)+rn*(zf(j)-zf(j-1)) |
|---|
| 191 | |
|---|
| 192 | write(*,*)'downward prob,detrain at',fdown/ftotal,j-1 |
|---|
| 193 | |
|---|
| 194 | endif ! downdraft |
|---|
| 195 | |
|---|
| 196 | else ! in /OUT |
|---|
| 197 | !! outside cloud |
|---|
| 198 | !! displace particle according to compensating subsidence |
|---|
| 199 | w_sub =(abs(umf(kk))-abs(dmf(kk)))/rho1d(kk)/dx/dx |
|---|
| 200 | zp = zp - w_sub*delt |
|---|
| 201 | write(*,*)'subsidence== distance(m)',w_sub*delt, & |
|---|
| 202 | 'rn,ftotal=',rn,ftotal |
|---|
| 203 | endif ! OUT |
|---|
| 204 | ENDIF !! prob_mix + forward |
|---|
| 205 | !!!!!!!!! END OF FORWARD treatment !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 206 | |
|---|
| 207 | !!!!!!!! ------------------------------------------------------------- |
|---|
| 208 | !!!!!!!! backward |
|---|
| 209 | !!!!!!!! reverse all processes, e.g., up --> down, down --> up, entrainment -> detrainment |
|---|
| 210 | !!!!!!!! detrainment --> entrainment,,,,,,, |
|---|
| 211 | !!!!!!!! ------------------------------------------------------------- |
|---|
| 212 | |
|---|
| 213 | IF (mix_option .eq. prob_mix .and. ldirect .lt. 0) then !! prob_mix + backward |
|---|
| 214 | fdown = abs(udr(kk))*delt/totalmass |
|---|
| 215 | fup = abs(ddr(kk))*delt/totalmass |
|---|
| 216 | ftotal=fup+fdown ! |
|---|
| 217 | ! ftotal=fup ! ignore downdraft |
|---|
| 218 | if (ftotal .le. 1e-10) goto 89 |
|---|
| 219 | |
|---|
| 220 | ! pick a random number to see if particle is inside cloud |
|---|
| 221 | rn = ran3(88) |
|---|
| 222 | if (rn .le. ftotal) then ! in |
|---|
| 223 | rn1=ran3(881) |
|---|
| 224 | if (rn1 .lt. fup/ftotal) then ! updraft |
|---|
| 225 | ! parcel will move "upward" till all particles inside being detrained |
|---|
| 226 | |
|---|
| 227 | do j=kk,nuvz-1 |
|---|
| 228 | rn=ran3(j) |
|---|
| 229 | if (dmf(j) .eq. 0.0) goto 981 |
|---|
| 230 | fde = abs(der(j+1)/dmf(j)) ! "detrainment" probability at level j+1 |
|---|
| 231 | if (rn .lt. fde .or. zh(j+1) .ge. downtop) goto 981 ! being detrtained from air parcel |
|---|
| 232 | ! or reach cloud top |
|---|
| 233 | enddo |
|---|
| 234 | 981 zp=zf(j+1) + rn*(zf(j+2)-zf(j+1)) ! partice is detrained between zf(j+1) and j+2 |
|---|
| 235 | |
|---|
| 236 | else !up/down |
|---|
| 237 | ! move downward till detrained or get to ground |
|---|
| 238 | do j=kk,2,-1 |
|---|
| 239 | rn=ran3(j) |
|---|
| 240 | if (umf(j) .eq. 0.0) goto 1021 |
|---|
| 241 | fde = abs(uer(j-1)/umf(j)) ! "detrainment" at level j-1 |
|---|
| 242 | if(rn .lt. fde .or. j .eq. 2) goto 1021 |
|---|
| 243 | enddo |
|---|
| 244 | 1021 zp=zf(j-1)+rn*(zf(j)-zf(j-1)) |
|---|
| 245 | |
|---|
| 246 | endif ! downdraft |
|---|
| 247 | |
|---|
| 248 | else ! in /OUT |
|---|
| 249 | !! outside cloud |
|---|
| 250 | !! displace particle according to compensating subsidence |
|---|
| 251 | w_sub =(abs(dmf(kk))-abs(umf(kk)))/rho1d(kk)/dx/dx |
|---|
| 252 | zp = zp - w_sub*delt |
|---|
| 253 | endif ! OUT |
|---|
| 254 | ENDIF !! prob_mix + backward |
|---|
| 255 | |
|---|
| 256 | !-------------- END OF BACKWARD treatment ------------------ |
|---|
| 257 | |
|---|
| 258 | |
|---|
| 259 | |
|---|
| 260 | |
|---|
| 261 | |
|---|
| 262 | !+++++++++ |
|---|
| 263 | ! END OF calculating particle position |
|---|
| 264 | |
|---|
| 265 | |
|---|
| 266 | 89 continue |
|---|
| 267 | |
|---|
| 268 | return |
|---|
| 269 | end subroutine redist_kf |
|---|
| 270 | |
|---|
| 271 | |
|---|
| 272 | |
|---|
