Module module_chem_plumerise_scalar
! USE module_plumerise1
  use module_model_constants
  use module_zero_plumegen_coms
  real,parameter :: rgas=r_d
  real,parameter :: cpor=1./rcp
  real,parameter :: p00=p1000mb
!  real, external:: esat_pr
CONTAINS
subroutine plumerise(m1,m2,m3,ia,iz,ja,jz,firesize,mean_fct   &
                    ,nspecies,eburn_in,eburn_out      &
                    ,up,vp,wp,theta,pp,dn0,rv,zt_rams,zm_rams)
! use module_zero_plumegen_coms, only : ucon,vcon,wcon,thtcon,rvcon,picon,tmpcon&
!                          ,dncon,prcon,zcon,zzcon,scon
  
  implicit none


  integer :: ng,m1,m2,m3,ia,iz,ja,jz,ibcon,mynum,i,j,k,iveg_ag,&
            imm,k1,k2,ixx,ispc,nspecies

  integer :: ncall = 0 
  integer :: kmt
 real,dimension(m1,nspecies), intent(out) :: eburn_out
 real,dimension(nspecies), intent(in) :: eburn_in

  real,    dimension(m1,m2,m3) :: up, vp, wp,theta,pp,dn0,rv
				 
  real,    dimension(m1)       :: zt_rams,zm_rams

  real                         :: burnt_area,STD_burnt_area,dz_flam,rhodzi,dzi
  real,    dimension(2)        :: ztopmax

  real                         :: q_smold_kgm2

! From plumerise1.F routine
  integer, parameter :: nveg_agreg      = 4
  integer, parameter :: tropical_forest = 1
  integer, parameter :: boreal_forest   = 2
  integer, parameter :: savannah        = 3

  integer, parameter :: grassland       = 4
  real, dimension(nveg_agreg) :: firesize,mean_fct

  INTEGER, PARAMETER :: wind_eff = 1


  !Fator de conversao de unidades
  !!fcu=1.	  !=> kg [gas/part] /kg [ar]
  !!fcu =1.e+12   !=> ng [gas/part] /kg [ar]
  !!real,parameter :: fcu =1.e+6 !=> mg [gas/part] /kg [ar] 
  !----------------------------------------------------------------------
  !               indexacao para o array "plume(k,i,j)" 
  ! k 
  ! 1   => area media (m^2) dos focos  em biomas floresta dentro do gribox i,j 
  ! 2   => area media (m^2) dos focos  em biomas savana dentro do gribox i,j
  ! 3   => area media (m^2) dos focos  em biomas pastagem dentro do gribox i,j
  ! 4   => desvio padrao da area media (m^2) dos focos : floresta
  ! 5   => desvio padrao da area media (m^2) dos focos : savana
  ! 6   => desvio padrao da area media (m^2) dos focos : pastagem
  ! 7 a 9 =>  sem uso
  !10(=k_CO_smold) => parte da emissao total de CO correspondente a fase smoldering
  !11, 12 e 13 => este array guarda a relacao entre
  !               qCO( flaming, floresta) e a quantidade total emitida 
  !               na fase smoldering, isto e;
  !               qCO( flaming, floresta) =  plume(11,i,j)*plume(10,i,j)
  !               qCO( flaming, savana  ) =  plume(12,i,j)*plume(10,i,j)
  !               qCO( flaming, pastagem) =  plume(13,i,j)*plume(10,i,j)
  !20(=k_PM25_smold),21,22 e 23 o mesmo para PM25               
  !
  !24-n1 =>  sem uso
  !----------------------------------------------------------------------
! print *,' Plumerise_scalar 1',ncall
  if (ncall == 0) then
    ncall = 1
    call zero_plumegen_coms
  endif
! print *,' Plumerise_scalar 1',ncall
  
    
! print *,' Plumerise_scalar 2',m1
  j=1
  i=1
