!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** subroutine redist (ipart,ktop,ipconv) !************************************************************************** ! Do the redistribution of particles due to convection ! This subroutine is called for each particle which is assigned ! a new vertical position randomly, based on the convective redistribution ! matrix !************************************************************************** ! Petra Seibert, Feb 2001, Apr 2001, May 2001, Jan 2002, Nov 2002 and ! Andreas Frank, Nov 2002 ! Caroline Forster: November 2004 - February 2005 use par_mod use com_mod use conv_mod implicit none real,parameter :: const=r_air/ga integer :: ipart, ktop,ipconv integer :: k, kz, levnew, levold real,save :: uvzlev(nuvzmax) real :: wsub(nuvzmax) real :: totlevmass, wsubpart real :: temp_levold,temp_levold1 real :: sub_levold,sub_levold1 real :: pint, pold, rn, tv, tvold, dlevfrac real :: ew,ran3, ztold,ffraction real :: tv1, tv2, dlogp, dz, dz1, dz2 integer :: iseed = -88 ! ipart ... number of particle to be treated ipconv=1 ! determine height of the eta half-levels (uvzlev) ! do that only once for each grid column ! i.e. when ktop.eq.1 !************************************************************** if (ktop .le. 1) then tvold=tt2conv*(1.+0.378*ew(td2conv)/psconv) pold=psconv uvzlev(1)=0. pint = phconv(2) ! determine next virtual temperatures tv1 = tconv(1)*(1.+0.608*qconv(1)) tv2 = tconv(2)*(1.+0.608*qconv(2)) ! interpolate virtual temperature to half-level tv = tv1 + (tv2-tv1)*(pconv(1)-phconv(2))/(pconv(1)-pconv(2)) if (abs(tv-tvold).gt.0.2) then uvzlev(2) = uvzlev(1) + & const*log(pold/pint)* & (tv-tvold)/log(tv/tvold) else uvzlev(2) = uvzlev(1)+ & const*log(pold/pint)*tv endif tvold=tv tv1=tv2 pold=pint ! integrate profile (calculation of height agl of eta layers) as required do kz = 3, nconvtop+1 ! note that variables defined in calcmatrix.f (pconv,tconv,qconv) ! start at the first real ECMWF model level whereas kz and ! thus uvzlev(kz) starts at the surface. uvzlev is defined at the ! half-levels (between the tconv, qconv etc. values !) ! Thus, uvzlev(kz) is the lower boundary of the tconv(kz) cell. pint = phconv(kz) ! determine next virtual temperatures tv2 = tconv(kz)*(1.+0.608*qconv(kz)) ! interpolate virtual temperature to half-level tv = tv1 + (tv2-tv1)*(pconv(kz-1)-phconv(kz))/ & (pconv(kz-1)-pconv(kz)) if (abs(tv-tvold).gt.0.2) then uvzlev(kz) = uvzlev(kz-1) + & const*log(pold/pint)* & (tv-tvold)/log(tv/tvold) else uvzlev(kz) = uvzlev(kz-1)+ & const*log(pold/pint)*tv endif tvold=tv tv1=tv2 pold=pint end do ktop = 2 endif ! (if ktop .le. 1) then ! determine vertical grid position of particle in the eta system !**************************************************************** ztold = ztra1(abs(ipart)) ! find old particle grid position do kz = 2, nconvtop if (uvzlev(kz) .ge. ztold ) then levold = kz-1 goto 30 endif end do ! Particle is above the potentially convective domain. Skip it. goto 90 30 continue ! now redistribute particles !**************************** ! Choose a random number and find corresponding level of destination ! Random numbers to be evenly distributed in [0,1] rn = ran3(iseed) ! initialize levnew levnew = levold ffraction = 0. totlevmass=dpr(levold)/ga do k = 1,nconvtop ! for backward runs use the transposed matrix if (ldirect.eq.1) then ffraction=ffraction+fmassfrac(levold,k) & /totlevmass else ffraction=ffraction+fmassfrac(k,levold) & /totlevmass endif if (rn.le.ffraction) then levnew=k ! avoid division by zero or a too small number ! if division by zero or a too small number happens the ! particle is assigned to the center of the grid cell if (ffraction.gt.1.e-20) then if (ldirect.eq.1) then dlevfrac = (ffraction-rn) / fmassfrac(levold,k) * totlevmass else dlevfrac = (ffraction-rn) / fmassfrac(k,levold) * totlevmass endif else dlevfrac = 0.5 endif goto 40 endif end do 40 continue ! now assign new position to particle if (levnew.le.nconvtop) then if (levnew.eq.levold) then ztra1(abs(ipart)) = ztold else dlogp = (1.-dlevfrac)* & (log(phconv(levnew+1))-log(phconv(levnew))) pint = log(phconv(levnew))+dlogp dz1 = pint - log(phconv(levnew)) dz2 = log(phconv(levnew+1)) - pint dz = dz1 + dz2 ztra1(abs(ipart)) = (uvzlev(levnew)*dz2+uvzlev(levnew+1)*dz1)/dz if (ztra1(abs(ipart)).lt.0.) & ztra1(abs(ipart))=-1.*ztra1(abs(ipart)) if (ipconv.gt.0) ipconv=-1 endif endif ! displace particle according to compensating subsidence ! this is done to those particles, that were not redistributed ! by the matrix !************************************************************** if (levnew.le.nconvtop.and.levnew.eq.levold) then ztold = ztra1(abs(ipart)) ! determine compensating vertical velocity at the levels ! above and below the particel position ! increase compensating subsidence by the fraction that ! is displaced by convection to this level if (levold.gt.1) then temp_levold = tconv(levold-1) + & (tconv(levold)-tconv(levold-1)) & *(pconv(levold-1)-phconv(levold))/ & (pconv(levold-1)-pconv(levold)) sub_levold = sub(levold)/(1.-sub(levold)/dpr(levold)*ga) wsub(levold)=-1.*sub_levold*r_air*temp_levold/(phconv(levold)) else wsub(levold)=0. endif temp_levold1 = tconv(levold) + & (tconv(levold+1)-tconv(levold)) & *(pconv(levold)-phconv(levold+1))/ & (pconv(levold)-pconv(levold+1)) sub_levold1 = sub(levold+1)/(1.-sub(levold+1)/dpr(levold+1)*ga) wsub(levold+1)=-1.*sub_levold1*r_air*temp_levold1/ & (phconv(levold+1)) ! interpolate wsub to the vertical particle position dz1 = ztold - uvzlev(levold) dz2 = uvzlev(levold+1) - ztold dz = dz1 + dz2 wsubpart = (dz2*wsub(levold)+dz1*wsub(levold+1))/dz ztra1(abs(ipart)) = ztold+wsubpart*real(lsynctime) if (ztra1(abs(ipart)).lt.0.) then ztra1(abs(ipart))=-1.*ztra1(abs(ipart)) endif endif !(levnew.le.nconvtop.and.levnew.eq.levold) ! Maximum altitude .5 meter below uppermost model level !******************************************************* 90 continue if (ztra1(abs(ipart)) .gt. height(nz)-0.5) & ztra1(abs(ipart)) = height(nz)-0.5 end subroutine redist