!
! Copyright (C) 1991-2004  ; All Rights Reserved ; Colorado State University
! Colorado State University Research Foundation ; ATMET, LLC
! 
! This file is free software; you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software 
! Foundation; either version 2 of the License, or (at your option) any later version.
! 
! This software is distributed in the hope that it will be useful, but WITHOUT ANY 
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 
! PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this 
! code; if not, write to the Free Software Foundation, Inc., 
! 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!======================================================================================

Subroutine eng_params ()

use mem_grid
use io_params

implicit none

!  Set some constants that were formerly defined in the namelist.
!  This set should be changed only for special tests or with code modification.

SSPCT   = 0.  ! Initial sound speed fraction (overwritten after t=0)
IMPL    = 1   ! Implicit flag for acoustic model  -  0=off, 1=on

return
END SUBROUTINE eng_params

!##############################################################################
Subroutine toptinit_user (n2,n3,ifm,topt,topzo)

use mem_grid

implicit none

integer :: i,j,n2,n3,ifm,flag_on
real :: hfwid,hfwid2,hgt
real, dimension(n2,n3) :: topt,topzo

!  This routine is the intended location for a user to customize TOPT,
!  the surface topography array.  It is called after all other types of
!  initialization of this field, so this routine has the last word.
!  By default the routine makes no change to the field.  The commented
!  lines below serve as a template for user-designed changes; the example
!  shown is the Witch of Agnesi mountain. Note that
!  this routine is called for each grid separately, so attention to the
!  current value of ngrid in this routine may be required by the user.

flag_on=1
if(flag_on==1)then
      topt(i,j)=16.

!      if(ifm.eq.1)then
!         hfwid=5000.
!         hgt=1000.
!         hfwid2=hfwid**2
!         do j=1,n3
!            do i=1,n2
!               topt(i,j)=hgt*hfwid2/(hfwid2+xtn(i,1)**2)
!               topzo(i,j) = 0.001 
!            enddo
!         enddo
!      elseif(ifm.eq.2)then
!      endif

endif

return
END SUBROUTINE toptinit_user

!##############################################################################
Subroutine sstinit_user (n2,n3,ifm,seatf)

implicit none

integer :: n2,n3,ifm,i,j,flag_on
real, dimension(n2,n3) :: seatf

!  This routine is the intended location for a user to customize the
!  SEATP and SEATF arrays.  It is called after all other types of 
!  initialization of these fields, so this routine has the last word.  
!  By default the routine makes no change to the fields.  The commented 
!  lines below serve as a template for user-designed changes.   Note that 
!  this routine is called for each grid separately, so attention to the 
!  current value of ngrid in this routine may be required by the user.

flag_on=0
if(flag_on==1)then

 if (ifm .eq. 1) then
    do j = 1,n3
       do i = 1,n2
          seatf(i,j) = 0.
       enddo
    enddo
 endif

endif

return
END SUBROUTINE sstinit_user

!##############################################################################
Subroutine ndviinit_user (n2,n3,npat,ifm,veg_ndvif)

implicit none

integer :: n2,n3,ipat,npat,ifm,i,j,flag_on
real, dimension(n2,n3,npat) :: veg_ndvif

!  This routine is the intended location for a user to customize the
!  NDVIP and NDVIF arrays.  It is called after all other types of 
!  initialization of these fields, so this routine has the last word.  
!  By default the routine makes no change to the fields.  The commented 
!  lines below serve as a template for user-designed changes.   Note that 
!  this routine is called for each grid separately, so attention to the 
!  current value of ngrid in this routine may be required by the user.

flag_on=0
if(flag_on==1)then

 if (ifm .eq. 1) then
    do j = 1,n3
       do i = 1,n2
          veg_ndvif(i,j,1) = 0.
          veg_ndvif(i,j,2) = 0.
          do ipat = 3,npat
             veg_ndvif(i,j,ipat) = 0.
          enddo
       enddo
    enddo
 endif

endif

return
END SUBROUTINE ndviinit_user

!##############################################################################
Subroutine sfcinit_file_user (n2,n3,nzg,npat,ifm  &
   ,patch_area,leaf_class,soil_text)

use rconstants

implicit none

integer :: n2,n3,nzg,npat,ifm,i,j,k,ipat,flag_on
real, dimension(nzg,n2,n3,npat) :: soil_text
real, dimension(n2,n3,npat) :: patch_area,leaf_class