! do j = ja,jz           ! loop em j
!   do i = ia,iz         ! loop em i
     
     !- if the max value of flaming is close to zero => there is not emission with
     !- plume rise => cycle

 	 do k = 1,m1
 	   ucon  (k)=up(k,i,j) 	      ! u wind
 	   vcon  (k)=vp(k,i,j) 	      ! v wind
 	   !wcon  (k)=wp(k,i,j) 	      ! w wind
 	   thtcon(k)=theta(k,i,j)	      ! pot temperature
 	   picon (k)=pp(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)
 	   !prcon (k)=(picon(k)/cp)**cpor*p00 ! pressure (Pa)
 	   rvcon (k)=rv(k,i,j)  	      ! water vapor mixing ratio
 	   zcon  (k)=zt_rams(k)               ! termod-point height
 	   zzcon (k)=zm_rams(k)               ! W-point height
!	   print*,'PL:',k,zcon(k),picon (k),thtcon(k),1000.*rvcon (k)
 	 enddo
         do ispc=1,nspecies
	      eburn_out(1,ispc) = eburn_in(ispc)
         enddo
 	
         !- get envinronmental state (temp, water vapor mix ratio, ...)
 	 call get_env_condition(1,m1,kmt,wind_eff)
 	    
         !- loop nos 4 biomas agregados com possivel queimada 
 	 do iveg_ag=1,nveg_agreg
!         print *,'iveg_ag = ',iveg_ag,mean_fct(iveg_ag)
  

         !- verifica a existencia de emissao flaming para um bioma especifico
    	 !orig: if( plume( k_CO_smold + iveg_ag ,i,j) < 1.e-6 ) cycle
           if(mean_fct(iveg_ag) < 1.e-6 ) cycle
 
	   ! burnt area and standard deviation
	   burnt_area    = firesize(iveg_ag)
	      
	   !not em use
	   !STD_burnt_area= plume(3+iveg_ag,i,j)
	    STD_burnt_area= 0.
	    
	   !- loop nos valores  minimo e maximo da taxa de calor
	   do imm=1,2
	   
!--------------------
                !ixx=iveg_ag*10 + imm
!               print*,'i j veg=',i, j, iveg_ag,imm
!--------------------             
	      
	     !- get fire properties (burned area, plume radius, heating rates ...)
 	     call get_fire_properties(imm,iveg_ag,burnt_area,STD_burnt_area)
     
             !------  generates the plume rise    ------
	     
	     !-- only one value for eflux of GRASSLAND
       !     if(iveg_ag == GRASSLAND .and. imm == 2) then 
	     if(iveg_ag == 4         .and. imm == 2) then 
 	        ztopmax(2)=ztopmax(1)
		ztopmax(1)=zzcon(1)
!               print *,'cycle',ztopmax(1),ztopmax(2)
                cycle
	     endif
	     
	     call makeplume (kmt,ztopmax(imm),ixx,imm)	      	      
	     	   
	   enddo ! enddo do loop em imm
	   
	   !- define o dominio vertical onde a emissao flaming ira ser colocada
	   call set_flam_vert(ztopmax,k1,k2,nkp,zzcon,W_VMD,VMD)
	   
	   !- espessura da camada vertical (compute the vertical layer thickness)
	   dz_flam=zzcon(k2+1)-zzcon(k1)
	   
	   !- distribui a emissao flaming entre os niveis k1 e k2	 
!          print *,'distribui, k1,k2,dz_flam',k1,k2,dz_flam
	   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)
	   
	    do ispc = 1, nspecies
	       
	      !- 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 = (1/dzt(2) *  dn0(2,i,j)	)*   &
	      ! 	     chem1_src_g(bburn,ispc,ng)%sc_src(2,i,j)	 
	      
	      !use this in case the emission src is tracer density
!       q_smold_kgm2 = ((zt_rams(2)-zt_rams(1))                 )*   &
!        	     eburn_in(ispc)
        q_smold_kgm2 = eburn_in(ispc)
	
	      ! units = already in ppbm,  don't need "fcu" factor 
	      eburn_out(k,ispc) = eburn_out(k,ispc) +&
	    						 mean_fct(iveg_ag)  *&
	        					 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
