!### 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 ! 17 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !--------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CONVECTION_INIT ! ! !DESCRIPTION: set %used flag of releveant meteo to True !\\ !\\ ! !INTERFACE: ! SUBROUTINE CONVECTION_INIT( status ) ! ! !USES: ! use dims , only : nregions use Meteo , only : Set use Meteodata , only : entu_dat, entd_dat, detu_dat, detd_dat use Meteodata , only : mfw_dat, omega_dat #ifndef without_wet_deposition use Meteodata , only : cp_dat #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Convection_Init' integer :: region ! 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 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CONVECTION_DONE ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE CONVECTION_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Convection_Done' ! nothing to do status = 0 END SUBROUTINE CONVECTION_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CONVEC ! ! !DESCRIPTION: ! !**** 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 !\\ !\\ ! !INTERFACE: ! SUBROUTINE CONVEC( region, status ) ! ! !USES: ! use dims, only : tref !, dx, dy, xref, yref use dims, only : itau, okdebug, kdebug, nconv, lmax_conv use dims, only : nregions use dims, only : revert use global_data, only : region_dat, conv_dat, mass_dat use meteodata, only : entu_dat, entd_dat, detu_dat, detd_dat use meteodata, only : cp_dat use meteodata, only : m_dat use chem_param, only : ra, ntracet #ifdef with_budgets #ifndef without_wet_deposition use budget_global, only : nzon_vg 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, plc_ipr_lwdep, plc_ipr_tcnvd, plc_kg_from_tm #endif use tm5_distgrid, only : dgrid, Get_DistGrid, reduce use partools, only : isRoot ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer,intent(out) :: status ! ! !REVISION HISTORY: ! 17 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/CONVEC' 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, i1, i2, j1, j2 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 ! for global wetdep budget of tracer #1 real :: g_sum_wet real, allocatable :: l_sum_wet(:,:) !====================== start ================================== 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;; !-------- local indices call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) #ifdef with_budgets #ifndef without_wet_deposition allocate(l_sum_wet(i1:i2,j1:j2)) l_sum_wet = 0.0 #endif #endif !-------- shortcuts m => m_dat(region)%data rm => mass_dat(region)%rm #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 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) & #endif !$OMP shared (region, ntracetloc, sceffdeep, dt, & !$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, li, ld, n, nzv, zxi, scavf, & !$OMP conv1, lbdcv, f, fu,amu, amd, zrm, zrxm, & !$OMP zrym, zwetrm, zcnvd, srm, srxm, srym, srzm, cp_value) do j = j1, j2 do i = i1, i2 #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 TRACER : do n = 1, ntracet ! 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 srxm(l)=conv1(l,l)*srxm(l) !no convection -> no multip. !WP-Emma van der Veen srym(l)=conv1(l,l)*srym(l) !no convection -> no multip. !WP-Emma van der Veen #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(n), 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 !WP-vdV! zrxm (l) = zrxm (l) + (conv1(l,k)-scavf) * srxm(k) !WP-vdV! 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 !WP-vdV! rxm(i,j,l,n)=zrxm(l) !WP-vdV! rym(i,j,l,n)=zrym(l) rxm(i,j,l,n)=srxm(l) !VDV rym(i,j,l,n)=srym(l) !VDv 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(n).gt.0.) then do l=1,lmax_conv nzv=nzon_vg(l) buddep_dat(region)%cp(i,j,nzv,n) = buddep_dat(region)%cp(i,j,nzv,n) + & zwetrm(l)/ra(n)*1.e3 !kg=> moles if ( n == 1 ) l_sum_wet(i,j) = l_sum_wet(i,j) + 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( region, plc_ipr_tcnvd, i, j, l, n, & zcnvd(l)*plc_kg_from_tm(n), status ) IF_NOTOK_RETURN(status=1) call PLC_AddValue( region, plc_ipr_lwdep, i, j, l, n, & zwetrm(l)*plc_kg_from_tm(n), status ) IF_NOTOK_RETURN(status=1) end do #endif end do TRACER !ntracet end do !i end do !j !$OMP END PARALLEL #ifndef without_wet_deposition #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 #endif !--- Cleanup 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 !EOC #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