!  This routine is the intended location for a user to customize the
!  PATCH_AREA, leaf_class, and SOIL_TEXT arrays.  It is called after all 
!  other types of initialization of these fields, so this routine has 
!  the last word.  By default the routine makes no change to the
!  fields.  The commented lines below serve as a template for user-designed
!  changes.   Note that this routine is called for each grid separately, so
!  attention to the current value of ngrid in this routine may be required
!  by the user.

flag_on=0
if(flag_on==1)then

 if (ifm .eq. 1) then
  do j = 1,n3
    do i = 1,n2

       patch_area(i,j,1) = 0.        ! patch 1
       leaf_class(i,j,1) = 0.        ! patch 1

       patch_area(i,j,2) = 0.        ! patch 2
       leaf_class(i,j,2) = 0.        ! patch 2

       do k = 1,nzg
          soil_text(k,i,j,1) = 0.    ! patch 1
          soil_text(k,i,j,2) = 0.    ! patch 2
       enddo

    enddo
  enddo

  do ipat = 3,npat
    do j = 1,n3
       do i = 1,n2

          patch_area(i,j,ipat) = 0.
          leaf_class(i,j,ipat) = 0.

          do k = 1,nzg
             soil_text(k,i,j,ipat) = 0.
          enddo

       enddo
    enddo
  enddo

 endif

endif

return
END SUBROUTINE sfcinit_file_user

!##############################################################################
Subroutine sfcinit_nofile_user (n1,n2,n3,mzg,mzs,npat,ifm  &
   ,theta,pi0,pp,rv,seatp,seatf                       &
   ,soil_water     ,soil_energy      ,soil_text       &
   ,sfcwater_mass  ,sfcwater_energy  ,sfcwater_depth  &
   ,veg_fracarea   ,veg_lai          ,veg_tai         &
   ,veg_rough      ,veg_height       ,veg_albedo      &
   ,patch_area     ,patch_rough      ,leaf_class      &
   ,soil_rough     ,sfcwater_nlev    ,stom_resist     &
   ,ground_rsat    ,ground_rvap      ,veg_water       &
   ,veg_temp       ,can_rvap         ,can_temp        & 
   ,veg_ndvip      ,veg_ndvic        ,veg_ndvif       &
   ,snow_mass      ,snow_depth       ,soil_moist_top  & 
   ,soil_moist_bot ,soil_temp_top    ,soil_temp_bot   &
   ,glat           ,glon             ,zot)

use mem_grid
use mem_leaf
use leaf_coms
use io_params
use rconstants
use mem_varinit

implicit none

logical :: there
character(len=7) :: cgrid
character(len=strl1) :: flnm
integer :: n1,n2,n3,mzg,mzs,npat,ifm,i,j,k,ipat &
 ,nveg,nsoil,initsurf,nc,runsoilingest,runsnowingest

real :: c1,airtemp,prsv,piv
real :: tsoil,fice !Saleeby: add these for frozen soil computation
real, dimension(n1,n2,n3) :: theta,pi0,pp,rv
real, dimension(n2,n3)    :: glat,glon,zot  &
                            ,seatp,seatf,snow_mass,snow_depth &
                            ,soil_moist_top,soil_moist_bot &
                            ,soil_temp_top,soil_temp_bot

real, dimension(mzg,n2,n3,npat) :: soil_water,soil_energy,soil_text
real, dimension(mzs,n2,n3,npat) :: sfcwater_mass,sfcwater_energy  &
                                  ,sfcwater_depth

real, dimension(n2,n3,npat) :: veg_fracarea ,veg_lai       ,veg_tai      &
                              ,veg_rough                                 &
                              ,veg_height   ,veg_albedo    ,patch_area   &
                              ,patch_rough  ,leaf_class   &
                              ,soil_rough   ,sfcwater_nlev ,stom_resist  &
                              ,ground_rsat  ,ground_rvap   ,veg_water    &
                              ,veg_temp     ,can_rvap      ,can_temp     &
                              ,veg_ndvip    ,veg_ndvic     ,veg_ndvif