!             if(ispc.eq. 1) print *,' Plume: ',k,eburn_out(k,ispc),mean_fct(iveg_ag),q_smold_kgm2,dzi
							 
							 

	      !srcCO(k,i,j)=   srcCO(k,i,j) + plume(k_CO_smold+iveg_ag,i,j)*&
	      ! 			    plume(k_CO_smold	    ,i,j)*&
	      ! 			    rhodzi*fcu
	    enddo
	   
	   enddo
	   
         enddo ! enddo do loop em iveg_ag
     !-----  
 	 !stop 333
 	 !endif  	
     !-----
!      enddo  ! loop em i
! enddo   ! loop em j
!stop 4544
end subroutine plumerise
!-------------------------------------------------------------------------

subroutine get_env_condition(k1,k2,kmt,wind_eff)

!se module_zero_plumegen_coms
!use rconstants
implicit none
integer :: k1,k2,k,kcon,klcl,kmt,nk,nkmid,i
real :: znz,themax,tlll,plll,rlll,zlll,dzdd,dzlll,tlcl,plcl,dzlcl,dummy
integer :: n_setgrid = 0 
integer :: wind_eff


if( n_setgrid == 0) then
  n_setgrid = 1
  call set_grid ! define vertical grid of plume model
                ! zt(k) =  thermo and water levels
                ! zm(k) =  dynamical levels 
endif

znz=zcon(k2)
do k=nkp,1,-1
  if(zt(k).lt.znz)go to 13
enddo
stop ' envir stop 12'
13 continue
!-srf-mb
kmt=min(k,nkp-1)

nk=k2-k1+1
!call htint(nk, wcon,zzcon,kmt,wpe,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)
!  print*,'ENV=',qvenv(k)*1000., te(k)-273.15,zt(k)
!-srf-mb
  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.
! 100 format(1x,I5,4f20.12)
!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 
end subroutine get_env_condition

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

subroutine set_grid()
!use module_zero_plumegen_coms  
implicit none
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,nkp,zzcon,W_VMD,VMD)

    REAL    , INTENT(IN)  :: ztopmax(2)
    INTEGER , INTENT(OUT) :: k1
    INTEGER , INTENT(OUT) :: k2

    ! plumegen_coms
    INTEGER , INTENT(IN)  :: nkp
    REAL    , INTENT(IN)  :: zzcon(nkp)

    INTEGER imm,k
    INTEGER, DIMENSION(2)  :: k_lim

    !- version 2
    REAL    , INTENT(IN)  :: W_VMD(nkp,2)
    REAL    , INTENT(OUT) ::   VMD(nkp,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,STD_burnt_area)
!use module_zero_plumegen_coms  
implicit none
integer ::  moist,  i,  icount,imm,iveg_ag
real::   bfract,  effload,  heat,  hinc ,burnt_area,STD_burnt_area,heat_fluxW
real,    dimension(2,4) :: heat_flux
INTEGER, parameter :: use_last = 0


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
!--------------------------------------------------------------------


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

!fluxo de calor para o bioma
heat_fluxW = heat_flux(imm,iveg_ag) * 1000. ! converte para W/m^2

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
!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"'

  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 
       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

return
end subroutine get_fire_properties
!-------------------------------------------------------------------------------
!
SUBROUTINE MAKEPLUME ( kmt,ztopmax,ixx,imm)  
!
! *********************************************************************
!
!    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)
!
!
!**********************************************************************
!**********************************************************************               
!use module_zero_plumegen_coms 
implicit none 
!logical :: endspace  
character (len=10) :: varn
integer ::  izprint, iconv,  itime, k, kk, kkmax, deltak,ilastprint,kmt &
           ,ixx,nrectotal,i_micro,n_sub_step
real ::  vc, g,  r,  cp,  eps,  &
         tmelt,  heatsubl,  heatfus,  heatcond, tfreeze, &
         ztopmax, wmax, rmaxtime, es, esat, heat,dt_save !ESAT_PR,
character (len=2) :: cixx	 
! Set threshold to be the same as dz=100., the constant grid spacing of plume grid model(meters) found in set_grid()
    REAL :: DELZ_THRESOLD = 100. 

    INTEGER     :: imm

