! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" #include "output.inc" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: WET_DEPOSITION ! ! !DESCRIPTION: holds convective precipitation (CP) and large scale ! precipitation (LSP) budgets variables. ! ! has all routines to deal with LSP, since it is done here. CP ! is done in convection. ! ! **** THIS VERSION DIFFERS FROM THE BASE IN TWO THRESHOLD VALUES **** ! ! (1) Scale factor relative to 100% scavenging (of HNO3) for scavenging ! type = 2 (ie variable factor, using henry solubility) is non-zero and ! computed if henry coeff > 1 (instead of 10 in base code). ! ! (2) Wet removal rates: in case of in-cloud scavenging, test on cloud ! cover has a threshold of 0.02 instead of 0.05 !\\ !\\ ! !INTERFACE: ! MODULE WET_DEPOSITION ! ! !USES: ! use GO, only : gol, goErr, goPr use dims, only : nregions use chem_param, only : ntracet use global_types, only : d3_data use tm5_distgrid, only : dgrid, Get_DistGrid, reduce, gather use ParTools, only : isRoot IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Wet_Deposition_Init, Wet_Deposition_Done public :: calc_cvsfac, calculate_lsp_scav, lspscav ! ! !PUBLIC DATA MEMBERS: ! real, public :: cvsfac(ntracet) = 0.0 ! scavenging efficiencies, used in convection real, public :: cp_scale = 0.5 ! default amount of rain (converted to m/s) with maximum CP removal on 1x1 degree #ifdef with_budgets real, dimension(nregions), public :: sum_wetdep ! global wet dep budget for tracer #1 (includes both LSP and CP; meaningful on root only) type, public :: buddep_data ! budgets in each column, split vertically... real,dimension(:,:,:,:),pointer :: lsp ! (i, j, nbud_vg, ntracet) computed here real,dimension(:,:,:,:),pointer :: cp ! (i, j, nbud_vg, ntracet) computed in convection end type buddep_data type(buddep_data), dimension(nregions), target, public :: buddep_dat ! ... for each region #endif ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'wet_deposition' ! ! Large scale scavenging coefficients [s-1] type(d3_data),dimension(nregions) :: rloss1 ! 1: incloud completely soluble gas !>>>TvN type(d3_data),dimension(nregions) :: rloss1_m7 ! as 1, but with ice as efficient as water ! now needed for M7 aerosols !<<ntot_ed) is used ! ! ! used from chem_param: ! ! nscav : selected species for scavenging ! nscav_index : index for scavenging: ! nscav_type : type of scavenging: ! 0 no scavenging ! 1 scavenging 100 % solubility assumed ! 2 scavenging henry solubility assumed ! 3 scavenging, bulk aerosol removal assumed ! 4 scavenging, special case for SO2 with aq phase diss. ! 5-11 scavenging, M7 aerosol removal ! !---------------------------------------------- ! acidity needed for explicit calculation of reactive removal of SO2. ! Parameterisation calculates reaction of SO2 with H2O2 and O3. ! Not yet implemented: for information Frank Dentener ! see routine wetS ! nacid : selected species for acidity ! nacid_index : indices of species for acidity : iso4,imsa,ihno3,inh3,inh4 ! nacid_eq : equivalents acid per mole ! integer,parameter :: nacid=5 ! integer,parameter,dimension(nacid) :: & ! nacid_index = (/iso4,imsa,ihno3,ino3_a,inh3,inh4/) ! integer,parameter,dimension(nacid) :: & ! nacid_eq = (/ 2, 1, 1, 1 , -1, -1/) !---------------------------------------------- ! ! !REVISION HISTORY: ! version March 2003, adapted for TM5 MK from KNMI version ! 29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: WET_DEPOSITION_INIT ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE WET_DEPOSITION_INIT( status ) ! ! !USES: ! use GO, only : TrcFile, Init, Done, ReadRc use dims, only : lmax_conv use global_data, only : rcfile use meteodata, only : Set, temper_dat, lwc_dat, iwc_dat, cc_dat, lsp_dat #ifdef with_budgets use budget_global, only : nbud_vg #endif use chem_param, only : ntrace ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Wet_Deposition_Init' integer :: region, imr,jmr,lmr, i1, i2, j1, j2 type(TrcFile) :: rcF ! --- begin ------------------------------------ ! select meteo to be used do region = 1, nregions call Set( temper_dat(region), status, used=.true. ) call Set( lwc_dat(region), status, used=.true. ) call Set( iwc_dat(region), status, used=.true. ) call Set( cc_dat(region), status, used=.true. ) call Set( lsp_dat(region), status, used=.true. ) end do ! allocate do region = 1, nregions lmr = lmax_conv CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate(rloss1(region)%d3(i1:i2, j1:j2, lmr)) !>>>TvN allocate(rloss1_m7(region)%d3(i1:i2, j1:j2, lmr)) !<<HSO3(-)<--->SO3(2-) ! for convective clouds T=273K and pH=5 is chosen ! 5-11) M7 modes cvsfac set equal to the scavenging parameters from Stier et al. (2005) ! for convective in-cloud scavenging !\\ !\\ ! !INTERFACE: ! SUBROUTINE CALC_CVSFAC ! ! !USES: ! use chem_param, only : nscav, nscav_index, nscav_type use chem_param, only : henry, ntlow, ntemp use chem_param, only : names ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC integer :: iscav,n,k real :: rtl, heff cvsfac=0.0 do iscav=1,nscav n=nscav_index(iscav) ! skip dummy tracers .. if ( n < 0 ) cycle ! fill cvsfac given scavenging type: select case(nscav_type(iscav)) case(0) cvsfac(n) = 0.0 !>>>TvN !case(1,3) case(1) !<<< TvN cvsfac(n) = 1.0 case(2) rtl=8.3148e-8*273. ! phase factor: ratio of aq. phase conc. to total conc. ! assuming LWC of 1e-6 k=nint(273.-float(ntlow)) k = min(max(1,k),ntemp) if ( henry(n,k) > 1. ) then cvsfac(n) = rtl*henry(n,k)/(1.+rtl*henry(n,k)) else cvsfac(n) = 0.0 end if !>>>TvN case(3) cvsfac(n) = 0.99 !<<>>TvN case(5) ! soluble nu !cvsfac(n) = 1.0 cvsfac(n) = 0.2 case(6) ! soluble ai !cvsfac(n) = 1.0 cvsfac(n) = 0.6 case(7) ! soluble ac !cvsfac(n) = 1.0 cvsfac(n) = 0.99 case(8) ! soluble co !cvsfac(n) = 1.0 cvsfac(n) = 0.99 case(9) ! insoluble ai !cvsfac(n) = 1.0 cvsfac(n) = 0.2 case(10) ! insoluble ac !cvsfac(n) = 1.0 cvsfac(n) = 0.4 case(11) ! insoluble co !cvsfac(n) = 1.0 cvsfac(n) = 0.4 !<< 0.0 ) then write (gol,'(" calc_cvsfac: Scavenging factor species ",i3,x,a,"; factor : ",e10.3)') & n, names(n), cvsfac(n); call goPr end if end do END SUBROUTINE CALC_CVSFAC !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: LSPSCAV ! ! !DESCRIPTION: Calculation of wet removal by large scale precipitation of ! soluble tracers ! ! Remove tracers with previously calculated rainout rate [s-1] ! separately for in- and below cloud scavenging and only for the ! cloud covered fraction of the gridcell ! ! Reference: ! Langner and Rodhe (1990) ! Junge (1959) !\\ !\\ ! !INTERFACE: ! SUBROUTINE LSPSCAV( region ) ! ! !USES: ! use binas, only : rgas use dims, only : im, jm, lm, nchem use dims, only : tref, lmax_conv, dy use chem_param, only : ntrace, ntracet, henry, ntlow, ra, ntemp use chem_param, only : xmhno3 !, ihno3 use chem_param, only : nscav, nscav_index, nscav_type #ifdef with_m7 use chem_param, only : inus_n, iais_n, iacs_n, icos_n, iaii_n, iaci_n, icoi_n #endif use meteodata, only : temper_dat, cc_dat use global_data, only : mass_dat #ifdef with_budgets use budget_global, only : nzon_vg #endif ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region ! ! !REVISION HISTORY: ! ! Programmed by Frank Dentener, August 1995; ! Ad Jeuken, KNMI, November 1998 ! Modifications Bas Henzing, KNMI, 2002 ! Adapted for TM5, Frank Dentener, JRC, 2002 ! Paralel, Maarten Krol, Jul 2003 ! 29 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/LSPSCAV' real :: dt_lagrangian ! liquid water content of raining cloud ! rgas (8.314 J/mol/K) ---> 0.08314 atm/(mol/l)/K ! (thesis Frank Dentener, p. 31) ! 1e-6 corresponds to 1 g/m3 dimensionless real,parameter :: rtl_fac=rgas/1e2*1e-6 ! assumed pH of rainwater real,parameter :: hplus = 1e-5 ! assumed fraction of in-cloud interstitial aerosol: real,parameter :: interstitial_fraction = 0.3 ! --- local ------------------------------ real,dimension(:,:,:,:), pointer :: rm #ifdef slopes real,dimension(:,:,:,:), pointer :: rxm, rym, rzm #ifdef secmom real,dimension(:,:,:,:), pointer :: rxxm, rxym, rxzm, ryym, ryzm, rzzm #endif #endif real,dimension(:,:,:), pointer :: t,cc real :: rtl,f,f1,f2,vol1,vol2,vol3,ahelp1,ahelp2 real :: incloud,belowcloud,newcov,c_old,corr_diff,fnchem integer :: n,iscav,i,j,k,itemp,nzone_v, nloc real :: ztr, dkso2, dkhso3, factor, heff ! oldcov: cloud cover in layer above, to calculate below-cloud scaveging. real,dimension(:,:),allocatable :: oldcov integer :: i1, i2, j1, j2, status ! for wetdep global budget of tracer #1 real :: g_sum_wet real, allocatable :: l_sum_wet(:,:) ! ! --- begin ------------------------------ ! !>>> TvN ! Tests have been performed at various resolutions ! using mixing times of 3, 6, 9, 12 and 24 hours. ! Simulations at 1x1 degrees are best reproduced ! by increasing the mixing time with 3 hours at 3x2 degrees ! and with 6 hours at 6x4 degrees. dt_lagrangian = 3600.0 * 3.0 ! value at 1x1 degrees or higher resolution if (dy == 2) dt_lagrangian = dt_lagrangian + 3600.0 * 3.0 if (dy == 4) dt_lagrangian = dt_lagrangian + 3600.0 * 6.0 !<<< TvN CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) #ifdef with_budgets allocate(l_sum_wet(i1:i2,j1:j2)) l_sum_wet = 0.0 #endif rm => mass_dat(region)%rm ! paralel over tracers #ifdef slopes rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm #ifdef secmom rxxm => mass_dat(region)%rxxm rxym => mass_dat(region)%rxym rxzm => mass_dat(region)%rxzm ryym => mass_dat(region)%ryym ryzm => mass_dat(region)%ryzm rzzm => mass_dat(region)%rzzm #endif #endif t => temper_dat(region)%data cc => cc_dat(region)%data !$OMP PARALLEL & !$OMP default (none) & #if defined (with_budgets) !$OMP shared (nzon_vg, buddep_dat) & #endif !$OMP shared (region, ier, isr, jer, jsr, tracer_loc, henry, & !$OMP tracer_active, fnchem, rm, rxm, rym, rzm, t, cc, & !$OMP nchem, rloss1, rloss1_m7, rloss2, rloss3, ra,nscav_type, & !$OMP nscav_index,tref,im ) & !$OMP private (i, j, k, n, nloc, iscav, jss, jee, rtl, itemp, & !$OMP corr_diff, belowcloud, incloud, f, f1, f2, newcov, & !$OMP vol1, vol2, vol3, c_old, nzone_v, ztr, dkso2, dkhso3, & !$OMP factor, heff, oldcov, l_sum_wet) fnchem=real(nchem/(2*tref(region))) ! allocate(oldcov(i1:i2,j1:j2)) do iscav=1,nscav ! n=nscav_index(iscav) ! tracer number in global count nloc = n ! tracer number in local count oldcov=0. ! ! assumption no stratiform precipitation above the maximum ! level of convection ! do k=lmax_conv,1,-1 ! top to bottom do j=j1,j2 do i=i1,i2 ! ! calculate, depending on solubility, scale factor relative ! to 100% scavenging (of HNO3) ! ! rtl - composite factor of liquid water content, rgas ! and liquid water content rtl=rtl_fac*t(i,j,k) ! multiplaction with Henry constant gives phase factor itemp=nint(t(i,j,k)-float(ntlow)) itemp = min(max(1,itemp),ntemp) ! corrected CMK dec2004 problems on ECMWF ! !corr_diff=sqrt(ra(ihno3)/ra(n)) corr_diff=sqrt(xmhno3/ra(n)) ! select case (nscav_type(iscav)) case(0) incloud = 0.0 belowcloud = 0.0 case(1) ! 100% solubility !AJS: note that old code used factor rtl ~ 1, !AJS: computed with huge henry factor ~ 1e7 : ! rtl = rtl*henry(n,itemp) / ( 1.0 + rtl*henry(n,itemp) ) ! near 1.0 ! incloud = rloss1(reion)%d3(i,j,k) * rtl incloud = rloss1(region)%d3(i,j,k) belowcloud = rloss2(region)%d3(i,j,k)*corr_diff case(2) ! henry solubility assumed rtl = rtl*henry(n,itemp) / ( 1.0 + rtl*henry(n,itemp) ) incloud = rloss1(region)%d3(i,j,k)*rtl belowcloud = rloss2(region)%d3(i,j,k)*rtl*corr_diff case(3) ! bulk aerosol incloud = rloss1(region)%d3(i,j,k)*(1.0 - interstitial_fraction) !>>>TvN ! Alternative would be to make the interstitial fraction for bulk aerosols ! consistent the values used for the M7 modes, ! which are taken from Bourgeois and Bey (JGR, 2011) ! and distinguish between warm, mixed and ice clouds !<<hso3m+hplus dkhso3 = 6.6e-8*exp(1510.*ztr) !hso3m<=>so3-- + hplus factor = 1.0 + dkso2/hplus + (dkso2*dkhso3)/(hplus**2) heff = factor*henry(n,itemp) rtl = rtl*heff/ ( 1.0 + rtl*heff ) incloud = rloss1(region)%d3(i,j,k)*rtl ! belowcloud = rloss2(region)%d3(i,j,k)*rtl*corr_diff !>>>TvN ! The in-cloud scavenging coefficients are defined as the fraction of the tracer ! in the cloudy part of the grid box that is embedded in the cloud liquid or ice water, ! i.e. the non-interstitial part. ! We distinguish between liquid, mixed and ice stratiform clouds (Stier et al., 2005), ! depending on the local temperature in the grid cell (Croft et al., ACP, 2010). ! The in-cloud scavenging coefficients depend on size and composition; ! revised values for the M7 modes were provided by Bourgeois and Bey (JGR, 2011). ! For mixed clouds, an alternative method was presented by Zhang et al. (ACP, 2012), ! which uses a continuous temperature dependency. ! Note that these in-cloud scavenging coefficients account for both nucleation scavenging ! and impaction scavenging (Croft et al., ACP, 2009; 2010). ! Thus, the below-cloud scavenging rates should only account for ! the impaction scavenging by precipitation coming from clouds above the current level. ! ! Estimates for below-cloud scavenging coefficients can be derived ! from Fig. 2 of Dana and Hales (AE, 1975). ! For estimating these values from the figure, I used aerodynamic radii of ! 0.007, 0.07, and 0.7 micron as the boundaries of the M7 modes ! (corresponding to a particle density of about 1800 g/cm^3). ! As in Stier et al. (2005), we do not distinguish between soluble and insoluble modes. ! Thus, dry particle radii can be used for estimating the scavenging coefficients from the figure ! (see also the mode boundaries applied in Fig. 2 in Croft et al., 2009). ! I thus arrive at the following rough estimates for below-cloud mass scavenging coefficients ! for the nucleation, aitken, accumulation, and coarse modes: ~0.01, 0.002, 0.01, and 1 mm^-1. ! These numbers are close to the estimates derived earlier from the same figure ! by Elisabetta Vignati, which were previously used: 0.005, 0.002, 0.008, and 1 mm^-1. ! ! However, both sets of estimates based on Dana and Hales are substantially higher ! than the values presented by Croft et al. (2009). ! From the curves presented in their Fig. 2 for the standard Marshall-Palmer rain distribution, ! rough estimates of the mass scavenging coefficients for the four size modes can be obtained. ! My estimates are 0.002, 0.0002, 0.03, and 0.7 mm^-1. ! Note that especially the value for the accumulation mode is very sensitive to the ! actual mean particle size, and hard to estimate from the figure. ! Since the mean droplet size of the Marshall-Palmer distribution depends on the rain intensity, ! these estimates are only valid for a rain rate of 1 mm/hr. ! For simplicity, we assume that the scavenging coefficients derived from the figure at 1 mm/hr ! can also be applied at other rain intensities. ! ! In the new implementation particle masses and numbers are scavenged at different rates. ! Rough estimates of the number scavenging coefficients for the four size modes ! can be obtained from the same figure in Croft et al. ! My estimates are 0.02, 0.001, 0.0003, and 0.3 mm^-1. ! Ideally, the below-cloud mass/number scavenging coefficients should be calculated ! using look-up tables to describe the dependence on median radius and precipitation rate, ! e.g. following the formulation/curves presented by Croft et al. ! case(5) ! soluble nu !belowcloud=0.1*rloss3(region)%d3(i,j,k) ! 0.1*0.05=0.005 mm^-1 #ifdef with_m7 if (n /= inus_n) then belowcloud=0.5*rloss3(region)%d3(i,j,k) ! 0.5*0.004 = 0.002 mm^-2 else #endif belowcloud=5.*rloss3(region)%d3(i,j,k) ! 5.*0.004 = 0.02 mm^-2 #ifdef with_m7 endif #endif !incloud=0.0 incloud=0.06*rloss1_m7(region)%d3(i,j,k) case(6) ! soluble ai !belowcloud=0.04*rloss3(region)%d3(i,j,k) ! 0.04*0.05=0.002 mm^-1 #ifdef with_m7 if (n /= iais_n) then belowcloud=0.05*rloss3(region)%d3(i,j,k) ! 0.05*0.004 = 0.0002 mm^-2 else #endif belowcloud=0.25*rloss3(region)%d3(i,j,k) ! 0.25*0.004 = 0.001 mm^-2 #ifdef with_m7 endif #endif !incloud=0.0 if (t(i,j,k).gt.273.15) then incloud=0.25*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case(7) ! soluble ac !belowcloud=0.16*rloss3(region)%d3(i,j,k) ! 0.16*0.05 = 0.008 mm^-1 #ifdef with_m7 if (n /= iacs_n) then belowcloud=7.5*rloss3(region)%d3(i,j,k) ! 7.5*0.004 = 0.03 mm^-1 else #endif belowcloud=0.075*rloss3(region)%d3(i,j,k) ! 0.075*0.004 = 0.0003 mm^-1 #ifdef with_m7 endif #endif !incloud=rloss1(region)%d3(i,j,k) if (t(i,j,k).gt.273.15) then incloud=0.85*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case(8) ! soluble co !belowcloud=20.*rloss3(region)%d3(i,j,k) ! 20*0.05 = 1 mm^-1 #ifdef with_m7 if (n /= icos_n) then belowcloud=175.*rloss3(region)%d3(i,j,k) ! 175*0.004 = 0.7 mm^-1 else #endif belowcloud=75.*rloss3(region)%d3(i,j,k) ! 75*0.004 = 0.3 mm^-1 #ifdef with_m7 endif #endif !incloud=rloss1(region)%d3(i,j,k) if (t(i,j,k).gt.273.15) then incloud=0.99*rloss1_m7(region)%d3(i,j,k) else if (t(i,j,k).gt.238.15) then incloud=0.75*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case(9) ! insoluble ai !belowcloud = 0.04*rloss3(region)%d3(i,j,k) #ifdef with_m7 if (n /= iaii_n) then belowcloud=0.05*rloss3(region)%d3(i,j,k) else #endif belowcloud=0.25*rloss3(region)%d3(i,j,k) #ifdef with_m7 endif #endif !incloud=0.0 if (t(i,j,k).gt.273.15) then incloud=0.2*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case(10) ! insoluble ac !belowcloud=0.16*rloss3(region)%d3(i,j,k) #ifdef with_m7 if (n /= iaci_n) then belowcloud=7.5*rloss3(region)%d3(i,j,k) else #endif belowcloud=0.075*rloss3(region)%d3(i,j,k) #ifdef with_m7 endif #endif !incloud=0.0 if (t(i,j,k).gt.273.15) then incloud=0.4*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case(11) ! insoluble co !belowcloud=20.*rloss3(region)%d3(i,j,k) #ifdef with_m7 if (n /= icoi_n) then belowcloud=175.*rloss3(region)%d3(i,j,k) else #endif belowcloud=75.*rloss3(region)%d3(i,j,k) #ifdef with_m7 endif #endif !incloud=0.0 if (t(i,j,k).gt.238.15) then incloud=0.4*rloss1_m7(region)%d3(i,j,k) else incloud=0.06*rloss1_m7(region)%d3(i,j,k) endif case default incloud = 0.0 belowcloud = 0.0 end select !if(incloud > 1e-4) then !print *, i,j,k,names(n),rtl, rloss1(region)%d3(i,j,k), rtl_fac !end if ! ! f1, f2 are the fractions of the tracermass that remain in the ! gridbox after scavenging. ! f1=exp(-dt_lagrangian*incloud) f2=exp(-dt_lagrangian*belowcloud) !CMK f1=exp(-fnchem*incloud) !CMK f2=exp(-fnchem*belowcloud) ! ! A grid box can be divided into three volumes ! 1) incloud scavenging (newcov) ! 2) below cloud scavenging (oldcov-newcov) ! 3) no in-cloud scavenging and no below cloud ! scavenging by precipitation (no removal) ! newcov=cc(i,j,k) vol1 = newcov !>>> TvN ! oldcov denotes the area fraction occupied by precipitating clouds above the current level. ! It is determined assuming maximum overlap of the cloudy fractions of the layers above (see below). ! The precipitation rate used for below-cloud scavenging correctly does not include the contribution ! from the current level. !<< 0.0 ) oldcov(i,j)=max(oldcov(i,j),cc(i,j,k)) end do !i end do !j end do !k end do !iscav deallocate(oldcov) !$OMP END PARALLEL #ifdef with_budgets call REDUCE( dgrid(region), l_sum_wet, 0, 'SUM', g_sum_wet, status) IF_NOTOK_RETURN(status=1) if ( isRoot ) sum_wetdep(region) = sum_wetdep(region) + g_sum_wet deallocate( l_sum_wet ) #endif nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #ifdef secmom nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) #endif #endif nullify(t) nullify(cc) END SUBROUTINE LSPSCAV !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CALCULATE_LSP_SCAV ! ! !DESCRIPTION: ! ! Calculate wet removal rates rloss1,rloss2,rloss3 (s-1) for ! stratiform precipitation from cloud-top to ground, ! distinguishing between below-cloud and in-cloud scavenging. ! ! Method: ! adapted from GJ Roelofs and Lelieveld, JGR, October 1995 ! ! fills array "rloss" should be called once new precipitation fields ! are available (MK: in trace_after_read) !\\ !\\ ! !INTERFACE: ! SUBROUTINE CALCULATE_LSP_SCAV( region ) ! ! !USES: ! use binas, only : grav, rgas, xmair use dims, only : im,jm,lm,lmax_conv use meteodata, only : temper_dat, lwc_dat, iwc_dat, cc_dat use meteodata, only : lsp_dat use global_data, only : emis_data use meteodata, only : phlb_dat use partools, only : isroot ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !REVISION HISTORY: ! Programmed by: Frank Dentener, IMAU, 1996 ! Ad Jeuken, KNMI, 1998 ! Modification, Bas Henzing, KNMI, 2002. ! Adapted for TM5, August 2002, Frank Dentener, JRC. ! And finally for the new version, MK, IMAU, March 2003 ! Parallel, MK Jul 2003 ! 25 Jan 2012 - P. Le Sager - adapted for mpi lat-lon domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC integer :: i1, i2, j1, j2 real,dimension(:,:,:),pointer :: lsp real,dimension(:,:,:),pointer :: t, lwc, iwc, cc, phlb real,parameter :: max_lwc=2.e-3 ! kg/m3 ! ! Microphysical parameters: ! ! rdrad2: square of raindroplet radius (20 microns) ! dghno3: in [cm2/s] ! dgair: in [cm2/s] ! rol: density of water in [kg/m3] ! roi: density of ice in [kg/m3] ! real,parameter :: rdrad2 = (2E-5)**2 real,parameter :: dghno3 = 0.136 real,parameter :: dgair = 0.133 real,parameter :: rol = 1000. real,parameter :: roi = 917. ! ! quantity used in calculation of Sherwood number ! ! bas henzing: bug repair real,parameter :: znsc=dgair/dghno3**(-3) ! bas henzing: should be: znsc=(dgair/dghno3)**(1./3.); ! znsc is now defined a real ! real,dimension(:),allocatable :: dzk real :: rflx,beta,fac,beta1,beta2,beta3,rlwc,rdrad,rn,ru !>>>TvN real :: fac_m7, beta1_m7 !<<< real :: press,aird,dgairx,dghno3x real :: znre,znsh,znsc,zkg,totw,sfz,sfz_no ! integer :: nfz,i,j,k,l,no_fz integer,dimension(:,:),allocatable :: kcltop real,dimension(:,:),allocatable :: oldcov,fz real,dimension(:,:,:),allocatable :: precip ! precipitation per level (kg/m2/s) real,dimension(:,:,:),allocatable :: prf ! precipitation formation rate. ! evaporation NOT YET AVAILABLE ! ! how much less efficient is tracer scavenged from ice ! cloud droplet compared to water cloud droplet. ! This should be tracer dependent. ! real,parameter :: ice_eff=0.2 real :: inc_rdf real,parameter :: scale_heigth= rgas/grav/xmair*1e3 ! about 29.3*T (m) ! ---------------- START ------------------------------------------------ t => temper_dat(region)%data lwc => lwc_dat(region)%data !mk: lm, and not lmax_conv iwc => iwc_dat(region)%data !mk: lm, and not lmax_conv cc => cc_dat(region)%data !mk: lm, and not lmax_conv phlb => phlb_dat(region)%data !mk: 1:lm+1 lsp => lsp_dat(region)%data CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( oldcov( i1:i2, j1:j2 ) ) allocate( fz( i1:i2, j1:j2 ) ) allocate( precip( i1:i2, j1:j2, lmax_conv) ) allocate( dzk( lmax_conv) ) allocate( kcltop( i1:i2, j1:j2 ) ) allocate( prf( i1:i2, j1:j2, lmax_conv) ) ! ! calculate precipitation formation rate prf ! call calfk ! ! initialize cloud top ! kcltop(:,:)=lmax_conv ! ! calculate and rescale precip ! sfz=0. nfz=0 precip(:,:,:)=0. ! do j= j1, j2 do i= i1, i2 ! ! Calculate precipitation intensity at the bottom of each layer ! do k=1,lmax_conv-1 ! thickness of layer, only correct in troposphere dzk(k)=scale_heigth*t(i,j,k)*alog(phlb(i,j,k)/(1.+phlb(i,j,k+1))) end do ! do k=lmax_conv-1,1,-1 precip(i,j,k)=precip(i,j,k+1)+prf(i,j,k)*dzk(k) !precip: kg/m2/s end do ! ! Rescale prf and precip such that these are consistent with ground lsp ! Note that lsp is defined as the total large-scale precipitation (rain+snow) ! produced by the cloud scheme. ! prf, precip and lsp thus all account for snow/ice ! no_fz = 0 ! for statistics !CMK was not initialised! sfz_no = 0.0 fz(i,j)=0. !cmk if (lsp(i,j) > 1.e-5) then old data came in mm/day if (lsp(i,j,1)*1e3*3600.*24. > 1.e-5) then !new data are in m/s if (precip(i,j,1) > 0.) then fz(i,j)=lsp(i,j,1)*1e3/precip(i,j,1) ! m/s ->kg/m2/s !new unit... !cmk fz(i,j)=lsp(i,j,1)/3600./24./precip(i,j,1) ! mm/day->kg/m2/s nfz=nfz+1 ! precipitation statistics ! (avoid 'strange' statistics by only few values) sfz=sfz+min(10.,fz(i,j)) else ! no precipitation calculated, but ECMWF fields contain rain value no_fz=no_fz+1 sfz_no=sfz_no+lsp(i,j,1)*1e3*86400. ! (in mm/day) end if else precip(i,j,:)=0. end if do k=1,lmax_conv precip(i,j,k)=precip(i,j,k)*fz(i,j) prf(i,j,k)=prf(i,j,k)*fz(i,j) end do !k end do ! i end do ! j ! !if(myid == root) then ! write(6,*) 'calculate_lsp_scav: region',region ! write(6,*) ' rainout: average scale factor for precipitation = ',sfz/nfz ! write(6,*) ' rainout: percentage of columns with precipitation = ', & ! 100.*nfz/real(im(region)*jm(region)),' %' ! if ( no_fz > 0 ) write(6,*) 'rainout: lsp and prf not consistent ', & ! no_fz,'average rainfall [mm/day]',sfz_no/no_fz !end if ! ! initialise ! oldcov=0. ! evap=0. rloss1(region)%d3=0. !>>>TvN rloss1_m7(region)%d3=0. !<<1E-15 ! if( (precip(i,j,k+1)>1e-15) .and. (k0.) )then ! ! rflx in [mm/hr] !MK? thought precip was in kg/m2/s ?? ! rlwc in [vol mixing ratio] ! rdrad in [cm] ! ru in [cm/s] (terminal velocity) ! znre = Reynolds number ! znsc = (Sherwood number)**(1./3.) ! zkg in [cm/s] = dghno3[cm^2/s]/rdrad[cm] ! rflx = precip(i,j,k+1)/oldcov(i,j)*3600. rlwc = 72.*rflx**0.88*1.e-9 rdrad = 0.1*0.3659*rflx**0.21 !******************************************* ! correction by Twan van Noije, Bas Henzing: ! ru = 9.58*(1.-exp(-(rdrad*10./0.885)**1.147)) ! the above equation gives an approximation to the terminal velocity in m/s !! ! a conversion factor to cm/s should therefore be included ! znre = 20.*rdrad*ru/dgair ! with ru in cm/s, the above expression overestimates the Reynolds number by a factor 10 ! the combined effect of having ru in m/s and the above expression for the Reynolds number ! is to underestimate the Reynolds number by a factor 10, resulting in an underestimation ! of the Sherwood number and an overestimation of the below-cloud scavenging ! Twan van Noije/Bas Henzing, 24-02-2004 !******************************************* ru = 100.*9.58*(1.-exp(-(rdrad*10./0.885)**1.147)) ! see Seinfeld (1986), p. 632 znre = 2.*rdrad*ru/dgair !WRONG ru = 9.58*(1.-exp(-(rdrad*10./0.885)**1.147)) !WRONG znre = 20.*rdrad*ru/dgair znsc = (dgair/dghno3)**(1./3.) znsh = 1.+0.3*(znre**0.5)*znsc zkg = dghno3/rdrad*znsh beta2 = 3.*zkg*rlwc/rdrad ! ! beta3 for below cloud scavenging of accumulation range aerosol ! (Dana and Hales, Atmos. Env. 1976, pp. 45-50 ! assuming a Whitby aerosol distrbution r=0.034 sigma=2; ! mass-mean radius r=0.14 microm; ! figure 2 gives a wash-out coefficient of 0.05 mm^-1 (raindepth) ! for other sigma and r look-up tables can be generated !>>>TvN ! A mass washout coefficient of 0.05 mm^-1 in Fig. 2 of Dana and Hales (1975) ! would correspond to a geometric mean/median particle radius of 0.14 micron. ! At a median radius of 0.034 micron and sigma=2, the value becomes ~0.002 mm^-1. ! However, as the aerodynamic radius is used in the plot, these values only apply ! to unit-density particles, with a density equal to that of water (1 g/cm^3). ! The aerodynamic radius equals the actual radius times ! the square root of ratio of the particle density over that for water. ! For the bulk inorganic aerosols, the particle density is about 1800 g/cm^3, ! so the aerodynamic radius is 1.34 times the actual radius, ! and a median aerodynamic radius of ~0.046 micron should be used. ! This gives a mass washout coefficient of ~0.004 mm^-1. ! Thus, based on the work by Dana and Hales, the value 0.05 mm^-1 seems too low, ! and a value ~0.004 mm^-1 would be more appropriate. ! ! mm-1*[mm hr-1] * [hr/s] => [s-1] ! !beta3= 0.05*rflx/3600. beta3= 0.004*rflx/3600. ! !<<1E-20.AND.precip(i,j,k+1)>1E-20) THEN ! rev1=max(0.,EVAP(i,j,k)/precip(i,j,k+1)) ! ! evaporation fraction ! rev(i,j,k)=MIN(1.,rev1) ! END IF ! !============================== ! in cloud scavenging (beta1) !============================== ! if (totw>1.e-9.and.prf(i,j,k)>0.and.cc(i,j,k)>0.02) then ! if(k>kcltop(i,j)) kcltop(i,j)=k !set new cloud top ! ! rlwc: convert mass mixing ratios to volume mixing ratios in cloud ! aird: air density in kg air / m^3 ! max_lwc: in kg H2O / m^3 ! prf: in kg H2O / m3 / s ! beta: in [s-1] = [vol mixing ratio]*[cm2/s]*1e-4/[m2] ! fac, beta1: in [s-1] ! aird=press/t(i,j,k)/rgas*xmair*1.e-3 rlwc=(lwc(i,j,k)/rol+iwc(i,j,k)/roi)*aird/cc(i,j,k) rlwc=min(max_lwc/rol,rlwc) ! !bas henzing: bug repair: beta=3.*rlwc*(dghno3*1e-4)/rdrad2 !bas henzing: should be: ! beta=(3.*rlwc*(dghno3*1e-4)/(2.*rdrad2))*znsh !bas henzing: reference: (Roelofs and Lelieveld, 1995) !fd beta=(3.*rlwc*(dghno3*1e-4)/(2.*rdrad2))*znsh ! fd 13/08/2002 use old defenition again (pers. comm Henzing) beta=3.*rlwc*(dghno3*1e-4)/rdrad2 ! inc_rdf=(iwc(i,j,k)/totw)*ice_eff+lwc(i,j,k)/totw fac=prf(i,j,k)*inc_rdf/(totw*aird) !>>>TvN ! In the calculation of fac (and thus fac_m7) currently the grid-cell average prf is used. ! Question: shouldn't we use the actual value in the cloudy part, i.e. prf/cc ? ! Note that scaling by 1/cc is also applied in the calculation of beta, ! and that for the below-cloud scavenging parameters beta2 and beta3 ! we also use the actual precipitation intensity in the precipitating fraction, ! i.e. precip/oldcov. !<<>>TvN ! The scavengy efficiency for in-cloud scavenging of M7 aerosols by ice in stratiform clouds ! is now reduced by applying a lower scavenging coefficient for mixed and ice clouds ! according to Bourgeois and Bey (JGR, 2011). fac_m7=prf(i,j,k)/(totw*aird) beta1_m7=1./(1./beta+1./fac_m7) !<<>>TvN rloss1_m7(region)%d3(i,j,k)=beta1_m7 ! as rloss1, but with ice as efficient as liquid water !<<1.e-5) then if ( lsp(i,j,1)*1e3*3600.*24. > 1.e-5 ) then !new data are in m/s ! determine cloud base k = 1 do if ( cc(i,j,k) >= 0.01 ) exit if ( k == lmax_conv ) exit k = k + 1 end do iclb = k ! determine cloud top k = lmax_conv do if ( cc(i,j,k) >= 0.01 ) exit if ( k == 1 ) exit k = k-1 end do iclt = k ! check for inconsistency in cloud cover fields if ( iclb >= iclt ) iclb=iclt if ( iclb < 1 ) iclb=1 !mvw: replace fixed iclb/t=14 by 120 hPa level (about 15km) ! if (iclb>14) iclb=14 ! if (iclt>14) iclt=14 icldtop=lvlpress(region,12000., 98400.) if ( iclb > icldtop ) iclb=icldtop if ( iclt > icldtop ) iclt=icldtop ! pr0=0. pcl=0. pci=0. rain_frac=0.00001 ! ! start at cloudtop do k=iclt,iclb,-1 zcc=max(0.01,cc(i,j,k)) ! ! precipitation formation from ice clouds ! ! pci in kg_H2O / (kg_air sec) ! if ( ( t(i,j,k) < 258.0 ) .and. ( k > 1 ) ) then ccp=max(0.01,cc(i,j,k+1)) ccm=max(0.01,cc(i,j,k-1)) qup=(iwc(i,j,k)/zcc+iwc(i,j,k+1)/ccp)*0.5 qdo=(iwc(i,j,k)/zcc+iwc(i,j,k-1)/ccm)*0.5 rup=(zrho(k)+zrho(k+1))*0.5 rdo=(zrho(k)+zrho(k-1))*0.5 vtup=alfa*(rup*qup)**beta vtdo=alfa*(rdo*qdo)**beta pci(k)=grav*(vtup*qup*rup-vtdo*qdo*rdo)/ & (phlb(i,j,k+1)-phlb(i,j,k)) pci(k)=max(pci(k),0.) end if ! ! precipitation formation from liquid clouds ! formulation ECMWF ! qcl=lwc(i,j,k)/zcc qcl=min(max_lwc/zrho(k),qcl) qcl=max(0.001e-3/zrho(k),qcl) ! ! pcl in kg_H2O / (kg_air sec) ! plocal=pr0/rain_frac f1=1.+b1*sqrt(plocal) f2=1.+b2*sqrt(max(0.,Tberg-t(i,j,k))) c0=cr1*f1*f2 qcr=qcrs/(f1*f2) pcl(k)=zcc*c0*qcl*(1.-exp(-(qcl/qcr)**2)) ! ! prec at top next layer in kg m-2 s-1 ! dzk=29.3*t(i,j,k)*alog(phlb(i,j,k)/(1.+phlb(i,j,k+1))) delta_prec=(pcl(k)+pci(k))*zrho(k)*dzk ahelp=(zcc*delta_prec+rain_frac*pr0)/(delta_prec+pr0) rain_frac=max(rain_frac,ahelp) pr0=pr0+(pcl(k)+pci(k))*zrho(k)*dzk ! ! liquid+ice precipitation formation rates in kg m-3 s-1 ! prf(i,j,k)= (pcl(k)+pci(k))*zrho(k) ! end do ! k end if ! lsp gt 1.e-5 end do ! i end do ! j deallocate(zrho) deallocate(pci) deallocate(pcl) ! end subroutine calfk END SUBROUTINE CALCULATE_LSP_SCAV !EOC END MODULE WET_DEPOSITION