!(Saleeby 4-8-2013): Variables for soil moisture ingest
integer :: ii,jj,minii,minjj,ksn
real :: disty,distx,dist,mindist,maxsoil
real, dimension(n2,n3) :: soil_moist_top_temp,soil_moist_bot_temp &
  ,soil_temp_top_temp,soil_temp_bot_temp

integer :: flag_on

! This routine fills the primary LEAF3 arrays for which standard RAMS
! data files do not exist with initial default values.  Many of the 
! initial values are horizontally homogeneous, although some depend on
! atmospheric conditions.  The default values assigned here may be 
! overridden by (1) specification from coarser grids or (2) specifying new 
! values in routine sfcinit_nofile_user in the file ruser.f.

!If this flag is turned on (==1) then proceed to override sfcinit_nofile
flag_on=0
if(flag_on==1)then 

c1 = .5 * cpi

! Time interpolation factor for updating SST

if (iupdsst == 0) then
   timefac_sst = 0.
else
   timefac_sst = (time - ssttime1(ifm)) / (ssttime2(ifm) - ssttime1(ifm))
endif

!--------------------------------------------------------------------------------
!Saleeby (4-8-2013): Ingest and check soil moisture, soil temperature, snow depth
!and snow water from varfiles. Also fill in water/land points when a disconnect
!exists between input dataset and RAMS grid. Soil moisture ingest will only
!work at initialization and not for history restarts, even if adding a grid.
!--------------------------------------------------------------------------------

!Set soil/snow ingest flags to zero
initsurf=0
runsoilingest=0
runsnowingest=0

!Check for initial start for soil/snow ingest
if(trim(runtype) == 'INITIAL' .or. trim(runtype) == 'ERROR') initsurf=1

!Check to see which varfiles are present since these contain the soil data
if(initial==2)then
 write(cgrid,'(a2,i1,a3)') '-g',ifm,'.h5'
 nc=len_trim(fnames_varf(nvarffl))
 flnm=fnames_varf(nvarffl)(1:nc-4)//trim(cgrid)
 inquire(file=trim(flnm),exist=there)
endif

!Set flag for soil ingest
if(isoildat==1 .and. initsurf==1 .and. initial==2 .and. there) runsoilingest=1
!Set flag for snow ingest
if(isnowdat==1 .and. initsurf==1 .and. initial==2 .and. there) runsnowingest=1

if(runsoilingest==1)then

 if(iprntstmt>=1 .and. print_msg)then
  print*,'-------------------------------------------------------'
  print*,'SETTING UP SOIL AND SNOW LEVEL VARIABLES ON GRID:',IFM
  print*,'-------------------------------------------------------'
 endif

 !Do some initial quality control
 do j = 1,n3
  do i = 1,n2
   if(soil_moist_top(i,j) > 0.49 .or. soil_moist_bot(i,j) > 0.49) then
    soil_moist_top(i,j) = .49
    soil_moist_bot(i,j) = .49
   endif
   if((soil_moist_top(i,j) > 0.0 .and. soil_moist_top(i,j) < 0.01) .or. & 
      (soil_moist_bot(i,j) > 0.0 .and. soil_moist_bot(i,j) < 0.01)) then
    soil_moist_top(i,j) = 0.0101
    soil_moist_bot(i,j) = 0.0101
   endif
   soil_moist_top_temp(i,j) = -9999.
   soil_moist_bot_temp(i,j) = -9999.
   soil_temp_top_temp(i,j) = -9999.
   soil_temp_bot_temp(i,j) = -9999.
  enddo
 enddo

 !Interpolate to what amounts to coastal point where data is null, but
 !RAMS says there is land that needs valid soil moisture values.
 do j = 1,n3
  do i = 1,n2
   if((soil_moist_top(i,j) < 0.01 .or. soil_moist_top(i,j) > 0.49 .or. &
       soil_moist_bot(i,j) < 0.01 .or. soil_moist_bot(i,j) > 0.49 ) .and. &
       patch_area(i,j,2) > 0.001 .and. & 
       nint(leaf_class(i,j,2)) >= 2 .and. nint(leaf_class(i,j,2)) <= 21) then
    mindist=9.99e12
    minii=i
    minjj=j
    do jj = 1,n3
     do ii = 1,n2
      if(soil_moist_top(ii,jj) > 0.01 .and. &
         nint(leaf_class(ii,jj,2)) >= 2 .and. nint(leaf_class(ii,jj,2)) <= 21) then
       disty=abs(glat(i,j)-glat(ii,jj))*111.12
       distx=abs(glon(i,j)-glon(ii,jj))*111.12*cos((glat(i,j)+glat(ii,jj))/2.0*3.14/180.)
       dist=sqrt(distx**2.0+disty**2.0)
       if(dist<mindist)then
        mindist=dist
        minii=ii
        minjj=jj
       endif
      endif
     enddo
    enddo
    ii=minii
    jj=minjj
    soil_moist_top_temp(i,j) = soil_moist_top(ii,jj)
    soil_moist_bot_temp(i,j) = soil_moist_bot(ii,jj)
    soil_temp_top_temp(i,j)  = soil_temp_top(ii,jj)
    soil_temp_bot_temp(i,j)  = soil_temp_bot(ii,jj)
   endif
  enddo
 enddo

 !Quality check our interpolation and look for bad points
 !This could indicate that varfile was not created.
 do j = 1,n3
  do i = 1,n2
   maxsoil=max(soil_moist_top_temp(i,j),soil_moist_top(i,j))
   if(patch_area(i,j,2) > 0.001 .and. (maxsoil >.50 .or. maxsoil <.01)) then
    print*,'WARNING: Some land areas have invalid soil moisture values:'
    print*,'Need varfiles for all grids at initial time for valid soil on all grids.'
    print*,ifm,i,j,soil_moist_top_temp(i,j),soil_moist_top(i,j)
    stop
   endif
  enddo
 enddo