!  real, external:: esat_pr!
!
! ******************* 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)  
!
tstpf = 2.0  	!- timestep factor
viscosity = 500.!- viscosity constant (original value: 0.001)

nrectotal=150
!
!*************** 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

!--- initialization
CALL INITIAL(kmt)  

!--- initial print fields:
izprint  = 0          ! if = 0 => no printout
if (izprint.ne.0) then
 write(cixx(1:2),'(i2.2)') ixx
 open(2, file = 'debug.'//cixx//'.dat')  
 open(19,file='plumegen9.'//cixx//'.gra',         &
     form='unformatted',access='direct',status='unknown',  &
     recl=4*nrectotal)  !PC   
!     recl=1*nrectotal) !sx6 e tupay
 call printout (izprint,nrectotal)
 ilastprint=2
endif     

! ******************* model evolution ******************************
rmaxtime = float(maxtime)
!
!print * ,' TIME=',time,' RMAXTIME=',rmaxtime
!print*,'======================================================='
 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))
                                
!-- 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

!-- bounday conditions (k=1)
    L=1
    call lbound()

!-- 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('SC',NM1)

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

!-- scalars entrainment, adiabatic
    call scl_misc(NM1)
    
!-- scalars dinamic entrainment
     call  scl_dyn_entrain(NM1,nkp,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,qv,qh,qi,qc,te,pe,qvenv)

!-- 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)
!-- microphysics
     do L=2,nm1-1
        WBAR    = 0.5*(W(L)+W(L-1))
        ES      = ESAT_PR (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  
     enddo
    enddo
    dt=dt_save
!
    101 continue
!
!-- W-viscosity for stability 
    call visc_W(nm1,deltak,kmt)

!-- update scalars
    call update_plumerise(nm1,'S')
    
    call hadvance_plumerise(1,nm1,dt,WC,WT,W,mintime) 

!-- Buoyancy
    call buoyancy_plumerise(NM1, T, TE, QV, QVENV, QH, QI, QC, WT, SCR1)
 
!-- Entrainment 
    call entrainment(NM1,W,WT,RADIUS,ALPHA)

!-- update W
    call update_plumerise(nm1,'W')

    call hadvance_plumerise(2,nm1,dt,WC,WT,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,qv,qh,qi,qc,te,pe,qvenv)
!---
       !- 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) then
      call printout (izprint,nrectotal)  
      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)  
 close (2)            
 close (19)            
endif

RETURN  
END SUBROUTINE MAKEPLUME
!-------------------------------------------------------------------------------
!
SUBROUTINE BURN(EFLUX, WATER)  
!	
!- calculates the energy flux and water content at lboundary
!use module_zero_plumegen_coms                               
!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)
real :: eflux,water
!
! 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 ()  
!
! ********** 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
!
!use module_zero_plumegen_coms  
implicit none
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,  esat, eflux, water,  pres, c1,  c2, f, zv,  denscor, xwater !,ESAT_PR
!  real, external:: esat_pr!

!            
QH (1) = QH (2)   !soak up hydrometeors
QI (1) = QI (2)              
QC (1) = 0.       !no cloud here
!
!
   CALL BURN (EFLUX, WATER)  
!
!  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  
!
RETURN  
END SUBROUTINE LBOUND
!-------------------------------------------------------------------------------
!
SUBROUTINE INITIAL ( kmt)  
!
! ************* SETS UP INITIAL CONDITIONS FOR THE PROBLEM ************
!use module_zero_plumegen_coms 
implicit none 
real, parameter :: tfreeze = 269.3
integer ::  isub,  k,  n1,  n2,  n3,  lbuoy,  itmp,  isubm1 ,kmt
real ::     xn1,  xi,  es,  esat!,ESAT_PR
!
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
  do k=2,N
     radius(k) = radius(k-1)+(6./5.)*alpha*(zt(k)-zt(k-1))
  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()

