
MODULE mod_chem_plumerise_scalar



  USE mem_chem1, ONLY:  &
       chem1_src_vars      ! Type
  
  USE mem_aer1, ONLY:   &
       aer1_vars           ! Type
  
  USE mem_plume_chem1, ONLY: &
                             plume_mean_vars,  &     ! Type
                             plume_fre_vars ,  &
            		     iflam_frac     ,  &
            		     imean_frp      ,  &
            		     istd_frp       ,  &
            		     imean_size     ,  &
            		     istd_size	      
!  USE mod_htint, ONLY: &
!       htint               ! Subroutine

!  USE mod_therm_lib, ONLY: &
!       esat                ! Function

 use node_mod, only:  &
       mynum,         &
       mchnum,        &
       master_num
     
  IMPLICIT NONE


  PRIVATE


  INTEGER, PARAMETER :: nkp       = 200
  INTEGER, PARAMETER :: ntime     = 200
  INTEGER            :: n_setgrid = 0

  REAL :: dz
  REAL :: zt (nkp)
  REAL :: zm (nkp)
  REAL :: dzt(nkp)
  REAL :: dzm(nkp)


  PUBLIC :: plumerise_driver   ! Subroutine



CONTAINS



  !----------------------------------------------------------------------------
  ! Plume rise model for vegetation fires (CPTEC/INPE 2005-2006)
  ! Refs.:
  ! Freitas, S. R., K. M. Longo, R. Chatfield, D. Latham, M. A. F. Silva Dias,
  !  M. O. Andreae, 
  ! E. Prins, J. C. Santos, R. Gielow and J. A. Carvalho Jr.: Including the 
  ! sub-grid scale 
  ! plume rise of vegetation fires in low resolution atmospheric transport 
  ! models. 
  !  Atmospheric Chemistry and Physics and discussion, 6, 2006.
  !-
  ! Freitas, S. R.; Longo, K. M.; M. Andreae. Impact of including the plume
  ! rise of vegetation 
  ! fires in numerical simulations of associated atmospheric pollutants. 
  ! Geophys. Res. Lett., 
  ! 33, L17808, doi:10.1029/2006GL026608, 2006. 
  !----------------------------------------------------------------------------

  SUBROUTINE plumerise_driver(mzp,mxp,myp,ia,iz,ja,jz,                        &
                              srctime1,chem_nspecies,spc_chem_alloc,src,      &
                              on,off,nmodes,aer_nspecies,spc_aer_alloc,       &
                              aer_bburn,nsrc,bburn,nveg_agreg,tropical_forest,&
                              boreal_forest,savannah,grassland,nzpmax,dtlt,   &
                              time,zt_rams,zm_rams,dzt_rams,g,cp,cpor,p00,    &
                              rgas,theta,pp,pi0,rv,up,vp,rtgt,lpw,pr_time,    &
                              aerosol,chem1_src_g,aer1_g,plume_mean_g,spc_aer_name,&
			      plume_fre_g,                                    &
			      plumerise_flag                                  )

    INTEGER,                INTENT(IN)    :: mzp
    INTEGER,                INTENT(IN)    :: mxp
    INTEGER,                INTENT(IN)    :: myp
    INTEGER,                INTENT(IN)    :: ia
    INTEGER,                INTENT(IN)    :: iz
    INTEGER,                INTENT(IN)    :: ja
    INTEGER,                INTENT(IN)    :: jz
    INTEGER,                INTENT(IN)    :: plumerise_flag
    ! grid_dims
    INTEGER ,               INTENT(IN)    :: nzpmax
    ! mem_basic
    REAL,                   INTENT(IN)    :: theta(:,:,:)
    REAL,                   INTENT(IN)    :: pp   (:,:,:)
    REAL,                   INTENT(IN)    :: pi0  (:,:,:)
    REAL,                   INTENT(IN)    :: rv   (:,:,:)
    REAL,                   INTENT(IN)    :: up   (:,:,:)!srf-AWE
    REAL,                   INTENT(IN)    :: vp   (:,:,:)!srf-AWE
    ! mem_grid
    REAL,                   INTENT(IN)    :: rtgt(:,:)
    INTEGER,                INTENT(IN)    :: lpw (:,:)
    REAL,                   INTENT(IN)    :: dtlt
    REAL,                   INTENT(IN)    :: time

    REAL,                   INTENT(IN)    :: zt_rams (:)
    REAL,                   INTENT(IN)    :: zm_rams (:)
    REAL,                   INTENT(IN)    :: dzt_rams(:)
    ! chem_sources
    REAL,                   INTENT(IN)    :: pr_time
    REAL,                   INTENT(IN)    :: srctime1
    ! chem1_list
    INTEGER,                INTENT(IN)    :: chem_nspecies
!!$    INTEGER,                INTENT(IN)    :: spc_chem_alloc(6,chem_nspecies)
    INTEGER,                INTENT(IN)    :: spc_chem_alloc(:,:)
    INTEGER,                INTENT(IN)    :: src
    INTEGER,                INTENT(IN)    :: on
    INTEGER,                INTENT(IN)    :: off
    ! aer1_list
    INTEGER,                INTENT(IN)    :: nmodes
    INTEGER,                INTENT(IN)    :: aer_nspecies
!!$    INTEGER,                INTENT(IN)    :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER,                INTENT(IN)    :: spc_aer_alloc(:,:,:)
    INTEGER,                INTENT(IN)    :: aer_bburn
    ! mem_chem1
    INTEGER,                INTENT(IN)    :: nsrc
!!$    TYPE (chem1_src_vars),  INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies)
    TYPE (chem1_src_vars),  INTENT(INOUT) :: chem1_src_g(:,:,:)
    INTEGER,                INTENT(IN)    :: bburn
    ! mem_aer1
!!$    TYPE (aer1_vars),       INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies)
    TYPE (aer1_vars),       INTENT(INOUT) :: aer1_g(:,:)
    INTEGER,                INTENT(IN)    :: aerosol
    ! mem_plume_chem1
    TYPE (plume_mean_vars), INTENT(IN)    :: plume_mean_g(:)
    TYPE (plume_fre_vars ), INTENT(INOUT) :: plume_fre_g (5)
    INTEGER,                INTENT(IN)    :: nveg_agreg
    INTEGER,                INTENT(IN)    :: tropical_forest
    INTEGER,                INTENT(IN)    :: boreal_forest
    INTEGER,                INTENT(IN)    :: savannah
    INTEGER,                INTENT(IN)    :: grassland
    ! rconstants
    REAL,                   INTENT(IN)    :: g
    REAL,                   INTENT(IN)    :: cp
    REAL,                   INTENT(IN)    :: cpor
    REAL,                   INTENT(IN)    :: p00
    REAL,                   INTENT(IN)    :: rgas
    CHARACTER(LEN=8) , INTENT(IN) :: spc_aer_name(nmodes,aer_nspecies)

    LOGICAL, PARAMETER :: IS2PRINT = .FALSE. !.TRUE.
    CHARACTER(LEN=8)   :: ctime
    CHARACTER(LEN=4)   :: cmynum
    INTEGER            :: nrec

    CHARACTER(LEN=*), PARAMETER :: h="**(plumerise_driver)**"
    CHARACTER(LEN=8) :: c0

    IF(MOD(time+.001,pr_time).LE.(dtlt).OR.time.LT..01 .OR. &
       ABS(time-srctime1).LT. 1.e-5 )THEN