endif !end if soil moisture ingest from varfile
!--------------------------------------------------------------------------------

!Initialize surface properties
do j = 1,n3
   do i = 1,n2
      piv = c1 * (pi0(1,i,j) + pi0(2,i,j)   &
                      + pp(1,i,j) + pp(2,i,j))
      airtemp = theta(2,i,j) * piv
      prsv = piv ** cpor * p00

      patch_rough(i,j,1) = 0.001
      can_temp(i,j,1) = airtemp
      can_rvap(i,j,1) = rv(2,i,j)

      soil_energy(mzg,i,j,1) = 334000.  &
         + 4186. * (seatp(i,j) + (seatf(i,j) - seatp(i,j))  &
         * timefac_sst - 273.15)

      !Loop over patches
      do ipat = 2,npat

         nveg = nint(leaf_class(i,j,ipat))

         soil_rough(i,j,ipat) = zrough
         patch_rough(i,j,ipat) = max(zrough,zot(i,j))
         veg_rough(i,j,ipat) = .13 * veg_ht(nveg)

         veg_height(i,j,ipat) = veg_ht(nveg)

         veg_albedo(i,j,ipat) = albv_green(nveg)
         stom_resist(i,j,ipat) = 1.e6

         veg_temp(i,j,ipat) = airtemp
         can_temp(i,j,ipat) = airtemp

         veg_water(i,j,ipat) = 0.
         can_rvap(i,j,ipat) = rv(2,i,j)

         !Loop over soil layers
         do k = 1,mzg

            nsoil = nint(soil_text(k,i,j,ipat))
            soil_water(k,i,j,ipat) = max(soilcp(nsoil),slmstr(k) * slmsts(nsoil))

            !Saleeby(4-9-2013):If using variable soil initialization from ingest data
            if(runsoilingest==1)then
             if(k==mzg) then
              soil_water(k,i,j,ipat) = max(soilcp(nsoil),soil_moist_top_temp(i,j) &
                                          ,soil_moist_top(i,j))
             else
              soil_water(k,i,j,ipat) = max(soilcp(nsoil),soil_moist_bot_temp(i,j) &
                                          ,soil_moist_bot(i,j))
             endif
            endif

            !For persistent wetlands (bogs, marshes, fens, swamps) and irrigated
            !crops, initialize with saturated soil. Currently, this corresponds to
            !leaf classes 16, 17, and 20.
            if (nint(leaf_class(i,j,ipat)) == 16 .or.  &
                nint(leaf_class(i,j,ipat)) == 17 .or.  &
                nint(leaf_class(i,j,ipat)) == 20) then
               soil_water(k,i,j,ipat) = slmsts(nsoil)
            endif

            !Limit soil water fraction maximum for soil type. Need to do this
            !since input soil water from gridded data can differ from RAMS due
            !to differences in soil classes and soil porosity.
            if(isfcl==2) &
             soil_water(k,i,j,ipat) = min(slmsts(nsoil),soil_water(k,i,j,ipat))

            !By default, initialize soil internal energy at a temperature equal
            !to airtemp + stgoff(k) and with all water assumed to be liquid.  If
            !the temperature is initially below 0C, this will immediately adjust
            !to soil at 0C with part ice.  In order to begin with partially or
            !totally frozen soil, reduce or remove the latent-heat-of-fusion
            !term (the one with the factor of 3.34) from soil_energy below. If
            !the soil is totally frozen and the temperature is below zero C, the
            !factor of 4.186 should be changed to 2.093 to reflect the reduced
            !heat capacity of ice compared to liquid. These changes may be
            !alternatively be done in routine sfcinit_user in ruser.f

            !Saleeby: Old Method sometimes allows surface runaway cooling!
            !soil_energy(k,i,j,ipat) = (airtemp - 273.15 + stgoff(k))  &
            !   * (slcpd(nsoil) + soil_water(k,i,j,ipat) * 4.186e6)  &
            !   + soil_water(k,i,j,ipat) * 3.34e8

            !*************************************************************
            !Saleeby: New Method allows for tsoil < 0C and set initial soil
            ! ice fraction (fice) to zero
            tsoil = (airtemp - 273.15 + stgoff(k))

            !Saleeby(4-9-2013):If we are ingesting soil temperature from varfiles
            if(runsoilingest==1)then
             if(k==mzg) then
              tsoil=max(tsoil,soil_temp_top_temp(i,j)-273.15,soil_temp_top(i,j)-273.15)
             else
              tsoil=max(tsoil,soil_temp_bot_temp(i,j)-273.15,soil_temp_bot(i,j)-273.15)
             endif
            endif

            fice  = 0.0
            !For soil temperature > 0C
            if (tsoil.gt.0.) then
              soil_energy(k,i,j,ipat) = tsoil  &
                 * (slcpd(nsoil) + soil_water(k,i,j,ipat) * 4.186e6)  &
                 + soil_water(k,i,j,ipat) * 3.34e8
            !For soil temperature <= 0C
            elseif (tsoil.le.0.) then
              !This soil ice fraction should be improved
              !Just a linear computation with fice=0.0 at 273K
              ! up to fice=1.0 at 272K
              fice = min(1.0,max(0.0, -1.0*tsoil*0.1 ))
              soil_energy(k,i,j,ipat) = tsoil  &
                 * (slcpd(nsoil) + &
                    (1.-fice)*soil_water(k,i,j,ipat) * 4.186e6 +  &
                         fice*soil_water(k,i,j,ipat) * 2.093e6)  &
                 + (1.-fice)*soil_water(k,i,j,ipat) * 3.34e8
            endif
            !*************************************************************

         enddo !end looping over soil layers

         !Determine number of surface water or snow levels and set values of
         !mass, depth, and energy for each level present.

         !Set number of inital surface water layers to zero
         sfcwater_nlev(i,j,ipat) = 0.

         !Loop over surface water / snow layers
         do k = 1,mzs

            sfcwater_mass(k,i,j,ipat) = 0.
            sfcwater_energy(k,i,j,ipat) = 0.
            sfcwater_depth(k,i,j,ipat) = 0.

            !For persistent wetlands (bogs, marshes, fens, swamps), initialize
            !with 10 cm water depth. This corresponds to leaf classes 17 and 20.
            if (nint(leaf_class(i,j,ipat)) == 17 .or.  &
                nint(leaf_class(i,j,ipat)) == 20) then
               if (k .eq. 1) then
                  sfcwater_mass(k,i,j,ipat) = 100.
                  sfcwater_energy(k,i,j,ipat) = (airtemp - 193.36) * 4186.
                  sfcwater_depth(k,i,j,ipat) = .1
               endif
            endif

            !Initialize surface water layer quantities
            if (runsnowingest==1) then
             if (snow_mass(i,j) > 0.) then
               if (k .eq. 1) then
                  sfcwater_mass(k,i,j,ipat)   = sfcwater_mass(k,i,j,ipat)   &
                     + snow_mass(i,j)
                  sfcwater_energy(k,i,j,ipat) = sfcwater_energy(k,i,j,ipat) &
                     + min(0., (airtemp - 273.15) * 2093.)
                  if(runsoilingest==1) then 
                    sfcwater_depth(k,i,j,ipat) = sfcwater_depth(k,i,j,ipat) & 
                        + snow_depth(i,j) ! From Varfiles
                  else
                    sfcwater_depth(k,i,j,ipat) = sfcwater_depth(k,i,j,ipat)  &
                        + snow_mass(i,j) * 5.e-3   ! 5x equivalent liquid depth
                  endif
               endif
             endif
            endif

            !Determine number of surface water layers
            if (sfcwater_mass(k,i,j,ipat) > 0.) sfcwater_nlev(i,j,ipat) = float(k)

         enddo !end loop over surface water layers

         !Call for all land surface options (isfcl)
         if (ipat >= 2) CALL ndvi (ifm  &
            ,leaf_g(ifm)%leaf_class (i,j,ipat)   &
            ,leaf_g(ifm)%veg_ndvip  (i,j,ipat)   &
            ,leaf_g(ifm)%veg_ndvic  (i,j,ipat)   &
            ,leaf_g(ifm)%veg_ndvif  (i,j,ipat)   )

         !Only call this for LEAF3. SIB will use its own values.
         if (ipat >= 2 .and. isfcl<=1)           &
            CALL veg (ifm                        &
            ,leaf_g(ifm)%leaf_class   (i,j,ipat) &
            ,leaf_g(ifm)%veg_fracarea (i,j,ipat) &
            ,leaf_g(ifm)%veg_lai      (i,j,ipat) &
            ,leaf_g(ifm)%veg_tai      (i,j,ipat) &
            ,leaf_g(ifm)%veg_rough    (i,j,ipat) &
            ,leaf_g(ifm)%veg_height   (i,j,ipat) &
            ,leaf_g(ifm)%veg_albedo   (i,j,ipat) &
            ,leaf_g(ifm)%veg_ndvic    (i,j,ipat) )

         ksn = nint(sfcwater_nlev(i,j,ipat))
         CALL grndvap (  &
             soil_energy(mzg,i,j,ipat),soil_water(mzg,i,j,ipat)  &
            ,soil_text  (mzg,i,j,ipat),sfcwater_energy(ksn,i,j,ipat)  &
            ,sfcwater_nlev(i,j,ipat),ground_rsat(i,j,ipat) &
            ,ground_rvap(i,j,ipat),can_rvap(i,j,ipat),prsv,i,j)

      enddo !end loop over patches
   enddo