RETURN  
END SUBROUTINE INITIAL
!-------------------------------------------------------------------------------
!
subroutine damp_grav_wave(ifrom,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
implicit none
integer nm1,ifrom,deltak
real dt
real, dimension(nm1) :: w,t,tt,qv,qh,qi,qc,te,pe,qvenv,dummy,zt,zm

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)
 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)
implicit none
real, dimension(nm1) :: var1,var2,vart,zt,zm
integer k,nfpt,kf,nm1,ifrom,deltak
real zmkf,ztop,distim,c1,c2,dt

!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.
 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)

implicit none
integer :: k,m1
real, dimension(m1) :: wc,wt,flxw,dzm,rho
real, dimension(m1) :: dn0 ! 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,dt,wc,wt,wp,mintime)

implicit none
integer :: k,iac
integer :: m1,mintime
real, dimension(m1) :: dummy, wc,wt,wp
real eps,dt
!     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,wt,dummy,iac,2.*dt,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,fa,af,iac,dtlp,epsu)
implicit none
integer :: npts,iac,m
real :: epsu,dtlp
real, dimension(*) :: ac,ap,fa,af

!     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
!   if (ngrid .eq. 1 .and. ipara .eq. 0) call cyclic(nzp,nxp,nyp,af,'T')
!   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)
implicit none
integer :: k,m1
real, parameter :: g = 9.8, eps = 0.622, gama = 0.5 ! mass virtual coeff.
real, dimension(m1) :: T, TE, QV, QVENV, QH, QI, QC, WT, scr1
real :: TV,TVE,QWTOTL,umgamai
real, parameter :: mu = 0.15 

!- orig
umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
                       ! das 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)
implicit none
integer :: k,m1
real, dimension(m1) :: w,wt,radius
REAL DMDTM,WBAR,RADIUS_BAR,umgamai,DYN_ENTR,ALPHA
real, parameter :: mu = 0.15 ,gama = 0.5 ! mass virtual coeff.

!- 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
                       ! das 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=0.5*(W(k)+W(k-1))           
      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)
      !print*,'W-ENTR=',k,w(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(varn,mzp)
!use module_zero_plumegen_coms
implicit none
integer :: mzp
character(len=*) :: varn
real :: dtlto2
integer :: k

!  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 advtndc_plumerise(mzp,T,scr1(1),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 advtndc_plumerise(mzp,QV,scr1(1),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 advtndc_plumerise(mzp,QC,scr1(1),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 advtndc_plumerise(mzp,QI,scr1(1),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 advtndc_plumerise(mzp,QH,scr1(1),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 advtndc_plumerise(mzp,VEL_P,scr1(1),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 advtndc_plumerise(mzp,rad_p,scr1(1),rad_t,dt)


   return
!
!- gas/particle advection tendency (SCT)
!    if(varn == 'SC')return
   scr1=SC
   call fa_zc_plumerise(mzp		    &
   	     	       ,SC	 ,scr1  (1)  &
   	     	       ,vt3dc (1) ,vt3df (1)  &
   	     	       ,vt3dg (1) ,vt3dk (1)  &
   	     	       ,vctr1,vctr2	     )
   
   call advtndc_plumerise(mzp,SC,scr1(1),SCT,dt)


return
end subroutine scl_advectc_plumerise
!-------------------------------------------------------------------------------
!
subroutine fa_zc_plumerise(m1,scp,scr1,vt3dc,vt3df,vt3dg,vt3dk,vctr1,vctr2)

implicit none
integer :: m1,k
real :: dfact
real, dimension(m1) :: scp,scr1,vt3dc,vt3df,vt3dg,vt3dk
real, dimension(m1) :: vctr1,vctr2

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)
implicit none
integer :: m1,k
real :: dtl,dtli
real, dimension(m1) :: scp,sca,sct

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
!use module_zero_plumegen_coms, only: nm1,wt,tt,qvt,qct,qht,qit,sct
 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)
!use module_zero_plumegen_coms
implicit none
real, parameter :: g = 9.81, cp=1004.
integer m1,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,nkp,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)
    implicit none

    INTEGER , INTENT(IN)    :: m1

    ! plumegen_coms
    INTEGER , INTENT(IN)    :: nkp
    REAL    , INTENT(INOUT) :: wbar 
    REAL    , INTENT(IN)    :: w(nkp)
    REAL    , INTENT(INOUT) :: adiabat 
    REAL    , INTENT(IN)    :: alpha
    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, 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,deltak,kmt)
!use module_zero_plumegen_coms
implicit none
integer m1,k,deltak,kmt,m2
real dz1t,dz1m,dz2t,dz2m,d2wdz,d2tdz  ,d2qvdz ,d2qhdz ,d2qcdz ,d2qidz ,d2scdz, &
 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,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)
!use module_zero_plumegen_coms
integer m1,k
character(len=*) :: varn
 
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)
!use module_zero_plumegen_coms
integer m1,k
real vtc, dfhz,dfiz,dz1
!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.
!
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)  
!use module_zero_plumegen_coms  
real, parameter :: tmelt = 273.3
integer, save :: nrec
data nrec/0/
integer :: ko,izprint,interval,nrectotal
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,F7.2,F7.2,F6.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  
!use module_zero_plumegen_coms  
!
                                        
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    !vapor to cloud,cloud to vapor  
!                             
CALL SUBLIMATE    !vapor to ice  
!                            
CALL GLACIATE     !rain to ice 
                           
