!### macro's ##################################################### ! #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" ! !################################################################# !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: convection ! ! !DESCRIPTION: methods to init and apply convection !\\ !\\ ! !INTERFACE: ! module convection ! ! !USES: ! use GO, only : gol, goPr, goErr implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Convection_Init, Convection_Done public :: Convec ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'convection' ! ! !REVISION HISTORY: ! 30 Mar 2010 - P. Le Sager - re-instate test on wet dep flag in convec ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ contains ! ================================================================ subroutine Convection_Init( status ) use dims , only : nregions use Meteo , only : Set use Meteo , only : entu_dat, entd_dat, detu_dat, detd_dat use Meteo , only : mfw_dat, omega_dat #ifndef without_wet_deposition use Meteo , only : cp_dat #endif ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Convection_Init' ! --- local ------------------------------------- integer :: region ! --- begin -------------------------------- ! loop over all (zoom) regions: do region = 1, nregions ! entrements/detrements call Set( entu_dat(region), status, used=.true. ) call Set( entd_dat(region), status, used=.true. ) call Set( detu_dat(region), status, used=.true. ) call Set( detd_dat(region), status, used=.true. ) #ifndef without_wet_deposition ! convective precip: call Set( cp_dat(region), status, used=.true. ) #endif ! omega is used when creating updr/downdr from raw data; ! computed from mfw : call Set( mfw_dat(region), status, used=.true. ) call Set( omega_dat(region), status, used=.true. ) end do ! ok status = 0 end subroutine Convection_Init ! *** subroutine Convection_Done( status ) ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Convection_Done' ! --- begin -------------------------------- ! nothing to be done ! ok status = 0 end subroutine Convection_Done ! ================================================================ !----------------------------------------------------------------------- ! !**** convec - mix tracers by convection v 8.5 ! ! programmed by mh mpi HH 1-oct-1991 ! modified by mh mpi HH 11-jun-1994 ! ! purpose ! ------- ! mix tracers vertically by convection ! ! interface ! --------- ! call convec ! ! method ! ------ ! the matrix conv is applied simultaneously to rm, rxm and rym. ! The diagonal elements of conv are applied to rzm. ! ! externals ! --------- ! subroutines: tstamp ! ! reference ! --------- ! see manual !----------------------------------------------------------------------- ! Following the method by Walter Guelle (1997) convective removal of ! tracers is added with some little changes: ! a vector cvsfac is added containing scavenging efficiencies (0-1) ! per tracer. ! aj, knmi may 1998 ! Maarten Krol, july 2000, implemented zoom version ! lmax_conv is now set as the maximum layer to which convection ! is taken into account. ! ! Second moments have been added: Bram Bregman, August 2004 !-------------------------------------------------------------- subroutine convec(region, status) use dims, only : dx, dy, xref, yref, tref use dims, only : isr, ier, jsr, jer use dims, only : im, jm, lm use dims, only : itau, okdebug, kdebug, nconv, lmax_conv use dims, only : nregions use dims, only : revert use global_data, only : mass_dat,region_dat use global_data, only : emis_data use global_data, only : conv_dat use meteo , only : entu_dat, entd_dat, detu_dat, detd_dat use meteo , only : cp_dat use meteodata , only : m_dat use chem_param, only : ra, ntracet #ifdef with_budgets use budget_global, only : nzon_vg, apply_budget_global #ifndef without_wet_deposition use wet_deposition,only : sum_wetdep, buddep_dat #endif #endif use toolbox, only : lvlpress #ifndef without_wet_deposition use wet_deposition,only : cvsfac, cp_scale #endif use datetime, only : tstamp #ifdef with_tendencies use tracer_data , only : PLC_AddValue_t, plc_ipr_lwdep, plc_ipr_tcnvd, plc_kg_from_tm #endif #ifdef MPI use mpi_const, only: my_real, & mpi_min,mpi_max,com_trac,ierr use mpi_const, only: pe_first_tracer #endif use ParTools, only : myid, root_t use ParTools, only : ntracetloc, ntracet_ar use omp_ParTools, only : my_omp_domain, my_omp_sum implicit none ! input/output integer,intent(in) :: region integer,intent(out) :: status ! 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 :: m real,dimension(:,:,:) , pointer :: entu, entd, detu, detd #ifndef without_diffusion real,dimension(:,:,:) , pointer :: dkg #endif integer,dimension(:,:) , pointer :: cloud_base, cloud_top, cloud_lfs, zoomed integer,dimension(3) :: mxloc ! arrays for convection matrix ! conv1: convection matrix, lbdcv: wet removal matrix real,dimension(lmax_conv,lmax_conv) :: conv1 real,dimension(lmax_conv,lmax_conv) :: lbdcv real,dimension(0:lmax_conv,lmax_conv) :: f real,dimension(0:lmax_conv,lmax_conv) :: fu real,dimension(0:lmax_conv) :: amu,amd real,dimension(lmax_conv) :: new_mass real,dimension(lmax_conv) :: zwetrm, zcnvd real,dimension(lmax_conv) :: zrm real,dimension(lmax_conv) :: srm #ifdef slopes real,dimension(lmax_conv) :: zrxm,zrym real,dimension(lmax_conv) :: srxm,srym,srzm #ifdef secmom real,dimension(lmax_conv) :: zrxxm,zrxym,zrxzm,zryym,zryzm real,dimension(lmax_conv) :: srxxm,srxym,srxzm,sryym,sryzm,srzzm #endif #endif integer :: n,l,k,j,i integer :: imr, jmr, lmr integer :: info,ld,li,kk, lshallowtop, nzv real :: scavf,dt,zxi,wetpct real :: minvalue, mxval,mnval,mxvald real :: minvalue_all, mxval_all,mnval_all,mxvald_all real :: sceffdeep, cp_value integer :: nglob, offsetn ! Ompenmp parameters integer :: jss,jee real :: l_sum_wet, g_sum_wet ! start imr=im(region) ; jmr=jm(region) ; lmr=lm(region) !WP! only PE's with nonzero ntracet if ( ntracetloc == 0 ) return offsetn = sum(ntracet_ar(0:myid-1)) sceffdeep = 1.0 ! rain scale on this resolution: e.g. on 6x4 this read 1x1/(6x4) ! thus 1 mm/hr becomes 1mm/day ! cp_scaler = cp_scale*xref(region)*yref(region)/(dx*dy) !CMK removed after test study: hardly resolution dependent;; m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #ifdef secmom rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t #endif #endif zoomed => region_dat(region)%zoomed entu => entu_dat(region)%data detu => detu_dat(region)%data entd => entd_dat(region)%data detd => detd_dat(region)%data cloud_top => conv_dat(region)%cloud_top cloud_base => conv_dat(region)%cloud_base cloud_lfs => conv_dat(region)%cloud_lfs #ifndef without_diffusion dkg => conv_dat(region)%dkg #endif if ( okdebug ) call tstamp(kdebug,itau,'convec ') if ( okdebug ) print *,'convec: Convection called (region)', region !timestep and position in packed arrays dt = nconv/(2.*tref(region)) if ( okdebug ) then minvalue = 2.0 mxvald = -1.0 mxval = 0.0 mnval = 2.0 end if ! start main loop over surface cells ! ! PLS : re-incorporate the test on wet dep, from older cy3base/src. ! It hides cvsfac and cp_scale when wet dep is off. !$OMP PARALLEL & !$OMP default (none) & #ifdef with_budgets !$OMP shared (nzon_vg, apply_budget_global) & #endif !$OMP shared (region, ntracetloc, sceffdeep, dt, offsetn, & !$OMP g_sum_wet, cloud_base, cloud_top, cloud_lfs, cp_dat, & !$OMP zoomed, rm, rxm, rym, rzm, m, entu, entd, detu, detd, & !$OMP dkg, revert, isr, ier, jsr, jer) & #ifndef without_wet_deposition !$OMP shared ( cvsfac,cp_scale ) & #ifdef with_budgets !$OMP shared ( buddep_dat,ra ) & #endif #endif !$OMP private (i, j, k,kk,l,n, jss, jee, li, ld, nglob, nzv, zxi, scavf, & !$OMP l_sum_wet, conv1, lbdcv, f, fu,amu, amd, zrm, zrxm, & !$OMP zrym, zwetrm, zcnvd, srm, srxm, srym, srzm, cp_value) call my_omp_domain (jsr(region), jer(region), jss, jee) l_sum_wet = 0.0 do j=jss,jee do i=isr(region),ier(region) #ifdef with_zoom if (zoomed(i,j) /= region) cycle #endif ! calculate convection matrix and directly apply to rm(i,j,*,1:ntracet) ! Set all amu's, amd's entu' etc to 0 for gridcells without clouds ! ! updraft - amu is positive upward ! f(:,:)=0. fu(:,:)=0. amu(:)=0. li = cloud_top(i,j) do k=1,li amu(k)=amu(k-1)+entu(i,j,k)-detu(i,j,k) if ( amu(k) > 0. ) then ! ! we limit zxi from below at 0. in case ! there are inconsistent en/detrainment ! rates. The case zxi>1 should not occur ! since eu and du are >=0. ! zxi=max(0.,1.-detu(i,j,k)/(amu(k-1)+entu(i,j,k))) else ! set massflux at upper boundary to strictly 0. amu(k)=0. zxi=0. end if do kk=1,k-1 fu(k,kk)=fu(k-1,kk)*zxi !all previous levels end do !cmk: note: the fu(k-1,k) values are zero (see manual Heimann) fu(k,k)=entu(i,j,k)/m(i,j,k)*zxi end do !k ! ! downdraft - amd is negative downward ! amd(:)=0. ld = cloud_lfs(i,j) do k=ld,2,-1 amd(k-1)=amd(k)-entd(i,j,k)+detd(i,j,k) if ( amd(k-1) < 0. ) then zxi=max(0.,1.+detd(i,j,k)/(amd(k)-entd(i,j,k))) else amd(k-1)=0. zxi=0. end if do kk=k+1,ld ! f(k-1,kk)=f(k,kk)*zxi !f is here only fd downdraftmatrix end do !kk f(k-1,k)=-entd(i,j,k)/m(i,j,k)*zxi !note the negatives for j end do !k ! ! add coefficients from up and downdraft together ! do k=1,lmax_conv-1 do kk=1,lmax_conv f(k,kk)=fu(k,kk)+f(k,kk) end do !kk ! add coefficient from subsidence !CMK SH f(k,k+1)=f(k,k+1)-(amu(k)+amd(k))/m(i,j,k+1) ! reported by Sander Houweling......... f(k,k+1)=f(k,k+1)-amu(k)/m(i,j,k+1) f(k,k )=f(k,k )-amd(k)/m(i,j,k ) ! #ifndef without_diffusion ! diffusion terms ! no diffusion in the stratosphere k>lmax_conv? f(k,k ) = f(k,k ) + dkg(i,j,k)/m(i,j,k ) f(k,k+1) = f(k,k+1) - dkg(i,j,k)/m(i,j,k+1) #endif end do ! >>> wet removal >>> ! initialize wet-removal lbdcv(:,:) = 0.0 #ifndef without_wet_deposition ! fill scaveging stuff if necessary ! ! WARNING: For matrices f(k,kk), fu(k,kk), and fd(k,kk), index k ! corresponds to the upper boundary of level k ! and index j to level j. ! For matrices g(i,j,k,kk) and lbdcv(i,j,k,kk), ! both indices k and j ! correspond to respective levels and not to boundaries. ! cp_value = cp_dat(region)%data(i,j,1) ! rain in m/s in region do k=2,lmax_conv do kk=1,k-1 !cmk lbdcv(k,kk)=sceffdeep*cp_eff(region)%surf(i,j)* & !cmk dt*(fu(k-1,kk) - fu(k,kk)) ! scale factor between 0 - 1, with 1 for cp > cp_scaler !CMK OLD lbdcv(k,kk)=sceffdeep*min(cp_value, cp_scaler)/cp_scaler * & lbdcv(k,kk)=sceffdeep*(1.0 - exp(-cp_value/cp_scale))* & dt*(fu(k-1,kk) - fu(k,kk)) enddo !kk enddo !k ! CMK end ! ! which ensures that the rate of scavenging is equal to dn/dt=-s Mu n ! ! TvN (8 July 2004) #endif ! <<< !---- WG end 21/04/98 ----- ! generate forward matrix do k=1,lmax_conv do kk=1,lmax_conv conv1(k,kk)=-dt*(f(k-1,kk)-f(k,kk)) end do !kk conv1(k,k) = conv1(k,k) + 1.0 end do !k call fastminv(conv1,lmax_conv) ! some checks...can be removed later! !if(okdebug) then ! minvalue = min(minvalue,minval(conv1)) ! do l=1,lmax_conv ! mxval = max(mxval,sum(conv1(:,l))) ! mnval = min(mnval,sum(conv1(:,l))) ! end do ! ! ! test deviations for a uniform mixing ratio of 1.....ppb ! ! do l=1,lmax_conv ! new_mass(l) = dot_product(conv1(l,:),m(i,j,1:lmax_conv)*1e-9) ! end do ! mxvald = max(mxvald, maxval(abs((m(i,j,1:lmax_conv)*1e-9-new_mass) ! /(m(i,j,1:lmax_conv)*1e-9)))) !end if !okdebug ! ! directly apply the matrix: ! copy i,j if(revert==-1) conv1 = transpose(conv1) ! loop over local tracers do n = 1, ntracetloc ! global tracer index: nglob=n+offsetn ! copy tracer concentrations for this column in srm/srxm/etc: do l=1,lmax_conv srm(l) = rm(i,j,l,n) #ifdef slopes srxm(l) = rxm(i,j,l,n) srym(l) = rym(i,j,l,n) srzm(l) = rzm(i,j,l,n) #ifdef secmom srxxm(l) = rxxm(i,j,l,n) srxym(l) = rxym(i,j,l,n) srxzm(l) = rxzm(i,j,l,n) sryym(l) = ryym(i,j,l,n) sryzm(l) = ryzm(i,j,l,n) srzzm(l) = rzzm(i,j,l,n) #endif #endif end do ! init budgets to zero: do l=1,lmax_conv zwetrm(l) = 0.0 zcnvd (l) = 0.0 end do ! init column of concentrations to zero: do l=1,lmax_conv zrm(l)=0. #ifdef slopes zrxm(l) = 0.0 zrym(l) = 0.0 #ifdef secmom zrxxm(l) = 0.0 zrxym(l) = 0.0 zrxzm(l) = 0.0 zryym(l) = 0.0 zryzm(l) = 0.0 #endif #endif end do ! apply convection matrix: do l=1,lmax_conv #ifdef slopes srzm(l)=conv1(l,l)*srzm(l) !rzm is directly calculated #ifdef secmom srzzm(l)=conv1(l,l)*srzzm(l) !rzzm is directly calculated #endif #endif do k=1,lmax_conv ! species dependent removal.... #ifndef without_wet_deposition scavf = min( lbdcv(l,k)*cvsfac(nglob), conv1(l,k) ) #else scavf = 0.0 #endif ! for budgets: zcnvd (l) = zcnvd (l) + conv1(l,k) * srm(k) zwetrm(l) = zwetrm(l) + scavf * srm(k) ! loss is positive ! store new concentrations in zrm/zrxm/etc ! using old concentrations in srm/srxm/etc zrm(l) = zrm(l) + (conv1(l,k)-scavf) * srm(k) #ifdef slopes zrxm (l) = zrxm (l) + (conv1(l,k)-scavf) * srxm(k) zrym (l) = zrym (l) + (conv1(l,k)-scavf) * srym(k) #ifdef secmom zrxxm(l) = zrxxm(l) + (conv1(l,k)-scavf) * srxxm(k) zrxym(l) = zrxym(l) + (conv1(l,k)-scavf) * srxym(k) zrxzm(l) = zrxzm(l) + (conv1(l,k)-scavf) * srxzm(k) zryym(l) = zryym(l) + (conv1(l,k)-scavf) * sryym(k) zryzm(l) = zryzm(l) + (conv1(l,k)-scavf) * sryzm(k) #endif #endif end do !k ! substract original tracer mass ! VH, TvN (2007-11-16) zcnvd(l)=zcnvd(l)-srm(l) end do !l ! restore concentrations: do l=1,lmax_conv rm(i,j,l,n)=zrm(l) #ifdef slopes rxm(i,j,l,n)=zrxm(l) rym(i,j,l,n)=zrym(l) rzm(i,j,l,n)=srzm(l) #ifdef secmom rxxm(i,j,l,n)=zrxxm(l) rxym(i,j,l,n)=zrxym(l) rxzm(i,j,l,n)=srxzm(l) ryym(i,j,l,n)=zryym(l) ryzm(i,j,l,n)=zryzm(l) rzzm(i,j,l,n)=srzzm(l) #endif #endif end do #ifdef with_budgets #ifndef without_wet_deposition if (cvsfac(nglob).gt.0.) then do l=1,lmax_conv nzv=nzon_vg(l) buddep_dat(region)%cp(i,j,nzv,nglob) = & buddep_dat(region)%cp(i,j,nzv,nglob) + & zwetrm(l)/ra(nglob)*1.e3 !kg=> moles if ( nglob == 1 ) l_sum_wet = l_sum_wet + & zwetrm(l) !in kg! end do !l end if #endif #endif #ifdef with_tendencies ! Add conv/diff and wet deposition budgets in kg to tendencies: do l = 1, lmax_conv call PLC_AddValue_t( region, plc_ipr_tcnvd, i, j, l, nglob, & zcnvd(l)*plc_kg_from_tm(nglob), status ) IF_NOTOK_RETURN(status=1) call PLC_AddValue_t( region, plc_ipr_lwdep, i, j, l, nglob, & zwetrm(l)*plc_kg_from_tm(nglob), status ) IF_NOTOK_RETURN(status=1) end do #endif end do !ntracet end do !i end do !j #if defined (with_budgets) call my_omp_sum (l_sum_wet, g_sum_wet) #endif !$OMP END PARALLEL #ifndef without_wet_deposition #if defined (with_budgets) sum_wetdep(region) = sum_wetdep(region) + g_sum_wet #endif #endif !if (okdebug) then ! ! call mpi_allreduce(minvalue,minvalue_all,1,my_real, & ! mpi_min,com_trac,ierr) ! call mpi_allreduce(mxval,mxval_all,1,my_real,mpi_max,com_trac,ierr) ! call mpi_allreduce(mnval,mnval_all,1,my_real,mpi_min,com_trac,ierr) ! call mpi_allreduce(mxvald,mxvald_all,1,my_real,mpi_max,com_trac,ierr) ! call barrier_t ! if (myid==root_t) then ! print* ,'-----convection matrix information----------------------' ! print* ,'Minimum value of convection matrix :',minvalue_all ! print *,'Max/Min along the columns (?/=1) :',mxval_all,mnval_all ! !these should be one/zero to ensure mass-conservation... ! print *,'Maximum fractional deviation from uniform mixing ', & ! 'ratio column:',mxvald_all ! print* ,'------convection matrix information---------------------' ! end if ! call barrier_t ! if (mxvald_all > 0.1) & ! call exitus(1, ' Convection matrix appears wrong') !end if nullify(m) 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(zoomed) nullify(entu) nullify(detu) nullify(entd) nullify(detd) nullify(cloud_top) nullify(cloud_base) nullify(cloud_lfs) #ifndef without_diffusion nullify(dkg) #endif ! ok status = 0 end subroutine convec #ifdef with_lapack !----------------------------------------------------------------------- ! compute inverse matrix using LU decompositon and forward-backward ! substitution without pivoting (for easy vectorization) ! ! LAPACK version, WP, Jan-2006 !----------------------------------------------------------------------- subroutine fastminv(a,lm) use toolbox, only: escape_tm ! input/output integer,intent(in) :: lm real,dimension(lm,lm),intent(inout) :: a ! local integer :: istat real,dimension(lm) :: work real,dimension(lm) :: ipiv call DGETRF( lm, lm, a, lm, ipiv, istat ) if(istat.ne.0) call escape_tm('fastminv: Failed to make LU decomposition ') call DGETRI( lm, a, lm, ipiv, work, lm, istat ) if(istat.ne.0) call escape_tm('fastminv: Failed to make inverse matrix') end subroutine fastminv #else !----------------------------------------------------------------------- ! compute inverse matrix using LU decompositon and forward-backward ! substitution without pivoting (for easy vectorization) ! ! Source: EMR and NR, sec. edition (p. 35ff). ! ! mh, {Sat, Jun 4, 1994} 16:27:15 !----------------------------------------------------------------------- subroutine fastminv(a,lm) ! input/output integer,intent(in) :: lm real,dimension(lm,lm),intent(inout) :: a ! local real,dimension(lm) :: y real,dimension(lm) :: b real,dimension(lm,lm) :: bb integer :: i,j,k do j=1,lm do i=1,j bb(i,j)=a(i,j) do k=1,i-1 bb(i,j)=bb(i,j)-bb(i,k)*bb(k,j) end do end do do i=j+1,lm bb(i,j)=a(i,j) do k=1,j-1 bb(i,j)=bb(i,j)-bb(i,k)*bb(k,j) end do bb(i,j)=bb(i,j)/bb(j,j) end do end do do k=1,lm do i=1,lm b(i)=0. end do b(k)=1. do i=1,lm y(i)=b(i) do j=1,i-1 y(i)=y(i)-bb(i,j)*y(j) end do end do do i=lm,1,-1 a(i,k)=y(i) do j=i+1,lm a(i,k)=a(i,k)-bb(i,j)*a(j,k) end do a(i,k)=a(i,k)/bb(i,i) end do end do end subroutine fastminv #endif end module convection