enddo

endif

return
END SUBROUTINE sfcinit_nofile_user

!##############################################################################
Subroutine bubble (m1,m2,m3,thp,rtp)

use micphys
use mem_grid

implicit none

integer :: m1,m2,m3,i,j,k
real, dimension(m1,m2,m3) :: thp,rtp
real R,rand_num

integer bubtemp
real acetmp1,acetmp2,acetmp3,acetmp4,acetmp5
real acetmp6,acetmp7,acetmp8,acetmp9
real bubctrx,bubctry,bubctrz
real bubradx,bubrady,bubradz
real, dimension(:,:), allocatable :: bub_rand_nums

if(ibubble==1) then
 print*,''
 print*,'INITIALIZING Square RAMSIN WARM BUBBLE HERE'
 print*,'On grid number=',IBUBGRD
 print*,'Bubble from I=',IBDXIA,'TO',IBDXIZ
 print*,'Bubble from J=',IBDYJA,'TO',IBDYJZ
 print*,'Bubble from K=',IBDZK1,'TO',IBDZK2
 print*,'THP perturbation=',BTHP
 print*,'RTP perturbation=',BRTP
 print*,''
 do j = IBDYJA,IBDYJZ
   do i = IBDXIA,IBDXIZ
     do k = IBDZK1,IBDZK2
            thp(k,i,j) = thp(k,i,j) + BTHP
            rtp(k,i,j) = rtp(k,i,j) * (1.0 + BRTP)
     enddo
   enddo
 enddo