CALL MELT         !ice to rain
!         
!if(ak1 > 0. .or. ak2 > 0.) &
CALL CONVERT () !(auto)conversion and accretion 
!CALL CONVERT2 () !(auto)conversion and accretion 
!

RETURN  
END SUBROUTINE WATERBAL
! *********************************************************************
SUBROUTINE EVAPORATE  
!
!- evaporates cloud,rain and ice to saturation
!
!use module_zero_plumegen_coms  
implicit none
!
!     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 ()  
!
!- ACCRETION AND AUTOCONVERSION
!
!use module_zero_plumegen_coms  
!
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, bc1,   bc2,  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
!
!**********************************************************************
!
SUBROUTINE CONVERT2 ()  
!use module_zero_plumegen_coms  
implicit none
LOGICAL  AEROSOL
parameter(AEROSOL=.true.)
!
real, parameter :: TNULL=273.16, LAT=2.5008E6 &
                  ,EPSI=0.622 ,DB=1. ,NB=1500. !ALPHA=0.2 
real :: KA,KEINS,KZWEI,KDREI,VT	
real :: A,B,C,D, CON,ACCRETE,total 
      
real Y(6),ROH
      
A=0.
B=0.
Y(1) = T(L)
Y(4) = W(L)
y(2) = QC(L)
y(3) = QH(L)
Y(5) = RADIUS(L)
ROH =  RHO(L)*1.e-3 ! dens (MKS) ??


! autoconversion

KA = 0.0005 
IF( Y(1) .LT. 258.15 )THEN
!   KEINS=0.00075
    KEINS=0.0009 
    KZWEI=0.0052
    KDREI=15.39
ELSE
    KEINS=0.0015
    KZWEI=0.00696
    KDREI=11.58
ENDIF
      
!   ROH=PE/RD/TE
VT=-KDREI* (Y(3)/ROH)**0.125

 
IF (Y(4).GT.0.0 ) THEN
 IF (AEROSOL) THEN
   A = 1/y(4)  *  y(2)*y(2)*1000./( 60. *( 5. + 0.0366*NB/(y(2)*1000.*DB) )  )
 ELSE
   IF (y(2).GT.(KA*ROH)) THEN
   !print*,'1',y(2),KA*ROH
   A = KEINS/y(4) *(y(2) - KA*ROH )
   ENDIF
 ENDIF
ELSE
   A = 0.0
ENDIF

! accretion

IF(y(4).GT.0.0) THEN
   B = KZWEI/(y(4) - VT) * MAX(0.,y(2)) *   &
       MAX(0.001,ROH)**(-0.875)*(MAX(0.,y(3)))**(0.875)
ELSE
   B = 0.0
