*
L. Bakels 2021: This module contains dry deposition related subroutines * *
!***************************************************************************** ! * ! L. Bakels 2021: This module contains dry deposition related subroutines * ! * !***************************************************************************** module wetdepo_mod use point_mod use par_mod use com_mod use particle_mod implicit none contains subroutine wetdepo(itime,ltsample,loutnext) ! i i i !***************************************************************************** ! * ! Calculation of wet deposition using the concept of scavenging coefficients.* ! For lack of detailed information, washout and rainout are jointly treated. * ! It is assumed that precipitation does not occur uniformly within the whole * ! grid cell, but that only a fraction of the grid cell experiences rainfall. * ! This fraction is parameterized from total cloud cover and rates of large * ! scale and convective precipitation. * ! * ! Author: A. Stohl * ! * ! 1 December 1996 * ! * ! Correction by Petra Seibert, Sept 2002: * ! use centred precipitation data for integration * ! Code may not be correct for decay of deposition! * ! * ! 2021 Andreas Plach: - moved backward wet depo. calc. here from timemanager * ! - bugfix in-cloud scavenging * ! * ! PS, AP 2021: followed up on some variable renaming and * ! corrected get_wetscav subroutine parameter list * !***************************************************************************** ! * ! Variables: * ! ix,jy indices of output grid cell for each particle * ! itime [s] actual simulation time [s] * ! jpart particle index * ! ldeltat [s] interval since radioactive decay was computed * ! loutnext [s] time for which gridded deposition is next output * ! loutstep [s] interval at which gridded deposition is output * ! ltsample [s] interval over which mass is deposited * ! wetdeposit mass that is wet deposited * ! wetgrid accumulated deposited mass on output grid * ! wetscav scavenging coefficient * ! * ! Constants: * ! * !***************************************************************************** #ifdef _OPENMP use omp_lib #endif use unc_mod implicit none integer :: i,jpart,itime,ltsample,loutnext,ldeltat integer :: itage,nage,inage,ithread,thread integer :: ks, kp, stat real :: gridfract,wetscav,restmass real,allocatable,dimension(:) :: wettmp real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled ! Compute interval since radioactive decay of deposited mass was computed !************************************************************************ if (itime.le.loutnext) then ldeltat=itime-(loutnext-loutstep) else ! first half of next interval ldeltat=itime-loutnext endif ! Loop over all particles !************************ #ifdef _OPENMP call omp_set_num_threads(numthreads_grid) #endif !$OMP PARALLEL PRIVATE(jpart,itage,nage,inage,ks,kp,thread,wetscav,wettmp, & !$OMP restmass, gridfract) #if (defined _OPENMP) thread = OMP_GET_THREAD_NUM() ! Starts with 0 #else thread = 0 #endif allocate( wettmp(nspec),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate wettmp inside of OMP loop' !$OMP DO do i=1,count%alive jpart=count%ialive(i) ! Determine age class of the particle - nage is used for the kernel !****************************************************************** itage=abs(itime-part(jpart)%tstart) nage=1 if (lagespectra.eq.1) then do inage=1,nageclass nage=inage if (itage.lt.lage(nage)) exit end do endif do ks=1,nspec ! loop over species if (.not. WETDEPSPEC(ks)) cycle !************************************************** ! CALCULATE DEPOSITION !************************************************** call get_wetscav(itime,jpart,ks,gridfract,wetscav) if (WETBKDEP) then if ((xscav_frac1(jpart,ks).lt.-0.1)) then ! particle marked as starting particle if (wetscav.gt.0.) then xscav_frac1(jpart,ks)=wetscav*(zpoint2(part(jpart)%npoint)-& zpoint1(part(jpart)%npoint))*gridfract else mass(jpart,ks)=0. xscav_frac1(jpart,ks)=0. endif endif endif if (wetscav.gt.0.) then wettmp(ks)=mass(jpart,ks)* & (1.-exp(-wetscav*abs(ltsample)))*gridfract ! wet deposition else ! if no scavenging wettmp(ks)=0. endif wetdeposit(jpart,ks)=wetdeposit(jpart,ks)+wettmp(ks) restmass = mass(jpart,ks)-wettmp(ks) if (ioutputforeachrelease.eq.1) then kp=part(jpart)%npoint else kp=1 endif if (restmass .gt. smallnum) then mass(jpart,ks)=restmass else mass(jpart,ks)=0. endif ! Correct deposited mass to the last time step when radioactive decay of ! gridded deposited mass was calculated if (decay(ks).gt.0.) then wettmp(ks)=wettmp(ks)*exp(abs(ldeltat)*decay(ks)) endif end do ! loop over species !************************************************************************ ! Sabine Eckhardt, June 2008 create deposition only for forward runs ! Add the wet deposition to accumulated amount on output grid ! and nested output grid if ((ldirect.eq.1).and.(iout.ne.0)) then !OMP reduction necessary for wetgridunc call wetdepokernel(part(jpart)%nclass,wettmp,real(part(jpart)%xlon), & real(part(jpart)%ylat),nage,kp,thread+1) if (nested_output.eq.1) call wetdepokernel_nest(part(jpart)%nclass, & wettmp,real(part(jpart)%xlon),real(part(jpart)%ylat),nage,kp,thread+1) endif end do ! all particles !$OMP END DO deallocate(wettmp) !$OMP END PARALLEL #ifdef _OPENMP call omp_set_num_threads(numthreads) #endif #ifdef _OPENMP if ((ldirect.eq.1).and.(iout.ne.0)) then do ithread=1,numthreads_grid wetgridunc(:,:,:,:,:,:)=wetgridunc(:,:,:,:,:,:)+gridunc_omp(:,:,1,:,:,:,:,ithread) gridunc_omp(:,:,1,:,:,:,:,ithread)=0. end do if (nested_output.eq.1) then do ithread=1,numthreads_grid wetgriduncn(:,:,:,:,:,:)=wetgriduncn(:,:,:,:,:,:)+griduncn_omp(:,:,1,:,:,:,:,ithread) griduncn_omp(:,:,1,:,:,:,:,ithread)=0. end do endif endif #endif end subroutine wetdepo subroutine get_wetscav(itime,jpart,ks,gridfract,wetscav) ! i i i i i o o !***************************************************************************** ! * ! Calculation of wet deposition using the concept of scavenging coefficients.* ! For lack of detailed information, washout and rainout are jointly treated. * ! It is assumed that precipitation does not occur uniformly within the whole * ! grid cell, but that only a fraction of the grid cell experiences rainfall. * ! This fraction is parameterized from total cloud cover and rates of large * ! scale and convective precipitation. * ! * ! Author: A. Stohl * ! * ! 1 December 1996 * ! * ! Correction by Petra Seibert, Sept 2002: * ! use centred precipitation data for integration * ! Code may not be correct for decay of deposition! * ! * ! ZHG, for v10: use below-cloud scavenging according to Laakso et al (2003) * ! and Kyro et al (2009) as described in Grytte et al (2017, GMD) * ! * ! PS, AP 04/2019: - put back temporal interpolation of rain, from v10.01 * ! - tansferred BCSCHEME parameters to par_mod.f90 * ! - added call to rain interpolation subroutine with new * ! interpolation for rain and all cloud params * ! - cleaned up scavenging determination algorithm * ! - added new below-cloud scavenging scheme * ! * !***************************************************************************** ! * ! Variables: * ! cc [0-1] total cloud cover * ! convp [mm/h] convective precipitation rate * ! grfraction [0-1] fraction of grid, for which precipitation occurs * ! ix,jy indices of output grid cell for each particle * ! itime [s] actual simulation time [s] * ! jpart particle index * ! lfr, cfr area fraction covered by precipitation for large scale * ! and convective precipitation (dependent on prec. rate) * ! loutnext [s] time for which gridded deposition is next output * ! loutstep [s] interval at which gridded deposition is output * ! lsp [mm/h] large-scale precipitation rate * ! ltsample [s] interval over which mass is deposited * ! prec [mm/h] precipitation rate in subgrid, where precipitation occurs* ! wetgrid accumulated deposited mass on output grid * ! wetscav scavenging coefficient * ! * ! Constants: * ! * !***************************************************************************** use interpol_mod use windfields_mod #ifdef ETA use coord_ecmwf_mod #endif implicit none integer,intent(in) :: jpart,itime,ks real,intent(out) :: gridfract,wetscav integer :: hz,interp_time, n,i,j,kz integer(kind=1) :: clouds_v integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count integer :: indcloud integer :: icbot,ictop real :: t_particle, si, cl, cle ! in cloud scavenging real :: lsp,convp,cc,prec !save lfr,cfr real :: xts,yts real, parameter :: precsub = 0.01 ! minimum precip rate (mm/h) real :: f logical :: readclouds_this_nest wetscav=0. ! Interpolate large-scale precipitation, convective precipitation, ! total cloud cover, particle temperature, cloud water content, ! cloud bottom and top !******************************************************************** interp_time=itime xts=real(part(jpart)%xlon) yts=real(part(jpart)%ylat) ! Determine which nesting level to be used !***************************************** readclouds_this_nest=.false. call find_ngrid(xts,yts) if ( (ngrid.gt.0) ) then if (lcw_nest(ngrid)) readclouds_this_nest=.true. endif ! If point at border of grid -> small displacement into grid !*********************************************************** if (ngrid.le.0) then if (xts.ge.real(nx-1)) xts=real(nx-1)-0.00001 if (yts.ge.real(ny-1)) yts=real(ny-1)-0.00001 else if (xts.ge.real(nx-1)) xts=real(nx-1)-0.00001 if (yts.ge.real(ny-1)) yts=real(ny-1)-0.00001 endif call find_grid_indices(xts,yts) call find_grid_distances(xts,yts) #ifdef ETA call update_zeta_to_z(itime,jpart) call find_z_level_eta_uv(real(part(jpart)%zeta)) kz=induv #else call find_z_level_meters(real(part(jpart)%z)) kz=indz #endif ! Interpolate cloud information call interpol_rain(itime,kz,lsp,convp,cc,t_particle,cl,icbot,ictop,icmv) ! cc = total cloud cover ! cl = ctwc ! If total precipitation is less than precsub=0.01 mm/h - no scavenging ! Note: original PS version (in order avoid step at 0.01) !----------------------------------------------------------------------- prec = lsp+convp if (prec .le. precsub) then return endif ! Remove the minimum 0.01 from the large scale precipitation (lsp) and ! convective precipitation (convp). Why 0.01??? f = (prec-precsub)/prec lsp = f*lsp convp = f*convp if (abs(memtime(1)-interp_time) .lt. abs(memtime(2)-interp_time)) then n=memind(1) else n=memind(2) endif ! if particle is above the clouds no scavenging is done !------------------------------------------------------ ! PS: part of 2011/2012 fix ! NOTE this is just for z coordinate ! Reverse sign for eta #ifdef ETA if (part(jpart)%zeta .gt. real(ictop)/eta_convert) then if (part(jpart)%zeta .le. real(icbot)/eta_convert) then #else if (part(jpart)%z .le. real(ictop)) then if (part(jpart)%z .gt. real(icbot)) then #endif indcloud = 2 ! in-cloud else indcloud = 1 ! below-cloud endif elseif (ictop .eq. icmv) then indcloud = 0 ! no cloud found, use old scheme else return ! above cloud endif ! 1) Parameterization of the the area fraction of the grid cell where the ! precipitation occurs: the absolute limit is the total cloud cover, but ! for low precipitation rates, an even smaller fraction of the grid cell ! is used. Large scale precipitation occurs over larger areas than ! convective precipitation. !************************************************************************** call get_gridfract(lsp,convp,cc,gridfract) ! 2) Computation of precipitation rate in sub-grid cell !****************************************************** prec=(lsp+convp)/gridfract ! 3) Computation of scavenging coefficients for all species !********************************************************** !------------------------------------------------------- if (indcloud .eq. 0) then ! NO CLOUD FOUND !------------------------------------------------------- ! Note: more complex formulation using H or particle diametre ! may be introduced later wetscav=wet_a*prec**wet_b !------------------------------------------------------- else if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING !------------------------------------------------------- if (dquer(ks).le.0. .and. & weta_gas(ks).gt.0. .and. wetb_gas(ks).gt.0.) then ! gas: if positive below-cloud parameters (A or B), and dquer<=0 call get_wetscav_belowcld_gas(ks,prec,wetscav) else if (dquer(ks).gt.0. .and. & (crain_aero(ks).gt.0. .or. csnow_aero(ks).gt.0.)) then ! aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0 if (t_particle .ge. 273. .and. crain_aero(ks).gt.0.) then ! Rain: call get_wetscav_belowcld_aerosol_rain(ks,prec,wetscav) else if (t_particle .lt. 273. .and. csnow_aero(ks).gt.0.) then ! Snow: call get_wetscav_belowcld_aerosol_snow(ks,prec,wetscav) !else ???????? endif endif ! gas or particle ! positive below-cloud scavenging parameters given in Species file ! end below cloud scavening !--------------------------------------------------------- elseif (indcloud .eq. 2) then ! IN CLOUD SCAVENGING !--------------------------------------------------------- ! if negative coefficients (turned off) set to zero for use in equation if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0. if (in_aero(ks).lt.0.) in_aero(ks)=0. if (dquer(ks).gt.0) then !aerosol call get_wetscav_incld_aerosol(ks,gridfract,prec,cl,cc,t_particle,wetscav) else !gas call get_wetscav_incld_gas(ks,gridfract,prec,cl,cc,t_particle,wetscav) endif !--------------------------------------------------------- endif ! incloud !--------------------------------------------------------- end subroutine get_wetscav subroutine get_wetscav_belowcld_gas(ks,prec,wetscav) implicit none integer,intent(in) :: ks ! Species index real, intent(in) :: prec ! precipitation in sub-grid cell real, intent(out) :: wetscav ! scavenging coefficient ! Sum number of particles below the cloud icnt_belowcld(ks)=icnt_belowcld(ks)+1 !weta_gas and wetb_gas are set in the SPECIES file wetscav=weta_gas(ks)*prec**wetb_gas(ks) end subroutine get_wetscav_belowcld_gas subroutine get_wetscav_belowcld_aerosol_rain(ks,prec,wetscav) implicit none integer,intent(in) :: ks ! Species index real, intent(in) :: prec ! precipitation in sub-grid cell real, intent(out) :: wetscav ! scavenging coefficient real :: dquer_m, ldquer real :: wetscavlim, logAd, B ! Sum number of particles below the cloud icnt_belowcld(ks)=icnt_belowcld(ks)+1 ! NIK 17.02.2015: local conversion particle diameter from um ! (SPECIES file, readspecies) to meter dquer_m = dquer(ks)*1.e-6 ! The parameterizations used by HG scheme are valid only for d < 10 um ! Therefore, locally the diametre is clipped. However, settling and dry dep ! are still calculated with the original diametre ! TODO check whether warning is written by readrelease for d > 10 um dquer_m = min( 10.e-6, dquer_m ) ! PS note that solid precip is usually expected for T<Tf, not T<273 ! also, if freezing/melting point is desired, why not 273.2? ! AT parameterization after WANG ET AL 2014 ! unit of dquer is in um ! unit of precip is in mm/h ! Wang et al. 2014: eq 6+7 ldquer = log10(dquer(ks)) if (dquer(ks) .le. 2.) then logAd = bclr_a(1) + & bclr_a(2) * ldquer + & bclr_a(3) * ldquer**2. + & bclr_a(4) * ldquer**3. B = bclr_c(1) + bclr_c(2)*ldquer else ! dquer .gt. 2. logAd = bclr_b(1) + & bclr_b(2) * ldquer + & bclr_b(3) * ldquer**2. + & bclr_b(4) * ldquer**3. + & bclr_b(5) * ldquer**4. + & bclr_b(6) * ldquer**5. + & bclr_b(7) * ldquer**6. B = bclr_e(1) + & bclr_e(2) * ldquer + & bclr_e(3) * ldquer**2. + & bclr_e(4) * ldquer**3. + & bclr_e(5) * ldquer**4. + & bclr_e(6) * ldquer**5. + & bclr_e(7) * ldquer**6. endif ! dquer ! Wang et al. 2014: eq. 4 wetscav = 10**(logAd+B*log10(prec)) end subroutine get_wetscav_belowcld_aerosol_rain subroutine get_wetscav_belowcld_aerosol_snow(ks,prec,wetscav) implicit none integer,intent(in) :: ks ! Species index real, intent(in) :: prec ! precipitation in sub-grid cell real, intent(out) :: wetscav ! scavenging coefficient real :: dquer_m, ldquer real :: wetscavlim, logAd, B ! Sum number of particles below the cloud icnt_belowcld(ks)=icnt_belowcld(ks)+1 ! NIK 17.02.2015: local conversion particle diameter from um ! (SPECIES file, readspecies) to meter dquer_m = dquer(ks)*1.e-6 ! The parameterizations used by HG scheme are valid only for d < 10 um ! Therefore, locally the diametre is clipped. However, settling and dry dep ! are still calculated with the original diametre ! TODO check whether warning is written by readrelease for d > 10 um dquer_m = min( 10.e-6, dquer_m ) ! AT parameterization after WANG ET AL 2014 ! unit of dquer is in um ! unit of precip is in mm/h ldquer = log10(dquer(ks)) ! Wang et al. 2014: eq. 8+9 if (dquer(ks) .le. 1.44) then logAd = bcls_a(1) + & bcls_a(2) * ldquer + & bcls_a(3) * ldquer**2. + & bcls_a(4) * ldquer**3. + & bcls_a(5) * ldquer**4. + & bcls_a(6) * ldquer**5. + & bcls_a(7) * ldquer**6. B = bcls_c(1) + & bcls_c(2) * ldquer + & bcls_c(3) * ldquer**2. + & bcls_c(4) * ldquer**3. + & bcls_c(5) * ldquer**4. + & bcls_c(6) * ldquer**5. + & bcls_c(7) * ldquer**6. else ! dquer .gt. 1.44 logAd = bcls_b(1) + & bcls_b(2) * ldquer + & bcls_b(3) * ldquer**2. + & bcls_b(4) * ldquer**3. + & bcls_b(5) * ldquer**4. + & bcls_b(6) * ldquer**5. + & bcls_b(7) * ldquer**6. B = bcls_e(1) + & bcls_e(2) * ldquer + & bcls_e(3) * ldquer**2. + & bcls_e(4) * ldquer**3. + & bcls_e(5) * ldquer**4. + & bcls_e(6) * ldquer**5. + & bcls_e(7) * ldquer**6. endif ! dquer ! Wang et al. 2014: eq. 4 wetscav = 10**(logAd+B*log10(prec)) end subroutine get_wetscav_belowcld_aerosol_snow subroutine get_wetscav_incld_aerosol(ks,gridfract,prec,cl,cc,t_particle,wetscav) implicit none integer,intent(in) :: ks real,intent(in) :: gridfract real,intent(in) :: t_particle ! temperature real, intent(in) :: prec,cc ! precipitation in sub-grid cell real, intent(inout) :: cl ! scavenging coefficient real,intent(out) :: wetscav real :: frac_act,si ! NIK 13 may 2015: do only if in-cloud aerosol scavenging parameters > 0 ! (all defined in SPECIES) if (ccn_aero(ks)+in_aero(ks) .le. 0.) return icnt_incld(ks)=icnt_incld(ks)+1 ! Compute cloud liquid and ice water call get_cloud_liquid(gridfract,prec,cc,cl) ! Compute actived fraction based on the temperature (rain vs. snow ) call get_activated_frac(ks,t_particle, frac_act) !ZHG Use the activated fraction and the liqid water to calculate the washout si= frac_act/cl ! scavenging coefficient based on Hertel et al 1995 - ! using si (S_i in paper) for both gas and aerosol ! wetscav = ratio_incloud*si*prec/3.6E6/cloud_thickness ! cloud_thickness cancels out since cl is computed without wetscav=ratio_incloud*si* prec/3.6E6 ! mm/h -> m/s end subroutine get_wetscav_incld_aerosol subroutine get_wetscav_incld_gas(ks,gridfract,prec,cl,cc,t_particle,wetscav) implicit none integer,intent(in) :: ks real,intent(in) :: gridfract real,intent(in) :: t_particle ! temperature real, intent(in) :: prec,cc ! precipitation in sub-grid cell real, intent(inout) :: cl ! scavenging coefficient real,intent(out) :: wetscav real :: frac_act,si,cle ! NIK 13 may 2015: do only if in-cloud aerosol scavenging parameters > 0 ! and Henry's constant > 0 (all defined in SPECIES) if (ccn_aero(ks)+in_aero(ks)+henry(ks) .le. 0.) return icnt_incld(ks)=icnt_incld(ks)+1 ! Compute cloud liquid and ice water call get_cloud_liquid(gridfract,prec,cc,cl) ! Compute actived fraction based on the temperature (rain vs. snow ) call get_activated_frac(ks,t_particle, frac_act) !ZHG Use the activated fraction and the liqid water to calculate the washout cle=(1.-cl)/(henry(ks)*(r_air/3500.)*t_particle) + cl si=1./cle ! scavenging coefficient based on Hertel et al 1995 - ! using si (S_i in paper) for both gas and aerosol ! wetscav = ratio_incloud*si*prec/3.6E6/cloud_thickness ! cloud_thickness cancels out since cl is computed without wetscav=ratio_incloud*si* prec/3.6E6 ! mm/h -> m/s end subroutine get_wetscav_incld_gas subroutine get_activated_frac(ks,t_particle,frac_act) implicit none integer, intent(in) :: ks real, intent(in) :: t_particle real, intent(out) :: frac_act real :: frac_liq, frac_ice ! AT use of correct Kelvin temperature for T_ice; after ECMWF if (t_particle .le. 250.16) then ! ice frac_liq=0. frac_ice=1. ! AT use of correct Kelvin temperature for T_liquid; after ECMWF elseif (t_particle .ge. 273.16) then ! liquid frac_liq=1. frac_ice=0. else ! mixed cloud ! Use exact value for the melting point and ice threshold as in ECMWF/IFS frac_liq= ((t_particle-250.16)/(273.16-250.16))**2. frac_ice = max(0., 1. - frac_liq) endif frac_act = frac_liq*ccn_aero(ks) + frac_ice*in_aero(ks) end subroutine get_activated_frac subroutine get_gridfract(lsp,convp,cc,gridfract) ! 1) Parameterization of the the area fraction of the grid cell where the ! precipitation occurs: the absolute limit is the total cloud cover, but ! for low precipitation rates, an even smaller fraction of the grid cell ! is used. Large scale precipitation occurs over larger areas than ! convective precipitation. !************************************************************************** implicit none real,intent(in) :: lsp,cc,convp real, intent(out) :: gridfract ! Where do these numbers come from? real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/) real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /) integer :: i, j if (.not. lgridfraction) then gridfract=1.0 return endif if (lsp.gt.20.) then i=5 else if (lsp.gt.8.) then i=4 else if (lsp.gt.3.) then i=3 else if (lsp.gt.1.) then i=2 else i=1 endif if (convp.gt.20.) then j=5 else if (convp.gt.8.) then j=4 else if (convp.gt.3.) then j=3 else if (convp.gt.1.) then j=2 else j=1 endif ! In the future, we may differentiate the gridfract for lsp and convp ! for now they are still mixed gridfract=max( 0.05, cc* (lsp*lfr(i) + convp*cfr(j)) / (lsp+convp)) end subroutine get_gridfract subroutine get_cloud_liquid(gridfract,prec,cc,cl) use interpol_mod implicit none real, intent(in) :: prec,cc,gridfract ! precipitation in sub-grid cell real, intent(inout) :: cl ! scavenging coefficient !ZHG 2015 use cloud liquid & ice water (CLWC+CIWC) from ECMWF ! MC -- Integrated water content: ! CTWC = SUM(CLWC * rho_water * cloud_thickness) [kg/kg * kg/m3 * m] ! -> Average water content: cl = CTWC/rho_water/cloud_thickness [m3(water)/m3(cloud)] ! Note that cloud_thickness is not included, since this will cancel out when ! computing the wetscavenging: Wetscav=ratio_incloud*S_i*(prec/3.6E6)/cloud_thickness ! S_i=1/cl ! Mother grid if (ngrid.le.0) then if (lcw) then if (lgridfraction) then cl=cl/rho_water*(gridfract/cc) else cl=cl/rho_water ! Grythe et al. eq (1), no cloud_thickness since this will cancel out later endif ! cl = cl*(gridfract/cc) ! A.Plach 2021 cl should not become too small ! cl=max(0.2*prec**0.36, cl*(gridfract/cc)) else ! no cloud water available, use parameterisation for cloud water [m2/m3] cl=0.2*prec**0.36 !max(0.2*prec**0.36, cl*(gridfract/cc)) ! ZHG updated parameterization to better reproduce the values from ECMWF ! cl = 1.E6*2E-7*prec**0.36 ! SEC ECMWF new, is also suitable for GFS ! SEC test: ! cl=1.E6*1.E-7*prec**0.3 ! SEC GFS new ! cl= 2.E-7*prec**0.36 ! Andreas ! cl= 1.6E-6*prec**0.36 ! Henrik endif else if (lcw_nest(ngrid)) then if (lgridfraction) then cl=cl/rho_water*(gridfract/cc) else cl=cl/rho_water ! Grythe et al. eq (1), no cloud_thickness since this will cancel out later endif !cl = cl*(gridfract/cc) else cl=0.2*prec**0.36 endif endif end subroutine get_cloud_liquid subroutine wetdepokernel(nunc,deposit,x,y,nage,kp,thread) ! i i i i i !***************************************************************************** ! * ! Attribution of the deposition from an individual particle to the * ! deposition fields using a uniform kernel with bandwidths dxout and dyout.* ! * ! Author: A. Stohl * ! * ! 26 December 1996 * ! * !***************************************************************************** ! * ! Variables: * ! * ! nunc uncertainty class of the respective particle * ! nage age class of the respective particle * ! deposit amount (kg) to be deposited * ! * !***************************************************************************** ! Changes: ! eso 10/2016: Added option to disregard kernel ! !***************************************************************************** use unc_mod implicit none integer,intent(in) :: thread real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp xl=(x*dx+xoutshift)/dxout yl=(y*dy+youtshift)/dyout ix=int(xl) jy=int(yl) ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! If no kernel is used, direct attribution to grid cell !****************************************************** if (.not.lusekerneloutput) then do ks=1,nspec if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & (jy.le.numygrid-1)) then #ifdef _OPENMP gridunc_omp(ix,jy,1,ks,kp,nunc,nage,thread)= & gridunc_omp(ix,jy,1,ks,kp,nunc,nage,thread)+deposit(ks) #else wetgridunc(ix,jy,ks,kp,nunc,nage)= & wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks) #endif end if end do else ! use kernel ! Determine mass fractions for four grid points !********************************************** do ks=1,nspec if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & (jy.le.numygrid-1)) then w=wx*wy #ifdef _OPENMP gridunc_omp(ix,jy,1,ks,kp,nunc,nage,thread)= & gridunc_omp(ix,jy,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgridunc(ix,jy,ks,kp,nunc,nage)= & wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. & (jyp.le.numygrid-1)) then w=(1.-wx)*(1.-wy) #ifdef _OPENMP gridunc_omp(ixp,jyp,1,ks,kp,nunc,nage,thread)= & gridunc_omp(ixp,jyp,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgridunc(ixp,jyp,ks,kp,nunc,nage)= & wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. & (jy.le.numygrid-1)) then w=(1.-wx)*wy #ifdef _OPENMP gridunc_omp(ixp,jy,1,ks,kp,nunc,nage,thread)= & gridunc_omp(ixp,jy,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgridunc(ixp,jy,ks,kp,nunc,nage)= & wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. & (jyp.le.numygrid-1)) then w=wx*(1.-wy) #ifdef _OPENMP gridunc_omp(ix,jyp,1,ks,kp,nunc,nage,thread)= & gridunc_omp(ix,jyp,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgridunc(ix,jyp,ks,kp,nunc,nage)= & wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w #endif endif end do end if end subroutine wetdepokernel subroutine wetdepokernel_nest(nunc,deposit,x,y,nage,kp,thread) ! i i i i i i !***************************************************************************** ! * ! Attribution of the deposition from an individual particle to the * ! nested deposition fields using a uniform kernel with bandwidths * ! dxoutn and dyoutn. * ! * ! Author: A. Stohl * ! * ! 26 December 1996 * ! * ! 2 September 2004: Adaptation from wetdepokernel. * ! * ! * !***************************************************************************** ! * ! Variables: * ! * ! nunc uncertainty class of the respective particle * ! nage age class of the respective particle * ! deposit amount (kg) to be deposited * ! * !***************************************************************************** use unc_mod implicit none integer,intent(in) :: thread real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage xl=(x*dx+xoutshiftn)/dxoutn yl=(y*dy+youtshiftn)/dyoutn ! old: ! ix=int(xl) ! jy=int(yl) ! ESO: for xl,yl in range <-.5,-1> we get ix,jy=0 and thus negative ! wx,wy as the int function rounds upwards for negative numbers. ! Either use the floor function, or (perhaps more correctly?) use "(xl.gt.-0.5)" ! in place of "(ix.ge.0)" and similar for the upper boundary. ! new: ix=floor(xl) jy=floor(yl) ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! Determine mass fractions for four grid points !********************************************** do ks=1,nspec if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & (jy.le.numygridn-1)) then w=wx*wy #ifdef _OPENMP griduncn_omp(ix,jy,1,ks,kp,nunc,nage,thread)= & griduncn_omp(ix,jy,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgriduncn(ix,jy,ks,kp,nunc,nage)= & wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. & (jyp.le.numygridn-1)) then w=(1.-wx)*(1.-wy) #ifdef _OPENMP griduncn_omp(ixp,jyp,1,ks,kp,nunc,nage,thread)= & griduncn_omp(ixp,jyp,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= & wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. & (jy.le.numygridn-1)) then w=(1.-wx)*wy #ifdef _OPENMP griduncn_omp(ixp,jy,1,ks,kp,nunc,nage,thread)= & griduncn_omp(ixp,jy,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgriduncn(ixp,jy,ks,kp,nunc,nage)= & wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w #endif endif if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. & (jyp.le.numygridn-1)) then w=wx*(1.-wy) #ifdef _OPENMP griduncn_omp(ix,jyp,1,ks,kp,nunc,nage,thread)= & griduncn_omp(ix,jyp,1,ks,kp,nunc,nage,thread)+deposit(ks)*w #else wetgriduncn(ix,jyp,ks,kp,nunc,nage)= & wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w #endif endif end do end subroutine wetdepokernel_nest subroutine writeprecip(itime,imem) !***************************************************************************** ! * ! This routine produces a file containing total precipitation for each * ! releases point. * ! * ! Author: S. Eckhardt * ! 7 Mai 2017 * !***************************************************************************** use point_mod use par_mod use com_mod use date_mod use windfields_mod implicit none integer :: jjjjmmdd,ihmmss,itime,i real(kind=dp) :: jul character :: adate*8,atime*6 integer :: ix,jy,imem real :: xp1,yp1 if (itime.eq.0) then open(unitprecip,file=path(2)(1:length(2))//'wetscav_precip.txt', & form='formatted',err=998) else open(unitprecip,file=path(2)(1:length(2))//'wetscav_precip.txt', & ACCESS='APPEND',form='formatted',err=998) endif jul=bdate+real(itime,kind=dp)/86400._dp call caldate(jul,jjjjmmdd,ihmmss) write(adate,'(i8.8)') jjjjmmdd write(atime,'(i6.6)') ihmmss do i=1,numpoint xp1=xpoint1(i)*dx+xlon0 !lat, long (real) coord yp1=ypoint1(i)*dy+ylat0 !lat, long (real) coord ix=int((xpoint1(i)+xpoint2(i))*0.5) jy=int((ypoint1(i)+ypoint2(i))*0.5) write(unitprecip,*) jjjjmmdd, ihmmss,xp1,yp1, & sum(lsprec(ix,jy,1,:,imem)),sum(convprec(ix,jy,1,:,imem)) !time is the same as in the ECMWF windfield ! units mm/h, valid for the time given in the windfield enddo close(unitprecip) return 998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### ' write(*,*) ' #### '//path(2)(1:length(2))//'header_txt'//' #### ' write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### ' write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### ' write(*,*) ' #### THE PROGRAM AGAIN. #### ' error stop end subroutine writeprecip end module wetdepo_mod