elseif(ibubble==2) then
 print*,'INITIALIZING Cosine-Squared RAMSIN WARM BUBBLE HERE'
 print*,'On grid number=',IBUBGRD
 print*,'Use of I,J coordinates for Cosine-Squared bubble is'
 print*,'only accurate and valid if the domain POLE POINT'
 print*,'is the same at the CENTRAL LAT and LON'
 !Set up X location of bubble center relative to grid center
 bubctrx = deltax * ( (IBDXIA+IBDXIZ)/2.0 - NNXP(1)/2.0 )
 !Set up Y location of bubble center relative to grid center
 bubctry = deltax * ( (IBDYJA+IBDYJZ)/2.0 - NNYP(1)/2.0 )
 !Set up Z location of bubble center
 print*,'IBDXIA,IBDXIZ',IBDXIA,IBDXIZ
 print*,'IBDYJA,IBDYJZ',IBDYJA,IBDYJZ
 print*,'X-center,Y-center:',bubctrx,bubctry
 bubtemp=int(IBDZK1+IBDZK2)/2.0
 bubctrz=ZMN(bubtemp,1)
 !Set up length, width, and depth of bubble
 bubradx=(IBDXIZ-IBDXIA) * deltax * 0.5
 bubrady=(IBDYJZ-IBDYJA) * deltax * 0.5
 bubradz=(ZMN(IBDZK2,1)-ZMN(IBDZK1,1)) * 0.5
 !Set up gaussian bubble
 acetmp8=atan(1.0)*4.0/2.0 ! pi/2
 do j=1,m3
  do i=1,m2
   do k=1,m1
    acetmp1=(XMN(i,1)+XMN(i+1,1))*0.5
    acetmp2=(YMN(j,1)+YMN(j+1,1))*0.5
    acetmp3=(ZMN(k,1)+ZMN(k+1,1))*0.5
    acetmp4=(acetmp1-bubctrx)/bubradx
    acetmp4=acetmp4**2
    acetmp5=(acetmp2-bubctry)/bubrady
    acetmp5=acetmp5**2
    acetmp6=(acetmp3-bubctrz)/bubradz
    acetmp6=acetmp6**2
    if(jdim==0)acetmp5=acetmp4 !For 2D simulation let Y-dir = X-dir
    acetmp7=acetmp4+acetmp5+acetmp6
    acetmp7=SQRT(acetmp7)
    if(acetmp7.ge.1.0)then
     acetmp9=0.0
    else
     acetmp9=(COS(acetmp8*acetmp7))**2
    endif
    thp(k,i,j)=thp(k,i,j)+(BTHP*acetmp9) !pt pert
    rtp(k,i,j)=rtp(k,i,j)*(1.0+(acetmp9*BRTP)) !mix ratio pert
   enddo
  enddo
 enddo