!!$       WRITE(c0,"(f8.1)") time
!!$       CALL MsgOne(h,' -----------------------------------------------------')
!!$       CALL MsgOne(h,' plumerise tendencies updated')
!!$       CALL MsgOne(h,' time = '//TRIM(ADJUSTL(c0))//' seconds')
!!$       CALL MsgOne(h,' -----------------------------------------------------')      

       CALL plumerise(mzp,mxp,myp,ia,iz,ja,jz, &
                      theta,pp,pi0,rv,up,vp,rtgt,lpw,zt_rams,zm_rams, &
                      nzpmax,dzt_rams,chem_nspecies,spc_chem_alloc,   &
                      src,on,off,nmodes,aer_nspecies,spc_aer_alloc,   &
                      aer_bburn,nsrc,chem1_src_g,bburn,aer1_g,aerosol,&
                      plume_mean_g,nveg_agreg,tropical_forest,        &
                      boreal_forest,savannah,grassland,g,cp,cpor,     &
                      p00,rgas,spc_aer_name,                          &          
		      plume_fre_g,                                    &
                      plumerise_flag                                       )

    ENDIF

  END SUBROUTINE plumerise_driver
  !-------------------------------------------------------------------------


  SUBROUTINE plumerise(m1,m2,m3,ia,iz,ja,jz, &
                       theta,pp,pi0,rv,up,vp,rtgt,lpw,zt_rams,zm_rams,     &
                       nzpmax,dzt_rams,chem_nspecies,spc_chem_alloc,src,   &
                       on,off,nmodes,aer_nspecies,spc_aer_alloc,aer_bburn, &
                       nsrc,chem1_src_g,bburn,aer1_g,aerosol,plume_mean_g, &
                       nveg_agreg,tropical_forest,boreal_forest,savannah,  &
                       grassland,g,cp,cpor,p00,rgas,spc_aer_name,          &         
		       plume_fre_g,                                        &
                       plumerise_flag                                           )

    INTEGER,               INTENT(IN)     :: m1
    INTEGER,               INTENT(IN)     :: m2
    INTEGER,               INTENT(IN)     :: m3
    INTEGER,               INTENT(IN)     :: ia
    INTEGER,               INTENT(IN)     :: iz
    INTEGER,               INTENT(IN)     :: ja
    INTEGER,               INTENT(IN)     :: jz
    INTEGER,               INTENT(IN)     :: plumerise_flag
    REAL,                  INTENT(IN)	  :: theta(:,:,:)
    REAL,                  INTENT(IN)	  :: pp   (:,:,:)
    REAL,                  INTENT(IN)	  :: pi0  (:,:,:)
    REAL,                  INTENT(IN)	  :: rv   (:,:,:)
    REAL,                  INTENT(IN)	  :: up   (:,:,:)!srf-AWE
    REAL,                  INTENT(IN)	  :: vp   (:,:,:)!srf-AWE
    REAL,                  INTENT(IN)	  :: rtgt(:,:)
    INTEGER,               INTENT(IN)	  :: lpw (:,:)
    REAL,                  INTENT(IN)     :: zt_rams(:)
    REAL,                  INTENT(IN)     :: zm_rams(:)
    ! grid_dims
    INTEGER,               INTENT(IN)     :: nzpmax
    ! mem_grid
!!$    REAL,                  INTENT(IN)     :: dzt_rams(nzpmax,m2,m3)
    REAL,                  INTENT(IN)     :: dzt_rams(:)
    ! chem1_list
    INTEGER,               INTENT(IN)     :: chem_nspecies
!!$    INTEGER,               INTENT(IN)     :: spc_chem_alloc(6,chem_nspecies)
    INTEGER,               INTENT(IN)     :: spc_chem_alloc(:,:)
    INTEGER,               INTENT(IN)     :: src
    INTEGER,               INTENT(IN)     :: on
    INTEGER,               INTENT(IN)     :: off
    ! aer1_list
    INTEGER,               INTENT(IN)     :: nmodes
    INTEGER,               INTENT(IN)     :: aer_nspecies
!!$    INTEGER,               INTENT(IN)     :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER,               INTENT(IN)     :: spc_aer_alloc(:,:,:)
    INTEGER,               INTENT(IN)     :: aer_bburn
    ! mem_chem1
    INTEGER,               INTENT(IN)     :: nsrc
!!$    TYPE (chem1_src_vars), INTENT(INOUT)  :: chem1_src_g(2,nsrc,chem_nspecies)
    TYPE (chem1_src_vars), INTENT(INOUT)  :: chem1_src_g(:,:,:)
    INTEGER,               INTENT(IN)     :: bburn

    ! mem_aer1
!!$    TYPE (aer1_vars),       INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies)
    TYPE (aer1_vars),       INTENT(INOUT) :: aer1_g(:,:)
    INTEGER,                INTENT(IN)    :: aerosol
    ! mem_plume_chem1
    TYPE (plume_mean_vars), INTENT(IN)    :: plume_mean_g(:)
    TYPE (plume_fre_vars) , INTENT(INOUT) :: plume_fre_g(5)
    INTEGER,                INTENT(IN)    :: nveg_agreg
    INTEGER,                INTENT(IN)    :: tropical_forest
    INTEGER,                INTENT(IN)    :: boreal_forest
    INTEGER,                INTENT(IN)    :: savannah
    INTEGER,                INTENT(IN)    :: grassland
    ! rconstants
    REAL,                   INTENT(IN)    :: g
    REAL,                   INTENT(IN)    :: cp
    REAL,                   INTENT(IN)    :: cpor
    REAL,                   INTENT(IN)    :: p00
    REAL,                   INTENT(IN)    :: rgas

    CHARACTER(LEN=8) , INTENT(IN) :: spc_aer_name(nmodes,aer_nspecies)

    INTEGER            :: i, j, k, ixx, iveg_ag, imm, k1, k2, ispc, imode, iloop
    INTEGER            :: kmt
    REAL               :: burnt_area, STD_burnt_area, dz_flam, dzi, FRP, convert_smold_to_flam
    REAL               :: ztopmax(2)
    REAL               :: W_VMD(nkp,2), VMD(nkp,2)
    REAL               :: q_smold_kgm2
    INTEGER            :: it1, it2
    INTEGER, PARAMETER :: use_last = 0
    INTEGER, PARAMETER :: izprint  = 0 ! if = 0 => no printout
    INTEGER, PARAMETER :: wind_eff = 1  

    REAL    :: area           
    INTEGER :: maxtime        
    REAL    :: alpha          
    REAL    :: rsurf          
    REAL    :: fmoist         
    REAL    :: tdur           
    REAL    :: zsurf          
    REAL    :: heating(ntime) 
    REAL    :: thtcon(nkp)    
    REAL    :: rvcon (nkp)    
    REAL    :: picon (nkp)    
    REAL    :: zcon  (nkp)    
    REAL    :: zzcon (nkp)    
    REAL    :: qv    (nkp)    
    REAL    :: qh    (nkp)    
    REAL    :: qi    (nkp)    
    REAL    :: qc    (nkp)    
    REAL    :: txs   (nkp)    
    REAL    :: cvi   (nkp)    
    REAL    :: vth   (nkp)    
    REAL    :: w     (nkp)    
    REAL    :: t     (nkp)    
    REAL    :: qsat  (nkp)    
    REAL    :: rho   (nkp)    
    REAL    :: radius(nkp)    
    REAL    :: visc  (nkp)    
    REAL    :: wc    (nkp)    
    REAL    :: wt    (nkp)    
    REAL    :: est   (nkp)    
    REAL    :: pe    (nkp)    
    REAL    :: te    (nkp)    
    REAL    :: tt    (nkp)    
    REAL    :: qvenv (nkp)    
    REAL    :: cvh   (nkp)    
    REAL    :: vti   (nkp)    
    REAL    :: pke   (nkp)    
    REAL    :: the   (nkp)    
    REAL    :: thve  (nkp)    
    REAL    :: dne   (nkp)    
    REAL    :: sc    (nkp)     
    REAL    :: sct   (nkp)     
    REAL    :: ucon  (nkp)    
    REAL    :: vcon  (nkp)    
    REAL    :: upe   (nkp)
    REAL    :: vpe   (nkp)
    REAL    :: vel_e (nkp)   
    REAL    :: vel_p (nkp)
    REAL    :: rad_p (nkp)
    REAL    :: vel_t (nkp)
    REAL    :: rad_t (nkp)   

    !- for biomass burn, only one memory allocation (only one time level)
    it1=1

    !- zera o termo fonte associado `as emissoes com plumerise (k>2)	
    !- chemistry section
    DO ispc = 1,chem_nspecies
       IF(spc_chem_alloc(src,ispc) == off) CYCLE
       chem1_src_g(it1,bburn,ispc)%sc_src(3:m1,:,:) = 0.
    ENDDO

    !- aerosol section
    IF(aerosol == 1) THEN 
       DO imode=1,nmodes
          !- only for biomass burning aerosols
          IF(spc_aer_alloc(src,imode,aer_bburn)== off ) CYCLE    
          aer1_g(imode,aer_bburn)%sc_src(3:m1,:,:) = 0.
       ENDDO
    ELSEIF(aerosol == 2) THEN ! for MATRIX and only for boc_ocar and boc_bcar species
		
       DO ispc = 1,aer_nspecies
	 DO imode=1,nmodes
          IF(spc_aer_alloc(src,imode,ispc) == off ) CYCLE

	  !-only for bburn aerosols) 
	  IF(spc_aer_name(imode,ispc)=="boc_bcar" .or. &
	     spc_aer_name(imode,ispc)=="boc_ocar" )	THEN
            aer1_g(imode,ispc)%sc_src(3:m1,:,:) = 0.
	  ENDIF
         ENDDO
       ENDDO	  
    ENDIF

    DO j = ja,jz
       DO i = ia,iz
          convert_smold_to_flam=0.0
	  
          !- if the max value of flaming is close to zero => there is not
          !  emission with plume rise => cycle
	  IF(PLUMERISE_flag == 1 ) THEN 
	    IF(plume_mean_g(tropical_forest)%flam_frac(i,j) + &
	       plume_mean_g(boreal_forest  )%flam_frac(i,j) + &
	       plume_mean_g(savannah       )%flam_frac(i,j) + &
	       plume_mean_g(grassland      )%flam_frac(i,j) < 1.e-6 ) CYCLE
            !- loop over the four types of aggregate biomes with fires
	    iloop = nveg_agreg
	       
	  ELSEIF(PLUMERISE_flag == 2 ) THEN
	     IF(plume_fre_g(iflam_frac)%pvar(i,j) < 1.e-6 ) CYCLE
             iloop = 1
	  
	  ENDIF
          DO k = 1,m1
             ucon  (k)=up(k,i,j)		! u wind (m/s)
             vcon  (k)=vp(k,i,j)		! v wind (m/s)
            !wcon  (k)=wp(k,i,j)		! w wind (m/s)
             thtcon(k)=theta(k,i,j)		! pot temperature (K)
             picon (k)=(pp(k,i,j)+pi0(k,i,j))   ! exner function
            !tmpcon(k)=thtcon(k)*picon(k)/cp    ! temperature (K)
            !dncon (k)=dn0(k,i,j)		! dry air density (basic state) (kg/m3)
            !prcon (k)=(picon(k)/cp)**cpor*p00  ! pressure (Pa)
             rvcon (k)=rv(k,i,j)		! water vapor mixing ratio (kg/kg)

             zcon  (k)=zt_rams(k) *rtgt(i,j)    ! termod-point height (m)
             zzcon (k)=zm_rams(k) *rtgt(i,j)    ! W-point height      (m)

	     zsurf = 0.
	     
          ENDDO

          !- get envinronmental state (temp, water vapor mix ratio, ...)
          CALL get_env_condition(lpw(i,j),m1-1,kmt,qvenv,pke,the,te,pe,thve,dne,upe,vpe,vel_e, &
	                         zcon,zt,thtcon,rvcon,picon,ucon,vcon,zm,dzm,dzt,dz,zsurf,     &
                                 wind_eff,g,cp,cpor,p00,rgas)

          !- only for testing (using an external sounding)
          !CALL get_env_condition_sound(lpw(i,j),m1-1,kmt,qvenv,pke,the,te,pe,thve,dne,upe,vpe,vel_e &
	  !			 ,zcon,zt,thtcon,rvcon,picon,ucon,vcon,zm,dzm,dzt,dz,zsurf,n_setgrid,wind_eff)


          !- loop over the four types of aggregate biomes with fires
          DO iveg_ag=1,iloop 

             !- verifica a existencia de emissao flaming para um bioma especifico
             IF(PLUMERISE_flag == 1 ) THEN 
	       IF(plume_mean_g(iveg_ag)%flam_frac(i,j) < 1.e-6 ) CYCLE 
	       
	       !-burnt area and standard deviation
               !burnt_area   = 50.*1e4           !m^2, only for testing with a sounding
                burnt_area   = plume_mean_g(iveg_ag)%fire_size(i,j)              
               
	        STD_burnt_area= 0.!not em use
	       
	       !-number to calculate the flaming emission from the amount emitted
	       !-during the smoldering phase
		convert_smold_to_flam=plume_mean_g(iveg_ag)%flam_frac(i,j)
             
	     ELSEIF(PLUMERISE_flag == 2 ) THEN 
	       !-number to calculate the emission during the flaming pahse 
	       !-from the amount emitted during the smoldering phase
	
		convert_smold_to_flam=plume_fre_g(iflam_frac)%pvar(i,j)
               !
	       !- check if there is only one fire in a given grid box (=> std =0.)
	       if(plume_fre_g(istd_frp )%pvar(i,j) < 1.0e-6) then
	       
	       !- if yes, we will set it as a 20% of the mean frp as a gross estimation
	       !- of the retrieval uncertainty by the sensors.
	       !- (we are not taking care about the fire size retrieval)
	          !print*,"xx=",I,J,plume_fre_g(istd_frp )%pvar(i,j),0.2*plume_fre_g(imean_frp )%pvar(i,j)
	          plume_fre_g(istd_frp )%pvar(i,j)=0.2*plume_fre_g(imean_frp )%pvar(i,j)
	       endif	       	       
	     ENDIF

             !- loop over the minimum and maximum heat fluxes/FRP
             DO imm=1,2
                !--------------------
                !ixx=iveg_ag*10 + imm
                ixx=0+imm
                !--------------------             
                
		IF(PLUMERISE_flag == 2 ) THEN
		  if(imm==1 ) then
                    !for imm = 1 => lower injection height
                    burnt_area = max(1.0e4,plume_fre_g(imean_size)%pvar(i,j) - 0.5*plume_fre_g(istd_size)%pvar(i,j))
		    FRP        = max(1000.,plume_fre_g(imean_frp )%pvar(i,j) - 0.5*plume_fre_g(istd_frp )%pvar(i,j))
                 
		  elseif(imm==2 ) then
                    !for imm = 2 => higher injection height
                    burnt_area = max(1.0e4,plume_fre_g(imean_size)%pvar(i,j) + 0.5*plume_fre_g(istd_size)%pvar(i,j))
		    FRP        = max(1000.,plume_fre_g(imean_frp )%pvar(i,j) + 0.5*plume_fre_g(istd_frp )%pvar(i,j))
                  endif
		ENDIF
		! print*,"imm",imm,plume_fre_g(imean_size)%pvar(i,j), plume_fre_g(istd_size)%pvar(i,j),&
		!                    plume_fre_g(imean_frp )%pvar(i,j), plume_fre_g(istd_frp )%pvar(i,j)
                !
		!-just for testing...
		!burnt_area =  burnt_area*30.
                !FRP	    =  FRP*600.    
		!-----------
		!print*,"size(ha)/frp=",i,j,imm, burnt_area/1.e+4, 0.88*1000.*FRP/(burnt_area) ; call flush(6)
        	!write(*,100)"size(ha)/frp(kW/m2)",i,j,imm, burnt_area/1.e+4, 0.88*1.e-3*FRP/(burnt_area) 
                !100 format (1x,A20,1x,3I4,1x,2F15.4)
		    
                !- get fire properties (burned area, plume radius, heating rates ...)
                CALL get_fire_properties(imm,iveg_ag,burnt_area,use_last, &
                                         heating,area,maxtime,alpha,      &
                                         rsurf,fmoist,tdur,               &
					 plumerise_flag,FRP)

                !------  generates the plume rise    ------

		IF(PLUMERISE_flag == 1 ) THEN
                  !-- only one value for eflux of GRASSLAND
                  IF(iveg_ag == GRASSLAND .AND. imm == 2) THEN 
                     ztopmax(2)=ztopmax(1)
                     ztopmax(1)=zzcon(1)
                     CYCLE
                  ENDIF
                ENDIF

                CALL makeplume(i,j,kmt,ztopmax(imm),ixx,imm,use_last,w,           &
                               t,pe,qsat,est,rho,zt,te,txs,wc,dz,maxtime,zm,  &
                               wt,dzm,tt,qv,qh,qi,qc,qvenv,radius,alpha,visc, &
                               cvi,vth,cvh,vti,dzt,sc,area,rsurf,tdur,heating,&
                               fmoist,vel_e,vel_p,rad_p,vel_t,rad_t,W_VMD,    &
                               izprint)     	      

             ENDDO ! enddo of the loop imm
             
             !- defines the vertical layer where the flaming phase emission will
             !  be released in the 3d atmospheric model
             CALL set_flam_vert(ztopmax,k1,k2,zzcon,W_VMD,VMD)

             !- thickness of the vertical layer
             dz_flam=zzcon(k2)-zzcon(k1-1)

             !print*,"   z1/z2=",k1,k2,zzcon(k1-1),zzcon(k2)!,convert_smold_to_flam ; call flush(6)
             !print*,"=================================================" ; call flush(6)

             !- emission during flamming phase is evenly distributed between levels k1 and k2 
             DO k=k1,k2
                !use this in case the emission src is already in mixing ratio
                !rhodzi= 1./(dn0(k,i,j) * dz_flam)
                
		!use this in case the emission src is tracer density
                dzi= 1./( dz_flam)

                !- chemistry section 
                DO ispc = 1, chem_nspecies
                   IF(spc_chem_alloc(src,ispc) == off ) CYCLE

                   !- get back the smoldering emission in kg/m2  (actually in
                   !   1e-9 kg/m2) 

                   !use this in case the emission src is already in mixing 
                   ! ratio
                   !q_smold_kgm2 = (rtgt(i,j)/dzt_rams(2) *  dn0(2,i,j)	)*   &
                   ! 	         chem1_src_g(it1,bburn,ispc,ng)%sc_src(2,i,j)

                   !use this in case the emission src is tracer density
                   q_smold_kgm2 = (rtgt(i,j)/dzt_rams(2)                   )*   &
                                   chem1_src_g(it1,bburn,ispc)%sc_src(2,i,j)	 

                   ! units = already in ppbm,  don't need "fcu" factor 
                   chem1_src_g(it1,bburn,ispc)%sc_src(k,i,j) = chem1_src_g(it1,bburn,ispc)%sc_src(k,i,j) +&
                                                               convert_smold_to_flam                     *&
                                                              !plume_mean_g(iveg_ag)%flam_frac(i,j)      *&
                                                               q_smold_kgm2                              *&
                                                               dzi    !use this in case the emission src is tracer density
                                                              !rhodzi !use this in case the emission src is already in mixing ratio
                   
                ENDDO

                !- aerosol section ( only for biomass burning aerosols )
                IF(aerosol == 1) THEN 

                   DO imode=1,nmodes

                      IF(spc_aer_alloc(src,imode,aer_bburn)== off ) CYCLE    

                      !- get back the smoldering emission in kg/m2  
                      !  (actually in 1e-9 kg/m2) 

                      !use this in case the emission src is already in 
                      !mixing ratio
                      !q_smold_kgm2 = (rtgt(i,j)/dzt_rams(2) * dn0(2,i,j) )*  &
                      !	               aer1_g(imode,aer_bburn)%sc_src(2,i,j)    

                      !use this in case the emission src is tracer density

                      q_smold_kgm2 = (rtgt(i,j)/dzt_rams(2)                 )*   &
                                      aer1_g(imode,aer_bburn)%sc_src(2,i,j)	 


                      ! units = already in ppbm,  don't need "fcu" factor 

                      aer1_g(imode,aer_bburn)%sc_src(k,i,j) = aer1_g(imode,aer_bburn)%sc_src(k,i,j) +&
                                                              convert_smold_to_flam                 *&
                                                             !plume_mean_g(iveg_ag)%flam_frac(i,j)  *&
                                                              q_smold_kgm2                          *&
                                                              dzi    !use this in case the emission src is tracer density
                                                             !rhodzi !use this in case the emission src is already in mixing ratio

                   ENDDO
                
		ELSEIF(aerosol == 2) THEN ! for MATRIX and only for boc_ocar and boc_bcar species
		
		   DO ispc = 1,aer_nspecies
	              DO imode=1,nmodes
                       IF(spc_aer_alloc(src,imode,ispc) == off ) CYCLE

		       !-only for bburn aerosols - MATRIX MECH=1 
	               IF(spc_aer_name(imode,ispc)=="boc_bcar" .or. &
	                  spc_aer_name(imode,ispc)=="boc_ocar" )    THEN
                        
                          q_smold_kgm2 = (rtgt(i,j)/dzt_rams(2))*   &
                                         aer1_g(imode,ispc)%sc_src(2,i,j)	 
                          
			  ! units = already in ppbm,  don't need "fcu" factor 
                          aer1_g(imode,ispc)%sc_src(k,i,j)= aer1_g(imode,ispc)%sc_src(k,i,j)     +&
                                                           !plume_mean_g(iveg_ag)%flam_frac(i,j) *&
                                                            convert_smold_to_flam                *&
                                                            q_smold_kgm2 * &
                                                            dzi     !use this in case the emission src is tracer density
                                                            !rhodzi !use this in case the emission src is already in mixing ratio
 
		       ENDIF
                      ENDDO
                   ENDDO
		ENDIF

             ENDDO

          ENDDO ! enddo do loop em iveg_ag
       ENDDO  ! loop em i
    ENDDO   ! loop em j
!stop "plumerise 400"
  END SUBROUTINE plumerise
  !-------------------------------------------------------------------------

  SUBROUTINE get_env_condition(k1,k2,kmt,qvenv,pke,the,te,pe,thve,dne,   &
                               upe,vpe,vel_e,zcon,zt,thtcon,rvcon,picon, &
                               ucon,vcon,zm,dzm,dzt,dz,zsurf,            &! n_setgrid,    &
                               wind_eff,g,cp,cpor,p00,rgas)

    INTEGER, INTENT(IN)    :: k1
    INTEGER, INTENT(IN)    :: k2
    INTEGER, INTENT(INOUT) :: kmt 
    INTEGER, INTENT(IN)    :: wind_eff
    ! plumegen_coms
    REAL,    INTENT(INOUT) :: dz
    REAL,    INTENT(IN)    :: zsurf
!   INTEGER, INTENT(INOUT) :: n_setgrid
!!$    REAL,    INTENT(INOUT) :: qvenv (nkp)
!!$    REAL,    INTENT(INOUT) :: pke   (nkp)
!!$    REAL,    INTENT(INOUT) :: the   (nkp)
!!$    REAL,    INTENT(INOUT) :: te    (nkp)
!!$    REAL,    INTENT(INOUT) :: pe    (nkp)
!!$    REAL,    INTENT(INOUT) :: thve  (nkp)
!!$    REAL,    INTENT(INOUT) :: dne   (nkp)  
!!$    REAL,    INTENT(INOUT) :: upe   (nkp)
!!$    REAL,    INTENT(INOUT) :: vpe   (nkp)
!!$    REAL,    INTENT(INOUT) :: vel_e (nkp)
    REAL,    INTENT(INOUT) :: zcon  (nkp)
!!$    REAL,    INTENT(INOUT) :: zt    (nkp)
!!$    REAL,    INTENT(INOUT) :: thtcon(nkp)
!!$    REAL,    INTENT(INOUT) :: rvcon (nkp)
!!$    REAL,    INTENT(INOUT) :: picon (nkp)
!!$    REAL,    INTENT(INOUT) :: ucon  (nkp)
!!$    REAL,    INTENT(INOUT) :: vcon  (nkp)
!!$    REAL,    INTENT(INOUT) :: zm    (nkp)
!!$    REAL,    INTENT(INOUT) :: dzm   (nkp)
!!$    REAL,    INTENT(INOUT) :: dzt   (nkp)
    REAL,    INTENT(INOUT) :: qvenv (:)
    REAL,    INTENT(INOUT) :: pke   (:)
    REAL,    INTENT(INOUT) :: the   (:)
    REAL,    INTENT(INOUT) :: te    (:)
    REAL,    INTENT(INOUT) :: pe    (:)
    REAL,    INTENT(INOUT) :: thve  (:)
    REAL,    INTENT(INOUT) :: dne   (:)  
    REAL,    INTENT(INOUT) :: upe   (:)
    REAL,    INTENT(INOUT) :: vpe   (:)
    REAL,    INTENT(INOUT) :: vel_e (:)
!!$    REAL,    INTENT(INOUT) :: zcon  (:)
    REAL,    INTENT(INOUT) :: zt    (:)
    REAL,    INTENT(INOUT) :: thtcon(:)
    REAL,    INTENT(INOUT) :: rvcon (:)
    REAL,    INTENT(INOUT) :: picon (:)
    REAL,    INTENT(INOUT) :: ucon  (:)
    REAL,    INTENT(INOUT) :: vcon  (:)
    REAL,    INTENT(INOUT) :: zm    (:)
    REAL,    INTENT(INOUT) :: dzm   (:)
    REAL,    INTENT(INOUT) :: dzt   (:)
    ! rconstants
    REAL,    INTENT(IN)    :: g
    REAL,    INTENT(IN)    :: cp
    REAL,    INTENT(IN)    :: cpor
    REAL,    INTENT(IN)    :: p00
    REAL,    INTENT(IN)    :: rgas


    CHARACTER(LEN=*), PARAMETER :: h="**(get_env_condition)**"

    INTEGER :: k,nk
    REAL    :: znz
    LOGICAL :: found
    
    IF( n_setgrid == 0) THEN
       n_setgrid = 1
       ! define vertical grid of plume model
       CALL set_grid(zt,zm,dzm,dzt,dz,zsurf) 
       ! zt(k) =  thermo and water levels
       ! zm(k) =  dynamical levels 
    ENDIF
    znz=zcon(k2)
    
    
    found = .FALSE.
    DO k=nkp,1,-1
       IF(zt(k).LT.znz) THEN
          found = .TRUE.
          exit
       END IF
    ENDDO
    
    IF (.not. found) THEN
       STOP ' ERROR: Envir stop 12 - chem_plumerise_scalar'
!!$       CALL FatalError(h//" Envir stop 12")
    END IF

!--(DMK-BRAMS-5-INI)------------------------------------------------------
    kmt=min(k,nkp-1)
!--(DMK-BRAMS-5-OLD)------------------------------------------------------
!    kmt=k
!--(DMK-BRAMS-5-FIM)------------------------------------------------------

    nk=k2-k1+1

    IF (nk > nkp) THEN
       STOP ' ERROR: nk > nkp - chem_plumerise_scalar'
    END IF

!!$    CALL htint(nk,  ucon,zcon(k1:nk+k1-1),kmt,upe  ,zt)
!!$    CALL htint(nk,  vcon,zcon(k1:nk+k1-1),kmt,vpe  ,zt)
!!$    CALL htint(nk,thtcon,zcon(k1:nk+k1-1),kmt,the  ,zt)
!!$    CALL htint(nk, rvcon,zcon(k1:nk+k1-1),kmt,qvenv,zt)
    CALL htint(nk,  ucon,zcon,kmt,upe  ,zt)
    CALL htint(nk,  vcon,zcon,kmt,vpe  ,zt)
    CALL htint(nk,thtcon,zcon,kmt,the  ,zt)
    CALL htint(nk, rvcon,zcon,kmt,qvenv,zt)
    DO k=1,kmt
       qvenv(k)=MAX(qvenv(k),1e-8)
    ENDDO

    pke(1)=picon(1)
    DO k=1,kmt
       thve(k)=the(k)*(1.+.61*qvenv(k)) ! virtual pot temperature
    ENDDO
    DO k=2,kmt
       pke(k)=pke(k-1)-g*2.*(zt(k)-zt(k-1))  & ! exner function
            /(thve(k)+thve(k-1))
    ENDDO
    DO k=1,kmt
       te(k)  = the(k)*pke(k)/cp         ! temperature (K) 
       pe(k)  = (pke(k)/cp)**cpor*p00    ! pressure (Pa)
       dne(k)= pe(k)/(rgas*te(k)*(1.+.61*qvenv(k))) !  dry air density (kg/m3)
       
       vel_e(k) = SQRT(upe(k)**2+vpe(k)**2)         !-env wind (m/s)
       !print*,'k,vel_e(k),te(k)=',vel_e(k),te(k)
    ENDDO
    
    !-ewe - env wind effect
    IF(wind_eff < 1)  vel_e(1:kmt) = 0.
    

    !-use este para gerar o RAMS.out
    ! ------- print environment state
    !print*,'k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000'
    !do k=1,kmt
    ! write(*,100)  k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000.
    !enddo
    !stop 333


    !--------- nao eh necessario este calculo
    !do k=1,kmt
    !  call thetae(pe(k),te(k),qvenv(k),thee(k))
    !enddo


    !--------- converte press de Pa para kPa para uso modelo de plumerise

    DO k=1,kmt
       pe(k) = pe(k)*1.e-3
    ENDDO

    RETURN 


    !para testes    ----------------------
    !
    !      open(13,file='prn.out')
    !      do 112, i=1,nkp
    !!  112 read (13,1000) zt(i),pe(i),te(i),rhe(i),dne(i),qvenv(i)
    !  112 read (13,1000) dummy,pe(i),te(i),rhe(i),dne(i),qvenv(i)
    !1000 format(1x,6F15.4)
    !      close(13)
    !
    !te(:)=te(:)+273.15
    !--------------------------------------

  END SUBROUTINE get_env_condition
  !-------------------------------------------------------------------------

  SUBROUTINE set_grid(zt,zm,dzm,dzt,dz,zsurf)

    ! plumegen_coms
    REAL, INTENT(INOUT) :: dz ! (DMK) alterado (OUT) para (INOUT)
    REAL, INTENT(IN)    :: zsurf
    REAL, INTENT(INOUT) :: zt (:)
    REAL, INTENT(INOUT) :: zm (:)
    REAL, INTENT(INOUT) :: dzm(:)
    REAL, INTENT(INOUT) :: dzt(:)

    INTEGER :: k, mzp

    dz=100. ! set constant grid spacing of plume grid model(meters)

    mzp=nkp
    zt(1) = zsurf
    zm(1) = zsurf
    zt(2) = zt(1) + 0.5*dz
    zm(2) = zm(1) + dz
    DO k=3,mzp
       zt(k) = zt(k-1) + dz ! thermo and water levels
       zm(k) = zm(k-1) + dz ! dynamical levels	
    ENDDO
    !print*,zsurf
    !Print*,zt(:)
    DO k = 1,mzp-1
       dzm(k) = 1. / (zt(k+1) - zt(k))
    ENDDO
    dzm(mzp)=dzm(mzp-1)

    DO k = 2,mzp
       dzt(k) = 1. / (zm(k) - zm(k-1))
    ENDDO
    dzt(1) = dzt(2) * dzt(2) / dzt(3)

    !   dzm(1) = 0.5/dz
    !   dzm(2:mzp) = 1./dz
    RETURN
  END SUBROUTINE set_grid
  !-------------------------------------------------------------------------

  SUBROUTINE set_flam_vert(ztopmax,k1,k2,zzcon,W_VMD,VMD)

!!$    REAL,    INTENT(IN)  :: ztopmax(2)
    REAL,    INTENT(IN)  :: ztopmax(:)
    INTEGER, INTENT(OUT) :: k1
    INTEGER, INTENT(OUT) :: k2
    ! plumegen_coms
!!$    REAL,    INTENT(IN)  :: zzcon(nkp)
    REAL,    INTENT(IN)  :: zzcon(:)
    !- version 2
!!$    REAL,    INTENT(IN)  :: W_VMD(nkp,2)
!!$    REAL,    INTENT(OUT) :: VMD(nkp,2)
    REAL,    INTENT(IN)  :: W_VMD(:,:)
    REAL,    INTENT(OUT) :: VMD(:,:)


    INTEGER :: imm,k
    INTEGER :: k_lim(2)
    REAL    :: w_thresold, xxx
    INTEGER :: k_initial, k_final, ko, kk4, kl

    !- version 1
    DO imm=1,2
       ! checar 
       !    do k=1,m1-1
       DO k=1,nkp-1
          IF(zzcon(k) > ztopmax(imm) ) EXIT
       ENDDO
       k_lim(imm) = k
    ENDDO
    k1=MAX(3,k_lim(1))
    k2=MAX(3,k_lim(2))

    IF(k2 < k1) THEN
       !print*,'1: ztopmax k=',ztopmax(1), k1
       !print*,'2: ztopmax k=',ztopmax(2), k2
       k2=k1
       !stop 1234
    ENDIF
    
    !- version 2    
    !- vertical mass distribution
    !- 
    w_thresold = 1.
    DO imm=1,2

    
       VMD(1:nkp,imm)= 0.
       xxx=0.
       k_initial= 0
       k_final  = 0
    
       !- define range of the upper detrainemnt layer
       DO ko=nkp-10,2,-1
     
        IF(w_vmd(ko,imm) < w_thresold) CYCLE
     
        IF(k_final==0) k_final=ko
     
        IF(w_vmd(ko,imm)-1. > w_vmd(ko-1,imm)) THEN
          k_initial=ko
          EXIT
        ENDIF
      
       ENDDO
       !- if there is a non zero depth layer, make the mass vertical 
       !  distribution 
       IF(k_final > 0 .AND. k_initial > 0) THEN 
       
           k_initial=INT((k_final+k_initial)*0.5)
       
           !- parabolic vertical distribution between k_initial and k_final
           kk4 = k_final-k_initial+2
           DO ko=1,kk4-1
               kl=ko+k_initial-1
               VMD(kl,imm) = 6.* float(ko)/float(kk4)**2 * (1. - float(ko)/float(kk4))
           ENDDO
	   IF(SUM(VMD(1:NKP,imm)) .NE. 1.) THEN
 	       xxx= ( 1.- SUM(VMD(1:NKP,imm)) )/float(k_final-k_initial+1)
 	       DO ko=k_initial,k_final
 	         VMD(ko,imm) = VMD(ko,imm)+ xxx !- values between 0 and 1.
 	       ENDDO
               ! print*,'new mass=',sum(mass)*100.,xxx
               !pause
           ENDIF
        ENDIF !k_final > 0 .and. k_initial > 

    ENDDO
    
  END SUBROUTINE set_flam_vert

  !-------------------------------------------------------------------------

  SUBROUTINE get_fire_properties(imm,iveg_ag,burnt_area,use_last,              &
                                 heating,area,maxtime,alpha,rsurf,fmoist,tdur, &
				 plumerise_flag,FRP)

    INTEGER, INTENT(IN)    :: imm
    INTEGER, INTENT(IN)    :: iveg_ag
    REAL,    INTENT(IN)    :: burnt_area,FRP
    INTEGER, INTENT(IN)    :: use_last
    INTEGER, INTENT(IN)    :: plumerise_flag
    ! plumegen_coms
    REAL,    INTENT(INOUT) :: heating(:)        
    REAL,    INTENT(INOUT) :: area              
    INTEGER, INTENT(INOUT) :: maxtime           
    REAL,    INTENT(OUT)   :: alpha
    REAL,    INTENT(OUT)   :: rsurf
    REAL,    INTENT(OUT)   :: fmoist
    REAL,    INTENT(OUT)   :: tdur


    CHARACTER(LEN=*), PARAMETER :: h="**(get_fire_properties)**"

    INTEGER :: moist, i, icount
    REAL    :: bfract, effload, heat, hinc, heat_fluxW
    REAL    :: heat_flux(2,4)

    INTEGER :: mdur           ! (DMK) scratch, not used
    REAL    :: bload          ! (DMK) scratch, not used

    DATA heat_flux/  &
         !---------------------------------------------------------------------
         !  heat flux      !IGBP Land Cover	    ! 
         ! min  ! max      !Legend and		    ! reference
         !    kW/m^2       !description  	    ! 
         !--------------------------------------------------------------------
         30.0,	 80.0,   &! Tropical Forest         ! igbp 2 & 4
         30.0,   80.0,   &! Boreal forest           ! igbp 1 & 3
         4.4,	 23.0,   &! cerrado/woody savanna   | igbp  5 thru 9
         3.3,	  3.3    /! Grassland/cropland      ! igbp 10 thru 17
         !--------------------------------------------------------------------
    real, parameter :: beta = 0.88  !ref.: Paugam et al., 2015
   !real, parameter :: beta = 5.0   !ref.: Wooster et al., 2005


    !-- fire at the surface
    !area = 20.e+4   ! area of burn, m^2
    area = burnt_area! area of burn, m^2

    IF ( PLUMERISE_flag == 1) THEN
    	!fluxo de calor para o bioma
	heat_fluxW = heat_flux(imm,iveg_ag) * 1000. ! converte para W/m^2
	
    ELSEIF ( PLUMERISE_flag == 2) THEN
	! "beta" factor converts FRP to convective energy
        heat_fluxW = beta*(FRP/area)/0.55 ! in W/m^2

    ENDIF

    mdur = 53        ! duration of burn, minutes
    bload = 10.      ! total loading, kg/m**2 
    moist = 10       ! fuel moisture, %. average fuel moisture,percent dry
    maxtime =mdur+2  ! model time, min
    !maxtime =mdur-1  ! model time, min

    !heat = 21.e6    !- joules per kg of fuel consumed                   
    !heat = 15.5e6   !joules/kg - cerrado
    heat = 19.3e6    !joules/kg - floresta em alta floresta (mt)
    !alpha = 0.1      !- entrainment constant
    alpha = 0.05      !- entrainment constant

    !-------------------- printout ----------------------------------------

    !!WRITE ( * ,  * ) ' SURFACE =', ZSURF, 'M', '  LCL =', ZBASE, 'M'  
    !
    !PRINT*,'======================================================='
    !print * , ' FIRE BOUNDARY CONDITION   :'  
    !print * , ' DURATION OF BURN, MINUTES =',MDUR  
    !print * , ' AREA OF BURN, HA	      =',AREA*1.e-4
    !print * , ' HEAT FLUX, kW/m^2	      =',heat_fluxW*1.e-3
    !print * , ' TOTAL LOADING, KG/M**2    =',BLOAD  
    !print * , ' FUEL MOISTURE, %	      =',MOIST !average fuel moisture,percent dry
    !print * , ' MODEL TIME, MIN.	      =',MAXTIME  
    !
    !
    !
    ! ******************** fix up inputs *********************************
    !

    !IF (MOD (MAXTIME, 2) .NE.0) MAXTIME = MAXTIME+1  !make maxtime even

    MAXTIME = MAXTIME * 60  ! and put in seconds
    !
    RSURF = SQRT (AREA / 3.14159) !- entrainment surface radius (m)

    FMOIST   = MOIST / 100.       !- fuel moisture fraction
    !
    !
    ! calculate the energy flux and water content at lboundary.
    ! fills heating() on a minute basis. could ask for a file at this po
    ! in the program. whatever is input has to be adjusted to a one
    ! minute timescale.
    !

    DO I = 1, ntime         !- make sure of energy release
       HEATING (I) = 0.0001  !- avoid possible divide by 0
    ENDDO
    !                                  
    TDUR = MDUR * 60.       !- number of seconds in the burn

    bfract = 1.             !- combustion factor

    EFFLOAD = BLOAD * BFRACT  !- patchy burning

    !     spread the burning evenly over the interval
    !     except for the first few minutes for stability
    ICOUNT = 1  
    !
    IF(MDUR > NTIME) &
         STOP 'Increase time duration (ntime) in min - see file "plumerise_mod.f90"'
!!$         CALL FatalError(h//" Increase time duration (ntime) in min - see file chem_plumerise_scalar.f90")

    DO WHILE (ICOUNT.LE.MDUR)                             
       !  HEATING (ICOUNT) = HEAT * EFFLOAD / TDUR  ! W/m**2 
       !  HEATING (ICOUNT) = 80000.  * 0.55         ! W/m**2 

       HEATING (ICOUNT) = heat_fluxW  * 0.55     ! W/m**2 (0.55 converte para energia convectiva)
       ICOUNT = ICOUNT + 1  
    ENDDO
    !     ramp for 5 minutes
    IF(use_last /= 1) THEN

       HINC = HEATING (1) / 4.  
       HEATING (1) = 0.1  
       HEATING (2) = HINC  
       HEATING (3) = 2. * HINC  
       HEATING (4) = 3. * HINC  
    ELSE
       IF(imm==1) THEN
          HINC = HEATING (1) / 4.  
          HEATING (1) = 0.1  
          HEATING (2) = HINC  
          HEATING (3) = 2. * HINC  
          HEATING (4) = 3. * HINC 
       ELSE 
          stop "check units and use heat_fluxW"
          HINC = (HEATING (1) - heat_flux(imm-1,iveg_ag) * 1000. *0.55)/ 4.
          HEATING (1) = heat_flux(imm-1,iveg_ag) * 1000. *0.55 + 0.1  
          HEATING (2) = HEATING (1)+ HINC  
          HEATING (3) = HEATING (2)+ HINC  
          HEATING (4) = HEATING (3)+ HINC 
       ENDIF
    ENDIF
    !srf-24jan2007 - end!<<<<<<<<<<<<<<<<<<<<02082007

    RETURN
  END SUBROUTINE get_fire_properties

  !----------------------------------------------------------------------------
  !
  SUBROUTINE MAKEPLUME(i,j,kmt,ztopmax,ixx,imm,use_last,w,t,pe,qsat, &
                       est,rho,zt,te,txs,wc,dz,maxtime,zm,wt,dzm,tt,qv,qh, &
                       qi,qc,qvenv,radius,alpha,visc,cvi,vth,cvh,vti,dzt,  &
                       sc,area,rsurf,tdur,heating,fmoist,vel_e,vel_p,rad_p,&
                       vel_t,rad_t,W_VMD,izprint)  
    !
    ! *********************************************************************
    !
    !    EQUATION SOURCE--Kessler Met.Monograph No. 32 V.10 (K)
    !    Alan Weinstein, JAS V.27 pp 246-255. (W),
    !    Ogura and Takahashi, Monthly Weather Review V.99,pp895-911 (OT)
    !    Roger Pielke,Mesoscale Meteorological Modeling,Academic Press,1984
    !    Originally developed by: Don Latham (USFS)
    !
    !
    ! ************************ VARIABLE ID ********************************
    !
    !     DT=COMPUTING TIME INCREMENT (SEC)
    !     DZ=VERTICAL INCREMENT (M)
    !     LBASE=LEVEL ,CLOUD BASE
    !
    !     CONSTANTS:
    !       G = GRAVITATIONAL ACCELERATION 9.80796 (M/SEC/SEC).
    !       R = DRY AIR GAS CONSTANT (287.04E6 JOULE/KG/DEG K)
    !       CP = SPECIFIC HT. (1004 JOULE/KG/DEG K)
    !       HEATCOND = HEAT OF CONDENSATION (2.5E6 JOULE/KG)
    !       HEATFUS = HEAT OF FUSION (3.336E5 JOULE/KG)
    !       HEATSUBL = HEAT OF SUBLIMATION (2.83396E6 JOULE/KG)
    !       EPS = RATIO OF MOL.WT. OF WATER VAPOR TO THAT OF DRY AIR (0.622)
    !       DES = DIFFERENCE BETWEEN VAPOR PRESSURE OVER WATER AND ICE (MB)
    !       TFREEZE = FREEZING TEMPERATURE (K)
    !
    !
    !     PARCEL VALUES:
    !       T = TEMPERATURE (K)
    !       TXS = TEMPERATURE EXCESS (K)
    !       QH = HYDROMETEOR WATER CONTENT (G/G DRY AIR)
    !       QHI = HYDROMETEOR ICE CONTENT (G/G DRY AIR)
    !       QC = WATER CONTENT (G/G DRY AIR)
    !       QVAP = WATER VAPOR MIXING RATIO (G/G DRY AIR)
    !       QSAT = SATURATION MIXING RATIO (G/G DRY AIR)
    !       RHO = DRY AIR DENSITY (G/M**3) MASSES = RHO*Q'S IN G/M**3
    !       ES = SATURATION VAPOR PRESSURE (kPa)
    !
    !     ENVIRONMENT VALUES:
    !       TE = TEMPERATURE (K)
    !       PE = PRESSURE (kPa)
    !       QVENV = WATER VAPOR (G/G)
    !       RHE = RELATIVE HUMIDITY FRACTION (e/esat)
    !       DNE = dry air density (kg/m^3)
    !
    !     HEAT VALUES:
    !       HEATING = HEAT OUTPUT OF FIRE (WATTS/M**2)
    !       MDUR = DURATION OF BURN, MINUTES
    !
    !       W = VERTICAL VELOCITY (M/S)
    !       RADIUS=ENTRAINMENT RADIUS (FCN OF Z)
    !	RSURF = ENTRAINMENT RADIUS AT GROUND (SIMPLE PLUME, TURNER)
    !	ALPHA = ENTRAINMENT CONSTANT
    !       MAXTIME = TERMINATION TIME (MIN)
    !
    !
    !**********************************************************************
    !**********************************************************************

    INTEGER, INTENT(IN)    :: kmt,i,j
    REAL,    INTENT(INOUT) :: ztopmax  ! (DMK) alterado (OUT) para (INOUT)
    INTEGER, INTENT(IN)    :: ixx
    INTEGER, INTENT(IN)    :: imm
    INTEGER, INTENT(IN)    :: use_last
    INTEGER, INTENT(IN)    :: izprint
    ! plumegen_coms
    REAL,    INTENT(IN)    :: dz
    INTEGER, INTENT(IN)    :: maxtime
    REAL,    INTENT(IN)    :: alpha
    REAL,    INTENT(IN)    :: area
    REAL,    INTENT(IN)    :: rsurf
    REAL,    INTENT(IN)    :: tdur
    REAL,    INTENT(IN)    :: fmoist
!!$    REAL,    INTENT(IN)    :: heating(ntime)
!!$    REAL,    INTENT(INOUT) :: w     (nkp)
!!$    REAL,    INTENT(INOUT) :: t     (nkp)
!!$    REAL,    INTENT(IN)    :: pe    (nkp)
!!$    REAL,    INTENT(INOUT) :: qsat  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: est   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: rho   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: zt    (nkp)
!!$    REAL,    INTENT(IN)    :: te    (nkp)
!!$    REAL,    INTENT(INOUT) :: txs   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: wc    (nkp)
!!$    REAL,    INTENT(IN)    :: zm    (nkp)
!!$    REAL,    INTENT(INOUT) :: wt    (nkp)
!!$    REAL,    INTENT(IN)    :: dzm   (nkp)
!!$    REAL,    INTENT(INOUT) :: tt    (nkp)
!!$    REAL,    INTENT(INOUT) :: qv    (nkp)
!!$    REAL,    INTENT(INOUT) :: qh    (nkp)
!!$    REAL,    INTENT(INOUT) :: qi    (nkp)
!!$    REAL,    INTENT(INOUT) :: qc    (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv (nkp)
!!$    REAL,    INTENT(INOUT) :: radius(nkp)
!!$    REAL,    INTENT(INOUT) :: visc  (nkp)
!!$    REAL,    INTENT(INOUT) :: cvi   (nkp)
!!$    REAL,    INTENT(INOUT) :: vth   (nkp)
!!$    REAL,    INTENT(INOUT) :: cvh   (nkp)
!!$    REAL,    INTENT(INOUT) :: vti   (nkp)
!!$    REAL,    INTENT(IN)    :: dzt   (nkp)
!!$    REAL,    INTENT(IN)    :: sc    (nkp)
!!$!srf-awe    
!!$    REAL,    INTENT(IN)    :: vel_e (nkp)
!!$    REAL,    INTENT(INOUT) :: vel_p (nkp)
!!$    REAL,    INTENT(INOUT) :: rad_p (nkp)
!!$    REAL,    INTENT(INOUT) :: vel_t (nkp)
!!$    REAL,    INTENT(INOUT) :: rad_t (nkp)
!!$    REAL,    INTENT(OUT)   :: W_VMD (nkp,2)
!!$!srf-awe    
    REAL,    INTENT(IN)    :: heating(:)
    REAL,    INTENT(INOUT) :: w     (:)
    REAL,    INTENT(INOUT) :: t     (:)
    REAL,    INTENT(IN)    :: pe    (:)
    REAL,    INTENT(INOUT) :: qsat  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: est   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: rho   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: zt    (:)
    REAL,    INTENT(IN)    :: te    (:)
    REAL,    INTENT(INOUT) :: txs   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: wc    (:)
    REAL,    INTENT(IN)    :: zm    (:)
    REAL,    INTENT(INOUT) :: wt    (:)
    REAL,    INTENT(IN)    :: dzm   (:)
    REAL,    INTENT(INOUT) :: tt    (:)
    REAL,    INTENT(INOUT) :: qv    (:)
    REAL,    INTENT(INOUT) :: qh    (:)
    REAL,    INTENT(INOUT) :: qi    (:)
    REAL,    INTENT(INOUT) :: qc    (:)
    REAL,    INTENT(IN)    :: qvenv (:)
    REAL,    INTENT(INOUT) :: radius(:)
    REAL,    INTENT(INOUT) :: visc  (:)
    REAL,    INTENT(INOUT) :: cvi   (:)
    REAL,    INTENT(INOUT) :: vth   (:)
    REAL,    INTENT(INOUT) :: cvh   (:)
    REAL,    INTENT(INOUT) :: vti   (:)
    REAL,    INTENT(IN)    :: dzt   (:)
    REAL,    INTENT(IN)    :: sc    (:)
!srf-awe    
    REAL,    INTENT(IN)    :: vel_e (:)
    REAL,    INTENT(INOUT) :: vel_p (:)
    REAL,    INTENT(INOUT) :: rad_p (:)
    REAL,    INTENT(INOUT) :: vel_t (:)
    REAL,    INTENT(INOUT) :: rad_t (:)
    REAL,    INTENT(OUT)   :: W_VMD (:,:)
!srf-awe    
    
 
    !logical :: endspace  
    INTEGER :: k, kk, kkmax, deltak, ilastprint, &
               nrectotal, i_micro, n_sub_step, nrec

    REAL    :: vc, g,  r,  cp,  eps,  &
               tmelt,  heatsubl,  heatfus,  heatcond, tfreeze, &
               wmax, rmaxtime, es,dt_save


    REAL, EXTERNAL :: esat

    REAL    :: DELZ_THRESOLD 

    CHARACTER (len=2) :: cixx	   
    
    !
    !
    ! ******************* SOME CONSTANTS **********************************
    !
    !      XNO=10.0E06 median volume diameter raindrop (K table 4)
    !      VC = 38.3/(XNO**.125) mean volume fallspeed eqn. (K)
    !
    PARAMETER (vc = 5.107387)  
    PARAMETER (g = 9.80796, r = 287.04, cp = 1004., eps = 0.622,  tmelt = 273.3)
    PARAMETER (heatsubl = 2.834e6, heatfus = 3.34e5, heatcond = 2.501e6)
    PARAMETER (tfreeze = 269.3)  
    !

    REAL    :: viscosity     ! (DMK) scratch
    REAL    :: tstpf         ! (DMK) scratch
    INTEGER :: mintime       ! (DMK) scratch
    REAL    :: ztop          ! (DMK) scratch
    REAL    :: dqsdz         ! (DMK) scratch
    REAL    :: time          ! (DMK) scratch
    INTEGER :: nm1           ! (DMK) scratch
    INTEGER :: l             ! (DMK) scratch
    REAL    :: adiabat       ! (DMK) scratch
    REAL    :: vhrel         ! (DMK) scratch
    REAL    :: wbar          ! (DMK) scratch
    REAL    :: dt            ! (DMK) scratch
    REAL    :: ztop_(ntime)  ! (DMK) scratch
    REAL    :: scr1(nkp)     ! (DMK) scratch
    REAL    :: qvt (nkp)     ! (DMK) scratch
    REAL    :: qct (nkp)     ! (DMK) scratch
    REAL    :: qht (nkp)     ! (DMK) scratch
    REAL    :: qit (nkp)     ! (DMK) scratch
  
    tstpf = 2.0     !- timestep factor
    viscosity = 100.!- viscosity constant (original value: 0.001)

    nrectotal=150
    nrec = 0
    !
    !*************** PROBLEM SETUP AND INITIAL CONDITIONS *****************
    mintime = 1  
    ztopmax = 0. 
    ztop    = 0. 
    time    = 0.  
    dt      = 1.
    wmax    = 1. 
    kkmax   = 10
    deltaK  = 20
    ilastprint=0
    L       = 1   ! L initialization
    wbar    = 0.
    W_VMD(:,:) = 0.
    
    !--- initialization
    ! CALL INITIAL(kmt)   
    DELZ_THRESOLD = 1.*dz
    IF(imm==1 .AND. use_last==1)  &
         CALL INITIAL(kmt,txs,w,t,wc,wt,qv,vth,vti,qh,qi,qc,est,qsat,rho, &
                      radius,visc,te,qvenv,rsurf,alpha,zt,viscosity,pe,   &
                      area,dt,l,wbar,dqsdz,cvi,tdur,heating,    &
                      mintime,fmoist,time,VEL_P,rad_p)
    IF(use_last/=1) &
         CALL INITIAL(kmt,txs,w,t,wc,wt,qv,vth,vti,qh,qi,qc,est,qsat,rho, &
                      radius,visc,te,qvenv,rsurf,alpha,zt,viscosity,pe,   &
                      area,dt,l,wbar,dqsdz,cvi,tdur,heating,    &
                      mintime,fmoist,time,VEL_P,rad_p)  
    !
    !--- initial print fields:
    IF (izprint.NE.0) THEN
       WRITE(cixx(1:2),'(i2.2)') ixx
       OPEN(2, file = 'debug.'//cixx//'.dat')  
       OPEN(19,file='plumegen.'//cixx//'.gra',         &
            form='unformatted',access='direct',status='unknown',  &
            recl=4*nrectotal)  !PC   
       !     recl=1*nrectotal) !sx6 e tupay
       CALL printout(izprint,nrectotal,mintime,dt,time,ztop, &
                     pe,t,te,qv,qsat,qc,qh, &
                     qi,zt,w,vth,sc,qvenv,nrec)
       ilastprint=2
    ENDIF

    ! ******************* model evolution ******************************
    rmaxtime = float(maxtime)
    !
    DO WHILE (TIME.LE.RMAXTIME)  !beginning of time loop

       !   do itime=1,120

       !-- set model top integration
       nm1 = MIN(kmt, kkmax + deltak)

       !-- set timestep
       !dt = (zm(2)-zm(1)) / (tstpf * wmax)  
       dt = MIN(5.,(zm(2)-zm(1)) / (tstpf * wmax))
       IF(dt < 0.25 .or. maxval(T(:))<0.) THEN
          ztop = 0.
          ztop_(mintime) = ztop
          ztopmax = 0.
          kkmax   = 1 
	  print*,"PRM_UNST=",mynum,i,j,dt
	  call flush(6)
          EXIT
       ENDIF

       !-- elapsed time, sec
       time = time+dt 
       !-- elapsed time, minutes                                      
       mintime = 1 + INT (time) / 60     
       wmax = 1.  !no zeroes allowed.
       !************************** BEGIN SPACE LOOP **************************

       !-- zerout all model tendencies
       CALL tend0_plumerise(nm1,wt,tt,qvt,qct,qht,qit,vel_t,rad_t)

       !-- bounday conditions (k=1)
       L=1
       wbar=W(1)
       CALL LBOUND(qh,qi,qc,w,t,wc,vth,vti,txs,visc,rho,qv,est,qsat,pe, &
                   alpha,area,rsurf,te,viscosity,dt,qvenv,l,wbar,dqsdz,cvi, &
                   tdur,heating, mintime,fmoist,time,VEL_P,rad_p)

       !-- dynamics for the level k>1 
       !-- W advection 
       !   call vel_advectc_plumerise(NM1,WC,WT,DNE,DZM)
       CALL vel_advectc_plumerise(NM1,WC,WT,RHO,DZM)

       !-- scalars advection 1
       CALL scl_advectc_plumerise(NM1,scr1,dt,w,wc,rho,dzm,zt,zm,dzt, &
                                  t,tt,qv,qvt,qc,qct,qi,qit,qh,qht,vel_p, &
                                  vel_t,rad_p,rad_t)

       !-- scalars advection 2
       !call scl_advectc_plumerise2('SC',NM1)

       !-- scalars entrainment, adiabatic
       CALL scl_misc(NM1,wbar,w,adiabat,alpha,radius,tt,t,te,qvt, &
                     qv,qvenv,qct,qc,qht,qh,qit,qi,vel_e,vel_p,vel_t, &
                     rad_p,rad_t)

       !-- scalars dinamic entrainment
       CALL  scl_dyn_entrain(NM1,wbar,w,adiabat,alpha,radius,tt,t,   &
                             te,qvt,qv,qvenv,qct,qc,qht,qh,qit,qi,vel_e, &
                             vel_p,vel_t,rad_p,rad_t)


       !-- gravity wave damping using Rayleigh friction layer fot T
       CALL damp_grav_wave(1,nm1,deltak,dt,zt,zm,w,t,tt,te)

       !-- microphysics
       !   goto 101 ! bypass microphysics
       dt_save=dt
       n_sub_step=3
       dt=dt/float(n_sub_step)

       DO i_micro=1,n_sub_step
          !-- sedim ?
          CALL fallpart(NM1,rho,vth,vhrel,w,cvh,vti,cvi,qh,qi,zm,qht,qit)
          !-- microphysics
          DO L=2,nm1-1
             WBAR    = 0.5*(W(L)+W(L-1))
             ES      = 0.1*ESAT (T(L))            !BLOB SATURATION VAPOR PRESSURE, EM KPA
             QSAT(L) = (EPS * ES) / (PE(L) - ES)  !BLOB SATURATION LWC G/G DRY AIR
             EST (L) = ES  
             RHO (L) = 3483.8 * PE (L) / T (L) ! AIR PARCEL DENSITY , G/M**3
             !srf18jun2005
             !	IF (W(L) .ge. 0.) DQSDZ = (QSAT(L  ) - QSAT(L-1)) / (ZT(L  ) -ZT(L-1))
             !	IF (W(L) .lt. 0.) DQSDZ = (QSAT(L+1) - QSAT(L  )) / (ZT(L+1) -ZT(L  ))
             IF (W(L) .GE. 0.) THEN 
                DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1 )-ZT(L-1))
             ELSE
                DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1) -ZT(L-1))
             ENDIF

             CALL waterbal(qc,l,qh,qi,qv,t,qsat,wbar,dqsdz,dt,rho,est,cvi)
          ENDDO
       ENDDO
       dt=dt_save
       !
101    CONTINUE
       !
       !-- W-viscosity for stability 
       CALL visc_W(nm1,kmt,zt,visc,zm,w,t,qv,qh,qc,qi,wt,tt,qvt,qct,qht,&
                   qit,vel_p,vel_t,rad_p,rad_t)

       !-- update scalars
       CALL update_plumerise(nm1,'S',wt,dt,tt,qvt,qct,qht,w,t,qv,qc,qh, &
                             qit,qi,vel_p,vel_t,rad_p,rad_t)

       CALL hadvance_plumerise(1,nm1,WC,W,mintime) 

       !-- Buoyancy
       CALL buoyancy_plumerise(nm1, t, te, qv, qvenv, qh, qi, qc, wt, scr1)

       !-- Entrainment 
       CALL entrainment(nm1,w,wt,radius,alpha,vel_p,vel_e)

       !-- update W
       CALL update_plumerise(nm1,'W',wt,dt,tt,qvt,qct,qht,w,t,qv,qc, &
                             qh,qit,qi, vel_p,vel_t,rad_p,rad_t)

       CALL hadvance_plumerise(2,nm1,WC,W,mintime) 


       !-- misc
       DO k=2,nm1
          !    pe esta em kpa  - esat do rams esta em mbar = 100 Pa = 0.1 kpa
          es       = 0.1*esat (t(k)) !blob saturation vapor pressure, em kPa
          !    rotina do plumegen calcula em kPa
          !    es    = esat_pr (t(k))  !blob saturation vapor pressure, em kPa
          qsat(k) = (eps * es) / (pe(k) - es)  !blob saturation lwc g/g dry air
          est (k) = es  
          txs (k) = t(k) - te(k)
          rho (k) = 3483.8 * pe (k) / t (k) ! air parcel density , g/m**3
          ! no pressure diff with radius

          IF((ABS(wc(k))).GT.wmax) wmax = ABS(wc(k)) ! keep wmax largest w
       ENDDO

       ! Gravity wave damping using Rayleigh friction layer for W
       CALL damp_grav_wave(2,nm1,deltak,dt,zt,zm,w,t,tt,te)
       !---


       !- update radius
       DO k=2,nm1
        radius(k) = rad_p(k)
       ENDDO

       !-- try to find the plume top (above surface height)
       kk = 1
       DO WHILE (w (kk) .GT. 1.)  
          kk = kk + 1  
          ztop =  zm(kk) 
          !print*,'W=',w (kk)
       ENDDO
       !
       ztop_(mintime) = ztop
       ztopmax = MAX (ztop, ztopmax) 
       kkmax   = MAX (kk  , kkmax  ) 
       !print * ,'ztopmax=', mintime,'mn ',ztop_(mintime), ztopmax

       !
       ! if the solution is going to a stationary phase, exit
       IF(mintime > 10) THEN                 
          !   if(mintime > 20) then                     
          !    if( abs(ztop_(mintime)-ztop_(mintime-10)) < DZ ) exit   
          IF( ABS(ztop_(mintime)-ztop_(mintime-10)) < DELZ_THRESOLD) THEN 
	   	  
	   	  !- determine W parameter to determine the VMD
           	  DO k=2,nm1
	   	   W_VMD(k,imm) = w(k)
           	  ENDDO
	   	  EXIT ! finish the integration
	   ENDIF  
       ENDIF

       IF(ilastprint == mintime .AND. izprint.NE.0) THEN
          CALL printout(izprint,nrectotal,mintime,dt,time,ztop,pe,t,te, &
                        qv,qsat,qc,qh,qi,zt,w,vth,sc,qvenv,nrec)
          ilastprint = mintime+1
       ENDIF


    ENDDO   !do next timestep

    !print * ,' ztopmax=',ztopmax,'m',mintime,'mn '
    !print*,'======================================================='
    !
    !the last printout
    IF (izprint.NE.0) THEN
       CALL printout(izprint,nrectotal,mintime,dt,time,ztop, &
                     pe,t,te,qv,qsat,qc,qh,qi,zt,w,vth,sc,qvenv, &
                     nrec)
       CLOSE (2)            
       CLOSE (19)            
    ENDIF

    RETURN  
  END SUBROUTINE MAKEPLUME

  !----------------------------------------------------------------------------
  !
  SUBROUTINE BURN(EFLUX,WATER,tdur,time,heating,mintime,dt,fmoist)  
    !	
    !- calculates the energy flux and water content at lboundary

    REAL,    INTENT(OUT) :: EFLUX
    REAL,    INTENT(OUT) :: WATER
    ! plumegen_coms
    REAL,    INTENT(IN)  :: tdur
    REAL,    INTENT(IN)  :: time
    INTEGER, INTENT(IN)  :: mintime
    REAL,    INTENT(IN)  :: dt
    REAL,    INTENT(IN)  :: fmoist
!!$    REAL,    INTENT(IN)  :: heating(ntime)
    REAL,    INTENT(IN)  :: heating(:)


    !real, parameter :: HEAT = 21.E6 !Joules/kg
    !real, parameter :: HEAT = 15.5E6 !Joules/kg - cerrado
    REAL, PARAMETER :: HEAT = 19.3E6 !Joules/kg - floresta em Alta Floresta (MT)
    !
    ! The emission factor for water is 0.5. The water produced, in kg,
    ! is then  fuel mass*0.5 + (moist/100)*mass per square meter.
    ! The fire burns for DT out of TDUR seconds, the total amount of
    ! fuel burned is AREA*BLOAD*(DT/TDUR) kg. this amount of fuel is
    ! considered to be spread over area AREA and so the mass burned per
    ! unit area is BLOAD*(DT/TDUR), and the rate is BLOAD/TDUR.
    !        
    IF (TIME.GT.TDUR) THEN !is the burn over?   
       EFLUX = 0.000001    !prevent a potential divide by zero
       WATER = 0.  
       RETURN  
    ELSE  
       !                                                   
       EFLUX = HEATING (MINTIME)                          ! Watts/m**2
       !  WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST)       ! kg/m**2 
       WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST) /0.55 ! kg/m**2 
       WATER = WATER * 1000.                              ! g/m**2
       !
       !        print*,'BURN:',time,EFLUX/1.e+9
    ENDIF
    !
    RETURN  
  END SUBROUTINE BURN
  !----------------------------------------------------------------------------
  !
  SUBROUTINE LBOUND(qh,qi,qc,w,t,wc,vth,vti,txs,visc,rho,qv,est,qsat,pe, &
                    alpha,area,rsurf,te,viscosity,dt,qvenv,l,wbar,dqsdz,cvi, &
                    tdur,heating, mintime,fmoist,time,VEL_P,rad_p)  

    !
    ! ********** BOUNDARY CONDITIONS AT ZSURF FOR PLUME AND CLOUD ********
    !
    ! source of equations: J.S. Turner Buoyancy Effects in Fluids
    !                      Cambridge U.P. 1973 p.172,
    !                      G.A. Briggs Plume Rise, USAtomic Energy Commissio
    !                      TID-25075, 1969, P.28
    !
    ! fundamentally a point source below ground. at surface, this produces
    ! a velocity w(1) and temperature T(1) which vary with time. There is
    ! also a water load which will first saturate, then remainder go into
    ! QC(1).
    ! EFLUX = energy flux at ground,watt/m**2 for the last DT
    !

    ! plumegen_coms
    REAL,    INTENT(IN)    :: alpha
    REAL,    INTENT(IN)    :: area
    REAL,    INTENT(IN)    :: rsurf
    REAL,    INTENT(IN)    :: viscosity
    REAL,    INTENT(IN)    :: dt
    INTEGER, INTENT(IN)    :: l
    REAL,    INTENT(IN)    :: wbar
    REAL,    INTENT(IN)    :: dqsdz
    REAL,    INTENT(IN)    :: tdur
    INTEGER, INTENT(IN)    :: mintime
    REAL,    INTENT(IN)    :: fmoist
    REAL,    INTENT(IN)    :: time
!!$    REAL,    INTENT(IN)    :: heating(ntime)
!!$    REAL,    INTENT(INOUT) :: qh   (nkp)
!!$    REAL,    INTENT(INOUT) :: qi   (nkp)
!!$    REAL,    INTENT(INOUT) :: qc   (nkp)
!!$    REAL,    INTENT(INOUT) :: w    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: t    (nkp)
!!$    REAL,    INTENT(INOUT) :: wc   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: vth  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: vti  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: txs  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: visc (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: rho  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qv   (nkp)
!!$    REAL,    INTENT(INOUT) :: est  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qsat (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: pe   (nkp)
!!$    REAL,    INTENT(IN)    :: te   (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv(nkp)
!!$    REAL,    INTENT(IN)    :: cvi  (nkp)
!!$    REAL,    INTENT(INOUT) :: VEL_P(nkp) 
!!$    REAL,    INTENT(INOUT) :: rad_p(nkp) 
    REAL,    INTENT(IN)    :: heating(:)
    REAL,    INTENT(INOUT) :: qh   (:)
    REAL,    INTENT(INOUT) :: qi   (:)
    REAL,    INTENT(INOUT) :: qc   (:)
    REAL,    INTENT(INOUT) :: w    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: t    (:)
    REAL,    INTENT(INOUT) :: wc   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: vth  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: vti  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: txs  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: visc (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: rho  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qv   (:)
    REAL,    INTENT(INOUT) :: est  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qsat (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: pe   (:)
    REAL,    INTENT(IN)    :: te   (:)
    REAL,    INTENT(IN)    :: qvenv(:)
    REAL,    INTENT(IN)    :: cvi  (:)
    REAL,    INTENT(INOUT) :: VEL_P(:) 
    REAL,    INTENT(INOUT) :: rad_p(:) 


    REAL, PARAMETER :: g = 9.80796, r = 287.04, cp = 1004.6, eps = 0.622,tmelt = 273.3
    REAL, PARAMETER :: tfreeze = 269.3, pi = 3.14159, e1 = 1./3., e2 = 5./3.
    REAL            :: es,  eflux, water,  pres, c1,  c2, f, zv,  denscor, xwater

    REAL, EXTERNAL :: esat

    !            
    QH (1) = QH (2)   !soak up hydrometeors
    QI (1) = QI (2)              
    QC (1) = 0.       !no cloud here
    !
    !
    CALL BURN (EFLUX,WATER,tdur,time,heating,mintime,dt,fmoist)  
    !
    !  calculate parameters at boundary from a virtual buoyancy point source
    !
    PRES = PE (1) * 1000.   !need pressure in N/m**2

    C1 = 5. / (6. * ALPHA)  !alpha is entrainment constant

    C2 = 0.9 * ALPHA  

    F = EFLUX / (PRES * CP * PI)  

    F = G * R * F * AREA  !buoyancy flux

    ZV = C1 * RSURF  !virtual boundary height

    W (1) = C1 * ( (C2 * F) **E1) / ZV**E1  !boundary velocity

    DENSCOR = C1 * F / G / (C2 * F) **E1 / ZV**E2   !density correction

    T (1) = TE (1) / (1. - DENSCOR)    !temperature of virtual plume at zsurf

    !
    WC(1) = W(1)

    VEL_P(1) = 0.
    rad_p(1) = rsurf
    !SC(1) = SCE(1)+F/1000.*dt  ! gas/particle (g/g)

    ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    !     match dw/dz,dt/dz at the boundary. F is conserved.
    !
    !WBAR = W (1) * (1. - 1. / (6. * ZV) )  
    !ADVW = WBAR * W (1) / (3. * ZV)  
    !ADVT = WBAR * (5. / (3. * ZV) ) * (DENSCOR / (1. - DENSCOR) )  
    !ADVC = 0.  
    !ADVH = 0.  
    !ADVI = 0.  
    !ADIABAT = - WBAR * G / CP  
    VTH (1) = - 4.  
    VTI (1) = - 3.  
    TXS (1) = T (1) - TE (1)  

    VISC (1) = VISCOSITY  

    RHO (1) = 3483.8 * PE (1) / T (1)   !air density at level 1, g/m**3

    XWATER = WATER / (W (1) * DT * RHO (1) )   !firewater mixing ratio

    QV (1) = XWATER + QVENV (1)  !plus what's already there 


    !  PE esta em kPa  - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
    ES       = 0.1*ESAT (T(1)) !blob saturation vapor pressure, em kPa
    !  rotina do plumegen ja calcula em kPa
    !  ES       = ESAT_PR (T(1))  !blob saturation vapor pressure, em kPa

    EST  (1)  = ES                                  
    QSAT (1) = (EPS * ES) / (PE (1) - ES)   !blob saturation lwc g/g dry air

    IF (QV (1) .GT. QSAT (1) ) THEN  
       QC (1) = QV   (1) - QSAT (1) + QC (1)  !remainder goes into cloud drops
       QV (1) = QSAT (1)  
    ENDIF
    !
    CALL WATERBAL(qc,l,qh,qi,qv,t,qsat,wbar,dqsdz,dt,rho,est,cvi)
    !
    RETURN  
  END SUBROUTINE LBOUND

  !----------------------------------------------------------------------------
  !
  SUBROUTINE INITIAL (kmt,txs,w,t,wc,wt,qv,vth,vti,qh,qi,qc,est,qsat,rho,   &
                      radius,visc,te,qvenv,rsurf,alpha,zt,viscosity,pe, &
                      area,dt,l,wbar,dqsdz,cvi,tdur,heating,mintime,  &
                      fmoist,time,VEL_P,rad_p) 

    !
    ! ************* SETS UP INITIAL CONDITIONS FOR THE PROBLEM ************

    INTEGER , INTENT(IN) :: kmt

    ! plumegen_coms
    REAL,    INTENT(IN)    :: rsurf
    REAL,    INTENT(IN)    :: alpha
    REAL,    INTENT(IN)    :: viscosity
    REAL,    INTENT(IN)    :: area
    REAL,    INTENT(IN)    :: dt
    INTEGER, INTENT(IN)    :: l
    REAL,    INTENT(IN)    :: wbar
    REAL,    INTENT(IN)    :: dqsdz
    REAL,    INTENT(IN)    :: tdur
    INTEGER, INTENT(IN)    :: mintime
    REAL,    INTENT(IN)    :: fmoist
    REAL,    INTENT(IN)    :: time
!!$    REAL,    INTENT(IN)    :: heating(ntime)
!!$    REAL,    INTENT(INOUT) :: txs   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: w     (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: t     (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: wc    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: wt    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qv    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: vth   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: vti   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qh    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qi    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qc    (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: est   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: qsat  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: rho   (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: radius(nkp)
!!$    REAL,    INTENT(INOUT) :: visc  (nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: te    (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv (nkp)
!!$    REAL,    INTENT(IN)    :: zt    (nkp)
!!$    REAL,    INTENT(IN)    :: pe    (nkp)
!!$    REAL,    INTENT(IN)    :: cvi   (nkp)
!!$    REAL,    INTENT(INOUT) :: VEL_P (nkp) 
!!$    REAL,    INTENT(INOUT) :: rad_p (nkp) 
    REAL,    INTENT(IN)    :: heating(:)
    REAL,    INTENT(INOUT) :: txs   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: w     (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: t     (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: wc    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: wt    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qv    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: vth   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: vti   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qh    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qi    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qc    (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: est   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: qsat  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: rho   (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: radius(:)
    REAL,    INTENT(INOUT) :: visc  (:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: te    (:)
    REAL,    INTENT(IN)    :: qvenv (:)
    REAL,    INTENT(IN)    :: zt    (:)
    REAL,    INTENT(IN)    :: pe    (:)
    REAL,    INTENT(IN)    :: cvi   (:)
    REAL,    INTENT(INOUT) :: VEL_P (:) 
    REAL,    INTENT(INOUT) :: rad_p (:) 


    REAL, PARAMETER :: tfreeze = 269.3
    INTEGER         :: k
    REAL            :: es

    INTEGER         :: n ! (DMK) scratch

    REAL, EXTERNAL :: esat

    !

    N=kmt
    ! initialize temperature structure,to the end of equal spaced sounding,
    DO k = 1, N			  
       TXS (k) = 0.0  
       W (k) = 0.0             
       T (k) = TE(k)       !blob set to environment		  
       WC(k) = 0.0
       WT(k) = 0.0
       QV(k) = QVENV (k)   !blob set to environment             
       VTH(k) = 0.		!initial rain velocity = 0	   
       VTI(k) = 0.		!initial ice  velocity = 0	              
       QH(k) = 0.		!no rain			     
       QI(k) = 0.		!no ice 			     
       QC(k) = 0.		!no cloud drops	                     
       !  PE esta em kPa  - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
       ES       = 0.1*ESAT (T(k)) !blob saturation vapor pressure, em kPa
       !  rotina do plumegen calcula em kPa
       !  ES       = ESAT_PR (T(k))  !blob saturation vapor pressure, em kPa
       EST  (k) = ES  
       QSAT (k) = (.622 * ES) / (PE (k) - ES) !saturation lwc g/g
       RHO  (k) = 3483.8 * PE (k) / T (k) 	!dry air density g/m**3    

       VEL_P(k) = 0.
       rad_p(k) = 0.
       
    ENDDO

    ! Initialize the entrainment radius, Turner-style plume
    radius(1) = rsurf
    rad_p(1)  = rsurf
    DO k=2,N
       radius(k) = radius(k-1)+(6./5.)*alpha*(zt(k)-zt(k-1))
       rad_p(k)  = radius(k)
   ENDDO

    !  Initialize the viscosity
    VISC (1) = VISCOSITY
    DO k=2,N
     !  VISC (k) = VISCOSITY!max(1.e-3,visc(k-1) - 1.* VISCOSITY/float(nkp))
       VISC (k) = MAX(1.e-3,visc(k-1) - 1.* VISCOSITY/float(nkp))
    ENDDO
    !--   Initialize gas/concentration
    !DO k =10,20
    !   SC(k) = 20.
    !ENDDO
    !stop 333

    CALL LBOUND(qh,qi,qc,w,t,wc,vth,vti,txs,visc,rho,qv,est,qsat,pe, &
                alpha,area,rsurf,te,viscosity,dt,qvenv,l,wbar,dqsdz,cvi, &
                tdur,heating,mintime,fmoist,time,VEL_P,rad_p)
    RETURN  
  END SUBROUTINE INITIAL

  !----------------------------------------------------------------------------
  !
  SUBROUTINE damp_grav_wave(ifrom,nm1,deltak,dt,zt,zm,w,t,tt,te)

    INTEGER, INTENT(IN)    :: ifrom
    INTEGER, INTENT(IN)    :: nm1
    INTEGER, INTENT(IN)    :: deltak
    REAL,    INTENT(IN)    :: dt
!!$    REAL,    INTENT(IN)    :: zt(nm1)
!!$    REAL,    INTENT(IN)    :: zm(nm1)
!!$    REAL,    INTENT(INOUT) :: w (nm1)
!!$    REAL,    INTENT(INOUT) :: t (nm1)
!!$    REAL,    INTENT(INOUT) :: tt(nm1)
!!$    REAL,    INTENT(IN)    :: te(nm1)
    REAL,    INTENT(IN)    :: zt(:)
    REAL,    INTENT(IN)    :: zm(:)
    REAL,    INTENT(INOUT) :: w (:)
    REAL,    INTENT(INOUT) :: t (:)
    REAL,    INTENT(INOUT) :: tt(:)
    REAL,    INTENT(IN)    :: te(:)

    REAL :: dummy(nm1)

    IF(ifrom==1) THEN
       CALL friction(ifrom,nm1,deltak,dt,zt,zm,t,tt,te)
       !call friction(ifrom,nm1,dt,zt,zm,qv,qvt,qvenv)
       ! call friction(ifrom,nm1,deltak,dt,zt,zm,vel_p,vel_t,vel_e)
       RETURN
    ENDIF

    dummy(:) = 0.
    IF(ifrom==2) CALL friction(ifrom,nm1,deltak,dt,zt,zm,w,dummy ,dummy)
    !call friction(ifrom,nm1,dt,zt,zm,qi,qit ,dummy)
    !call friction(ifrom,nm1,dt,zt,zm,qh,qht ,dummy)
    !call friction(ifrom,nm1,dt,zt,zm,qc,qct ,dummy)
    RETURN
  END SUBROUTINE damp_grav_wave

  !----------------------------------------------------------------------------
  !
  SUBROUTINE friction(ifrom,nm1,deltak,dt,zt,zm,var1,vart,var2)

    INTEGER, INTENT(IN)    :: ifrom
    INTEGER, INTENT(IN)    :: nm1
    INTEGER, INTENT(IN)    :: deltak
    REAL,    INTENT(IN)    :: dt
    REAL,    INTENT(IN)    :: zt  (nm1)
    REAL,    INTENT(IN)    :: zm  (nm1)
    REAL,    INTENT(INOUT) :: var1(nm1)
    REAL,    INTENT(INOUT) :: vart(nm1)
    REAL,    INTENT(IN)    :: var2(nm1)


    INTEGER :: k, kf
    REAL    :: zmkf, ztop, distim, c1, c2

    !nfpt=50
    !kf = nm1 - nfpt
    !kf = nm1 - INT(deltak/2)
    kf = nm1 - INT(deltak)

    zmkf = zm(kf) !old: float(kf )*dz
    ztop = zm(nm1)
    !distim = min(4.*dt,200.)
    !distim = 60. ! orig
    distim = MIN(3.*dt,60.)

    c1 = 1. / (distim * (ztop - zmkf))
    c2 = dt * c1

    IF(ifrom == 1) THEN  
       DO k = nm1,2,-1
          IF (zt(k) .LE. zmkf) CYCLE
          vart(k) = vart(k)   + c1 * (zt(k) - zmkf)*(var2(k) - var1(k))
       ENDDO
    ELSEIF(ifrom == 2) THEN
       DO k = nm1,2,-1
          IF (zt(k) .LE. zmkf) CYCLE
          var1(k) =  var1(k) + c2 * (zt(k) - zmkf)*(var2(k) - var1(k))
       ENDDO
    ENDIF
    RETURN
  END SUBROUTINE friction

  !----------------------------------------------------------------------------
  !
  SUBROUTINE vel_advectc_plumerise(m1,wc,wt,rho,dzm)

    INTEGER, INTENT(IN)    :: m1
!!$    REAL,    INTENT(IN)    :: wc (m1)
!!$    REAL,    INTENT(INOUT) :: wt (m1)
!!$    REAL,    INTENT(IN)    :: rho(m1)
!!$    REAL,    INTENT(IN)    :: dzm(m1)
    REAL,    INTENT(IN)    :: wc (:)
    REAL,    INTENT(INOUT) :: wt (:)
    REAL,    INTENT(IN)    :: rho(:)
    REAL,    INTENT(IN)    :: dzm(:)

    INTEGER :: k
    REAL    :: flxw(m1)
    REAL    :: dn0(m1) ! var local
    REAL    :: c1z

    !dzm(:)= 1./dz

    dn0(1:m1)=rho(1:m1)*1.e-3 ! converte de cgs para mks

    flxw(1) = wc(1) * dn0(1) 

    DO k = 2,m1-1
       flxw(k) = wc(k) * .5 * (dn0(k) + dn0(k+1))
    ENDDO

    ! Compute advection contribution to W tendency

    c1z = .5 

    DO k = 2,m1-2

       wt(k) = wt(k)  &
            + c1z * dzm(k) / (dn0(k) + dn0(k+1)) *     (   &
            (flxw(k) + flxw(k-1))  * (wc(k) + wc(k-1))   &
            - (flxw(k) + flxw(k+1))  * (wc(k) + wc(k+1))   &
            + (flxw(k+1) - flxw(k-1)) * 2.* wc(k)       )

    ENDDO

    RETURN
  END SUBROUTINE vel_advectc_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE hadvance_plumerise(iac,m1,wc,wp,mintime)

    INTEGER, INTENT(IN)    :: iac
    INTEGER, INTENT(IN)    :: m1
!!$    REAL,    INTENT(INOUT) :: wc(m1)
!!$    REAL,    INTENT(INOUT) :: wp(m1)
    REAL,    INTENT(INOUT) :: wc(:)
    REAL,    INTENT(INOUT) :: wp(:)
    INTEGER, INTENT(IN)    :: mintime


    REAL :: dummy(m1)
    REAL :: eps

    !     It is here that the Asselin filter is applied.  For the velocities
    !     and pressure, this must be done in two stages, the first when
    !     IAC=1 and the second when IAC=2.


    eps = .2
    IF(mintime == 1) eps=0.5

    !     For both IAC=1 and IAC=2, call PREDICT for U, V, W, and P.
    !
    CALL predict_plumerise(m1,wc,wp,dummy,iac,eps)
    !print*,'mintime',mintime,eps
    !do k=1,m1
    !   print*,'W-HAD',k,wc(k),wp(k),wt(k)
    !enddo
    RETURN
  END SUBROUTINE hadvance_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE predict_plumerise(npts,ac,ap,af,iac,epsu)

    INTEGER, INTENT(IN)    :: npts
!!$    REAL,    INTENT(INOUT) :: ac(npts)
!!$    REAL,    INTENT(INOUT) :: ap(npts)
!!$    REAL,    INTENT(INOUT) :: af(npts) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: ac(:)
    REAL,    INTENT(INOUT) :: ap(:)
    REAL,    INTENT(INOUT) :: af(:) ! (DMK) alterado (OUT) para (INOUT)
    INTEGER, INTENT(IN)    :: iac
    REAL,    INTENT(IN)    :: epsu


    INTEGER :: m


    !     For IAC=3, this routine moves the arrays AC and AP forward by
    !     1 time level by adding in the prescribed tendency. It also
    !     applies the Asselin filter given by:

    !              {AC} = AC + EPS * (AP - 2 * AC + AF)

    !     where AP,AC,AF are the past, current and future time levels of A.
    !     All IAC=1 does is to perform the {AC} calculation without the AF
    !     term present.  IAC=2 completes the calculation of {AC} by adding
    !     the AF term only, and advances AC by filling it with input AP
    !     values which were already updated in ACOUSTC.
    !

    IF (iac .EQ. 1) THEN
       DO m = 1,npts
          ac(m) = ac(m) + epsu * (ap(m) - 2. * ac(m))
       ENDDO
       RETURN
    ELSEIF (iac .EQ. 2) THEN
       DO m = 1,npts
          af(m) = ap(m)
          ap(m) = ac(m) + epsu * af(m)
       ENDDO
       !elseif (iac .eq. 3) then
       !   do m = 1,npts
       !      af(m) = ap(m) + dtlp * fa(m)
       !   enddo
       !   do m = 1,npts
       !      ap(m) = ac(m) + epsu * (ap(m) - 2. * ac(m) + af(m))
       !   enddo
    ENDIF

    DO m = 1,npts
       ac(m) = af(m)
    ENDDO
    RETURN
  END SUBROUTINE predict_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE  buoyancy_plumerise(m1, T, TE, QV, QVENV, QH, QI, QC, WT, scr1)

    INTEGER, INTENT(IN)    :: m1
!!$    REAL,    INTENT(IN)    :: T    (m1)
!!$    REAL,    INTENT(IN)    :: TE   (m1)
!!$    REAL,    INTENT(IN)    :: QV   (m1)
!!$    REAL,    INTENT(IN)    :: QVENV(m1)
!!$    REAL,    INTENT(IN)    :: QH   (m1)
!!$    REAL,    INTENT(IN)    :: QI   (m1)
!!$    REAL,    INTENT(IN)    :: QC   (m1)
!!$    REAL,    INTENT(INOUT) :: wt   (m1)
!!$    REAL,    INTENT(INOUT) :: scr1 (m1) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: T    (:)
    REAL,    INTENT(IN)    :: TE   (:)
    REAL,    INTENT(IN)    :: QV   (:)
    REAL,    INTENT(IN)    :: QVENV(:)
    REAL,    INTENT(IN)    :: QH   (:)
    REAL,    INTENT(IN)    :: QI   (:)
    REAL,    INTENT(IN)    :: QC   (:)
    REAL,    INTENT(INOUT) :: wt   (:)
    REAL,    INTENT(INOUT) :: scr1 (:) ! (DMK) alterado (OUT) para (INOUT)


    REAL, PARAMETER :: g = 9.8, eps = 0.622, gama = 0.5 ! mass virtual coeff.
    REAL, PARAMETER :: mu = 0.15 

    INTEGER         :: k
    REAL            :: TV, TVE, QWTOTL, umgamai

    !- orig
    umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
                           ! pertubacoes nao-hidrostaticas no campo de pressao

    !- new                 ! Siesbema et al, 2004
    !umgamai = 1./(1.-2.*mu)

    DO k = 2,m1-1

       TV =   T(k) * (1. + (QV(k)   /EPS))/(1. + QV(k)   )  !blob virtual temp.                                        	   
       TVE = TE(k) * (1. + (QVENV(k)/EPS))/(1. + QVENV(k))  !and environment

       QWTOTL = QH(k) + QI(k) + QC(k)                       ! QWTOTL*G is drag
       !- orig
       !scr1(k)= G*( umgamai*(  TV - TVE) / TVE   - QWTOTL) 
       scr1(k)= G*  umgamai*( (TV - TVE) / TVE   - QWTOTL) 

       !if(k .lt. 10)print*,'BT',k,TV,TVE,TVE,QWTOTL
    ENDDO

    DO k = 2,m1-2
       wt(k) = wt(k)+0.5*(scr1(k)+scr1(k+1))
       !   print*,'W-BUO',k,wt(k),scr1(k),scr1(k+1)
    ENDDO

  END SUBROUTINE  buoyancy_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE ENTRAINMENT(m1,w,wt,radius,ALPHA,vel_p,vel_e)

    INTEGER, INTENT(IN)    :: m1
    REAL,    INTENT(IN)    :: alpha
!!$    REAL,    INTENT(IN)    :: w     (m1)
!!$    REAL,    INTENT(INOUT) :: wt    (m1)
!!$    REAL,    INTENT(IN)    :: radius(m1)
!!$    REAL,    INTENT(IN)    :: vel_p (m1)
!!$    REAL,    INTENT(IN)    :: vel_e (m1)
    REAL,    INTENT(IN)    :: w     (:)
    REAL,    INTENT(INOUT) :: wt    (:)
    REAL,    INTENT(IN)    :: radius(:)
    REAL,    INTENT(IN)    :: vel_p (:)
    REAL,    INTENT(IN)    :: vel_e (:)


    REAL, PARAMETER :: mu = 0.15 ,gama = 0.5 ! mass virtual coeff.

    INTEGER         :: k
    REAL            :: DMDTM, WBAR, RADIUS_BAR, umgamai, DYN_ENTR

    !- new - Siesbema et al, 2004
    !umgamai = 1./(1.-2.*mu)

    !- orig
    !umgamai = 1
    umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
                           ! pertubacoes nao-hidrostaticas no campo de pressao

    !
    !-- ALPHA/RADIUS(L) = (1/M)DM/DZ  (W 14a)
    DO k=2,m1-1

       !-- for W: WBAR is only W(k)
       WBAR=W(k)          
       RADIUS_BAR = 0.5*(RADIUS(k) + RADIUS(k-1))
       ! orig
       !DMDTM =           2. * ALPHA * ABS (WBAR) / RADIUS_BAR  != (1/M)DM/DT
       DMDTM  = umgamai * 2. * ALPHA * ABS (WBAR) / RADIUS_BAR  != (1/M)DM/DT
       
       !--  DMDTM*W(L) entrainment,
       wt(k) = wt(k)  - DMDTM*ABS (WBAR)

       !if(VEL_P (k) - VEL_E (k) > 0.) cycle

       !-   dynamic entrainment
       DYN_ENTR =  (2./3.1416)*0.5*ABS (VEL_P(k)-VEL_E(k)+VEL_P(k-1)-&
            VEL_E(k-1)) /RADIUS_BAR

       wt(k) = wt(k)  - DYN_ENTR*ABS (WBAR)
       
       !- entraiment acceleration for output only
       !dwdt_entr(k) =  - DMDTM*ABS (WBAR)- DYN_ENTR*ABS (WBAR)

    ENDDO
  END SUBROUTINE  ENTRAINMENT

  !----------------------------------------------------------------------------
  !
  SUBROUTINE scl_advectc_plumerise(mzp,scr1,dt,w,wc,rho,dzm,zt,zm,dzt, &
                                   t,tt,qv,qvt,qc,qct,qi,qit,qh,qht,vel_p, &
                                   vel_t,rad_p,rad_t)

    INTEGER, INTENT(IN) :: mzp
    ! plumegen_coms
    REAL, INTENT(IN)    :: dt
!!$    REAL, INTENT(INOUT) :: scr1 (nkp)
!!$    REAL, INTENT(IN)    :: w    (nkp)
!!$    REAL, INTENT(IN)    :: wc   (nkp)
!!$    REAL, INTENT(IN)    :: rho  (nkp)
!!$    REAL, INTENT(IN)    :: dzm  (nkp)
!!$    REAL, INTENT(IN)    :: zt   (nkp)
!!$    REAL, INTENT(IN)    :: zm   (nkp)
!!$    REAL, INTENT(IN)    :: dzt  (nkp)
!!$    REAL, INTENT(IN)    :: t    (nkp)
!!$    REAL, INTENT(INOUT) :: tt   (nkp)
!!$    REAL, INTENT(IN)    :: qv   (nkp)
!!$    REAL, INTENT(INOUT) :: qvt  (nkp)
!!$    REAL, INTENT(IN)    :: qc   (nkp)
!!$    REAL, INTENT(INOUT) :: qct  (nkp)
!!$    REAL, INTENT(IN)    :: qi   (nkp)
!!$    REAL, INTENT(INOUT) :: qit  (nkp)
!!$    REAL, INTENT(IN)    :: qh   (nkp)
!!$    REAL, INTENT(INOUT) :: qht  (nkp)
!!$    REAL, INTENT(IN)    :: vel_p(nkp)
!!$    REAL, INTENT(INOUT) :: vel_t(nkp)
!!$    REAL, INTENT(IN)    :: rad_p(nkp)
!!$    REAL, INTENT(INOUT) :: rad_t(nkp)
    REAL, INTENT(INOUT) :: scr1 (:)
    REAL, INTENT(IN)    :: w    (:)
    REAL, INTENT(IN)    :: wc   (:)
    REAL, INTENT(IN)    :: rho  (:)
    REAL, INTENT(IN)    :: dzm  (:)
    REAL, INTENT(IN)    :: zt   (:)
    REAL, INTENT(IN)    :: zm   (:)
    REAL, INTENT(IN)    :: dzt  (:)
    REAL, INTENT(IN)    :: t    (:)
    REAL, INTENT(INOUT) :: tt   (:)
    REAL, INTENT(IN)    :: qv   (:)
    REAL, INTENT(INOUT) :: qvt  (:)
    REAL, INTENT(IN)    :: qc   (:)
    REAL, INTENT(INOUT) :: qct  (:)
    REAL, INTENT(IN)    :: qi   (:)
    REAL, INTENT(INOUT) :: qit  (:)
    REAL, INTENT(IN)    :: qh   (:)
    REAL, INTENT(INOUT) :: qht  (:)
    REAL, INTENT(IN)    :: vel_p(:)
    REAL, INTENT(INOUT) :: vel_t(:)
    REAL, INTENT(IN)    :: rad_p(:)
    REAL, INTENT(INOUT) :: rad_t(:)


    REAL    :: dtlto2
    INTEGER :: k
    REAL    :: vt3dc(nkp) ! (DMK) scratch
    REAL    :: vt3df(nkp) ! (DMK) scratch
    REAL    :: vt3dg(nkp) ! (DMK) scratch
    REAL    :: vt3dk(nkp) ! (DMK) scratch
    REAL    :: vctr1(nkp) ! (DMK) scratch
    REAL    :: vctr2(nkp) ! (DMK) scratch

    !  wp => w
    !- Advect  scalars
    dtlto2   = .5 * dt
    !  vt3dc(1) =      (w(1) + wc(1)) * dtlto2 * dne(1)
    vt3dc(1) =      (w(1) + wc(1)) * dtlto2 * rho(1)*1.e-3!converte de CGS p/ MKS
    vt3df(1) = .5 * (w(1) + wc(1)) * dtlto2 * dzm(1)

    DO k = 2,mzp
       !     vt3dc(k) =  (w(k) + wc(k)) * dtlto2 *.5 * (dne(k) + dne(k+1))
       vt3dc(k) =  (w(k) + wc(k)) * dtlto2 *.5 * (rho(k) + rho(k+1))*1.e-3
       vt3df(k) =  (w(k) + wc(k)) * dtlto2 *.5 *  dzm(k)
       !print*,'vt3df-vt3dc',k,vt3dc(k),vt3df(k)
    ENDDO


    !-srf-24082005
    !  do k = 1,mzp-1
    DO k = 1,mzp
       vctr1(k) = (zt(k+1) - zm(k)) * dzm(k)
       vctr2(k) = (zm(k)   - zt(k)) * dzm(k)
       !    vt3dk(k) = dzt(k) / dne(k)
       vt3dk(k) = dzt(k) /(rho(k)*1.e-3)
       !print*,'VT3dk',k,dzt(k) , dne(k)
    ENDDO

    !      scalarp => scalar_tab(n,ngrid)%var_p
    !      scalart => scalar_tab(n,ngrid)%var_t

    !- temp advection tendency (TT)
    scr1=T
!!$    CALL fa_zc_plumerise(mzp,T,scr1(1),vt3dc(1),vt3df(1),vt3dg(1), &
!!$                         vt3dk(1),vctr1,vctr2)
    CALL fa_zc_plumerise(mzp,T,scr1,vt3dc,vt3df,vt3dg, &
                         vt3dk,vctr1,vctr2)

!!$    CALL advtndc_plumerise(mzp,T,scr1(1),TT,dt)
    CALL advtndc_plumerise(mzp,T,scr1,TT,dt)

    !- water vapor advection tendency (QVT)
    scr1=QV
!!$    CALL fa_zc_plumerise(mzp,QV,scr1(1),vt3dc(1),vt3df(1),vt3dg(1), &
!!$                         vt3dk(1),vctr1,vctr2)
    CALL fa_zc_plumerise(mzp,QV,scr1,vt3dc,vt3df,vt3dg, &
                         vt3dk,vctr1,vctr2)

!!$    CALL advtndc_plumerise(mzp,QV,scr1(1),QVT,dt)
    CALL advtndc_plumerise(mzp,QV,scr1,QVT,dt)

    !- liquid advection tendency (QCT)
    scr1=QC
!!$    CALL fa_zc_plumerise(mzp,QC,scr1(1),vt3dc(1),vt3df(1),vt3dg(1), &
!!$                         vt3dk(1),vctr1,vctr2)
    CALL fa_zc_plumerise(mzp,QC,scr1,vt3dc,vt3df,vt3dg, &
                         vt3dk,vctr1,vctr2)

!!$    CALL advtndc_plumerise(mzp,QC,scr1(1),QCT,dt)
    CALL advtndc_plumerise(mzp,QC,scr1,QCT,dt)

    !- ice advection tendency (QIT)
    scr1=QI
!!$    CALL fa_zc_plumerise(mzp,QI,scr1(1),vt3dc(1),vt3df(1),vt3dg(1), &
!!$                         vt3dk(1),vctr1,vctr2)
    CALL fa_zc_plumerise(mzp,QI,scr1,vt3dc,vt3df,vt3dg, &
                         vt3dk,vctr1,vctr2)

!!$    CALL advtndc_plumerise(mzp,QI,scr1(1),QIT,dt)
    CALL advtndc_plumerise(mzp,QI,scr1,QIT,dt)

    !- hail/rain advection tendency (QHT)
    !   if(ak1 > 0. .or. ak2 > 0.) then

    scr1=QH
!!$    CALL fa_zc_plumerise(mzp,QH,scr1(1),vt3dc(1),vt3df(1),vt3dg(1), &
!!$                         vt3dk(1),vctr1,vctr2)
    CALL fa_zc_plumerise(mzp,QH,scr1,vt3dc,vt3df,vt3dg, &
                         vt3dk,vctr1,vctr2)

!!$    CALL advtndc_plumerise(mzp,QH,scr1(1),QHT,dt)
    CALL advtndc_plumerise(mzp,QH,scr1,QHT,dt)
    !   endif

    !- horizontal wind advection tendency (VEL_T)
    scr1=VEL_P
!!$    CALL fa_zc_plumerise(mzp		       &
!!$    			,VEL_P     ,scr1  (1)  &
!!$    			,vt3dc (1) ,vt3df (1)  &
!!$    			,vt3dg (1) ,vt3dk (1)  &
!!$    			,vctr1,vctr2	     )
    CALL fa_zc_plumerise(mzp		  &
    			,VEL_P, scr1  &
    			,vt3dc, vt3df &
    			,vt3dg, vt3dk &
    			,vctr1, vctr2	  )

!!$    CALL advtndc_plumerise(mzp,VEL_P,scr1(1),VEL_T,dt)
    CALL advtndc_plumerise(mzp,VEL_P,scr1,VEL_T,dt)

    !- vertical radius transport

    scr1=rad_p
!!$    CALL fa_zc_plumerise(mzp                   &
!!$             	        ,rad_p     ,scr1  (1)  &
!!$             	        ,vt3dc (1) ,vt3df (1)  &
!!$             	        ,vt3dg (1) ,vt3dk (1)  &
!!$             	        ,vctr1,vctr2	       )
    CALL fa_zc_plumerise(mzp           &
             	        ,rad_p, scr1   &
             	        ,vt3dc, vt3df  &
             	        ,vt3dg, vt3dk  &
             	        ,vctr1, vctr2	   )

!!$    CALL advtndc_plumerise(mzp,rad_p,scr1(1),rad_t,dt)
    CALL advtndc_plumerise(mzp,rad_p,scr1,rad_t,dt)


    RETURN

  END SUBROUTINE scl_advectc_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE fa_zc_plumerise(m1,scp,scr1,vt3dc,vt3df,vt3dg,vt3dk,vctr1,vctr2)

    INTEGER, INTENT(IN)    :: m1
!!$    REAL,    INTENT(IN)    :: scp  (m1)
!!$    REAL,    INTENT(INOUT) :: scr1 (m1)
!!$    REAL,    INTENT(IN)    :: vt3dc(m1)
!!$    REAL,    INTENT(IN)    :: vt3df(m1)
!!$    REAL,    INTENT(INOUT) :: vt3dg(m1) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: vt3dk(m1)
!!$    REAL,    INTENT(IN)    :: vctr1(m1)
!!$    REAL,    INTENT(IN)    :: vctr2(m1)
    REAL,    INTENT(IN)    :: scp  (:)
    REAL,    INTENT(INOUT) :: scr1 (:)
    REAL,    INTENT(IN)    :: vt3dc(:)
    REAL,    INTENT(IN)    :: vt3df(:)
    REAL,    INTENT(INOUT) :: vt3dg(:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: vt3dk(:)
    REAL,    INTENT(IN)    :: vctr1(:)
    REAL,    INTENT(IN)    :: vctr2(:)


    INTEGER :: k
    REAL    :: dfact

    dfact = .5

    ! Compute scalar flux VT3DG
    DO k = 1,m1-1
       vt3dg(k) = vt3dc(k)                   &
            * (vctr1(k) * scr1(k)        &
            +  vctr2(k) * scr1(k+1)      &
            +  vt3df(k) * (scr1(k) - scr1(k+1)))
    ENDDO

    ! Modify fluxes to retain positive-definiteness on scalar quantities.
    !    If a flux will remove 1/2 quantity during a timestep,
    !    reduce to first order flux. This will remain positive-definite
    !    under the assumption that ABS(CFL(i)) + ABS(CFL(i-1)) < 1.0 if
    !    both fluxes are evacuating the box.

    DO k = 1,m1-1
       IF (vt3dc(k) .GT. 0.) THEN
          IF (vt3dg(k) * vt3dk(k)    .GT. dfact * scr1(k)) THEN
             vt3dg(k) = vt3dc(k) * scr1(k)
          ENDIF
       ELSEIF (vt3dc(k) .LT. 0.) THEN
          IF (-vt3dg(k) * vt3dk(k+1) .GT. dfact * scr1(k+1)) THEN
             vt3dg(k) = vt3dc(k) * scr1(k+1)
          ENDIF
       ENDIF

    ENDDO

    ! Compute flux divergence

    DO k = 2,m1-1
       scr1(k) = scr1(k)  &
            + vt3dk(k) * ( vt3dg(k-1) - vt3dg(k) &
            + scp  (k) * ( vt3dc(k)   - vt3dc(k-1)))
    ENDDO
    RETURN
  END SUBROUTINE fa_zc_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE advtndc_plumerise(m1,scp,sca,sct,dtl)

    INTEGER, INTENT(IN)    :: m1
    REAL,    INTENT(IN)    :: dtl
!!$    REAL,    INTENT(IN)    :: scp(m1)
!!$    REAL,    INTENT(IN)    :: sca(m1)
!!$    REAL,    INTENT(INOUT) :: sct(m1)
    REAL,    INTENT(IN)    :: scp(:)
    REAL,    INTENT(IN)    :: sca(:)
    REAL,    INTENT(INOUT) :: sct(:)


    INTEGER :: k
    REAL    :: dtli

    dtli = 1. / dtl
    DO k = 2,m1-1
       sct(k) = sct(k) + (sca(k)-scp(k)) * dtli
    ENDDO
    RETURN
  END SUBROUTINE advtndc_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE tend0_plumerise(nm1,wt,tt,qvt,qct,qht,qit,vel_t,rad_t)

    ! plumegen_coms
    INTEGER, INTENT(IN)    :: nm1
!!$    REAL,    INTENT(INOUT) :: wt   (nkp)  
!!$    REAL,    INTENT(INOUT) :: tt   (nkp)  
!!$    REAL,    INTENT(INOUT) :: qvt  (nkp) 
!!$    REAL,    INTENT(INOUT) :: qct  (nkp) 
!!$    REAL,    INTENT(INOUT) :: qht  (nkp) 
!!$    REAL,    INTENT(INOUT) :: qit  (nkp) 
!!$    REAL,    INTENT(INOUT) :: vel_t(nkp) 
!!$    REAL,    INTENT(INOUT) :: rad_t(nkp) 
    REAL,    INTENT(INOUT) :: wt   (:)  
    REAL,    INTENT(INOUT) :: tt   (:)  
    REAL,    INTENT(INOUT) :: qvt  (:) 
    REAL,    INTENT(INOUT) :: qct  (:) 
    REAL,    INTENT(INOUT) :: qht  (:) 
    REAL,    INTENT(INOUT) :: qit  (:) 
    REAL,    INTENT(INOUT) :: vel_t(:) 
    REAL,    INTENT(INOUT) :: rad_t(:) 


    wt   (1:nm1)  = 0.
    tt   (1:nm1)  = 0.
    qvt  (1:nm1)  = 0.
    qct  (1:nm1)  = 0.
    qht  (1:nm1)  = 0.
    qit  (1:nm1)  = 0.
    vel_t(1:nm1)  = 0.
    rad_t(1:nm1)  = 0.
    !sct (1:nm1)  = 0.
  END SUBROUTINE tend0_plumerise

  !     ****************************************************************
  SUBROUTINE scl_misc(m1,wbar,w,adiabat,alpha,radius,tt,t,te,qvt,qv, &
                      qvenv,qct,qc,qht,qh,qit,qi,vel_e,vel_p,vel_t,rad_p,&
                      rad_t)

    INTEGER, INTENT(IN)    :: m1
    ! plumegen_coms
    REAL,    INTENT(INOUT) :: wbar    ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: adiabat ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: alpha
!!$    REAL,    INTENT(IN)    :: w     (nkp)
!!$    REAL,    INTENT(IN)    :: radius(nkp)
!!$    REAL,    INTENT(INOUT) :: tt    (nkp)
!!$    REAL,    INTENT(IN)    :: t     (nkp)
!!$    REAL,    INTENT(IN)    :: te    (nkp)
!!$    REAL,    INTENT(INOUT) :: qvt   (nkp)
!!$    REAL,    INTENT(IN)    :: qv    (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv (nkp)
!!$    REAL,    INTENT(INOUT) :: qct   (nkp)
!!$    REAL,    INTENT(IN)    :: qc    (nkp)
!!$    REAL,    INTENT(INOUT) :: qht   (nkp)
!!$    REAL,    INTENT(IN)    :: qh    (nkp)
!!$    REAL,    INTENT(INOUT) :: qit   (nkp)
!!$    REAL,    INTENT(IN)    :: qi    (nkp)
!!$    REAL,    INTENT(IN)    :: vel_e (nkp)
!!$    REAL,    INTENT(IN)    :: vel_p (nkp)
!!$    REAL,    INTENT(INOUT) :: vel_t (nkp)
!!$    REAL,    INTENT(INOUT) :: rad_T (nkp)
!!$    REAL,    INTENT(IN)    :: rad_p (nkp)
    REAL,    INTENT(IN)    :: w     (:)
    REAL,    INTENT(IN)    :: radius(:)
    REAL,    INTENT(INOUT) :: tt    (:)
    REAL,    INTENT(IN)    :: t     (:)
    REAL,    INTENT(IN)    :: te    (:)
    REAL,    INTENT(INOUT) :: qvt   (:)
    REAL,    INTENT(IN)    :: qv    (:)
    REAL,    INTENT(IN)    :: qvenv (:)
    REAL,    INTENT(INOUT) :: qct   (:)
    REAL,    INTENT(IN)    :: qc    (:)
    REAL,    INTENT(INOUT) :: qht   (:)
    REAL,    INTENT(IN)    :: qh    (:)
    REAL,    INTENT(INOUT) :: qit   (:)
    REAL,    INTENT(IN)    :: qi    (:)
    REAL,    INTENT(IN)    :: vel_e (:)
    REAL,    INTENT(IN)    :: vel_p (:)
    REAL,    INTENT(INOUT) :: vel_t (:)
    REAL,    INTENT(INOUT) :: rad_T (:)
    REAL,    INTENT(IN)    :: rad_p (:)


    REAL, PARAMETER :: g = 9.81, cp=1004.
    INTEGER         :: k
    REAL            :: dmdtm

    DO k=2,m1-1
       WBAR    = 0.5*(W(k)+W(k-1))  
       !-- dry adiabat
       ADIABAT = - WBAR * G / CP 
       !      
       !-- entrainment     
       DMDTM = 2. * ALPHA * ABS (WBAR) / RADIUS (k)  != (1/M)DM/DT

       !-- tendency temperature = adv + adiab + entrainment
       TT(k) = TT(K) + ADIABAT - DMDTM * ( T  (k) -    TE (k) ) 

       !-- tendency water vapor = adv  + entrainment
       QVT(K) = QVT(K)         - DMDTM * ( QV (k) - QVENV (k) )

       QCT(K) = QCT(K)	      - DMDTM * ( QC (k)  )
       QHT(K) = QHT(K)	      - DMDTM * ( QH (k)  )
       QIT(K) = QIT(K)	      - DMDTM * ( QI (k)  )


      !-- tendency horizontal speed = adv  + entrainment
      VEL_T(K) = VEL_T(K)     - DMDTM * ( VEL_P (k) - VEL_E (k) )

      !-- tendency horizontal speed = adv  + entrainment
      rad_t(K) = rad_t(K)     + 0.5*DMDTM*(6./5.)*RADIUS (k)

       !-- tendency gas/particle = adv  + entrainment
       !      SCT(K) = SCT(K)         - DMDTM * ( SC (k) -   SCE (k) )

    ENDDO
  END SUBROUTINE scl_misc

!     ****************************************************************
  SUBROUTINE scl_dyn_entrain(m1,wbar,w,adiabat,alpha,radius,tt,t, &
                             te,qvt,qv,qvenv,qct,qc,qht,qh,qit,qi,    &
                             vel_e,vel_p,vel_t,rad_p,rad_t)

    INTEGER, INTENT(IN)    :: m1
    ! plumegen_coms
    REAL,    INTENT(INOUT) :: wbar 
    REAL,    INTENT(INOUT) :: adiabat 
    REAL,    INTENT(IN)    :: alpha
!!$    REAL,    INTENT(IN)    :: w     (nkp)
!!$    REAL,    INTENT(IN)    :: radius(nkp)
!!$    REAL,    INTENT(INOUT) :: tt    (nkp)
!!$    REAL,    INTENT(IN)    :: t     (nkp)
!!$    REAL,    INTENT(IN)    :: te    (nkp)
!!$    REAL,    INTENT(INOUT) :: qvt   (nkp)
!!$    REAL,    INTENT(IN)    :: qv    (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv (nkp)
!!$    REAL,    INTENT(INOUT) :: qct   (nkp)
!!$    REAL,    INTENT(IN)    :: qc    (nkp)
!!$    REAL,    INTENT(INOUT) :: qht   (nkp)
!!$    REAL,    INTENT(IN)    :: qh    (nkp)
!!$    REAL,    INTENT(INOUT) :: qit   (nkp)
!!$    REAL,    INTENT(IN)    :: qi    (nkp)
!!$    REAL,    INTENT(IN)    :: vel_e (nkp)
!!$    REAL,    INTENT(IN)    :: vel_p (nkp)
!!$    REAL,    INTENT(INOUT) :: vel_t (nkp)
!!$    REAL,    INTENT(INOUT) :: rad_T (nkp)
!!$    REAL,    INTENT(IN)    :: rad_p (nkp)
    REAL,    INTENT(IN)    :: w     (:)
    REAL,    INTENT(IN)    :: radius(:)
    REAL,    INTENT(INOUT) :: tt    (:)
    REAL,    INTENT(IN)    :: t     (:)
    REAL,    INTENT(IN)    :: te    (:)
    REAL,    INTENT(INOUT) :: qvt   (:)
    REAL,    INTENT(IN)    :: qv    (:)
    REAL,    INTENT(IN)    :: qvenv (:)
    REAL,    INTENT(INOUT) :: qct   (:)
    REAL,    INTENT(IN)    :: qc    (:)
    REAL,    INTENT(INOUT) :: qht   (:)
    REAL,    INTENT(IN)    :: qh    (:)
    REAL,    INTENT(INOUT) :: qit   (:)
    REAL,    INTENT(IN)    :: qi    (:)
    REAL,    INTENT(IN)    :: vel_e (:)
    REAL,    INTENT(IN)    :: vel_p (:)
    REAL,    INTENT(INOUT) :: vel_t (:)
    REAL,    INTENT(INOUT) :: rad_T (:)
    REAL,    INTENT(IN)    :: rad_p (:)


    REAL, PARAMETER :: g = 9.81, cp=1004., pi=3.1416
    INTEGER         :: k
    REAL            :: dmdtm

    DO k=2,m1-1
      !      
      !-- tendency horizontal radius from dyn entrainment
     	   !rad_t(K) = rad_t(K)   +	(vel_e(k)-vel_p(k)) /pi
     	    rad_t(K) = rad_t(K)   + ABS((vel_e(k)-vel_p(k)))/pi
      
      !-- entrainment	  
     	   !DMDTM = (2./3.1416)  *     (VEL_E (k) - VEL_P (k)) / RADIUS (k)  
     	    DMDTM = (2./3.1416)  *  ABS(VEL_E (k) - VEL_P (k)) / RADIUS (k)  
      
      !-- tendency horizontal speed  from dyn entrainment
     	    VEL_T(K) = VEL_T(K)     - DMDTM * ( VEL_P (k) - VEL_E (k) )
      
      !     if(VEL_P (k) - VEL_E (k) > 0.) cycle
      
      !-- tendency temperature  from dyn entrainment
     	    TT(k) = TT(K)	    - DMDTM * ( T (k) - TE  (k) ) 
      
      !-- tendency water vapor  from dyn entrainment
   	    QVT(K) = QVT(K)	    - DMDTM * ( QV (k) - QVENV (k) )
      
     	    QCT(K) = QCT(K)	    - DMDTM * ( QC (k)  )
     	    QHT(K) = QHT(K)	    - DMDTM * ( QH (k)  )
     	    QIT(K) = QIT(K)	    - DMDTM * ( QI (k)  )
      
      !-- tendency gas/particle  from dyn entrainment
      !	 SCT(K) = SCT(K)	 - DMDTM * ( SC (k) - SCE (k) )
    
    ENDDO
   END SUBROUTINE scl_dyn_entrain
  !     ****************************************************************

  SUBROUTINE  visc_W(m1,kmt,zt,visc,zm,w,t,qv,qh,qc,qi,wt,tt,qvt, &
                     qct,qht,qit,vel_p,vel_t,rad_p,rad_t)

    INTEGER, INTENT(IN)    :: m1
    INTEGER, INTENT(IN)    :: kmt
    ! plumegen_coms
!!$    REAL,    INTENT(IN)    :: zt   (nkp)
!!$    REAL,    INTENT(IN)    :: visc (nkp)
!!$    REAL,    INTENT(IN)    :: zm   (nkp)
!!$    REAL,    INTENT(IN)    :: w    (nkp)
!!$    REAL,    INTENT(IN)    :: t    (nkp)
!!$    REAL,    INTENT(IN)    :: qv   (nkp)
!!$    REAL,    INTENT(IN)    :: qh   (nkp)
!!$    REAL,    INTENT(IN)    :: qc   (nkp)
!!$    REAL,    INTENT(IN)    :: qi   (nkp)
!!$    REAL,    INTENT(INOUT) :: wt   (nkp)
!!$    REAL,    INTENT(INOUT) :: tt   (nkp)
!!$    REAL,    INTENT(INOUT) :: qvt  (nkp)
!!$    REAL,    INTENT(INOUT) :: qct  (nkp)
!!$    REAL,    INTENT(INOUT) :: qht  (nkp)
!!$    REAL,    INTENT(INOUT) :: qit  (nkp)  
!!$    REAL,    INTENT(IN)    :: vel_p(nkp)
!!$    REAL,    INTENT(INOUT) :: vel_t(nkp)
!!$    REAL,    INTENT(INOUT) :: rad_T(nkp)
!!$    REAL,    INTENT(IN)    :: rad_p(nkp)
    REAL,    INTENT(IN)    :: zt   (:)
    REAL,    INTENT(IN)    :: visc (:)
    REAL,    INTENT(IN)    :: zm   (:)
    REAL,    INTENT(IN)    :: w    (:)
    REAL,    INTENT(IN)    :: t    (:)
    REAL,    INTENT(IN)    :: qv   (:)
    REAL,    INTENT(IN)    :: qh   (:)
    REAL,    INTENT(IN)    :: qc   (:)
    REAL,    INTENT(IN)    :: qi   (:)
    REAL,    INTENT(INOUT) :: wt   (:)
    REAL,    INTENT(INOUT) :: tt   (:)
    REAL,    INTENT(INOUT) :: qvt  (:)
    REAL,    INTENT(INOUT) :: qct  (:)
    REAL,    INTENT(INOUT) :: qht  (:)
    REAL,    INTENT(INOUT) :: qit  (:)  
    REAL,    INTENT(IN)    :: vel_p(:)
    REAL,    INTENT(INOUT) :: vel_t(:)
    REAL,    INTENT(INOUT) :: rad_T(:)
    REAL,    INTENT(IN)    :: rad_p(:)


    INTEGER :: k, m2
    REAL    :: dz1t, dz1m, dz2t, dz2m, d2wdz, d2tdz, &
               d2qvdz, d2qhdz, d2qcdz, d2qidz,  &
               d2vel_pdz, d2rad_dz

    !srf--- 17/08/2005
    !m2=min(m1+deltak,kmt)
    !m2=MIN(m1,kmt)

    !do k=2,m1-1
    DO k=2,MIN(m1,kmt-1)!m2-1
       DZ1T   = 0.5*(ZT(K+1)-ZT(K-1))
       DZ2T   = VISC (k) / (DZ1T * DZ1T)  
       DZ1M   = 0.5*(ZM(K+1)-ZM(K-1))
       DZ2M   = VISC (k) / (DZ1M * DZ1M)  
       D2WDZ  = (W  (k + 1) - 2 * W  (k) + W  (k - 1) ) * DZ2M  
       D2TDZ  = (T  (k + 1) - 2 * T  (k) + T  (k - 1) ) * DZ2T  
       D2QVDZ = (QV (k + 1) - 2 * QV (k) + QV (k - 1) ) * DZ2T  
       D2QHDZ = (QH (k + 1) - 2 * QH (k) + QH (k - 1) ) * DZ2T 
       D2QCDZ = (QC (k + 1) - 2 * QC (k) + QC (k - 1) ) * DZ2T  
       D2QIDZ = (QI (k + 1) - 2 * QI (k) + QI (k - 1) ) * DZ2T  
       !D2SCDZ = (SC (k + 1) - 2 * SC (k) + SC (k - 1) ) * DZ2T 
       
       d2vel_pdz=(vel_P  (k + 1) - 2 * vel_P  (k) + vel_P  (k - 1) ) * DZ2T
       d2rad_dz =(rad_p  (k + 1) - 2 * rad_p  (k) + rad_p  (k - 1) ) * DZ2T

       WT(k) =   WT(k) + D2WDZ 
       TT(k) =   TT(k) + D2TDZ                          
       QVT(k) =  QVT(k) + D2QVDZ 
       QCT(k) =  QCT(k) + D2QCDZ 
       QHT(k) =  QHT(k) + D2QHDZ 
       QIT(k) =  QIT(k) + D2QIDZ     

       vel_t(k) =   vel_t(k) + d2vel_pdz
       rad_t(k) =   rad_t(k) + d2rad_dz
       
       !SCT(k) =  SCT(k) + D2SCDZ
       !print*,'W-VISC=',k,D2WDZ
    ENDDO

  END SUBROUTINE visc_W

  !     ****************************************************************

  SUBROUTINE update_plumerise(m1,varn,wt,dt,tt,qvt,qct,qht,w,t,qv,qc, &
                              qh,qit,qi, vel_p,vel_t,rad_p,rad_t)

    INTEGER ,         INTENT(IN)    :: m1
    CHARACTER(len=*), INTENT(IN)    :: varn
    ! plumegen_coms
    REAL,             INTENT(IN)    :: dt
!!$    REAL,             INTENT(IN)    :: wt   (nkp)
!!$    REAL,             INTENT(IN)    :: tt   (nkp)
!!$    REAL,             INTENT(IN)    :: qvt  (nkp)
!!$    REAL,             INTENT(IN)    :: qct  (nkp)
!!$    REAL,             INTENT(IN)    :: qht  (nkp)
!!$    REAL,             INTENT(INOUT) :: w    (nkp)
!!$    REAL,             INTENT(INOUT) :: t    (nkp)
!!$    REAL,             INTENT(INOUT) :: qv   (nkp)
!!$    REAL,             INTENT(INOUT) :: qc   (nkp)
!!$    REAL,             INTENT(INOUT) :: qh   (nkp)
!!$    REAL,             INTENT(IN)    :: qit  (nkp) ! (DMK) alterado (INOUT) (IN)
!!$    REAL,             INTENT(INOUT) :: qi   (nkp)
!!$    REAL,             INTENT(INOUT) :: vel_p(nkp)
!!$    REAL,             INTENT(IN)    :: vel_t(nkp)
!!$    REAL,             INTENT(IN)    :: rad_T(nkp)
!!$    REAL,             INTENT(INOUT) :: rad_p(nkp)
    REAL,             INTENT(IN)    :: wt   (:)
    REAL,             INTENT(IN)    :: tt   (:)
    REAL,             INTENT(IN)    :: qvt  (:)
    REAL,             INTENT(IN)    :: qct  (:)
    REAL,             INTENT(IN)    :: qht  (:)
    REAL,             INTENT(INOUT) :: w    (:)
    REAL,             INTENT(INOUT) :: t    (:)
    REAL,             INTENT(INOUT) :: qv   (:)
    REAL,             INTENT(INOUT) :: qc   (:)
    REAL,             INTENT(INOUT) :: qh   (:)
    REAL,             INTENT(IN)    :: qit  (:) ! (DMK) alterado (INOUT) (IN)
    REAL,             INTENT(INOUT) :: qi   (:)
    REAL,             INTENT(INOUT) :: vel_p(:)
    REAL,             INTENT(IN)    :: vel_t(:)
    REAL,             INTENT(IN)    :: rad_T(:)
    REAL,             INTENT(INOUT) :: rad_p(:)


    INTEGER :: k

    IF(varn == 'W') THEN

       DO k=2,m1-1
          W(k) =  W(k) +  WT(k) * DT  
       ENDDO
       RETURN

    ELSE 
       DO k=2,m1-1
          T(k) =  T(k) +  TT(k) * DT  

          QV(k) = QV(k) + QVT(k) * DT  

          QC(k) = QC(k) + QCT(k) * DT !cloud drops travel with air 
          QH(k) = QH(k) + QHT(k) * DT  
          QI(k) = QI(k) + QIT(k) * DT 
          ! SC(k) = SC(k) + SCT(k) * DT 

          !srf---18jun2005  
          QV(k) = MAX(0., QV(k))
          QC(k) = MAX(0., QC(k))
          QH(k) = MAX(0., QH(k))
          QI(k) = MAX(0., QI(k))
	  
	  
	  VEL_P(k) =  VEL_P(k) + VEL_T(k) * DT  

          rad_p(k) =  rad_p(k) + rad_t(k) * DT  

          ! SC(k) = max(0., SC(k))

       ENDDO
    ENDIF
  END SUBROUTINE update_plumerise

  !----------------------------------------------------------------------------
  !
  SUBROUTINE fallpart(m1,rho,vth,vhrel,w,cvh,vti,cvi,qh,qi,zm,qht,qit)

    INTEGER, INTENT(IN)    :: m1
    ! plumegen_coms
    REAL,    INTENT(INOUT) :: vhrel   ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: rho(nkp)
!!$    REAL,    INTENT(INOUT) :: vth(nkp)
!!$    REAL,    INTENT(IN)    :: w  (nkp)
!!$    REAL,    INTENT(INOUT) :: cvh(nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(INOUT) :: vti(nkp)
!!$    REAL,    INTENT(INOUT) :: cvi(nkp) ! (DMK) alterado (OUT) para (INOUT)
!!$    REAL,    INTENT(IN)    :: qh (nkp)
!!$    REAL,    INTENT(IN)    :: qi (nkp)
!!$    REAL,    INTENT(IN)    :: zm (nkp)
!!$    REAL,    INTENT(INOUT) :: qht(nkp)
!!$    REAL,    INTENT(INOUT) :: qit(nkp)
    REAL,    INTENT(IN)    :: rho(:)
    REAL,    INTENT(INOUT) :: vth(:)
    REAL,    INTENT(IN)    :: w  (:)
    REAL,    INTENT(INOUT) :: cvh(:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(INOUT) :: vti(:)
    REAL,    INTENT(INOUT) :: cvi(:) ! (DMK) alterado (OUT) para (INOUT)
    REAL,    INTENT(IN)    :: qh (:)
    REAL,    INTENT(IN)    :: qi (:)
    REAL,    INTENT(IN)    :: zm (:)
    REAL,    INTENT(INOUT) :: qht(:)
    REAL,    INTENT(INOUT) :: qit(:)


    !srf==================================
    !   verificar se o gradiente esta correto 
    !  
    !srf==================================
    !
    !     XNO=1.E7  [m**-4] median volume diameter raindrop,Kessler
    !     VC = 38.3/(XNO**.125), median volume fallspeed eqn., Kessler
    !     for ice, see (OT18), use F0=0.75 per argument there. rho*q
    !     values are in g/m**3, velocities in m/s

    REAL, PARAMETER :: VCONST = 5.107387, EPS = 0.622, F0 = 0.75  
    REAL, PARAMETER :: G = 9.81, CP = 1004.

    INTEGER         :: k
    REAL            :: vtc, dfhz, dfiz, dz1
    REAL            :: virel ! (DMK) scratch

    !
    DO k=2,m1-1

       VTC = VCONST * RHO (k) **.125   ! median volume fallspeed (KTable4)

       !  hydrometeor assembly velocity calculations (K Table4)
       !  VTH(k)=-VTC*QH(k)**.125  !median volume fallspeed, water            
       VTH (k) = - 4.	    !small variation with qh

       VHREL = W (k) + VTH (k)  !relative to surrounding cloud

       !  rain ventilation coefficient for evaporation
       CVH(k) = 1.6 + 0.57E-3 * (ABS (VHREL) ) **1.5  
       !
       !  VTI(k)=-VTC*F0*QI(k)**.125    !median volume fallspeed,ice          
       VTI (k) = - 3.                !small variation with qi

       VIREL = W (k) + VTI (k)       !relative to surrounding cloud
       !
       !  ice ventilation coefficient for sublimation
       CVI(k) = 1.6 + 0.57E-3 * (ABS (VIREL) ) **1.5 / F0  
       !
       !
       IF (VHREL.GE.0.0) THEN  
          DFHZ=QH(k)*(RHO(k  )*VTH(k  )-RHO(k-1)*VTH(k-1))/RHO(k-1)
       ELSE  
          DFHZ=QH(k)*(RHO(k+1)*VTH(k+1)-RHO(k  )*VTH(k  ))/RHO(k)
       ENDIF
       !
       !
       IF (VIREL.GE.0.0) THEN  
          DFIZ=QI(k)*(RHO(k  )*VTI(k  )-RHO(k-1)*VTI(k-1))/RHO(k-1)
       ELSE   
          DFIZ=QI(k)*(RHO(k+1)*VTI(k+1)-RHO(k  )*VTI(k  ))/RHO(k)
       ENDIF

       DZ1=ZM(K)-ZM(K-1)

       qht(k) = qht(k) - DFHZ / DZ1 !hydrometeors don't

       qit(k) = qit(k) - DFIZ / DZ1  !nor does ice? hail, what about

    ENDDO

  END SUBROUTINE fallpart



  !----------------------------------------------------------------------------
  !
  SUBROUTINE printout(izprint,nrectotal,mintime,dt,time,ztop,pe,t,te,qv, &
                      qsat,qc,qh,qi,zt,w,vth,sc,qvenv,nrec)  

    INTEGER, INTENT(IN)    :: izprint
    INTEGER, INTENT(IN)    :: nrectotal
    ! plumegen_coms
    INTEGER, INTENT(IN)    :: mintime
    REAL,    INTENT(IN)    :: dt
    REAL,    INTENT(IN)    :: time
    REAL,    INTENT(IN)    :: ztop
!!$    REAL,    INTENT(IN)    :: pe   (nkp)
!!$    REAL,    INTENT(IN)    :: t    (nkp)
!!$    REAL,    INTENT(IN)    :: te   (nkp)
!!$    REAL,    INTENT(IN)    :: qv   (nkp)
!!$    REAL,    INTENT(IN)    :: qsat (nkp)
!!$    REAL,    INTENT(IN)    :: qc   (nkp)
!!$    REAL,    INTENT(IN)    :: qh   (nkp)
!!$    REAL,    INTENT(IN)    :: qi   (nkp)
!!$    REAL,    INTENT(IN)    :: zt   (nkp)
!!$    REAL,    INTENT(IN)    :: w    (nkp)
!!$    REAL,    INTENT(IN)    :: vth  (nkp)
!!$    REAL,    INTENT(IN)    :: sc   (nkp)
!!$    REAL,    INTENT(IN)    :: qvenv(nkp)
    REAL,    INTENT(IN)    :: pe   (:)
    REAL,    INTENT(IN)    :: t    (:)
    REAL,    INTENT(IN)    :: te   (:)
    REAL,    INTENT(IN)    :: qv   (:)
    REAL,    INTENT(IN)    :: qsat (:)
    REAL,    INTENT(IN)    :: qc   (:)
    REAL,    INTENT(IN)    :: qh   (:)
    REAL,    INTENT(IN)    :: qi   (:)
    REAL,    INTENT(IN)    :: zt   (:)
    REAL,    INTENT(IN)    :: w    (:)
    REAL,    INTENT(IN)    :: vth  (:)
    REAL,    INTENT(IN)    :: sc   (:)
    REAL,    INTENT(IN)    :: qvenv(:)
    ! save attribute
    INTEGER, INTENT(INOUT) :: nrec


    REAL, PARAMETER :: tmelt = 273.3
    INTEGER         :: ko, interval
    REAL            :: pea, btmp, etmp, vap1, vap2, gpkc, gpkh, gpki, deficit

    interval = 1              !debug time interval,min

    !
    IF (IZPRINT.EQ.0) RETURN  

    IF(MINTIME == 1) nrec = 0
    !
    WRITE (2, 430) MINTIME, DT, TIME  
    WRITE (2, 431) ZTOP  
    WRITE (2, 380)  
    !
    ! do the print
    !
    DO 390 KO = 1, nrectotal, interval  

       PEA = PE (KO) * 10.       !pressure is stored in decibars(kPa),print in mb;
       BTMP = T (KO) - TMELT     !temps in Celsius
       ETMP = T (KO) - TE (KO)   !temperature excess
       VAP1 = QV (KO)   * 1000.  !printout in g/kg for all water,
       VAP2 = QSAT (KO) * 1000.  !vapor (internal storage is in g/g)
       GPKC = QC (KO)   * 1000.  !cloud water
       GPKH = QH (KO)   * 1000.  !raindrops
       GPKI = QI (KO)   * 1000.  !ice particles 
       DEFICIT = VAP2 - VAP1     !vapor deficit
       !
       WRITE (2, 400) zt(KO)/1000., PEA, W (KO), BTMP, ETMP, VAP1, &
            VAP2, GPKC, GPKH, GPKI, VTH (KO), SC(KO)
       !
       !
       !                                    !end of printout

390    CONTINUE		  

       nrec=nrec+1
       WRITE (19,rec=nrec) (W (KO), KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (T (KO), KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (TE(KO), KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (QV(KO)*1000., KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (QC(KO)*1000., KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (QH(KO)*1000., KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (QI(KO)*1000., KO=1,nrectotal)
       nrec=nrec+1
       !   write (19,rec=nrec) (SC(KO), KO=1,nrectotal)
       WRITE (19,rec=nrec) (QSAT(KO)*1000., KO=1,nrectotal)
       nrec=nrec+1
       WRITE (19,rec=nrec) (QVENV(KO)*1000., KO=1,nrectotal)

       PRINT*,'ntimes=',nrec/(9)
       !
       RETURN  
       !
       ! ************** FORMATS *********************************************
       !
380    FORMAT(/,' Z(KM) P(MB) W(MPS) T(C)  T-TE   VAP   SAT   QC    QH' &
            '     QI    VTH(MPS) SCAL'/)
       !
400    FORMAT(1H , F4.1,F8.2,F8.2,F7.1,6F6.2,F7.2,1X,F6.2)  
       !
430    FORMAT(1H ,//I5,' MINUTES       DT= ',F6.2,' SECONDS   TIME= ' &
            ,F8.2,' SECONDS')
431    FORMAT(' ZTOP= ',F10.2)  
       !
     END SUBROUTINE printout
     !
     ! *********************************************************************

     SUBROUTINE WATERBAL(qc,l,qh,qi,qv,t,qsat,wbar,dqsdz,dt,rho,est,cvi)

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL,    INTENT(IN)    :: wbar
       REAL,    INTENT(IN)    :: dqsdz
       REAL,    INTENT(IN)    :: dt
!!$       REAL,    INTENT(INOUT) :: qc  (nkp)
!!$       REAL,    INTENT(INOUT) :: qh  (nkp)
!!$       REAL,    INTENT(INOUT) :: qi  (nkp)
!!$       REAL,    INTENT(INOUT) :: qv  (nkp)
!!$       REAL,    INTENT(INOUT) :: t   (nkp)
!!$       REAL,    INTENT(IN)    :: qsat(nkp)
!!$       REAL,    INTENT(IN)    :: rho (nkp)
!!$       REAL,    INTENT(IN)    :: est (nkp)
!!$       REAL,    INTENT(IN)    :: cvi (nkp)
       REAL,    INTENT(INOUT) :: qc  (:)
       REAL,    INTENT(INOUT) :: qh  (:)
       REAL,    INTENT(INOUT) :: qi  (:)
       REAL,    INTENT(INOUT) :: qv  (:)
       REAL,    INTENT(INOUT) :: t   (:)
       REAL,    INTENT(IN)    :: qsat(:)
       REAL,    INTENT(IN)    :: rho (:)
       REAL,    INTENT(IN)    :: est (:)
       REAL,    INTENT(IN)    :: cvi (:)
       !

       IF (QC (L) .LE.1.0E-10) QC (L) = 0.  !DEFEAT UNDERFLOW PROBLEM
       IF (QH (L) .LE.1.0E-10) QH (L) = 0.  
       IF (QI (L) .LE.1.0E-10) QI (L) = 0.  
       !
       CALL EVAPORATE(qc,qh,qi,qv,t,qsat,l,wbar,dqsdz,dt,rho,est,cvi)  !vapor to cloud,cloud to vapor  
       !                             
       CALL SUBLIMATE(t,l,qv,qsat,rho,qi,est,dt)  !vapor to ice  
       !                            
       CALL GLACIATE(qh,l,qv,qsat,t,dt,qi) !rain to ice 

       CALL MELT(qi,l,t,dt,rho,cvi,qh)  !ice to rain
       !         
       !if(ak1 > 0. .or. ak2 > 0.) &
       CALL CONVERT(t,l,qc,rho,qh,dt) !(auto)conversion and accretion 
       !CALL CONVERT2 () !(auto)conversion and accretion 
       !

       RETURN  
     END SUBROUTINE WATERBAL

     ! *********************************************************************
     SUBROUTINE EVAPORATE(qc,qh,qi,qv,t,qsat,l,wbar,dqsdz,dt,rho,est,cvi)  
       !
       !- evaporates cloud,rain and ice to saturation
       !

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL,    INTENT(IN)    :: wbar
       REAL,    INTENT(IN)    :: dqsdz
       REAL,    INTENT(IN)    :: dt
!!$       REAL,    INTENT(INOUT) :: qc  (nkp)
!!$       REAL,    INTENT(INOUT) :: qh  (nkp)
!!$       REAL,    INTENT(INOUT) :: qi  (nkp)
!!$       REAL,    INTENT(INOUT) :: qv  (nkp)
!!$       REAL,    INTENT(INOUT) :: t   (nkp)
!!$       REAL,    INTENT(IN)    :: qsat(nkp)
!!$       REAL,    INTENT(IN)    :: rho (nkp)
!!$       REAL,    INTENT(IN)    :: est (nkp)
!!$       REAL,    INTENT(IN)    :: cvi (nkp)
       REAL,    INTENT(INOUT) :: qc  (:)
       REAL,    INTENT(INOUT) :: qh  (:)
       REAL,    INTENT(INOUT) :: qi  (:)
       REAL,    INTENT(INOUT) :: qv  (:)
       REAL,    INTENT(INOUT) :: t   (:)
       REAL,    INTENT(IN)    :: qsat(:)
       REAL,    INTENT(IN)    :: rho (:)
       REAL,    INTENT(IN)    :: est (:)
       REAL,    INTENT(IN)    :: cvi (:)


       !
       !     XNO=10.0E06
       !     HERC = 1.93*1.E-6*XN035        !evaporation constant
       !
       REAL, PARAMETER :: HERC = 5.44E-4, CP = 1.004, HEATCOND = 2.5E3  
       REAL, PARAMETER :: HEATSUBL = 2834., TMELT = 273., TFREEZE = 269.3
       REAL, PARAMETER :: FRC = HEATCOND / CP, SRC = HEATSUBL / CP

       REAL            :: evhdt, evidt, evrate, evap, sd, &
                          quant, dividend, divisor, devidt

       !
       !
       SD = QSAT (L) - QV (L)  !vapor deficit
       IF (SD.EQ.0.0)  RETURN  
       !IF (abs(SD).lt.1.e-7)  RETURN  

       EVHDT = 0.  
       EVIDT = 0.  
       !evrate =0.; evap=0.; sd=0.0; quant=0.0; dividend=0.0; divisor=0.0; devidt=0.0

       EVRATE = ABS (WBAR * DQSDZ)   !evaporation rate (Kessler 8.32)
       EVAP = EVRATE * DT            !what we can get in DT

       IF (SD.LE.0.0) THEN  !     condense. SD is negative

          IF (EVAP.GE.ABS (SD) ) THEN    !we get it all

             QC (L) = QC  (L) - SD  !deficit,remember?
             QV (L) = QSAT(L)       !set the vapor to saturation  
             T  (L) = T   (L) - SD * FRC  !heat gained through condensation
             !per gram of dry air
             RETURN  

          ELSE  

             QC (L) = QC (L) + EVAP         !get what we can in DT 
             QV (L) = QV (L) - EVAP         !remove it from the vapor
             T  (L) = T  (L) + EVAP * FRC   !get some heat

             RETURN  

          ENDIF
          !
       ELSE                                !SD is positive, need some water
          !
          ! not saturated. saturate if possible. use everything in order
          ! cloud, rain, ice. SD is positive

          IF (EVAP.LE.QC (L) ) THEN        !enough cloud to last DT  
             !

             IF (SD.LE.EVAP) THEN          !enough time to saturate

                QC (L) = QC (L) - SD       !remove cloud                                          
                QV (L) = QSAT (L)          !saturate
                T (L) = T (L) - SD * FRC   !cool the parcel                                          
                RETURN  !done
                !

             ELSE   !not enough time

                SD = SD-EVAP               !use what there is
                QV (L) = QV (L) + EVAP     !add vapor
                T (L) = T (L) - EVAP * FRC !lose heat
                QC (L) = QC (L) - EVAP     !lose cloud
                !go on to rain.                                      
             ENDIF
             !
          ELSE                !not enough cloud to last DT
             !      
             IF (SD.LE.QC (L) ) THEN   !but there is enough to sat

                QV (L) = QSAT (L)  !use it
                QC (L) = QC (L) - SD  
                T  (L) = T (L) - SD * FRC  
                RETURN  

             ELSE            !not enough to sat
                SD = SD-QC (L)  
                QV (L) = QV (L) + QC (L)  
                T  (L) = T (L) - QC (L) * FRC         
                QC (L) = 0.0  !all gone

             ENDIF       !on to rain                           
          ENDIF          !finished with cloud
          !
          !  but still not saturated, so try to use some rain
          !  this is tricky, because we only have time DT to evaporate. 
          !  if there is enough rain, we can evaporate it for dt. ice can 
          !  also sublimate at the same time. there is a compromise 
          !  here.....use rain first, then
          !  ice. saturation may not be possible in one DT time.
          !  rain evaporation rate (W12),(OT25),(K Table 4). 
          !  evaporate rain first
          !  sd is still positive or we wouldn't be here.

          IF (QH (L) .LE.1.E-10) GOTO 33                                  

          !srf-25082005
          !  QUANT = ( QC (L)  + QV (L) - QSAT (L) ) * RHO (L)   !g/m**3
          QUANT = ( QSAT (L)- QC (L) - QV (L)   ) * RHO (L)   !g/m**3
          !
          EVHDT = (DT * HERC * (QUANT) * (QH (L) * RHO (L) ) **.65) / RHO (L)
          !             rain evaporation in time DT

          IF (EVHDT.LE.QH (L) ) THEN           !enough rain to last DT

             IF (SD.LE.EVHDT) THEN  		!enough time to saturate
                QH (L) = QH (L) - SD   	!remove rain	  
                QV (L) = QSAT (L)  		!saturate	  
                T (L) = T (L) - SD * FRC  	!cool the parcel	  

                RETURN  			!done
                !                       
             ELSE                               !not enough time
                SD = SD-EVHDT  		 !use what there is
                QV (L) = QV (L) + EVHDT  	 !add vapor
                T (L) = T (L) - EVHDT * FRC  	 !lose heat
                QH (L) = QH (L) - EVHDT  	 !lose rain

             ENDIF  				  !go on to ice.
             !                                    
          ELSE  !not enough rain to last DT
             !
             IF (SD.LE.QH (L) ) THEN             !but there is enough to sat
                QV (L) = QSAT (L)                !use it
                QH (L) = QH (L) - SD  
                T (L) = T (L) - SD * FRC  
                RETURN  
                !                            
             ELSE                              !not enough to sat
                SD = SD-QH (L)  
                QV (L) = QV (L) + QH (L)  
                T (L) = T (L) - QH (L) * FRC    
                QH (L) = 0.0                   !all gone

             ENDIF                             !on to ice
             !

          ENDIF                                !finished with rain
          !
          !
          !  now for ice
          !  equation from (OT); correction factors for units applied
          !
33        CONTINUE
          IF (QI (L) .LE.1.E-10) RETURN            !no ice there
          !
          DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (SD / QSAT (L) &
               - 1) * (QI (L) **0.525) * 1.13
          DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )  

          DEVIDT = - CVI(L) * DIVIDEND / DIVISOR   !rate of change

          EVIDT = DEVIDT * DT                      !what we could get
          !
          ! logic here is identical to rain. could get fancy and make 
          ! subroutine but duplication of code is easier. God bless the 
          ! screen editor.
          !

          IF (EVIDT.LE.QI (L) ) THEN             !enough ice to last DT
             !

             IF (SD.LE.EVIDT) THEN  		  !enough time to saturate
                QI (L) = QI (L) - SD   	  !remove ice
                QV (L) = QSAT (L)  		  !saturate
                T (L) = T (L) - SD * SRC  	  !cool the parcel

                RETURN  			  !done
                !

             ELSE                                !not enough time

                SD = SD-EVIDT  		  !use what there is
                QV (L) = QV (L) + EVIDT  	  !add vapor
                T (L) =  T (L) - EVIDT * SRC  	  !lose heat
                QI (L) = QI (L) - EVIDT  	  !lose ice

             ENDIF  				  !go on,unsatisfied
             !                                          
          ELSE                                   !not enough ice to last DT
             !                                         
             IF (SD.LE.QI (L) ) THEN             !but there is enough to sat

                QV (L) = QSAT (L)                !use it
                QI (L) = QI   (L) - SD  
                T (L) =  T   (L) - SD * SRC  

                RETURN  
                !
             ELSE                                 !not enough to sat
                SD = SD-QI (L)  
                QV (L) = QV (L) + QI (L)  
                T (L) = T (L) - QI (L) * SRC             
                QI (L) = 0.0                      !all gone

             ENDIF                                !on to better things
             !finished with ice
          ENDIF
          !                                 
       ENDIF                                   !finished with the SD decision
       !
       RETURN  
       !
     END SUBROUTINE EVAPORATE

     !
     ! *********************************************************************
     SUBROUTINE CONVERT (t,l,qc,rho,qh,dt)  
       !
       !- ACCRETION AND AUTOCONVERSION
       !

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL,    INTENT(IN)    :: dt
!!$       REAL,    INTENT(IN)    :: t  (nkp)
!!$       REAL,    INTENT(INOUT) :: qc (nkp)
!!$       REAL,    INTENT(IN)    :: rho(nkp)
!!$       REAL,    INTENT(INOUT) :: qh (nkp)
       REAL,    INTENT(IN)    :: t  (:)
       REAL,    INTENT(INOUT) :: qc (:)
       REAL,    INTENT(IN)    :: rho(:)
       REAL,    INTENT(INOUT) :: qh (:)


       REAL,    PARAMETER ::  AK1 = 0.001    !conversion rate constant
       REAL,    PARAMETER ::  AK2 = 0.0052   !collection (accretion) rate
       REAL,    PARAMETER ::  TH  = 0.5      !Kessler threshold
       INTEGER, PARAMETER ::  iconv = 1      !- Kessler conversion (=0)
       !real,   parameter :: ANBASE =  50.!*1.e+6   !Berry-number at cloud base #/m^3(maritime)
       REAL,    PARAMETER :: ANBASE =100000.!*1.e+6 !Berry-number at cloud base #/m^3(continental)
       !real,   parameter :: BDISP = 0.366          !Berry--size dispersion (maritime)
       REAL,    PARAMETER :: BDISP = 0.146          !Berry--size dispersion (continental)
       REAL,    PARAMETER :: TFREEZE = 269.3        !ice formation temperature
       !
       REAL               :: accrete, con, q, h, total


       IF (T (L)  .LE. TFREEZE) RETURN  !process not allowed above ice
       !
       IF (QC (L) .EQ. 0.     ) RETURN  

       ACCRETE = 0.  
       CON = 0.  
       Q = RHO (L) * QC (L)  
       H = RHO (L) * QH (L)  
       !
       !     selection rules
       !                         
       !            
       IF (QH (L) .GT. 0.     ) ACCRETE = AK2 * Q * (H**.875)  !accretion, Kessler
       !
       IF (ICONV.NE.0) THEN   !select Berry or Kessler
          !
          !old   BC1 = 120.  
          !old   BC2 = .0266 * ANBASE * 60.  
          !old   CON = BDISP * Q * Q * Q / (BC1 * Q * BDISP + BC2) 	  

          CON = Q*Q*Q*BDISP/(60.*(5.*Q*BDISP+0.0366*ANBASE))
          !
       ELSE  
          !                             
          !   CON = AK1 * (Q - TH)   !Kessler autoconversion rate
          !      
          !   IF (CON.LT.0.0) CON = 0.0   !havent reached threshold

          CON = MAX(0.,AK1 * (Q - TH)) ! versao otimizada
          !
       ENDIF
       !
       !
       TOTAL = (CON + ACCRETE) * DT / RHO (L)  

       !
       IF (TOTAL.LT.QC (L) ) THEN  
          !
          QC (L) = QC (L) - TOTAL  
          QH (L) = QH (L) + TOTAL    !no phase change involved
          RETURN  
          !
       ELSE  
          !              
          QH (L) = QH (L) + QC (L)    !uses all there is
          QC (L) = 0.0  
          !
       ENDIF
       !
       RETURN  
       !
     END SUBROUTINE CONVERT

     ! ice - effect on temperature
     !      TTD = 0.0 
     !      TTE = 0.0  
     !       CALL ICE(QSATW,QSATE,Y(1),Y(2),Y(3), &
     !               TTA,TTB,TTC,DZ,ROH,D,C,TTD,TTE)
     !       DYDX(1) = DYDX(1) + TTD  + TTE ! DT/DZ on Temp
     !
     !**********************************************************************
     !
     SUBROUTINE SUBLIMATE(t,l,qv,qsat,rho,qi,est,dt)  
       !
       ! ********************* VAPOR TO ICE (USE EQUATION OT22)***************

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL,    INTENT(IN)    :: dt
!!$       REAL,    INTENT(INOUT) :: t   (nkp)
!!$       REAL,    INTENT(INOUT) :: qv  (nkp)
!!$       REAL,    INTENT(IN)    :: qsat(nkp)
!!$       REAL,    INTENT(IN)    :: rho (nkp)
!!$       REAL,    INTENT(INOUT) :: qi  (nkp)
!!$       REAL,    INTENT(IN)    :: est (nkp)
       REAL,    INTENT(INOUT) :: t   (:)
       REAL,    INTENT(INOUT) :: qv  (:)
       REAL,    INTENT(IN)    :: qsat(:)
       REAL,    INTENT(IN)    :: rho (:)
       REAL,    INTENT(INOUT) :: qi  (:)
       REAL,    INTENT(IN)    :: est (:)


       !
       REAL, PARAMETER :: EPS = 0.622, HEATFUS = 334., HEATSUBL = 2834., CP = 1.004
       REAL, PARAMETER :: SRC = HEATSUBL / CP, FRC = HEATFUS / CP, TMELT = 273.3
       REAL, PARAMETER :: TFREEZE = 269.3

       REAL            :: dtsubh, dividend, divisor, subl
       !
       DTSUBH = 0.  
       !
       !selection criteria for sublimation
       IF (T (L)  .GT. TFREEZE  ) RETURN  
       IF (QV (L) .LE. QSAT (L) ) RETURN  
       !
       !     from (OT); correction factors for units applied
       !
       DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (QV (L) / QSAT (L) &
            - 1) * (QI (L) **0.525) * 1.13
       DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )  
       !

       DTSUBH = ABS (DIVIDEND / DIVISOR)   !sublimation rate
       SUBL = DTSUBH * DT                  !and amount possible
       !
       !     again check the possibilities
       !
       IF (SUBL.LT.QV (L) ) THEN  
          !
          QV (L) = QV (L) - SUBL             !lose vapor
          QI (L) = QI (L) + SUBL  	      !gain ice
          T (L) = T (L) + SUBL * SRC         !energy change, warms air

          !print*,'5',l,qi(l),SUBL

          RETURN  
          !
       ELSE  
          !                                     
          QI (L) = QV (L)                    !use what there is
          T  (L) = T (L) + QV (L) * SRC      !warm the air
          QV (L) = 0.0  
          !print*,'6',l,qi(l)
          !
       ENDIF
       !
       RETURN  
     END SUBROUTINE SUBLIMATE

     !
     ! *********************************************************************
     !
     SUBROUTINE GLACIATE(qh,l,qv,qsat,t,dt,qi)
       !
       ! *********************** CONVERSION OF RAIN TO ICE *******************
       !     uses equation OT 16, simplest. correction from W not applied, but
       !     vapor pressure differences are supplied.
       !

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL   , INTENT(IN)    :: dt
!!$       REAL   , INTENT(INOUT) :: qh  (nkp)
!!$       REAL   , INTENT(IN)    :: qv  (nkp)
!!$       REAL   , INTENT(IN)    :: qsat(nkp)
!!$       REAL   , INTENT(INOUT) :: t   (nkp)
!!$       REAL   , INTENT(INOUT) :: qi  (nkp)
       REAL   , INTENT(INOUT) :: qh  (:)
       REAL   , INTENT(IN)    :: qv  (:)
       REAL   , INTENT(IN)    :: qsat(:)
       REAL   , INTENT(INOUT) :: t   (:)
       REAL   , INTENT(INOUT) :: qi  (:)


       !
       REAL, PARAMETER :: HEATFUS = 334., CP = 1.004, EPS = 0.622, HEATSUBL = 2834.
       REAL, PARAMETER :: FRC = HEATFUS / CP, FRS = HEATSUBL / CP, TFREEZE =  269.3
       REAL, PARAMETER :: GLCONST = 0.025   !glaciation time constant, 1/sec
       REAL            :: dfrzh
       !

       DFRZH = 0.    !rate of mass gain in ice
       !
       !selection rules for glaciation
       IF (QH (L) .LE. 0.       ) RETURN  
       IF (QV (L) .LT. QSAT (L) ) RETURN                                        
       IF (T  (L) .GT. TFREEZE  ) RETURN  
       !
       !      NT=TMELT-T(L)
       !      IF (NT.GT.50) NT=50
       !

       DFRZH = DT * GLCONST * QH (L)    ! from OT(16)
       !
       IF (DFRZH.LT.QH (L) ) THEN  
          !
          QI (L) = QI (L) + DFRZH  
          QH (L) = QH (L) - DFRZH  
          T (L) = T (L) + FRC * DFRZH  !warms air

          !print*,'7',l,qi(l),DFRZH


          RETURN  
          !
       ELSE  
          !
          QI (L) = QI (L) + QH (L)  
          T  (L) = T  (L) + FRC * QH (L)  
          QH (L) = 0.0  

          !print*,'8',l,qi(l), QH (L)  
          !
       ENDIF
       !
       RETURN  
       !
     END SUBROUTINE GLACIATE

     !
     !
     ! *********************************************************************
     SUBROUTINE MELT(qi,l,t,dt,rho,cvi,qh)  
       !
       ! ******************* MAKES WATER OUT OF ICE **************************

       ! plumegen_coms
       INTEGER, INTENT(IN)    :: l
       REAL   , INTENT(IN)    :: dt
!!$       REAL   , INTENT(INOUT) :: qi (nkp)
!!$       REAL   , INTENT(INOUT) :: t  (nkp)
!!$       REAL   , INTENT(IN)    :: rho(nkp)
!!$       REAL   , INTENT(IN)    :: cvi(nkp)
!!$       REAL   , INTENT(INOUT) :: qh (nkp)
       REAL   , INTENT(INOUT) :: qi (:)
       REAL   , INTENT(INOUT) :: t  (:)
       REAL   , INTENT(IN)    :: rho(:)
       REAL   , INTENT(IN)    :: cvi(:)
       REAL   , INTENT(INOUT) :: qh (:)


       !                                              
       REAL, PARAMETER :: FRC = 332.27, TMELT = 273., F0 = 0.75   !ice velocity factor
       REAL            :: DTMELT
       !                                    
       DTMELT = 0.   !conversion,ice to rain
       !
       !selection rules
       IF (QI (L) .LE. 0.0  ) RETURN  
       IF (T (L)  .LT. TMELT) RETURN  
       !
       !OT(23,24)
       DTMELT = DT * (2.27 / RHO (L) ) * CVI(L) * (T (L) - TMELT) * ( (RHO(L)  &
            * QI (L) * 1.E-6) **0.525) * (F0** ( - 0.42) )
       !after Mason,1956
       !
       !     check the possibilities
       !
       IF (DTMELT.LT.QI (L) ) THEN  
          !
          QH (L) = QH (L) + DTMELT  
          QI (L) = QI (L) - DTMELT  
          T  (L) = T (L) - FRC * DTMELT     !cools air
          !print*,'9',l,qi(l),DTMELT


          RETURN  
          !
       ELSE  
          !
          QH (L) = QH (L) + QI (L)   !get all there is to get
          T  (L) = T (L) - FRC * QI (L)  
          QI (L) = 0.0  
          !print*,'10',l,qi(l)
          !
       ENDIF
       !
       RETURN  
       !
     END SUBROUTINE MELT

     !-------------------------------------------------------------------------
     FUNCTION ESAT_PR (TEM)  

       REAL, INTENT(IN) :: TEM


       !
       ! ******* Vapor Pressure  A.L. Buck JAM V.20 p.1527. (1981) ***********
       !
       REAL, PARAMETER :: CI1 = 6.1115, CI2 = 22.542, CI3 = 273.48
       REAL, PARAMETER :: CW1 = 6.1121, CW2 = 18.729, CW3 = 257.87, CW4 = 227.3
       REAL, PARAMETER :: TMELT = 273.3

       REAL            :: ESAT_PR, temc, esatm
       !
       !     formulae from Buck, A.L., JAM 20,1527-1532
       !     custom takes esat wrt water always. formula for h2o only
       !     good to -40C so:
       !
       !
       TEMC = TEM - TMELT  
       IF (TEMC.GT. - 40.0) THEN
          ESATM = CW1 * EXP ( ( (CW2 - (TEMC / CW4) ) * TEMC) / (TEMC + CW3))
          ESAT_PR = ESATM / 10.	!kPa			  
       ELSE
          ESATM = CI1 * EXP (CI2 * TEMC / (TEMC + CI3) )  !ice, millibars  
          ESAT_PR = ESATM / 10.	!kPa			  
       END IF

       RETURN  

     END FUNCTION ESAT_PR


   

END MODULE mod_chem_plumerise_scalar