ENDIF
   
   
      !PSATW=610.7*EXP( 17.25 *( Y(1) - TNULL )/( Y(1)-36. ) )
      !PSATE=610.7*EXP( 22.33 *( Y(1) - TNULL )/( Y(1)- 2. ) )

      !QSATW=EPSI*PSATW/( PE-(1.-EPSI)*PSATW )
      !QSATE=EPSI*PSATE/( PE-(1.-EPSI)*PSATE )
      
      !MU=2.*ALPHA/Y(5)

      !C = MU*( ROH*QSATW - ROH*QVE + y(2) )
      !D = ROH*LAT*QSATW*EPSI/Y1/Y1/RD *DYDX1

      
      !DYDX(2) = - A - B - C - D  ! d rc/dz
      !DYDX(3) = A + B            ! d rh/dz
 
 
      ! rc=rc+dydx(2)*dz
      ! rh=rh+dydx(3)*dz
 
CON      = A
ACCRETE  = B
 
TOTAL = (CON + ACCRETE) *(1/DZM(L)) /ROH     ! DT / RHO (L)  

!print*,'L=',L,total,QC(L),dzm(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 CONVERT2
! 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  
!
! ********************* VAPOR TO ICE (USE EQUATION OT22)***************
!use module_zero_plumegen_coms  
!
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  
!
! *********************** CONVERSION OF RAIN TO ICE *******************
!     uses equation OT 16, simplest. correction from W not applied, but
!     vapor pressure differences are supplied.
!
!use module_zero_plumegen_coms  
!
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  
!
! ******************* MAKES WATER OUT OF ICE **************************
!use module_zero_plumegen_coms  
!                                              
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

SUBROUTINE htint (nzz1, vctra, eleva, nzz2, vctrb, elevb)
  IMPLICIT NONE
  INTEGER, INTENT(IN ) :: nzz1
  INTEGER, INTENT(IN ) :: nzz2
  REAL,    INTENT(IN ) :: vctra(nzz1)
  REAL,    INTENT(OUT) :: vctrb(nzz2)
  REAL,    INTENT(IN ) :: eleva(nzz1)
  REAL,    INTENT(IN ) :: elevb(nzz2)

  INTEGER :: l
  INTEGER :: k
  INTEGER :: kk
  REAL    :: wt

  l=1

  DO k=1,nzz2
     DO
        IF ( (elevb(k) <  eleva(1)) .OR. &
             ((elevb(k) >= eleva(l)) .AND. (elevb(k) <= eleva(l+1))) ) THEN
           wt       = (elevb(k)-eleva(l))/(eleva(l+1)-eleva(l))
           vctrb(k) = vctra(l)+(vctra(l+1)-vctra(l))*wt
           EXIT
        ELSE IF ( elevb(k) >  eleva(nzz1))  THEN
           wt       = (elevb(k)-eleva(nzz1))/(eleva(nzz1-1)-eleva(nzz1))
           vctrb(k) = vctra(nzz1)+(vctra(nzz1-1)-vctra(nzz1))*wt
           EXIT
        END IF

        l=l+1
        IF(l == nzz1) THEN
           PRINT *,'htint:nzz1',nzz1
           DO kk=1,l
              PRINT*,'kk,eleva(kk),elevb(kk)',eleva(kk),elevb(kk)
           END DO
           STOP 'htint'
        END IF
     END DO
  END DO
END SUBROUTINE htint
!-----------------------------------------------------------------------------
FUNCTION ESAT_PR (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
real temc , tem,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) GOTO 230  
ESATM = CI1 * EXP (CI2 * TEMC / (TEMC + CI3) )  !ice, millibars  
ESAT_PR = ESATM / 10.	!kPa			  

RETURN  
!
230 ESATM = CW1 * EXP ( ( (CW2 - (TEMC / CW4) ) * TEMC) / (TEMC + CW3))
                          
ESAT_PR = ESATM / 10.	!kPa			  
RETURN  
END function ESAT_PR
!     ******************************************************************

!      ------------------------------------------------------------------------
END Module module_chem_plumerise_scalar