elseif(ibubble==3) then
 if(print_msg) then
  print*,'Activating random temperature perturbation with'
  print*,'amplitude 0.1K at z=2, linearly-decreasing to 0K over 500m'
  print*,'On grid number=',IBUBGRD
  print*,''
 endif

 do k=2,m1
   ! select levels for temp. pert. based on altitude
   if( zt(k) <= (500.+zt(2)) ) then ! only over lowest 500 m
     ! At this point bub_rand_nums covers the entire domain and is filled with
     ! random numbers (between -1 and 1) that were generated by only one
     ! process. Go through the subdomain and copy the corresponding section
     ! of bub_rand_nums when reading the random numbers.
     do j = 1, m3
       do i = 1, m2
         CALL random_number (rand_num)  ! random number between 0 and 1
         R = 1.0 - (2.0 * rand_num)     ! random number between -1 and 1
         ! Pert with amplitude of 0.1K, largest at level 2, decrease linearly over 500m
         thp(k,i,j)=thp(k,i,j) + R*0.1*(500.+zt(2)-zt(k))/500.
         if(k==2) thp(1,i,j)=thp(k,i,j) ! set level 1 to level 2
       enddo ! i
     enddo ! j

   endif ! in lowest 500 m
 enddo ! k

endif

return
END SUBROUTINE bubble

!##############################################################################
Subroutine conv_forcing (ut,vt)

use mem_grid
use mem_basic
use micphys
use node_mod
use mem_tend

implicit none
!   +------------------------------------------------------------------
!   ! This routine creates a convergenze zone for layer lifting
!   !   following Loftus et al. (2008, MWR), except convergence
!   !   zone is implemented as a tendency, similarly to 
!   !   Schumacher (2009, JAS)
!   +------------------------------------------------------------------

integer :: i,j,k,ii,jj
real :: fz,camp,uprime,vprime,wmax,bot,top
real, dimension(mzp,mxp,myp) :: ut,vt

!check for termination type
if (ctmax.lt.0.0) then
  if(vertvel(ngrid) >= abs(ctmax)) iconv=0 !Turn off convergence
  if(iconv > 0 .and. print_msg) &
    print*,'Run Convergence Forcing, ICONV=',iconv,' Stop at W >=',abs(ctmax)
else
  if(time >= ctmax) iconv=0 !Turn off convergence
  if(iconv > 0 .and. print_msg) &
    print*,'Run Convergence Forcing, ICONV=',iconv,' Stop at Time=',ctmax
endif

if(iconv == 0) return

!calculate 1/2 wavelength boundaries if ICVERT=2
if (icvert.eq.2) then
  bot=zt(ickcent)-czrad ! m
  top=czrad+zt(ickcent) ! m
endif

! calculate camp from specified namelist settings and divergence
! (Loftus et al. 2008 eqn 7)  units: m^2/s
camp = -0.5*cdivmax / ( 1./(cxrad**2.) + 1./(cyrad**2.) )

do j=1,myp
  do i=1,mxp
    do k=2,mzp

      !create vertical decrease with height dependence - max at k=2
      if (icvert.eq.1 .and. k.le.ickmax) then
        fz=1.-(2.-k)/(2.-ickmax)
      endif
      !create 1/2 vertical wavelength elevated convergence (Schumacher 2009)
      if (icvert.eq.2 .and. (zt(k).ge.bot .and. zt(k).le.top)) then
        fz=cos(0.5*3.14159*abs((zt(k)-zt(ickcent))/czrad))**2.
      elseif (icvert.eq.2) then
        fz=0.
      endif

      !Get absolute grid points for parallel (& sequential) computation
      ii = i+mi0(ngrid)
      jj = j+mj0(ngrid)

      !create u' and v'. Units: m/s. From Loftus et al. eqns. 5 & 6
      if (iconv.eq.1) then !Gaussian in x and y, U and V
        uprime = (-2.*camp*float(ii-icicent)*deltax/cxrad**2.)*&
          exp( -1.*(float(ii-icicent)*deltax/cxrad)**2. )*&
          exp( -1.*(float(jj-icjcent)*deltax/cyrad)**2. )*fz
        vprime = (-2.*camp*float(jj-icjcent)*deltax/cyrad**2.)*&
          exp( -1.*(float(ii-icicent)*deltax/cxrad)**2. )*&
          exp( -1.*(float(jj-icjcent)*deltax/cyrad)**2. )*fz
      elseif (iconv.eq.2) then !Gaussian in x and y, U only
        uprime = (-2.*camp*float(ii-icicent)*deltax/cxrad**2.)*&
          exp( -1.*(float(ii-icicent)*deltax/cxrad)**2. )*&
          exp( -1.*(float(jj-icjcent)*deltax/cyrad)**2. )*fz
        vprime = 0.
      elseif (iconv.eq.3) then !Gaussian in x and y, V only
        uprime = 0.
        vprime = (-2.*camp*float(jj-icjcent)*deltax/cyrad**2.)*&
          exp( -1.*(float(ii-icicent)*deltax/cxrad)**2. )*&
          exp( -1.*(float(jj-icjcent)*deltax/cyrad)**2. )*fz
      elseif (iconv.eq.4) then !Gaussian in x, U only
        uprime = (-2.*camp*float(ii-icicent)*deltax/cxrad**2.)*&
          exp( -1.*(float(ii-icicent)*deltax/cxrad)**2. )*fz
        vprime = 0.
      elseif (iconv.eq.5) then !Gaussian in y, V only
        uprime = 0.
        vprime = (-2.*camp*float(jj-icjcent)*deltax/cyrad**2.)*&
          exp( -1.*(float(jj-icjcent)*deltax/cyrad)**2. )*fz
      endif ! iconv check

      ! add tendency over timescale specified in namelist
      ut(k,i,j) = ut(k,i,j) + uprime/ctau
      vt(k,i,j) = vt(k,i,j) + vprime/ctau

    enddo ! k = 2:mzp
  enddo ! i loop
enddo ! j loop

return
END SUBROUTINE conv_forcing
