!------------------------------------------------------------------------
!************************************************************************
!------------------------------------------------------------------------
!
!  EVV: radial variation
!
!  ndays = 300 DEF run
!  double precision, parameter :: rain_event_min = 25.0
!
!  ndays = 100
!  double precision, parameter :: rain_event_min = 40.0
!
!  double precision, parameter :: fevap_low = 0.6
!  double precision, parameter :: fevap_add = 0.0
!
!  double precision, parameter :: fevap_add = 0.7
!
!  ndays = 300 DEF run
!   to see what div.pdf looks like (likely same as 50 days)
!
!  ndays = 100  Directory: DD=0.5
!    weird clustering but no MJO
!    shows need higher rates of downdraft formation at higher dp mass fluxes
!    double precision, parameter :: fevap_low = 0.5
!    double precision, parameter :: fevap_add = 0.0
!
!  ndays = 50 
!  double precision, parameter :: colrh_dp_half = 0.75   Directory: COLRH=0.75
!
!   ndays = 100  dist_lt = 650.*1000. SLOWS DOWN MJO
!   AVERAGE MJO speed =    4.7360420729234125
!
!  You can't do a simple change in the weight function to change the MJO speed, since if do
!  at the surface for both east and west directions, will shift both up and down motions, and
!  likely no net effect. Likely have to shift the eastward direction only to affect the up
!  motion only
!
!  Order of loops in spread_mass makes a difference (Sep 30 2025)
!  --------------------------------------------------------------
!  - These are the AA,BB,YY loops. The structure of the hov rain pattern is
!    the same, but the specific pattern depends on loop order. The calculation of the
!    weight function for all individual ii,jj values is not affected by the order,
!    but when these values are summed into weight_total the sum does slightly
!    depend on the loop order, and this affects rainfall variance after a few weeks.
!    It appears likely that the dependence of a sum on the order of addition occurs
!    when you are adding numbers of very different magnitude. It is likely better to
!    add the smallest first, or use quadruple precision. The origin of the bug is simply
!    that a computer can not add numbers of order 1 to 1.0E-12 in a different order and
!    get the same result every time. Likely know fully accurate way of doing at double
!    precision.
!
!  ****************************************************************************
!  *                                                                          *
!  *   -----  Changing spatial and temporal resolution  -----                 *
!  *   -changing temporal resolution doesn't change much before convection    *
!  *   starts. However, a bit unclear about what to do with dp_inc and dp_dec *
!  *   when change timestep.                                                  *
!  *   - changing from 2 deg to 1.5 degrees speeds things up a lot. It does   *
!  *     change the behaviour. Rain variability is similar but not the same.  *
!  *                                                                          *
!  ****************************************************************************
!
!
!   Processes:
!   ----------
!   CG: cumulus congestus
!   - cg prod and loss. prod depends on colrh,caplt,dyn production
!   CGE: CG evaporation
!   DP: Deep Convection
!    - depends on caput,colrh, prod and loss terms
!   AN: anvil rain
!    - when upper trop RH goes above a value
!   DD: downdrafts
!    - now generated by rain_an and rain_dp
!   DPE: Deep evaporation
!   DN: Dynamically Forced Ascent/Descent
!   SD: Radiative Subsidence SD
!   - Radiative Cooling: one value for each layer
!   SP: horizontal motion
!
!------------------------------------------------------------------------
module model
!------------------------------------------------------------------------
!
!   Use
!
!------------------------------------------------------------------------
use help
use netcdf
implicit none
!------------------------------------------------------------------------
!
!  SP: spread_mass parameters 
!
!------------------------------------------------------------------------
double precision, parameter :: pi = 3.1415
double precision, parameter :: dlat = 30.0
double precision, parameter :: weight_min = 1.0E-05
double precision, parameter :: dmdry_tolerance = 1.0E-04
double precision, parameter :: colwat_tolerance = 5.0E-06
!
integer, parameter :: j_weight_save_upper = 13 ! 12S
integer, parameter :: i_weight_save_upper = 80
integer, parameter :: j_weight_save_lower = 28 ! 12N
integer, parameter :: i_weight_save_lower = 200
!------------------------------------------------------------------------
!
!  Solving for the divergent circulation
!
!------------------------------------------------------------------------
integer, parameter :: niter_massflux = 50000
!------------------------------------------------------------------------
!
!  Module Integer Parameters
!
!------------------------------------------------------------------------
integer, parameter :: nlon = 240
integer, parameter :: nlat = 41
integer, parameter :: nlev = 3
integer, parameter :: nlon_mjo = 121  ! from ii = 60 t0 180
!------------------------------------------------------------------------
!  General Control Parameters
!  The only change I can see in increasing the timestep from 30 min to 60 min
!  is that drizzle gets noisier, as you would expect. Everything else stays the same.
!------------------------------------------------------------------------
!integer, parameter :: ndays = 1  ! need to reduce ndays_spinup = 25
!integer, parameter :: ndays = 50  ! need to reduce ndays_spinup = 25
!
!integer, parameter :: ndays = 30  ! 30*24 = 840
!integer, parameter :: ndays = 35  ! 35*24 = 840
!integer, parameter :: ndays = 50  ! 1200
integer, parameter :: ndays = 100
!integer, parameter :: ndays = 150
!integer, parameter :: ndays = 200  ! 200*24 = 4800
!integer, parameter :: ndays = 300  ! 300*24 = 7200
!
integer, parameter :: nsteps_day = 24
integer, parameter :: nsteps = ndays*nsteps_day
!------------------------------------------------------------------------
!  Options for processes
!------------------------------------------------------------------------
integer, parameter :: do_trade_wind = 0
integer, parameter :: do_flux = 1
integer, parameter :: do_an = 1
integer, parameter :: do_dp = 1
integer, parameter :: do_dr = 1
integer, parameter :: do_cg = 1
integer, parameter :: do_cg_evap = 1
integer, parameter :: do_dp_evap = 1
integer, parameter :: do_down = 1
integer, parameter :: do_spread = 1
integer, parameter :: do_rad_vert = 1
integer, parameter :: do_dyn_vert = 1
!------------------------------------------------------------------------
!  Pressure Parameters
!------------------------------------------------------------------------
double precision, parameter :: re = 6378000.  ! radius of the earth in m
double precision, parameter :: p1 = 100000. ! surface pressure in Pa
double precision, parameter :: p2 = 90000. ! p at top of BL
double precision, parameter :: p3 = 50000. ! p at top of LT
double precision, parameter :: p4 = 15000. ! tropopause p
double precision, parameter :: z1 = 500.
double precision, parameter :: z2 = 3500.
double precision, parameter :: z3 = 10000.
double precision, parameter :: sec_day = 3600.0000000*24.0000000
!------------------------------------------------------------------------
!   T/RH Initialization
!------------------------------------------------------------------------
double precision, parameter :: t_bl_init = 295.
double precision, parameter :: t_lt_init = 280.
double precision, parameter :: t_ut_init = 235.
double precision, parameter :: rh_bl_init = 0.8
double precision, parameter :: rh_lt_init = 0.6
double precision, parameter :: rh_ut_init = 0.8
!------------------------------------------------------------------------
!  SST parameters
!  - with sst_eq = 28 and sst_dist = 5000 km, sst = 19.6 at 30N.
!  - with sst_eq = 28 and sst_dist = 6000 km, sst = 21.8 at 30N.
!------------------------------------------------------------------------
!double precision, parameter :: sst_eq = 8.0  ! Have to increase for interactive SST
!double precision, parameter :: sst_back = 25.0 ! interactive
double precision, parameter :: sst_eq = 6.0  ! Have to increase for interactive SST
double precision, parameter :: sst_back = 23.
double precision, parameter :: sst_dist = 1500.*1000.
double precision, parameter :: sst_tscale = 5.0*sec_day ! (sst recovery time)
!double precision, parameter :: sst_mass = 4000.0  !(4000 kg = 4m deep) Interactive
double precision, parameter :: sst_mass = 900000000.  ! fixed sst
!------------------------------------------------------------------------
!  Flux Parameters
!------------------------------------------------------------------------
double precision, parameter :: a_evap = 1.0E-10
double precision, parameter :: a_heat = 0.5E-06
!------------------------------------------------------------------------
!
!  t_from_hm control parameters
!  - routine be a time sink if too many loops
!
!------------------------------------------------------------------------
integer, parameter :: niter_hm = 100   ! number of iterations (should not matter)
double precision, parameter :: dtt_start = 20.0
double precision, parameter :: hmerr_want = 1.0E-02 ! max allowed hmerr difference
double precision, parameter :: dtt_reduce = 0.05  ! reduce dttt each iteration
!------------------------------------------------------------------------
!
!  Trade Wind circulation
!  In 30 min, 5 m/s goes 9000 m or 9 km.
!  This routine does not conserve column water.
!
!------------------------------------------------------------------------
double precision, parameter :: u_trade_wind = 0.0
double precision, parameter :: dist_tw = 1000.*1000.
!------------------------------------------------------------------------
!
!  dpdx varies by about 2000/(100*100*1000) = 2/(100*100) = 2/10000 = 2E-04
!
!------------------------------------------------------------------------
double precision, parameter :: omega = 7.292E-05  ! inverse seconds
double precision, parameter :: frict = 5.0E-05 ! friction (inverse seconds) 
!------------------------------------------------------------------------
!
!  Module Integer Parameters: control diagnostics
!
!------------------------------------------------------------------------
integer, parameter :: nprint_sigma = 0
integer, parameter :: nprint_vert = 0
integer, parameter :: nprint_dp = 0
integer, parameter :: nprint_dn = 0
integer, parameter :: nprint_dyn = 0
integer, parameter :: nprint_pg = 0 !  printing for pga 
integer, parameter :: nprint_hmerr = 0
integer, parameter :: i_want = 1
integer, parameter :: j_want = 11
integer, parameter :: i_one = 114
integer, parameter :: j_one = 21
!------------------------------------------------------------------------
!
!  Module Real Variables Arrays
!
!------------------------------------------------------------------------
double precision :: pp(nlev),ph(nlev+1),pph(nlon,nlat,nlev+1),ppp(nlon,nlat,nlev)
double precision :: tt(nlon,nlat,nlev),rv(nlon,nlat,nlev),hm(nlon,nlat,nlev),zz(nlon,nlat,nlev)
double precision :: zz_mn(nlon,nlat,nlev),senflux_mn(nlon,nlat),watflux_mn(nlon,nlat)
double precision :: pt(nlon,nlat,nlev),zh(nlon,nlat,nlev+1)
double precision :: sigma1(nlon,nlat),sigma2(nlon,nlat)
double precision :: omega1(nlon,nlat),omega2(nlon,nlat)
double precision :: mextra(nlon,nlat,nlev)
double precision :: rh(nlon,nlat,nlev)
double precision :: mass_ut_lt(nlon,nlat)  ! Units are kg/(m2*tstep)
double precision :: mass_lt_bl(nlon,nlat)
double precision :: mfu(nlon_mjo,nlat,nlev)
double precision :: mfv(nlon_mjo,nlat,nlev)
double precision :: mconv(nlon,nlat,nlev)
double precision :: mconv_mn(nlon,nlat,nlev)
double precision :: sub_ut_lt(nlon,nlat)
double precision :: sub_lt_bl(nlon,nlat)
double precision :: speed_mjo(nsteps),speed_mjo_av
!------------------------------------------------------------------------
!  prev variables
!------------------------------------------------------------------------
double precision :: rv_prev(nlon,nlat,nlev),mass_ut_lt_prev(nlon,nlat),mass_lt_bl_prev(nlon,nlat)
double precision :: sub_ut_lt_prev(nlon,nlat),sub_lt_bl_prev(nlon,nlat),rain_tt_prev(nlon,nlat)
double precision :: colrh_prev(nlon,nlat),wat_bl_prev(nlon,nlat)
!------------------------------------------------------------------------
! mn
!------------------------------------------------------------------------
double precision :: sigma1_mn(nlon,nlat)
double precision :: sigma2_mn(nlon,nlat)
double precision :: sub_ut_lt_mn(nlon,nlat)
double precision :: sub_lt_bl_mn(nlon,nlat)
double precision :: cape_lt_mn(nlon,nlat)
double precision :: cape_ut_mn(nlon,nlat)
double precision :: f_dp_mn(nlon,nlat)
double precision :: f_cg_mn(nlon,nlat)
double precision :: f_dp(nlon,nlat)
double precision :: f_cg(nlon,nlat)
double precision :: f_caput_mn(nlon,nlat)
double precision :: f_caplt_mn(nlon,nlat)
double precision :: f_caput(nlon,nlat)
double precision :: f_caplt(nlon,nlat)
double precision :: f_colrh_dp_mn(nlon,nlat)
double precision :: f_colrh_cg_mn(nlon,nlat)
double precision :: f_colrh_dp(nlon,nlat)
double precision :: f_colrh_cg(nlon,nlat)
double precision :: mass_dp_mn(nlon,nlat)
double precision :: mass_cg_mn(nlon,nlat)
double precision :: wconv_bl_mn(nlon,nlat)
double precision :: wdelta_bl_mn(nlon,nlat)
double precision :: sst_mn(nlon,nlat)
double precision :: lon(nlon),lat(nlat)
double precision :: lon_dyn(nlon),lat_dyn(nlat)
double precision :: weight_save(nlon,nlat)
double precision :: weight_save_upper(nlon,nlat)
double precision :: weight_save_lower(nlon,nlat)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
double precision :: rain_an(nlon,nlat)
double precision :: rain_dr(nlon,nlat)
double precision :: rain_dp(nlon,nlat)
double precision :: rain_tt(nlon,nlat)  ! total rain
double precision :: rain_cg(nlon,nlat)
double precision :: evap_cg(nlon,nlat)
double precision :: evap_dp(nlon,nlat)
double precision :: down_dp(nlon,nlat)
double precision :: mass_dp(nlon,nlat)
double precision :: dp_prod(nlon,nlat)
double precision :: mass_dp_old(nlon,nlat)
double precision :: mass_cg(nlon,nlat)
double precision :: mass_dn(nlon,nlat)
double precision :: sst_tar(nlon,nlat)
double precision :: sst(nlon,nlat)
double precision :: cape_lt(nlon,nlat)
double precision :: cape_ut(nlon,nlat)
double precision :: colrh(nlon,nlat)
double precision :: wat_bl(nlon,nlat)
double precision :: wdelta_bl(nlon,nlat)
double precision :: wconv_bl(nlon,nlat)
double precision :: rvmass_bl(nlon,nlat)
!------------------------------------------------------------------------
!
!  Module
!
!------------------------------------------------------------------------
double precision :: f_caput_av,f_caplt_av,colwat_error,dm_extra,dm_final,dm_want,dmdry,qq1,qq2
double precision :: f_dp_av,f_cg_av,f_colrh_dp_av,f_colrh_cg_av,rain_an_av,rain_dp_av,rain_dr_av
double precision :: rain_cg_av,evap_cg_av,drymass_error,sst_av,tstep
double precision :: sub_ut_lt_av,sub_lt_bl_av,rad_cool,f_sub,evap_dp_av,down_dp_av
integer :: nstep,num_speed_mjo_av,num_speed_mjo
!------------------------------------------------------------------------
!
!   *******************************************************************
!   ********************* Horizontal Spreading  ***********************
!   *******************************************************************
!
!  value of zonal_amp affects direct of propagation but not speed
!  dist_lt = 400: MJO gets faster and pattern partially breaks down
!
!------------------------------------------------------------------------
!integer, parameter :: dist_ut = 1000.*1000.  ! DEF
integer, parameter :: dist_ut = 1100.*1000.  ! MJO but no apparent change in speed?
!
integer, parameter :: dist_lt = 500.*1000.  ! DEF
!integer, parameter :: dist_lt = 650.*1000.  ! SLOWS MJO, STRONGER LOBES (Nov 26)
!integer, parameter :: dist_lt = 700.*1000.  ! Do not get an MJO (Nov 26)
!
!integer, parameter :: dist_bl = 1200.*1000.  ! DEF
integer, parameter :: dist_bl = 1100.*1000.
!
double precision, parameter :: lat_rossby_ut = 10.0 !  DEF
!double precision, parameter :: lat_rossby_ut = 12.0  ! NO IMP
!
! lat_rossby_lt = 10 rain moves away from being centered at the equator.
! Increasing lat_rossby_lt also makes MJO bigger and slower
double precision, parameter :: lat_rossby_lt = 5.0 ! DEF
!double precision, parameter :: lat_rossby_lt = 6.0  ! starts separating rain 
!
double precision, parameter :: lat_rossby_bl = 10.0 ! DEF
!double precision, parameter :: lat_rossby_bl = 12.0  ! NO IMP
!
double precision, parameter :: f_move_eq = 0.6 ! DEF
!double precision, parameter :: f_move_eq = 0.4 ! I think noisy

! ZZZ
!double precision, parameter :: zonal_amp = 1.0 ! Nice heart pattern: symmetrical
double precision, parameter :: zonal_amp = 0.95  ! DEF
!double precision, parameter :: zonal_amp = 1.00
!
!double precision, parameter :: zonal_amp = 1.2
!double precision, parameter :: zonal_amp = 1.0
integer :: apply_zonal_amp(nlev) = (/1,1,1/)
double precision :: ros_radius(nlat,nlev)
double precision :: ros_radius_add = 0.20
!------------------------------------------------------------------------
!
!   *******************************************************************
!   ***********************  Radiative Subsidence  ********************
!   *******************************************************************
!
! SD Subsidence
! - modify radiative subsidence depending on local value of extra mass
!   (applied to both layers)
! - Try to make sub_scale as big as possible. Otherwise get too much subsidence
!   from upper troposphere at high rain: better that this get pushed outward
!
!  sub_scale = 120: thought this might slow down the MJO but doesn't
!
!------------------------------------------------------------------------
double precision, parameter :: sub_half = 0.0
double precision, parameter :: sub_scale = 80.0  ! DEF
double precision, parameter :: sub_low = 0.0
double precision, parameter :: sub_add = 2.0
!------------------------------------------------------------------------
!
!  Radiative Cooling
!
!------------------------------------------------------------------------
double precision, parameter :: rad_cool_ut = 0.8 ! K/day
double precision, parameter :: rad_cool_lt = 1.4 ! K/day
double precision, parameter :: rad_cool_bl = 1.3 ! K/day
!------------------------------------------------------------------------
!
!   *******************************************************************
!   ************ Dynamically Forced Ascent/Descent  *******************
!   *******************************************************************
!
!  DN: Dynamical/Wave Motion: subsidence/ascent timescale
!  - proportional to convergence at upper level minus divergence at low level
!  - should be divided by the stability
!  - Needed to create colrh variability to support deep convection
!
!   dmdry_move = f_dyn_mass_ut*(mconv(i,j,3)-mconv(i,j,2)) ! > 0. for down
!
!------------------------------------------------------------------------
double precision, parameter :: sigma_small = 5.0E-07
double precision, parameter :: sigma_scale = 3.0E-05
double precision, parameter :: dmdry_move_max = 1000.
double precision, parameter :: f_dyn_mass_lt = 0.24  ! DEF
double precision, parameter :: f_dyn_mass_ut = 0.24  ! DEF
!double precision, parameter :: f_dyn_mass_lt = 0.23  ! Messy MJO at start, likely west
!double precision, parameter :: f_dyn_mass_ut = 0.23  ! Messy MJO at start, likely west
!double precision, parameter :: f_dyn_mass_lt = 0.22  ! NO MJO propagation
!double precision, parameter :: f_dyn_mass_ut = 0.22  ! NO MJO propagation
!------------------------------------------------------------------------
!
!   *******************************************************************
!   ***********************  DEEP CONVECTION  *************************
!   *******************************************************************
!
! DPP
! DP Deep convection parameters
! If dmdry_cn_max = 100 and rv = 0.02, this is 2 kg/m2 is 30 min. 
! This corresponds to about 100 mm/day
! at rv = 0.02, dmdry_cn_max = 500.0 should produce 10 mm of rain.
!
!------------------------------------------------------------------------
double precision, parameter :: mass_thresh = 1.0E-09
double precision, parameter :: rain_thresh = 0.01
double precision, parameter :: dmdry_cn_max = 500.0
double precision, parameter :: dmdry_bl_fract = 0.3
double precision, parameter :: rh_max_detrain = 1.4
double precision, parameter :: cape_dp_min = 200.
!------------------------------------------------------------------------
!
!  DP Amplitude
!
!------------------------------------------------------------------------
!double precision, parameter :: dp_inc = 3.0
!double precision, parameter :: dp_dec = 5.0
double precision, parameter :: dp_inc = 6.0
double precision, parameter :: dp_dec = 10.0
!------------------------------------------------------------------------
!
! DP Deep convection colrh dependence
!
!------------------------------------------------------------------------
double precision, parameter :: colrh_dp_half = 0.80  ! DEF
!double precision, parameter :: colrh_dp_half = 0.75  ! MJO goes west
!double precision, parameter :: colrh_dp_half = 0.70 ! MJO does not propagate MJO bigger
double precision, parameter :: colrh_dp_scale = 0.02
double precision, parameter :: colrh_dp_low = 0.0
double precision, parameter :: colrh_dp_add = 1.0
!------------------------------------------------------------------------
!
! DP Deep convection caput dependence
! - this should help cut off convection from downdrafts and help move it.
!
!------------------------------------------------------------------------
double precision, parameter :: caput_half = 1200.0
double precision, parameter :: caput_scale = 300.0
double precision, parameter :: caput_low = 0.0
double precision, parameter :: caput_add = 1.0
!------------------------------------------------------------------------
!
!   *******************************************************************
!   *********************** Drizzle Rain   ******************************
!   *******************************************************************
!
!  DR: 
!  - drizzle is slightly noisier with a larger timestep of 1 hour.
!    could make less sensitive to timestep change by decreasing rh_dr_max
!    and increasing timescale.
!
!------------------------------------------------------------------------
!double precision, parameter :: rh_dr_max = 0.9
double precision, parameter :: rh_dr_max = 0.85   ! DEF
double precision, parameter :: tau_dr = 1.0*3600.
!------------------------------------------------------------------------
!
!   *******************************************************************
!   *********************** Anvil Rain   ******************************
!   *******************************************************************
!  AN
!------------------------------------------------------------------------
double precision, parameter :: rh_an_max = 0.9
double precision, parameter :: tau_an = 1.0*3600.
!------------------------------------------------------------------------
!
!   *******************************************************************
!   *********************** Deep Downdrafts  *************************
!   *******************************************************************
!
!  DD: Downdraft parameters
!------------------------------------------------------------------------
integer, parameter :: num_evap = 200
!double precision, parameter :: b_down = -0.05
double precision, parameter :: b_down = -0.01
double precision, parameter :: rv_evap_int = 0.0004
!double precision, parameter :: rh_dn_max = 1.00
double precision, parameter :: rh_dn_max = 0.85
!double precision, parameter :: fevap_max = 0.5
double precision, parameter :: fevap_max = 0.6
!------------------------------------------------------------------------
!
!  Max Fractional mass as a function of updraft mass
!  - 100 kg/m2 with rv = 0.025 has 2.5 kg/m2 water is about 2.5*48 = 125 mm/day
!  - Want fmass_dn_half to correspond to about 30 mm/day
!
!  Used for fractional evaporation and mass flux.
!
!------------------------------------------------------------------------
double precision, parameter :: fevap_half = 25.0  ! DEF
double precision, parameter :: fevap_scale = 8.0  ! DEF
double precision, parameter :: fevap_low = 0.0  ! DEF
double precision, parameter :: fevap_add = 0.8  ! DEF
!
!double precision, parameter :: fevap_add = 0.7  ! MJO with some irregularities
!double precision, parameter :: fevap_add = 0.6  ! MJO destroyed
!  DD=0.5
!double precision, parameter :: fevap_low = 0.6  ! No MJo
!double precision, parameter :: fevap_add = 0.0  ! stationary clumps
!
!double precision, parameter :: fevap_low = 0.7  ! No MJo
!double precision, parameter :: fevap_add = 0.0  ! stationary clumps
!
!double precision, parameter :: fevap_low = 0.8  ! Normal MJO
!double precision, parameter :: fevap_add = 0.0  ! 
!------------------------------------------------------------------------
!
!   DPE Deep evaporation
!
!------------------------------------------------------------------------
double precision, parameter :: rh_dp_evap_max = 0.80  ! DEF
!------------------------------------------------------------------------
!
!   DGE Deep evaporation
!   - fraction of rain available to evaporate as a function of ambient level 2 RH
!
!------------------------------------------------------------------------
double precision, parameter :: f_evap_dp_half = 0.75 
double precision, parameter :: f_evap_dp_scale = 0.05
double precision, parameter :: f_evap_dp_low = 0.5
double precision, parameter :: f_evap_dp_add = -0.5
!------------------------------------------------------------------------
!
!   *******************************************************************
!   ***********************  Cumulous Congestus ***********************
!   *******************************************************************
!
!  pt_add: Level 2 PT must be this much larger than level 1 PT for evap 
!   to occur. Try and avoid instability.
!  CGG
!------------------------------------------------------------------------
!double precision, parameter :: rh_cg_evap_max = 0.85
double precision, parameter :: rh_cg_evap_max = 0.80
double precision, parameter :: pt_add = 2.0
!------------------------------------------------------------------------
!
!   CGE Congestus evaporation
!   June 1: seems to be better not to set f_evap_cg_low to big (like 0.8).
!   Seems to generate negative sigma between BL and LT.
!------------------------------------------------------------------------
double precision, parameter :: f_evap_cg_half = 0.75
double precision, parameter :: f_evap_cg_scale = 0.05
double precision, parameter :: f_evap_cg_low = 0.5
double precision, parameter :: f_evap_cg_add = -0.5
!------------------------------------------------------------------------
!
!  CG
!
!------------------------------------------------------------------------
double precision, parameter :: cg_amp_d = 3.0
double precision, parameter :: cg_amp_c = 0.0
!------------------------------------------------------------------------
!
! CG Congestus colrh dependence
!
!------------------------------------------------------------------------
double precision, parameter :: colrh_cg_half = 0.65   ! def
double precision, parameter :: colrh_cg_scale = 0.05  ! def
double precision, parameter :: colrh_cg_low = 0.0
double precision, parameter :: colrh_cg_add = 1.0
!------------------------------------------------------------------------
!
! CG Congestus caplt dependence
!
!------------------------------------------------------------------------
!double precision, parameter :: caplt_half = 150.0  ! def
!double precision, parameter :: caplt_scale = 50.0  ! def
double precision, parameter :: caplt_half = 100.0
double precision, parameter :: caplt_scale = 25.0
double precision, parameter :: caplt_low = 0.0
double precision, parameter :: caplt_add = 1.0
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
integer, parameter :: i_prect = 1
integer, parameter :: i_t_bl = 2
integer, parameter :: i_t_lt = 3
integer, parameter :: i_t_ut = 4
integer, parameter :: i_cape_ut = 5
integer, parameter :: i_cape_lt = 6
integer, parameter :: i_mass_dp = 7
integer, parameter :: i_mass_dn = 8
integer, parameter :: i_mass_cg = 9
integer, parameter :: i_f_dp = 10
integer, parameter :: i_f_cg = 11
integer, parameter :: i_f_caput = 12
integer, parameter :: i_f_caplt = 13
integer, parameter :: i_f_colrh_dp = 14
integer, parameter :: i_f_colrh_cg = 15
integer, parameter :: i_colrh = 16
integer, parameter :: i_rh_bl = 17
integer, parameter :: i_rh_lt = 18
integer, parameter :: i_rh_ut = 19
integer, parameter :: i_mconv1 = 20
integer, parameter :: i_mconv2 = 21
integer, parameter :: i_mconv3 = 22
integer, parameter :: i_mass_lt_bl = 23
integer, parameter :: i_mass_ut_lt = 24
integer, parameter :: i_rain_an = 25
integer, parameter :: i_rain_cg = 26
integer, parameter :: i_zz_bl = 27
integer, parameter :: i_zz_lt = 28
integer, parameter :: i_zz_ut = 29
integer, parameter :: i_down_dp = 30
integer, parameter :: i_mconv_diff = 31
integer, parameter :: i_sub_lt_bl = 32
integer, parameter :: i_sub_ut_lt = 33
integer, parameter :: i_sum_lt_bl = 34
integer, parameter :: i_sum_ut_lt = 35
integer, parameter :: i_evap_dp = 36
integer, parameter :: i_rain_dr = 37   ! NEW
integer, parameter :: i_evap_cg = 38   ! NEW
integer, parameter :: i_mextra_ut = 39
integer, parameter :: i_mextra_lt = 40
integer, parameter :: i_mextra_bl = 41
integer, parameter :: i_colrh_diff = 42  ! NEW
integer, parameter :: i_rain_tt_diff = 43  ! NEW
integer, parameter :: i_mass_ut_lt_diff = 44  ! NEW
integer, parameter :: i_mass_lt_bl_diff = 45  ! NEW
integer, parameter :: i_sub_ut_lt_diff = 46  ! NEW
integer, parameter :: i_sub_lt_bl_diff = 47  ! NEW
integer, parameter :: i_rv_bl = 48  ! NEW
integer, parameter :: i_dp_prod = 49  ! NEW
integer, parameter :: i_mextra_diff = 50  ! NEW
integer, parameter :: i_mextra_all = 51  ! NEW
integer, parameter :: i_mconv_all = 52  ! NEW
!------------------------------------------------------------------------
!
!  Rain Events  EV
!  - MUST CHANGE nsteps_stuff and ns_mid if change timestep!!!
!
!------------------------------------------------------------------------
double precision, parameter :: rain_event_min = 25.0  ! DEF
!double precision, parameter :: rain_event_min = 40.0
integer, parameter :: nsteps_events = 100 ! min timesteps before start storing events
integer, parameter :: nsteps_stuff = 25 ! -12 hours to +12 hours; 1 hour timestep
integer, parameter :: ns_mid = 13  ! 1 hour timestep
integer, parameter :: nvar_stuff = 52
double precision :: stuff(nlon,nlat,nsteps_stuff,nvar_stuff)
double precision :: events(nsteps_stuff,nvar_stuff)
double precision :: cofm(nvar_stuff)
double precision, parameter :: cofm_lat = 30.0
character(len=10) :: name_var(nvar_stuff)
!------------------------------------------------------------------------
!
!  EVV
!
!------------------------------------------------------------------------
integer, parameter :: num_radial = 8
integer, parameter :: num_var = 5
double precision, parameter :: dist_max = 1000.
double precision rdist(num_radial),dist_int 
double precision var_rad(num_var,num_radial)
double precision num_rad(num_radial)
!------------------------------------------------------------------------
!
!  Divergence DIV
!
!------------------------------------------------------------------------
integer, parameter :: num_div = 18
double precision :: pp_div(num_div),div(nsteps_stuff,num_div),mm_ut,mm_lt,mm_bl
double precision :: div_ut,div_lt,div_bl,pph_1,pph_2,pph_3,pph_4
!------------------------------------------------------------------------
!
!  mjo_mean
!  mjo basically 80 degrees wide, or ~50 grid cells
!  mjo separated by 170 degrees, so can go 90 degrees either side so ndist can be 60
!  MMM
!
!------------------------------------------------------------------------
!integer, parameter :: ndist_speed = 35
!integer, parameter :: ndist_speed = 50
integer, parameter :: ndist_mjo = 60  ! DEF
integer, parameter :: j_lat = 21
!double precision, parameter :: rain_mjo_min = 20.0
double precision, parameter :: rain_mjo_min = 2.0
double precision :: prectt_smooth(nlon_mjo),lon_mjo(nlon_mjo)
double precision :: mjo_mean(nlon_mjo,nlat,nvar_stuff)
double precision :: temp(nlon_mjo,nlat)
double precision :: uu(nlon_mjo,nlat,nlev)
double precision :: vv(nlon_mjo,nlat,nlev)
double precision :: mjo_mean_lon_15S_15N(nlon_mjo,nvar_stuff)
double precision :: mjo_mean_lon_10S_10N(nlon_mjo,nvar_stuff)
double precision :: mjo_mean_lon_5S_5N(nlon_mjo,nvar_stuff)
double precision :: rain_big,yyy_smooth,mjo_sp,mass_dp_max
integer :: i_mjo,ndiff,num_mjo,i_lead
!------------------------------------------------------------------------
!
!***********************************************************************
!*****************  end module toy  ************************************
!***********************************************************************
!
!------------------------------------------------------------------------
end module model
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  program toy
!
!------------------------------------------------------------------------
program toy
use netcdf
use model
implicit none
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
integer :: ncid, status
integer :: TimDimID,LatDimID,LonDimID,LevDimID
!------------------------------------------------------------------------
!
!  Local Real Variables
!
!------------------------------------------------------------------------
double precision :: w_a,colwat_a,colwat_b,rv_b,dmdry_b,dm_have,pt_diff,sst_senflux,sst_watflux,convert
double precision :: x1,x2,xdiff,xxx,zzz,yyy,dx,dx_deg,dy,es,dz1,dz2,rvs,tvp,tdd,tdp,dlnp,dtd_lt,dtd,dtd_ut,rv_avail
double precision :: hmp,kmp,rtp,rvp,rip,rlp,tp,zp,rv_bl,z_bl,bowen,bowen_av,xxx1,xxx2,xxx3,xlon_mjo,xlon_mjo_prev
double precision :: colwat_sat,colwat,es_t,drv,dtt,es_sst,drv_av,dtt_av,ps_av,xlon_weight
double precision :: dmdry_av,colwat_av,colwat_av_prev,rh_bl_av,rh_lt_av,rh_ut_av
double precision :: wat_added,watflux_av,cape_lt_av,cape_ut_av,colrh_av,zz_bl,zz_lt,zz_ut
double precision :: zz_bl_av,zz_lt_av,zz_ut_av,dpp,dpp1,dpp2
double precision :: f_colrh,e_bl,rhp,tt_bl_av,tt_lt_av,tt_ut_av,wat_error
double precision :: rain_mmday,days
double precision :: dm_initial,day_max,dm_bl
double precision :: colwat_bf_dn,colwat_af_dn,mextra_lt_initial
double precision :: hm_bl,rvsp,dmdry_bl,rv_final,rv_initial,hm_initial,dmdry_final,dmdry_initial
double precision :: rt_final,rl_final,ri_final,hm_final,km_final,tp_final
double precision :: rv_lt,rv_evap,km_wat,bp,rvp_bl,tvp_bl,colwat_start,colwat_end
double precision :: rv_bl_initial,dmdry_bl_initial,dmdry_bl_final,rvmass_bl_initial,rv_bl_final,rv_down
double precision :: hm_bl_initial,hmmass_bl,hm_bl_final,hm_down,rvmass_bl_final
double precision :: rvmass_lt_initial,rvmass_lt_final,dmdry_lt_initial,dmdry_lt_final,rv_lt_initial
double precision :: dmdry_end,dmdry_start,dmdry_move,rvmass_initial
double precision :: dmdry_bl_start,dmdry_lt_start,dmdry_ut_start,dmdry_bl_tot,dmdry_lt_tot,dmdry_ut_tot
double precision :: dmdry_bl_end,dmdry_lt_end,dmdry_ut_end,drymass_bl_error,drymass_lt_error,drymass_ut_error
double precision :: mass_cg_av,mass_dp_av,mass_dn_av,mass_bl_lt_av,mass_lt_bl_av,mass_ut_lt_av
double precision :: mextra_bl_av,mextra_lt_av,mextra_ut_av,fff
double precision :: pp1,pp2,pp3,phh,tt1,pt1,p1_bot,p2_bot,p1_top,p2_top,fevap
double precision :: km_initial,pppp,rvmass_evap,rvmass_max,rvmass_have,t_final,t_initial,dmdry_error
double precision :: f_mass_dn_max,rh_max,rv_bl_av,rv_lt_av,rv_ut_av
!------------------------------------------------------------------------
!
!  Local Real Variables Arrays: averages
!
!------------------------------------------------------------------------
double precision :: rain_mn(nlon,nlat)
double precision :: rain_an_mn(nlon,nlat)
double precision :: rain_dp_mn(nlon,nlat)
double precision :: down_dp_mn(nlon,nlat)
double precision :: evap_dp_mn(nlon,nlat)
double precision :: rain_cg_mn(nlon,nlat)
double precision :: evap_cg_mn(nlon,nlat)
double precision :: colrh_mn(nlon,nlat)
double precision :: rh_bl_mn(nlon,nlat)
double precision :: rh_lt_mn(nlon,nlat)
double precision :: rh_ut_mn(nlon,nlat)
double precision :: t_bl_mn(nlon,nlat)
double precision :: t_lt_mn(nlon,nlat)
double precision :: t_ut_mn(nlon,nlat)
double precision :: mass_dn_mn(nlon,nlat)
double precision :: mass_ut_lt_mn(nlon,nlat)
double precision :: mass_lt_bl_mn(nlon,nlat)
double precision :: mextra_ut_mn(nlon,nlat)
double precision :: mextra_lt_mn(nlon,nlat)
double precision :: mextra_bl_mn(nlon,nlat)
double precision :: mdiff_lt(nlon,nlat)
double precision :: mdiff_ut(nlon,nlat)
double precision :: senflux(nlon,nlat)
double precision :: watflux(nlon,nlat)
!------------------------------------------------------------------------
!
!  Speed PDF
!
!------------------------------------------------------------------------
integer, parameter :: num_pdf = 50
double precision :: speed_pdf(num_pdf),speed_pdf_tot,speed1,speed2,pdf_int
double precision :: speed_pdf_grid(num_pdf)
double precision, parameter :: speed_pdf_min = -20.
double precision, parameter :: speed_pdf_max = 30.
!------------------------------------------------------------------------
!
! WKK
! can't store the output file in a local directory since too big.
!
!------------------------------------------------------------------------
character (len = *), parameter :: dir_temp = "/temp-1/folkins/toy/"
character (len = *), parameter :: FILE_NAME_PRECT_WK = dir_temp//"rain.WK.nc"
integer :: ncid_wk
integer :: ncid_ts
integer, parameter :: NLATS_WK = 21 ! at equator + 10 on either side +-15
integer, parameter :: NLONS_WK = nlon
integer, parameter :: NDIMS_WK = 3 ! rain has 3 dimensions: lat,lon,time
integer, parameter :: NDIMS_TS = 3 ! time has one dimension
!------------------------------------------------------------------------
!
!  WKK:
!  You can't use NRECS_WK = nsteps because prect_wk would then be too big.
!  Use a 3 hour interval.
!  You also want to exclude about 45 days at the beginning of the record, for
!  the model to spin up.
! 
!------------------------------------------------------------------------
integer, parameter :: ndays_spinup = 25
!integer, parameter :: ndays_spinup = 2
!integer, parameter :: ndays_spinup = 0
integer, parameter :: NRECS_WK = nsteps - ndays_spinup*nsteps_day  ! exclude first spinup days
integer, parameter :: NRECS_TS = nsteps - ndays_spinup*nsteps_day  ! exclude first spinup days
!------------------------------------------------------------------------
!
! WK Dimension Names
!
!------------------------------------------------------------------------
character (len = *), parameter :: LAT_NAME = "latitude"
character (len = *), parameter :: LON_NAME = "longitude"
character (len = *), parameter :: REC_NAME = "time"
integer :: lon_dimid, lat_dimid, rec_dimid
!------------------------------------------------------------------------
!
! WK
! The start and count arrays will tell the netCDF library where to
! write our data.
!
!------------------------------------------------------------------------
integer :: start(NDIMS_WK), count(NDIMS_WK)
integer :: start_ts(NDIMS_TS), count_ts(NDIMS_TS)
!------------------------------------------------------------------------
!
! WK
! These program variables hold the latitudes and longitudes.
!
!------------------------------------------------------------------------
real :: lats_wk(NLATS_WK), lons_wk(NLONS_WK), tim_wk(NRECS_WK)
integer :: lon_varid, lat_varid, rec_varid
!------------------------------------------------------------------------
!
!  WK variable array
!
!------------------------------------------------------------------------
character (len = *), parameter :: RAIN_NAME="rain"
integer :: rain_varid
integer :: dimids_wk(NDIMS_WK)
integer :: dimids_ts(NDIMS_TS)
!------------------------------------------------------------------------
!
!  WKK
!
!  Wheeler-Kiladis paper: Every 6 hours tropical values
!  - This is for the wk.ncl
!  - prect_wk is not written out at every timestep but only at the end, so
!    is a huge file.
!  Here NLONS_WK = 240 and NLATS_WK = 21, so each timestep is 5040 numbers
!  In /trmm/read.50 I had an array of 155 million numbers.
!  This would appear to allow 30754 timesteps, or 1280 days at 1 hour timesteps.
!  This would allow 3.5 years, so maybe I can store at 1 hour resolution.
!
!------------------------------------------------------------------------
real :: prect_wk(NLONS_WK,NLATS_WK,NRECS_WK)
!------------------------------------------------------------------------
!
!  Local Real Variables Arrays
!
!------------------------------------------------------------------------
double precision :: rho(nlon,nlat,nlev)  ! not defined
double precision :: rhomid1(nlon,nlat),rhomid2(nlon,nlat)
!------------------------------------------------------------------------
!
!  HOV variables
!
!------------------------------------------------------------------------
double precision :: hov_15(nlon,nsteps)
double precision :: hov_colrh(nlon,nsteps)
integer, parameter :: n_hov_start = 400
integer, parameter :: n_hov_end = 800
double precision, parameter :: lat_hov = 16.
!------------------------------------------------------------------------
!
!  Diagnostics - Rain pdf
!
!------------------------------------------------------------------------
integer, parameter :: num_rain_int = 100     !  number of rain grid values
double precision, parameter :: rain_seed = 0.01
double precision, parameter :: rain_mult = 1.1
double precision :: grid_pdf(num_rain_int)
double precision :: rain_pdf(num_rain_int)
double precision :: rain1, rain2, rr
!------------------------------------------------------------------------
!
!  Local Integer Variables
!
!------------------------------------------------------------------------
integer :: nnn,ii,i,j,k,jj,k1,k2,i1,i2,j1,j2,niter,ndir,num,nprint,ntot,nr,nr_want,ns,nv
integer :: num_events,num_wk,num_15,num_10,num_5,nd
!------------------------------------------------------------------------
!
!  ***********************************************************************
!  ***************  START PROGRAM  ***************************************
!  ***********************************************************************
!
!
!  Define dx,dy,lat,lon
!
!------------------------------------------------------------------------
tstep = 3600.0*24.0/float(nsteps_day)
dx = (2.0*pi*re/float(nlon))
dy = (2.0*pi*re/float(nlat-1))*(2.0*dlat/360.0)
!print*,'dx = ',dx
!print*,'dy = ',dy
do i = 1,nlon
  xxx = 360./float(nlon)
  lon(i) = float(i-1)*xxx
  lon_dyn(i) = lon(i) + 0.5*xxx
!  print*,'i lon(i) = ',i,lon(i)
end do
do j = 1,nlat
  xxx = 2.0*dlat/float(nlat-1)
  lat(j) = -dlat + float(j-1)*xxx
  lat_dyn(j) = lat(j) + 0.5*xxx
!  print*,'j lat(j) = ',j,lat(j)
end do
!------------------------------------------------------------------------
!
!  EVV
!
!------------------------------------------------------------------------
dist_int = dist_max/float(num_radial)
do nr = 1,num_radial
  rdist(nr) = float(nr-1)*dist_int
  num_rad(nr) = 0.
  do nv = 1,num_var
    var_rad(nv,nr) = 0.
  end do
end do
!------------------------------------------------------------------------
!
! Initialize name_var
!
!------------------------------------------------------------------------
name_var(1) = 'PRECT     '
name_var(2) = 'T_BL      '
name_var(3) = 'T_LT      '
name_var(4) = 'T_UT      '
name_var(5) = 'CAPE_UT   '
name_var(6) = 'CAPE_LT   '
name_var(7) = 'MASS_DP   '
name_var(8) = 'MASS_DN   '
name_var(9) = 'MASS_CG   '
name_var(10) = 'F_DP      '
name_var(11) = 'F_CG      '
name_var(12) = 'F_CAPUT   '
name_var(13) = 'F_CAPLT   '
name_var(14) = 'F_COLRH_DP'
name_var(15) = 'F_COLRH_CG'
name_var(16) = 'COLRH     '
name_var(17) = 'RH_BL     '
name_var(18) = 'RH_LT     '
name_var(19) = 'RH_UT     '
name_var(20) = 'MCONV1    '
name_var(21) = 'MCONV2    '
name_var(22) = 'MCONV3    '
name_var(23) = 'MASS_LT_BL'
name_var(24) = 'MASS_UT_LT'
name_var(25) = 'RAIN_AN   '
name_var(26) = 'RAIN_CG   '
name_var(27) = 'ZZ_BL     '
name_var(28) = 'ZZ_LT     '
name_var(29) = 'ZZ_UT     '
name_var(30) = 'DOWN_DP   '
name_var(31) = 'MCONV_DIFF'
name_var(32) = 'SUB_LT_BL '
name_var(33) = 'SUB_UT_LT '
name_var(34) = 'SUM_LT_BL '
name_var(35) = 'SUM_UT_LT '
name_var(36) = 'EVAP_DP   '
name_var(37) = 'RAIN_DR   '
name_var(38) = 'EVAP_CG   '
name_var(39) = 'MEXTRA_UT '
name_var(40) = 'MEXTRA_LT '
name_var(41) = 'MEXTRA_BL '
name_var(42) = 'COLRH_DIFF'
name_var(43) = 'RAIN_TT_DI'
name_var(44) = 'M_UT_LT_DI'
name_var(45) = 'M_LT_BL_DI'
name_var(46) = 'S_UT_LT_DI'
name_var(47) = 'S_LT_BL_DI'
name_var(48) = 'RV_BL     '
name_var(49) = 'RV_BL_DIFF'
name_var(50) = 'MEXTRA_DIF'
name_var(51) = 'MEXTRA_ALL'
!------------------------------------------------------------------------
!
!  Define grid_pdf
!
!------------------------------------------------------------------------
grid_pdf(1) = rain_seed
rain_pdf(1) = 0.
do nr = 2,num_rain_int
  grid_pdf(nr) = rain_mult*grid_pdf(nr-1)
  rain_pdf(nr) = 0.
end do
!------------------------------------------------------------------------
!
!  Initialize hov
!
!------------------------------------------------------------------------
do i = 1,nlon
do nstep = 1,nsteps
  hov_15(i,nstep) = 0.
  hov_colrh(i,nstep) = 0.
end do
end do
!------------------------------------------------------------------------
!
!  Initialize mean fields before start
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  wdelta_bl_mn(i,j) = 0.
  wconv_bl_mn(i,j) = 0.
  rain_mn(i,j) = 0.
  rain_an_mn(i,j) = 0.
  rain_dp_mn(i,j) = 0.
  down_dp_mn(i,j) = 0.
  evap_dp_mn(i,j) = 0.
  rain_cg_mn(i,j) = 0.
  evap_cg_mn(i,j) = 0.
  cape_lt_mn(i,j) = 0.
  cape_ut_mn(i,j) = 0.
  colrh_mn(i,j) = 0.
  t_bl_mn(i,j) = 0.
  t_lt_mn(i,j) = 0.
  t_ut_mn(i,j) = 0.
  rh_bl_mn(i,j) = 0.
  rh_lt_mn(i,j) = 0.
  rh_ut_mn(i,j) = 0.
  f_caput_mn(i,j) = 0.
  f_caplt_mn(i,j) = 0.
  f_colrh_dp_mn(i,j) = 0.
  f_colrh_cg_mn(i,j) = 0.
  f_colrh_dp(i,j) = 0.
  f_colrh_cg(i,j) = 0.
  senflux_mn(i,j) = 0.
  watflux_mn(i,j) = 0.
  sst_mn(i,j) = 0.
  sigma1_mn(i,j) = 0.
  sigma2_mn(i,j) = 0.
  f_dp_mn(i,j) = 0.
  f_cg_mn(i,j) = 0.
  mass_ut_lt_mn(i,j) = 0.
  mass_lt_bl_mn(i,j) = 0.
  sub_ut_lt_mn(i,j) = 0.
  sub_lt_bl_mn(i,j) = 0.
  mass_dp_mn(i,j) = 0.
  mass_cg_mn(i,j) = 0.
  mass_dn_mn(i,j) = 0.
  mextra_ut_mn(i,j) = 0.
  mextra_lt_mn(i,j) = 0.
  mextra_bl_mn(i,j) = 0.
  do k = 1,nlev
    mconv_mn(i,j,k) = 0.
    zz_mn(i,j,k) = 0.
  end do
end do
end do
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
speed_mjo_av = 0.
num_speed_mjo = 0
num_speed_mjo_av = 0
pdf_int = (speed_pdf_max - speed_pdf_min)/float(num_pdf)
speed_pdf_tot = 0.
do num = 1,num_pdf
  speed_pdf_grid(num) = speed_pdf_min + 0.5*pdf_int + float(num-1)*pdf_int
  speed_pdf(num) = 0.
end do
!------------------------------------------------------------------------
!
! Initialze ph and pp before timesteps
! These never change.
!
!------------------------------------------------------------------------
ph(1) = p1
ph(2) = p2
ph(3) = p3
ph(4) = p4
do k = 1,nlev
  pp(k) = 0.5*(ph(k) + ph(k+1))
end do
!------------------------------------------------------------------------
!
!  Initialize sst_tar Before timesteps
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dp_prod(i,j) = 0.0
  mass_dp(i,j) = 0.0
  mass_dp_old(i,j) = 0.0
  mass_cg(i,j) = 0.0
  xxx = abs((lat(j)*dy)/sst_dist)
  sst_tar(i,j) = sst_eq*exp(-xxx*xxx) + tkelvin + sst_back
  sst(i,j) = sst_tar(i,j)
  mass_lt_bl(i,j) = 0. 
  do k = 1,nlev
    mextra(i,j,k) = 0.
  end do
end do
end do
!------------------------------------------------------------------------
!
!  Define initial tt/rv/zz
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
do k = 1,nlev
!------------------------------------------------------------------------
!  Define tt
!------------------------------------------------------------------------
  if (k == 1) tt(i,j,k) = t_bl_init
  if (k == 2) tt(i,j,k) = t_lt_init
  if (k == 3) tt(i,j,k) = t_ut_init
!------------------------------------------------------------------------
!  Define rv
!------------------------------------------------------------------------
  rvs = rsat(tt(i,j,k),pp(k))
  if (k == 1) rv(i,j,k) = rh_bl_init*rvs
  if (k == 2) rv(i,j,k) = rh_lt_init*rvs
  if (k == 3) rv(i,j,k) = rh_ut_init*rvs
!------------------------------------------------------------------------
!  Check for negative rv
!------------------------------------------------------------------------
  if (rv(i,j,k) < 0.) then
    print*,'rv(i,j,k) = ',rv(i,j,k)
    stop 'rv negative at initialization stage'
  endif
!------------------------------------------------------------------------
!  Define zz
!------------------------------------------------------------------------
  if (k == 1) zz(i,j,k) = z1
  if (k == 2) zz(i,j,k) = z2
  if (k == 3) zz(i,j,k) = z3
!------------------------------------------------------------------------
!  Define hm from tt/rv/zz
!  From now on, hm treated as a fundamental variable
!------------------------------------------------------------------------
  rtp = rv(i,j,k)
  rip = 0.
  hm(i,j,k) = hmoist(tt(i,j,k),rtp,rv(i,j,k),rip,zz(i,j,k))
!------------------------------------------------------------------------
!  Check for bad hm
!------------------------------------------------------------------------
  if (hm(i,j,k) > 410000.) then
    print*,'-------- Problem with unrealistic hm during initialization ----'
    print*,'tt(i,j,k) = ',tt(i,j,k)
    print*,'rv(i,j,k) = ',rv(i,j,k)
    print*,'zz(i,j,k) = ',zz(i,j,k)
    print*,'i j k = ',i,j,k
    print*,'hm(i,j,k) = ',hm(i,j,k)
    stop '---------- bad hm ---'
  endif
!------------------------------------------------------------------------
!  Define ppp
!------------------------------------------------------------------------
  ppp(i,j,k) = pp(k)
  pph(i,j,k) = ph(k)
  if (k == 3) pph(i,j,k+1) = ph(k+1)
!------------------------------------------------------------------------
!  end loop over i,j,k
!------------------------------------------------------------------------
end do
end do
end do
print*,'Initialized T rv zz hm'
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dmdry = (ph(1) - ph(2))/(g*(1.0+rv(i,j,1))) + mextra(i,j,1)
  wat_bl_prev(i,j) = wat_bl(i,j)
  wat_bl(i,j) = rv(i,j,1)*dmdry
end do
end do
!------------------------------------------------------------------------
!
!  Redefine tt/zz from hm/rv/ppp/pph to get a consistent set
!
!------------------------------------------------------------------------
print*,'Calling height at start of program'
call height
print*,'finished call to height at start of program'
!------------------------------------------------------------------------
!
!  Initialize events
!
!------------------------------------------------------------------------
num_events = 0
do ns = 1,nsteps_stuff
do nv = 1,nvar_stuff
   events(ns,nv) = 0.
end do
end do
print*,'finished initialization'
!------------------------------------------------------------------------
!
!  MMM
!  Initialize
!
!------------------------------------------------------------------------
num_mjo = 0
do i = 1,nlon_mjo
do j = 1,nlat
do nv = 1,nvar_stuff
  mjo_mean(i,j,nv) = 0.
end do
end do
end do
do i = 1,nlon_mjo
  lon_mjo(i) = float(i-1)*(360./float(nlon))
end do
!------------------------------------------------------------------------
!
!  ***********   Start Loop over Timesteps  ******************8
!
!------------------------------------------------------------------------
do nstep = 1,nsteps
if (mod(nstep,1) == 0) print*,'Starting timestep = ',nstep
!if (mod(nstep,10) == 0) print*,'Starting timestep = ',nstep
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  Starting timestep nstep = ',nstep
    print*,'starting tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'starting tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'starting tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'starting pt(i_want,j_want,1) = ',pt(i_want,j_want,1)
    print*,'starting pt(i_want,j_want,2) = ',pt(i_want,j_want,2)
    print*,'starting pt(i_want,j_want,3) = ',pt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
!  Initialize every timestep
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
!  wconv_bl(i,j) = 0.
  rain_an(i,j) = 0.
  rain_dp(i,j) = 0.
  rain_dr(i,j) = 0.
  evap_dp(i,j) = 0.
  down_dp(i,j) = 0.
  rain_cg(i,j) = 0.
  evap_cg(i,j) = 0.
  mass_dn(i,j) = 0.
  f_dp(i,j) = 0.
  f_cg(i,j) = 0.
  f_caput(i,j) = 0.
  f_caplt(i,j) = 0.
  f_colrh_dp(i,j) = 0.
  f_colrh_cg(i,j) = 0.
  mass_dp_old(i,j) = mass_dp(i,j)
  do k = 1,nlev
    rv_prev(i,j,k) = rv(i,j,k)
  end do
end do
end do
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dmdry = (ph(1) - ph(2))/(g*(1.0+rv(i,j,1))) + mextra(i,j,1)
  wat_bl(i,j) = rv(i,j,1)*dmdry
  wdelta_bl(i,j) = wat_bl(i,j) - wat_bl_prev(i,j)
end do
end do
!------------------------------------------------------------------------
!
!  Define tt (in timestep)
!
!  - Regard hm/zz/rv as fundamental: get km
!  - Then find tt
!  - This is not entirely self consistent approach since the zz obtained below
!    from hydrostatic balance will not form a consistent set. However, error
!    in tt/zz from doing this likely small if timestep not too big.
!  - Likely some way to iteratively derive a self consistent profile.
!  - Assumes no background condensate
!
!------------------------------------------------------------------------
!print*,'Calling height at start of timestep'
call height
!print*,'finished call to height at start of timestep'
!------------------------------------------------------------------------
!
!  Define stability  (in timestep)
!  where the static stability is sigma_p = -(T/theta)*(dtheta/dp)
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  sigma between BL and LT
!------------------------------------------------------------------------
  pp1 = ppp(i,j,1)
  pp2 = ppp(i,j,2)
  phh = pph(i,j,2)
  tt1 = (pp1 - phh)*tt(i,j,1) + (phh - pp2)*tt(i,j,2)  ! mean T
  tt1 = tt1/(pp1 - pp2)
  pt1 = (pp1 - phh)*pt(i,j,1) + (phh - pp2)*pt(i,j,2)  ! mean PT
  pt1 = pt1/(pp1 - pp2)
  sigma1(i,j) = -(tt1/pt1)*((pt(i,j,2)-pt(i,j,1))/(pp2-pp1))
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (sigma1(i,j) /= sigma1(i,j)) then
     print*,'-------- sigma1 problem -----'
     print*,'tt1 = ',tt1
     print*,'pt1 = ',pt1
     print*,'pp2 - pp1 = ',pp2-pp1
     print*,'sigma1(i,j) = ',sigma1(i,j)
     stop 'sigma1(i,j) NaN after def'
  endif
!------------------------------------------------------------------------
! Check for negative sigma1
!------------------------------------------------------------------------
  if (sigma1(i,j) < 0.) then
    if (nprint_sigma == 1) then
    print*,'--------  WARNING sigma1 negative ------------'           
    print*,'i j = ',i,j
    print*,'evap_dp(i,j) = ',evap_dp(i,j)
    print*,'sigma1(i,j) = ',sigma1(i,j)
    print*,'tt(i,j,1) = ',tt(i,j,1)
    print*,'tt(i,j,2) = ',tt(i,j,2)
    print*,'pt(i,j,1) = ',pt(i,j,1)
    print*,'pt(i,j,2) = ',pt(i,j,2)
    print*,'pt(i,j,2)-pt(i,j,1) = ',pt(i,j,2)-pt(i,j,1)
    print*,'tt1 = ',tt1
    print*,'mean pt pt1 = ',pt1
    print*,'pp1 = ',pp1
    print*,'pp2 = ',pp2
    print*,'mextra(i,j,1) = ',mextra(i,j,1)
    print*,'mextra(i,j,2) = ',mextra(i,j,2)
    print*,'mextra(i,j,3) = ',mextra(i,j,3)
    endif 
    sigma1(i,j) = sigma_small
!    stop '------------- sigma1 negative -------------'
  endif
!------------------------------------------------------------------------
! sigma1 info
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'--------  sigma1 info ------------'           
    print*,'i j = ',i,j
    print*,'evap_dp(i,j) = ',evap_dp(i,j)
    print*,'sigma1(i,j) = ',sigma1(i,j)
    print*,'pt(i,j,2) = ',pt(i,j,2)
    print*,'pt(i,j,1) = ',pt(i,j,1)
    print*,'pt(i,j,2)-pt(i,j,1) = ',pt(i,j,2)-pt(i,j,1)
    print*,'tt1 = ',tt1
    print*,'mean pt pt1 = ',pt1
    print*,'pp1 = ',pp1
    print*,'pp2 = ',pp2
    print*,'mextra(i,j,1) = ',mextra(i,j,1)
    print*,'mextra(i,j,2) = ',mextra(i,j,2)
    print*,'mextra(i,j,3) = ',mextra(i,j,3)
    print*,'------------- end info -------------'
  endif
!------------------------------------------------------------------------
!  sigma2 between LT and UT
!------------------------------------------------------------------------
  pp1 = ppp(i,j,2) 
  pp2 = ppp(i,j,3)
  phh = pph(i,j,3)
  tt1 = (pp1 - phh)*tt(i,j,2) + (phh - pp2)*tt(i,j,3)  ! Mean layer T
  tt1 = tt1/(pp1 - pp2)
  pt1 = (pp1 - phh)*pt(i,j,2) + (phh - pp2)*pt(i,j,3)  ! Mean layer PT
  pt1 = pt1/(pp1 - pp2)
  sigma2(i,j) = -(tt1/pt1)*((pt(i,j,3)-pt(i,j,2))/(pp2-pp1))
  if (sigma2(i,j) /= sigma2(i,j)) stop 'sigma2(i,j) NaN after def'
!------------------------------------------------------------------------
! Check for negative sigma2
!------------------------------------------------------------------------
  if (sigma2(i,j) < 0.) then
    if (nprint_sigma == 1) then
    print*,'--------  WARNING sigma2 negative ------------'           
    print*,'sigma2(i,j) = ',sigma2(i,j)
    print*,'pt(i,j,3) = ',pt(i,j,3)
    print*,'pt(i,j,2) = ',pt(i,j,2)
    print*,'pt(i,j,3)-pt(i,j,2) = ',pt(i,j,3)-pt(i,j,2)
    print*,'tt1 = ',tt1
    print*,'pt1 = ',pt1
    print*,'pp1 = ',pp1
    print*,'pp2 = ',pp2
    print*,'mextra(i,j,1) = ',mextra(i,j,1)
    print*,'mextra(i,j,2) = ',mextra(i,j,2)
    print*,'mextra(i,j,3) = ',mextra(i,j,3)
    endif
    sigma2(i,j) = sigma_small
!    stop '------------- sigma2 negative -------------'
  endif
!------------------------------------------------------------------------
!  Calculate sigma
!   omega[Pa/s] = Q[J/s]/(cp*sigma_p)
!   qq is J/s
!   cpd is J/kg*K
!   convert Q to K/s 
!   multiply by cpd to get K/s*(J/(kg*K)) = J/(kg*s)
!   rad_cool_ut = K/day
!   [qq1] = (K/day)*(J/(kg*K))/(s/day))
!         = K*(J/(kg*K)/s
!         = (J/(kg)/s
!         = J/kg*s
!
!------------------------------------------------------------------------
   qq1 = (rad_cool_ut*cpd)/sec_day  !  K/day*(day/s*(J/(kg*K))
   qq2 = (rad_cool_lt*cpd)/sec_day  !  K/day*(day/s*(J/(kg*K))
   omega1(i,j) = qq1/(cpd*sigma1(i,j))  ! should be positive down  LT -> BL
   omega2(i,j) = qq2/(cpd*sigma2(i,j))  ! should be positive down  UT-> LT
!------------------------------------------------------------------------
!------------------------------------------------------------------------
   if (sigma1(i,j) /= sigma1(i,j)) stop 'sigma1(i,j) NaN'
   if (sigma2(i,j) /= sigma2(i,j)) stop 'sigma1(i,j) NaN'
   if (omega1(i,j) /= omega1(i,j)) stop 'omega1(i,j) NaN'
   if (omega2(i,j) /= omega2(i,j)) stop 'omega2(i,j) NaN'
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
     print*,'-------- omega calculation --- nstep = ',nstep
     print*,'sigma2(i,j) = ',sigma2(i,j)
     print*,'omega2(i,j) = ',omega2(i,j)
     print*,'qq2 = ',qq2
  endif
!  print*,'pp1 = ',pp1
!  print*,'pp2 = ',pp2
!  print*,'phh = ',phh
!  print*,'tt1 = ',tt1
!  print*,'pt1 = ',pt1
!  print*,'pt(i,j,1) = ',pt(i,j,1)
!  print*,'pt(i,j,2) = ',pt(i,j,2)
!  print*,'sigma1(i,j) = ',sigma1(i,j)
!   print*,'omega1(i,j) = ',omega1(i,j)
!  stop
!------------------------------------------------------------------------
!  end loops over i,j
!------------------------------------------------------------------------
end do
end do
!------------------------------------------------------------------------
!
!   Radiative cooling (in timestep)
!   - cool T and the redefine hm (using old zz)
!
!------------------------------------------------------------------------
!print*,'before rad_cool tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
do i = 1,nlon
do j = 1,nlat
do k = 1,nlev
  if (k == 1) rad_cool = rad_cool_bl
  if (k == 2) rad_cool = rad_cool_lt
  if (k == 3) rad_cool = rad_cool_ut
  tt(i,j,k) = tt(i,j,k) - ((rad_cool*tstep)/sec_day)
  rtp = rv(i,j,k)
  rip = 0.
  hm(i,j,k) = hmoist(tt(i,j,k),rtp,rv(i,j,k),rip,zz(i,j,k))
!------------------------------------------------------------------------
!  Check for bad hm
!------------------------------------------------------------------------
  if (hm(i,j,k) > 410000.) then
    print*,'-------- Problem with unrealistic hm ----'
    print*,'tt(i,j,k) = ',tt(i,j,k)
    print*,'rv(i,j,k) = ',rv(i,j,k)
    print*,'zz(i,j,k) = ',zz(i,j,k)
    print*,'i j k = ',i,j,k
    print*,'hm(i,j,k) = ',hm(i,j,k)
    stop '---------- bad hm ---'
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end do
end do
end do
!print*,'after rad_cool tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
!------------------------------------------------------------------------
!
!  Find new zz and hm (in timestep)
!
!------------------------------------------------------------------------
!call height
!------------------------------------------------------------------------
!
!  Define colwat before water fluxes (in timestep)
!
!------------------------------------------------------------------------
colwat_av = 0.
dmdry_av = 0.
dmdry_bl_tot = 0.
dmdry_lt_tot = 0.
dmdry_ut_tot = 0.
do i = 1,nlon
do j = 1,nlat
  colwat = 0.
  dmdry = 0.
  do k = 1,nlev
    dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
    dmdry = dmdry + mextra(i,j,k)
    colwat = colwat + rv(i,j,k)*dmdry
    dmdry_av = dmdry_av + dmdry
    if (k == 1) dmdry_bl_tot = dmdry_bl_tot + dmdry
    if (k == 2) dmdry_lt_tot = dmdry_lt_tot + dmdry
    if (k == 3) dmdry_ut_tot = dmdry_ut_tot + dmdry
  end do
  colwat_av = colwat_av + colwat
end do
end do
colwat_av = colwat_av/(float(nlat)*float(nlon))
dmdry_av = dmdry_av/(float(nlat)*float(nlon))
colwat_start = colwat_av
dmdry_start = dmdry_av
dmdry_bl_start = dmdry_bl_tot/(float(nlat)*float(nlon))
dmdry_lt_start = dmdry_lt_tot/(float(nlat)*float(nlon))
dmdry_ut_start = dmdry_ut_tot/(float(nlat)*float(nlon))
!------------------------------------------------------------------------
!
!   Heat/Moisture fluxes (in timestep)
!   - In the new expression for hmoist below, don't bother adjusting zz
!     due to new heat in BL. OK if zz a bit out of hydrosttaic balance.
!   - but do adjust surface pressure for increase in BL water mass
!
!  lv/cpd = 2500
!  So in order for Bowen ratio to be 0.05 = cpd*dtt/lv*drv
!  need dtt/drv = 0.05*(2500) = 125
!
!------------------------------------------------------------------------
dtt_av = 0.
drv_av = 0.
watflux_av = 0.
colwat_av_prev = colwat_av
bowen_av = 0.
if (do_flux == 1) then
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!   Calculate dmdy
!   - have to use initial rv here
!------------------------------------------------------------------------
  dmdry = (ph(1) - ph(2))/(g*(1.0+rv(i,j,1))) + mextra(i,j,1)
!------------------------------------------------------------------------
! Get drv
!------------------------------------------------------------------------
  es_sst = esat(sst(i,j))
  es_t = esat(tt(i,j,1))
  if (tt(i,j,1) > 350.) then
    print*,'surface temperature i j tt(i,j,1) = ',i,j,tt(i,j,1)
    stop 'surface T too hot before call to rsat in bowen ratio calculation'
  endif
  rvs = rsat(tt(i,j,1),pp(1))
  rhp = rv(i,j,1)/rvs
  e_bl = rhp*es_t
  drv = a_evap*(es_sst - e_bl)*tstep
  dtt = a_heat*(sst(i,j) - tt(i,j,1))*tstep
  tt(i,j,1) = tt(i,j,1) + dtt
  rv(i,j,1) = rv(i,j,1) + drv
  senflux(i,j) = cpd*dtt*dmdry/tstep  ! [J/(kg*K)]*K*(kg/m2) = J/m2
  watflux(i,j) = drv*latent(tt(i,j,1))*dmdry/tstep  ! (kg water)/(kg dry air)*(J/kg water)*(kg dry air/m2) = J/m2
!  print*,'senflux(i,j) = ',senflux(i,j)
!  print*,'watflux(i,j) = ',watflux(i,j)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'---------  In flux --------'
    print*,'sst(i,j) = ',sst(i,j)
    print*,'tt(i,j,1) = ',tt(i,j,1)
    print*,'dtt = ',dtt
    print*,'drv = ',drv
    print*,'rhp = ',rhp
    print*,'---------  end flux   ---'
  endif
!------------------------------------------------------------------------
!  bowen ratio for diagnostics
!------------------------------------------------------------------------
  if (abs(drv) < 1.0E-08) stop 'drv tto small problem'
  bowen = (cpd*dtt)/(latent(tt(i,j,1))*drv)
  bowen_av = bowen_av + bowen
!------------------------------------------------------------------------
!  Find new hm from updated tt/rv
!------------------------------------------------------------------------
  rtp = rv(i,j,1)
  rip = 0.
  hm(i,j,1) = hmoist(tt(i,j,1),rtp,rv(i,j,1),rip,zz(i,j,1))
!------------------------------------------------------------------------
!  Diagnostics
!------------------------------------------------------------------------
  dtt_av = dtt_av + dtt
  drv_av = drv_av + drv
  wat_added = drv*dmdry
  watflux_av = watflux_av + wat_added
!------------------------------------------------------------------------
! Adjust mextra
! - even though drymass of layer is unchanged, have more total mass since
!   more water vapor, so must move some dry mass to drymass_extra
!------------------------------------------------------------------------
  dm_bl = dmdry*(1.0+rv(i,j,1))  ! New total mass
  dm_want = (ph(1) - ph(2))/g  ! Available total mass
  dm_extra = dm_bl - dm_want
  mextra(i,j,1) = dm_extra/(1.0 + rv(i,j,1))
!------------------------------------------------------------------------
!  end loop over i,j
!------------------------------------------------------------------------
end do
end do
endif
!print*,'after surface fluxes tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
!------------------------------------------------------------------------
!
!  Update sst
!  say watflux = 400
!  sst_mass = 2000
!  cl = 4000
!  dsst = 1
!  dsst*tstep/sst_tscale = tstep*400/(4000*2000)
!  dsst/sst_tscale = 400/(4000*2000)
!  dsst/sst_tscale = 4/(40*2000)
!  dsst/sst_tscale = 1/(10*2000)
!  1/sst_tscale = 1/(10*2000)
!  sst_tscale = 10*2000 = 20000
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  sst(i,j) = sst(i,j) + (((sst_tar(i,j) - sst(i,j))*tstep)/sst_tscale)
  sst_watflux = -(tstep*watflux(i,j))/(cl*sst_mass)
  sst_senflux= -(tstep*senflux(i,j))/(cl*sst_mass)
  sst(i,j) = sst(i,j) + sst_watflux + sst_senflux
!  print*,'sst_watflux = ',sst_watflux
!  print*,'sst_senflux = ',sst_senflux
end do
end do
!------------------------------------------------------------------------
!
!  Define colwat after water fluxes
!
!------------------------------------------------------------------------
colwat_av = 0.
do i = 1,nlon
do j = 1,nlat
  colwat = 0.
  do k = 1,nlev
    dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
    dmdry = dmdry + mextra(i,j,k)
    colwat = colwat + rv(i,j,k)*dmdry
  end do
  colwat_av = colwat_av + colwat
  if ((i==1).and.(j==1)) colwat_a = colwat
end do
end do
colwat_av = colwat_av/(float(nlat)*float(nlon))
!------------------------------------------------------------------------
!
!  Determine some tropical mean values.
!
!------------------------------------------------------------------------
dtt_av = dtt_av/float(nlon*nlat)
drv_av = drv_av/float(nlon*nlat)
ps_av = ps_av/float(nlon*nlat)
watflux_av = watflux_av/float(nlon*nlat)
bowen_av = bowen_av/float(nlat*nlon)
!------------------------------------------------------------------------
!
!  Define new height and hm
!
!------------------------------------------------------------------------
!print*,'Calling height after fluxes'
call height
!print*,'Finished calling height after fluxes'
!------------------------------------------------------------------------
!
!  Calculate cape_ut and cape_lt (inside timesteps)
!
!------------------------------------------------------------------------
!print*,'Starting cape calculation'
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  First find cape_lt
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!  Calculate parcel T,rv,rl,ri at lt height
!  - impose conservation of hm
!  - allow rlp/rip formation
!------------------------------------------------------------------------
!  print*,'input hm(i,j,1) = ',hm(i,j,1)
!  print*,'input kmm = ',kmm
  hmp = hm(i,j,1)
  rvp = rv(i,j,1)
  zp = zz(i,j,2) ! New height of parcel
  rtp = rvp
  rlp = 0. ! No background condensate in BL
  rip = 0. ! No background condensate in BL
  kmp = hmp - (1. + rtp)*g*zp ! enthalpy at new height
  call t_from_km(kmp,rtp,ppp(i,j,2),rvp,rlp,rip,tp)  ! Do not allow saturation
!  print*,'tt(i,j,2) = ',tt(i,j,2)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (rvp < 0.) then
    print*,'rvp = ',rvp
    stop 'rvp negative after lt t_from_km call in cape_lt calculation'
  endif
!------------------------------------------------------------------------
!  Calculate contribution to cape
!------------------------------------------------------------------------
  tdd = tt(i,j,2)*(1. + epsi*rv(i,j,2))/(1. + rv(i,j,2)) ! grid cell Td
!  tdp = tp*(1. + epsi*rvp)/(1. + rtp)  ! parcel Td
  tdp = tp*(1. + epsi*rvp)/(1. + rvp)  ! parcel Td  assume rainout
  dtd_lt = tdp - tdd
  dlnp = log(ppp(i,j,1)) - log(ppp(i,j,2))
  cape_lt(i,j) = rdd*0.5*dtd_lt*dlnp ! multiply by half assume B = 0 at bl
!------------------------------------------------------------------------
!  Diagnostics
!------------------------------------------------------------------------
  nprint = 1
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'---------- cape_lt calculation --- nstep = ',nstep
    print*,'nstep = ',nstep
    print*,'grid cell tdd = ',tdd
    print*,'parcel tp = ',tp
    print*,'parcel tdp = ',tdp
    print*,'kmp = ',kmp
    print*,'rvp = ',rvp
    print*,'hmp = ',hmp
    print*,'dlnp = ',dlnp
    print*,'dtd_lt = tdp - tdd = ',tdp - tdd
    print*,'cape_lt = ',cape_lt(i,j)
    print*,'---------'
  endif
!------------------------------------------------------------------------
!  Calculate parcel properties in ut
!------------------------------------------------------------------------
  kmp = hmp - (1. + rtp)*g*zz(i,j,3)
!  print*,'input hm(i,j,1) = ',hm(i,j,1)
!  print*,'input kmm = ',kmm
  call t_from_km(kmp,rtp,ppp(i,j,3),rvp,rlp,rip,tp)
!  print*,'output rvp = ',rvp
!  print*,'output tp = ',tp
!  print*,'tt(i,j,3) = ',tt(i,j,3)
!  print*,'output rlp = ',rlp
!  print*,'Found cape_lt = ',cape_lt(i,j)
!------------------------------------------------------------------------
!  Calculate cape_ut
!------------------------------------------------------------------------
  tdd = tt(i,j,3)*(1. + epsi*rv(i,j,3))/(1. + rv(i,j,3)) ! grid cell Td
  tdp = tp*(1. + epsi*rvp)/(1. + rtp)  ! parcel Td
  dtd_ut = tdp - tdd
  dtd = 0.5*(dtd_ut + dtd_lt)
  dlnp = log(ppp(i,j,2)) - log(ppp(i,j,3))
  cape_ut(i,j) = cape_lt(i,j) + rdd*dtd*dlnp ! 
  if (cape_lt(i,j) < 0.) cape_lt(i,j) = 0.
  if (cape_ut(i,j) < 0.) cape_ut(i,j) = 0.
!  print*,'grid cell tdd = ',tdd
!  print*,'parcel tdp = ',tdp
!  print*,'cape_ut = ',cape_ut(i,j)
!  stop
end do
end do
!print*,'Defined cape_ut cape_lt'
!------------------------------------------------------------------------
!
!  Define colrh (after water fluxes)
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  colwat_sat = 0.
  colwat = 0.
  do k = 1,nlev
    dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
    dmdry = dmdry + mextra(i,j,k)
    rvs = rsat(tt(i,j,k),ppp(i,j,k))
    rh(i,j,k) = rv(i,j,k)/rvs
    if (rh(i,j,k) > 15.) then
      print*,'----- high rh i j k = ',i,j,k
      print*,'ppp(i,j,k) = ',ppp(i,j,k)
      print*,'rv(i,j,k) = ',rv(i,j,k)
      print*,'rh(i,j,k) = ',rh(i,j,k)
!      stop '---------  RH out of bounds ----'
    endif
    colwat = colwat + rv(i,j,k)*dmdry
    colwat_sat = colwat_sat + rvs*dmdry
  end do
  colrh_prev(i,j) = colrh(i,j)
  colrh(i,j) = colwat/colwat_sat
!  print*,'i j colrh(i,j) = ',i,j,colrh(i,j)
end do
end do
!print*,'Defined colrh after water fluxes'
!------------------------------------------------------------------------
!
!  Define tropical mean quantities
!  - cape_lt,cape_ut,colrh,etc
!
!------------------------------------------------------------------------
zz_bl_av = 0.
zz_lt_av = 0.
zz_ut_av = 0.
cape_lt_av = 0.
cape_ut_av = 0.
colrh_av = 0.
do i = 1,nlon
do j = 1,nlat
  cape_lt_av = cape_lt_av + cape_lt(i,j)
  cape_ut_av = cape_ut_av + cape_ut(i,j)
  colrh_av = colrh_av + colrh(i,j)
  zz_bl_av = zz_bl_av + zz(i,j,1)
  zz_lt_av = zz_lt_av + zz(i,j,2)
  zz_ut_av = zz_ut_av + zz(i,j,3)
end do
end do
cape_lt_av = cape_lt_av/float(nlon*nlat)
cape_ut_av = cape_ut_av/float(nlon*nlat)
colrh_av = colrh_av/float(nlon*nlat)
zz_bl_av = zz_bl_av/float(nlon*nlat)
zz_lt_av = zz_lt_av/float(nlon*nlat)
zz_ut_av = zz_ut_av/float(nlon*nlat)
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  Drizzle (DR)          ********************************
! *********************************************************************
!
!------------------------------------------------------------------------
if (do_dr == 1) call drizzle
!------------------------------------------------------------------------
!
!  Update heights before horizontal transport
!
!------------------------------------------------------------------------
call height
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  Horizontal Transport  ********************************
! *********************************************************************
!
!   Initialize mconv
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dmdry = (ph(1) - ph(2))/(g*(1.0 + rv(i,j,1))) + mextra(i,j,1)
  rvmass_bl(i,j) = rv(i,j,1)*dmdry
  do k = 1,nlev
    mconv(i,j,k) = 0.
  end do
end do
end do
!------------------------------------------------------------------------
!
!   Trade Wind Circulation
!   - do this first since CG conv needs mconv
!
!------------------------------------------------------------------------
if (do_trade_wind == 1) call trade_wind
!------------------------------------------------------------------------
!
!   Spread Mass Horizontally
!   - do this first since CG conv needs mconv
!
!------------------------------------------------------------------------
!stop 'calling spread_mass'
if (do_spread == 1) call spread_mass
!print*,'after spread tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
!stop 'after first call spread_mass'  ! no problem to here
!------------------------------------------------------------------------
!
!  Determine  BL water vapor convergence
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dmdry = (ph(1) - ph(2))/(g*(1.0 + rv(i,j,1))) + mextra(i,j,1)
  xxx = rv(i,j,1)*dmdry
  wconv_bl(i,j) = xxx - rvmass_bl(i,j) 
end do
end do
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  After spreading timestep nstep = ',nstep
    print*,'tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'pt(i_want,j_want,1) = ',pt(i_want,j_want,1)
    print*,'pt(i_want,j_want,2) = ',pt(i_want,j_want,2)
    print*,'pt(i_want,j_want,3) = ',pt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  DP and CG Convection (DC)  ********************************
! *********************************************************************
!
!------------------------------------------------------------------------
if (do_cg == 1) call convect(2)
print*,'Finished dp call to convect starting cg call'
if (do_dp == 1) call convect(3)
!print*,'after convect tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
print*,'Finished cg call to convect'
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  AN: Anvil Rain   *************************************
! *********************************************************************
!
!  Should put in subroutine
!
!------------------------------------------------------------------------
print*,'starting anvil rain'
if (do_an == 1) then
do i = 1,nlon
do j = 1,nlat
  rvsp = rsat(tt(i,j,3),ppp(i,j,3))
  rv_avail = rv(i,j,3) - rh_an_max*rvsp
  if (rv_avail > 0.) then
    dmdry = (ph(3) - ph(4))/(g*(1.0 + rv(i,j,3))) + mextra(i,j,3)
    rv_initial = rv(i,j,3)
    xxx = rain_dp(i,j)*(sec_day/tstep)
    xxx = xxx*(tstep/tau_an)
    if (xxx > 1.) xxx = 1.0
    xxx = 1.0
    rv_final = rv(i,j,3) - xxx*rv_avail
    ri_final = 0.
    km_initial = hm(i,j,3) - (1. + rv_initial)*g*zz(i,j,3)  
    km_wat = cl*tt(i,j,3)
    km_final = km_initial - km_wat*(rv_initial - rv_final)
    call t_from_km_allow_sat(km_final,rv_final,ppp(i,j,3),t_final)
    rv(i,j,3) = rv_final
    tt(i,j,3) = t_final
    rain_an(i,j) = rain_an(i,j) + (rv_initial - rv_final)*dmdry
    hm(i,j,3) = hmoist(t_final,rv_final,rv_final,ri_final,zz(i,j,3))
    dm_want = (ph(3) - ph(4))/g  
    dm_final = (1.0 + rv_final)*dmdry  ! total mass after rain removal
    dm_extra = dm_final - dm_want
    mextra(i,j,3) = dm_extra/(1.0 + rv_final)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
    if (rv_final < 0.) then
    print*,'rv_final = ',rv_final
    print*,'rv_initial = ',rv_initial
    stop '------rv_final negative in anvil --------'
    endif
!------------------------------------------------------------------------
!------------------------------------------------------------------------
  endif
end do
end do
endif
!print*,'after anvil tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
print*,'finished anvil'
!------------------------------------------------------------------------
!
! *********************************************************************
! *************    Downdrafts (DD)  ***********************************
! *********************************************************************
!
!  Should put in subroutine
!
!------------------------------------------------------------------------
if (do_down == 1) then
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  Calculate colwat before downdrafts: colwat_bf_dn
!------------------------------------------------------------------------
    colwat_bf_dn = 0.
    do k = 1,nlev
      dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
      colwat_bf_dn = colwat_bf_dn + rv(i,j,k)*(dmdry + mextra(i,j,k))
    end do
!------------------------------------------------------------------------
!  Downdraft diagnostics
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dn == 1)) then
    print*,'******** Starting downdraft diagnostics ***********'
    print*,'Updraft mass mass_dp(i,j) = ',mass_dp(i,j)
    print*,'Initial colwat_bf_dn = ',colwat_bf_dn 
    print*,'UT mextra(i,j,3) = ',mextra(i,j,3)
    print*,'LT mextra(i,j,2) = ',mextra(i,j,2)
    print*,'BL mextra(i,j,1) = ',mextra(i,j,1)
  endif
!------------------------------------------------------------------------
!   DD
!------------------------------------------------------------------------
  mextra_lt_initial = mextra(i,j,2)
  rv_lt_initial = rv(i,j,2)
  rv_bl_initial = rv(i,j,1)
  hm_bl_initial = hm(i,j,1)
  dmdry_bl_initial = (ph(1) - ph(2))/(g*(1.0+rv_bl_initial)) + mextra(i,j,1)
  rvmass_bl_initial = rv_bl_initial*dmdry_bl_initial
  dmdry_lt_initial = (ph(2) - ph(3))/(g*(1.0+rv_lt_initial)) + mextra(i,j,2)
  rvmass_lt_initial = rv_lt_initial*dmdry_lt_initial
!------------------------------------------------------------------------
!  Determine required evaporation rv_evap into a LT test parcel to reach the BL
!  with prescribed negative buoyancy
!------------------------------------------------------------------------
  if ((rh(i,j,1) < 1.0).and.(rh(i,j,2) < 1.0)) then
!  if ((rh(i,j,1) < 2.0).and.(rh(i,j,2) < 1.0)) then
  if ((rain_an(i,j)+rain_dp(i,j)) > rain_thresh) then
    bp = 1.0
    rhp = 0.5
    rv_lt = rv(i,j,2)
    do num = 1,num_evap
    if ((bp > b_down).and.(rhp < rh_dn_max)) then
      rv_evap = float(num)*rv_evap_int
      rvp = rv_lt + rv_evap
      km_wat = cl*tt(i,j,2)
      hmp = hm(i,j,2) + rv_evap*(km_wat + g*zz(i,j,3))
      z_bl = zz(i,j,1)
      kmp = hmp - (1. + rvp)*g*zz_bl  ! km of ddraft parcel moved to BL
      call t_from_km_allow_sat(kmp,rvp,ppp(i,j,1),tp)
      tvp = tp*(1. + epsi*rvp)/(1. + rvp)
      rvsp = rsat(tp,pp(1))
      rhp = rvp/rvsp
      rvp_bl = rv(i,j,1)
      tvp_bl = tt(i,j,1)*(1. + epsi*rvp_bl)/(1. + rvp_bl)
      bp = g*(tvp-tvp_bl)/tvp_bl
!     print*,'------------- num ----------',num
!     print*,'rv_evap = ',rv_evap
!     print*,'rhp = ',rhp
!     print*,'tp = ',tp
!     print*,'bp = ',bp
!     print*,'tvp_bl = ',tvp_bl
    endif
    end do
!------------------------------------------------------------------------
!  Make initial guess of downdraft mass
!  rv_evap is like a water mass per kg.
!------------------------------------------------------------------------
!    if (bp < b_down) then
    if (bp < 0.) then
!    if (bp < 0.1) then
      xxx = (rain_an(i,j)+rain_dp(i,j))*(sec_day/tstep)
      xxx = mass_dp(i,j)
!      fevap = sigmoidal(xxx,fevap_half,fevap_scale,fevap_low,fevap_add)
      fff = sigmoidal(xxx,fevap_half,fevap_scale,fevap_low,fevap_add)
      mass_dn(i,j) = fevap_max*(rain_an(i,j)+rain_dp(i,j)-evap_dp(i,j)-down_dp(i,j))/rv_evap
!      if (mass_dn(i,j) > f_down_max*mass_dp(i,j)) mass_dn(i,j) = f_down_max*mass_dp(i,j) 
      if (mass_dn(i,j) > fff*mass_dp(i,j)) mass_dn(i,j) = fff*mass_dp(i,j) 
    else
      mass_dn(i,j) = 0.
    endif
!------------------------------------------------------------------------
! Have final rv_evap,mass_dn
! Define evap
!------------------------------------------------------------------------
    hm_down = hmp
    rv_down = rv_lt + rv_evap
    down_dp(i,j) = rv_evap*mass_dn(i,j)
    down_dp_av = down_dp_av + down_dp(i,j)
!------------------------------------------------------------------------
!   Diagnostics
!------------------------------------------------------------------------
    if ((i == i_want).and.(j == j_want).and.(nprint_dn == 1)) then
      print*,'Final bp = ',bp
      print*,'rv_down = ',rv_down
      print*,'rv_lt = ',rv_lt
      print*,'rv_evap = ',rv_evap
      print*,'mass_dn(i,j) = ',mass_dn(i,j)
      print*,'hm_down = ',hm_down
      print*,'hm_bl_initial = ',hm_bl_initial
    endif
!------------------------------------------------------------------------
!  Remove downdraft parcel from LT
!  - LT rv/hm should not change
!------------------------------------------------------------------------
    mextra(i,j,2) = mextra(i,j,2) - mass_dn(i,j)
    dmdry_lt_final = dmdry_lt_initial - mass_dn(i,j)  ! need?
    rvmass_lt_final = dmdry_lt_final*rv_lt            ! need?
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
    if (dmdry_lt_final < 0.) then
      print*,' ------- Mass problem in Downdrafts ---'
      print*,'dmdry_lt_final = ',dmdry_lt_final
      print*,'dmdry_lt_initial = ',dmdry_lt_initial
      print*,'rain_an(i,j) in mm/day = ',rain_an(i,j)*(sec_day/tstep)
      print*,'mass_dp(i,j) = ',mass_dp(i,j)
      print*,'mass_dn(i,j) = ',mass_dn(i,j)
      print*,'i = ',i
      print*,'j = ',j
      stop '----------  dmdry_lt_final less than zero in Downdrafts ---'
    endif
!------------------------------------------------------------------------
!  Determine final BL rv value
!------------------------------------------------------------------------
    rv_bl_initial = rv(i,j,1)
    dmdry_bl_final = dmdry_bl_initial + mass_dn(i,j)
    rvmass_bl_final = dmdry_bl_initial*rv_bl_initial + mass_dn(i,j)*rv_down
    rv_bl_final = rvmass_bl_final/dmdry_bl_final
    rv(i,j,1) = rv_bl_final
!------------------------------------------------------------------------
!  Check for NaN
!------------------------------------------------------------------------
    if (rv(i,j,1) /= rv(i,j,1)) then 
      print*,'------------ NaN problem in downdrafts  ------------'
      print*,'f_mass_dn_max = ',f_mass_dn_max
      print*,'mass_dn(i,j) = ',mass_dn(i,j)
      print*,'rv_evap = ',rv_evap
      stop 'rv(i,j,1) NaN problem'
    endif
!------------------------------------------------------------------------
!  Determine final BL hm value
!------------------------------------------------------------------------
    hm_bl_initial = hm(i,j,1)
    hmmass_bl = dmdry_bl_initial*hm_bl_initial + mass_dn(i,j)*hm_down
    if (abs(dmdry_bl_final) < 1.0E-06) stop 'dmdry_bl_final problem'
    hm_bl_final = hmmass_bl/dmdry_bl_final
    hm(i,j,1) = hm_bl_final
!------------------------------------------------------------------------
!  Add downdraft parcel to LT by increasing surface pressure
!  - Use rv_down here and added downdraft mass
!------------------------------------------------------------------------
    dm_have = dmdry_bl_final*(1.0+rv_bl_final)
    dm_want = (ph(1) - ph(2))/g  ! Available total mass
    mextra(i,j,1) = (dm_have - dm_want)/(1.0+rv_bl_final)
!------------------------------------------------------------------------
!  Calculate colwat after downdrafts: colwat_af_dn
!------------------------------------------------------------------------
    colwat_af_dn = 0.
    do k = 1,nlev
      dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
      colwat_af_dn = colwat_af_dn + rv(i,j,k)*(dmdry + mextra(i,j,k))
    end do
!------------------------------------------------------------------------
!
!  Check for column water error (downdrafts)
!
!------------------------------------------------------------------------
colwat_error = colwat_af_dn - colwat_bf_dn - down_dp(i,j)
  if (abs(colwat_error) > colwat_tolerance) then
    print*,'-----------  colwat error after downdrafts  ----'
    print*,'colwat_error = ',colwat_error
    print*,'colwat_bf_dn = ',colwat_bf_dn
    print*,'colwat_af_dn = ',colwat_af_dn
    print*,'down_dp(i,j) = ',down_dp(i,j)
    print*,'Final bp = ',bp
    print*,'rv_down = ',rv_down
    print*,'rv_lt = ',rv_lt
    print*,'rv_evap = ',rv_evap
    print*,'mass_dn(i,j) = ',mass_dn(i,j)
    print*,'hm_down = ',hm_down
    print*,'hm_bl_initial = ',hm_bl_initial
    print*,'Initial LT mextra_lt_initial = ',mextra_lt_initial
    print*,'Final mextra(i,j,2) = ',mextra(i,j,2)
    print*,'rv_evap = ',rv_evap
    print*,'rv_bl_initial = ',rv_bl_initial
    print*,'rv_bl_final = ',rv_bl_final
    print*,'hm_bl_final = ',hm_bl_final
    print*,'Evaporation downdraft fraction= ',down_dp(i,j)/rain_an(i,j)
    print*,'num bp rv_evap = ',num,bp,rv_evap
    print*,'mass_dp(i,j) = ',mass_dp(i,j)
    print*,'colwat_af_dn = ',colwat_af_dn 
    print*,'rvmass_lt_initial = ',rvmass_lt_initial
    print*,'rvmass_lt_final = ',rvmass_lt_final
    print*,'Change in LT rvmass = ',rvmass_lt_final-rvmass_lt_initial
    print*,'rvmass_bl_initial = ',rvmass_bl_initial
    print*,'rvmass_bl_final = ',rvmass_bl_final
    print*,'Increase in BL rvmass = ',rvmass_bl_final-rvmass_bl_initial
    print*,'down_dp(i,j) = ',down_dp(i,j)
    print*,'Change in colwat due to downdrafts = ',colwat_af_dn - colwat_bf_dn
    print*,'Equals downdraft evaporation? down_dp(i,j) = ',down_dp(i,j)
    print*,'WATER ERROR = ',colwat_af_dn - colwat_bf_dn - down_dp(i,j)
    print*,'******** Ending downdraft diagnostics ***********'
    stop '----  colwat error after downdrafts ---'
  endif
!------------------------------------------------------------------------
!  endif for minimum rain and rh_bl < 1
!------------------------------------------------------------------------
  endif
  endif
end do
end do
endif
print*,'finished downdrafts'
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  After downdrafts timestep nstep = ',nstep
    print*,'tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!  print*,'after downdrafts tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  CGE: Cumulus Congestus/Deep Evaporation  *************
! *********************************************************************
!
!------------------------------------------------------------------------
if (do_cg_evap == 1) call evap(2)
if (do_dp_evap == 1) call evap(3)
!------------------------------------------------------------------------
!
! *********************************************************************
! *************    SD: Radiative Subsidence  **************************
! *********************************************************************
!
!  Initialize
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  sub_ut_lt_prev(i,j) = sub_ut_lt(i,j)
  sub_lt_bl_prev(i,j) = sub_lt_bl(i,j)
end do
end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
mextra_ut_av = 0.
mextra_lt_av = 0.
mextra_bl_av = 0.
do i = 1,nlon
do j = 1,nlat
  mextra_ut_av = mextra_ut_av + mextra(i,j,3)
  mextra_lt_av = mextra_lt_av + mextra(i,j,2)
  mextra_bl_av = mextra_bl_av + mextra(i,j,1)
end do
end do
mextra_ut_av = mextra_ut_av/(float(nlon)*float(nlat))
mextra_lt_av = mextra_lt_av/(float(nlon)*float(nlat))
mextra_bl_av = mextra_bl_av/(float(nlon)*float(nlat))
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (do_rad_vert == 1) then
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  Subsidence from UT -> LT
!------------------------------------------------------------------------
!   xxx = mextra(i,j,3)
   xxx = mextra(i,j,3) - mextra(i,j,2)
   f_sub = sigmoidal(xxx,sub_half,sub_scale,sub_low,sub_add)
   dmdry_move = f_sub*tstep*omega2(i,j)/g
   if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
   call move_vert(i,j,3,2,dmdry_move,1)
   sub_ut_lt(i,j) = dmdry_move
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'---------- SUBSIDENCE from UT -> LT ----- nstep = ',nstep
    print*,'i j = ',i,j
    print*,'dmdry_move = ',dmdry_move
    print*,'dmdry_move = f_sub*tstep*omega2(i,j)/g'
    print*,'sigmodial fn of xxx f_sub = ',f_sub
    print*,'omega2(i,j) = ',omega2(i,j)
    print*,'sigma2(i,j) = ',sigma2(i,j)
    print*,'mextra(i,j,3) - mextra(i,j,2)  = ',xxx
    print*,'mextra(i,j,3) = ',mextra(i,j,3)
    print*,'mextra(i,j,2) = ',mextra(i,j,2)
    print*,'--------end SD  ---'
  endif
!------------------------------------------------------------------------
!  Subsidence from LT -> BL
!------------------------------------------------------------------------
!xxx = mextra(i,j,2) ! Doing this destros the MJO
   xxx = mextra(i,j,2) - mextra(i,j,1)
   f_sub = sigmoidal(xxx,sub_half,sub_scale,sub_low,sub_add)
   dmdry_move = f_sub*tstep*omega1(i,j)/g
   if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
   call move_vert(i,j,2,1,dmdry_move,2)
   sub_lt_bl(i,j) = dmdry_move
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'---------- SUBSIDENCE from LT -> BL ----- nstep = ',nstep
    print*,'dmdry_move = ',dmdry_move
    print*,'omega1(i,j) = ',omega1(i,j)
    print*,'sigma1(i,j) = ',sigma1(i,j)
    print*,'sigmodial fn of xxx f_sub = ',f_sub
    print*,'mextra(i,j,2) - mextra(i,j,1)  = ',xxx
    print*,'mextra(i,j,2) = ',mextra(i,j,2)
    print*,'mextra(i,j,1) = ',mextra(i,j,1)
    print*,'--------end SD  ---'
  endif
!------------------------------------------------------------------------
!  end loop over i,j
!------------------------------------------------------------------------
end do
end do
endif
!print*,'after rad subsidence tt(i_one,j_one,3) = ',tt(i_one,j_one,3)
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  After subsidence timestep nstep = ',nstep
    print*,'tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'pt(i_want,j_want,1) = ',pt(i_want,j_want,1)
    print*,'pt(i_want,j_want,2) = ',pt(i_want,j_want,2)
    print*,'pt(i_want,j_want,3) = ',pt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
! *********************************************************************
! *************    DN: Dynamically Induced Motion  ********************
! *********************************************************************
!
!   Must be after Spreading
!
!  Initialialize
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  mass_ut_lt_prev(i,j) = mass_ut_lt(i,j)
  mass_lt_bl_prev(i,j) = mass_lt_bl(i,j)
end do
end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (do_dyn_vert == 1) then
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  DYN Ascent/Descent fom UT -> LT
!------------------------------------------------------------------------
   xxx = sigma_scale/sigma2(i,j)
   if (xxx > 1.) xxx = 1.0
   dmdry_move = f_dyn_mass_ut*(mconv(i,j,3)-mconv(i,j,2)) ! > 0. for down
!------------------------------------------------------------------------
!------------------------------------------------------------------------
   if ((i == i_want).and.(j == j_want).and.(nprint_dyn == 1)) then
       print*,'---------------  mass_ut_lt calc ----'
       print*,'dmdry_move = ',dmdry_move
       print*,'mconv(i,j,3) = ',mconv(i,j,3)
       print*,'mconv(i,j,2) = ',mconv(i,j,2)
   endif
!------------------------------------------------------------------------
!------------------------------------------------------------------------
   if (dmdry_move > 0.) then 
     if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
     call move_vert(i,j,3,2,dmdry_move,3) ! descent
     mass_ut_lt(i,j) = dmdry_move
   else
     dmdry_move = -dmdry_move
     if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
     call move_vert(i,j,2,3,dmdry_move,4) ! ascent
     mass_ut_lt(i,j) =  -dmdry_move
   endif
!------------------------------------------------------------------------
!  DYN Ascent/Descent from LT -> BL
!------------------------------------------------------------------------
   xxx = sigma_scale/sigma1(i,j)
   if (xxx > 1.) xxx = 1.0
   dmdry_move = f_dyn_mass_lt*(mconv(i,j,2)-mconv(i,j,1))
!dmdry_move = f_dyn_mass_lt*mconv(i,j,2) ! NO MJO is do this
   if (dmdry_move > 0.) then 
     if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
     call move_vert(i,j,2,1,dmdry_move,5) ! descent
     mass_lt_bl(i,j) = dmdry_move
   else
     dmdry_move = -dmdry_move
     if (dmdry_move > dmdry_move_max) dmdry_move = dmdry_move_max
     call move_vert(i,j,1,2,dmdry_move,6) ! ascent
     mass_lt_bl(i,j) = -dmdry_move
   endif
!------------------------------------------------------------------------
!------------------------------------------------------------------------
   if ((i == i_want).and.(j == j_want).and.(nprint_dyn == 1)) then
     print*,'mass_ut_lt(i,j) = ',mass_ut_lt(i,j)
     print*,'mass_lt_bl(i,j) = ',mass_lt_bl(i,j)
   endif
!------------------------------------------------------------------------
!  end loop over i,j
!------------------------------------------------------------------------
end do
end do
endif
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  At end of timestep nstep = ',nstep
    print*,'tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'pt(i_want,j_want,1) = ',pt(i_want,j_want,1)
    print*,'pt(i_want,j_want,2) = ',pt(i_want,j_want,2)
    print*,'pt(i_want,j_want,3) = ',pt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
! *********************************************************************
! *************   START DG: Diagnostics  *****************************
! *********************************************************************
!
!  DG: Tropical mean quantities  (In timestep loop)
!
!------------------------------------------------------------------------
rain_mmday = rain_dp_av*(sec_day/tstep)
!------------------------------------------------------------------------
!
!  Define rain_tt
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  rain_tt_prev(i,j) = rain_tt(i,j)
!
  rain_tt(i,j) = rain_an(i,j) + rain_dp(i,j) + rain_dr(i,j) + rain_cg(i,j) - &
                 evap_dp(i,j) - down_dp(i,j) - evap_cg(i,j)
end do
end do
!------------------------------------------------------------------------
!
!  MMM: MJO diagnostics
!  Define prectt_smooth
!  ii = 60 -> i = 1  set prectt_smooth = 0 since need 60 on either side
!  ii = 180 -> i = 121  allow pectt_smooth /= 0 since have 60
!  Largest x1 = 240 - 1 + 59 = 298
!  smallest x1 = 1 - 120 + 59 = -60
!  Suppose i = 2, corresponds to ii = 61
!  x1 should be zero when i = 2 and ii = 61
!  x1 = ii - i - 59 = 61 - 2 - 59 = 0
!
!------------------------------------------------------------------------
do i = 1,nlon_mjo   ! nlon_mjo = 120  from ii = 60 to 240
  if (i == 1) then
    prectt_smooth(i) = 0.  ! don't want MJO at ii = 60 since not 60 spaces
  else
    prectt_smooth(i) = 0.
    yyy_smooth = 0.
    do ii = 1,nlon  ! 1 to 240
    do jj = 1,nlat
      x1 = float(ii - i - 59)
      x1 = abs(x1) ! Better since want the longitude centre
      yyy = (1. - (x1/float(ndist_mjo)))
      if (yyy > 0.) then
        prectt_smooth(i) = prectt_smooth(i) + yyy*rain_tt(ii,jj)
        yyy_smooth = yyy_smooth + yyy
      endif
    end do
    end do
    if (abs(yyy_smooth) < 1.0E-06) stop 'yy_smooth problem'
    prectt_smooth(i) = prectt_smooth(i)/yyy_smooth
  endif
end do
!stop
!------------------------------------------------------------------------
!
!  MMM
!  Find rain_big,i_mjo
!  double precision, parameter :: rain_mjo_min = 20.0
!
!------------------------------------------------------------------------
i_mjo = 0
rain_big = 0.
do i = 1,nlon_mjo
!  print*,'i prectt_smooth(i)  = ',i,prectt_smooth(i)
  if (prectt_smooth(i) > rain_big) then
    rain_big = prectt_smooth(i)
    i_mjo = i
  endif
end do
print*,'rain_big*(sec_day/tstep) = ',rain_big*(sec_day/tstep)
!------------------------------------------------------------------------
!
!  MMM
!  Find current xlon_mjo
!  xxx    lat index
!   0   i_mjo - 0.5
!   1     i_mjo
!
!------------------------------------------------------------------------
xlon_mjo_prev = xlon_mjo
xlon_mjo = bad
nnn = ndays_spinup*nsteps_day
if (nstep > nnn) then
if (i_mjo /= 0) then
  xlon_mjo = 0.
  xlon_weight = 0.
  do ii = 1,nlon  ! 1 to 240
  do jj = 1,nlat
    x1 = float(ii - i_mjo - 59)  ! = 0 if i_mjo = 1, ii = 60
    x1 = abs(x1)
    yyy = (1. - (x1/float(ndist_mjo)))
    if (yyy > 0.) then
      xlon_mjo = xlon_mjo + lon(ii)*yyy*rain_tt(ii,jj)
      xlon_weight = xlon_weight + yyy*rain_tt(ii,jj)
    endif
  end do
  end do
  if (abs(xlon_weight) < 1.0E-06) stop 'xlon_weight problem'
  xlon_mjo = xlon_mjo/xlon_weight
endif
endif
!------------------------------------------------------------------------
!
!  Find speed_mjo from two xlon
!
!------------------------------------------------------------------------
if ((xlon_mjo_prev /= bad).and.(xlon_mjo /= bad)) then
  dx_deg = (2.0*pi*re/360.)
  xxx = ((xlon_mjo - xlon_mjo_prev)*dx_deg)/(tstep)
  if (abs(xxx) < 100.) then
    num_speed_mjo = num_speed_mjo + 1           
    num_speed_mjo_av = num_speed_mjo_av + 1           
    speed_mjo_av = speed_mjo_av + xxx
    speed_mjo(num_speed_mjo) = xxx
  endif
endif
!------------------------------------------------------------------------
!
!  Check if meet criteria for an MJO
!
!------------------------------------------------------------------------
rain_big = rain_big*(sec_day/tstep)
nnn = ndays_spinup*nsteps_day
if ((rain_big > rain_mjo_min).and.(nstep > nnn)) then
  num_mjo = num_mjo + 1
else
  i_mjo = 0
endif
print*,'i_mjo num_mjo = ',i_mjo,num_mjo
!------------------------------------------------------------------------
!
!  Initialize av variables In time step loop
!
!------------------------------------------------------------------------
sst_av = 0.
rain_an_av = 0.
rain_dr_av = 0.
rain_dp_av = 0.
rain_cg_av = 0.
evap_cg_av = 0.
down_dp_av = 0.
evap_dp_av = 0.
ps_av = 0.
rh_bl_av = 0.
rh_lt_av = 0.
rh_ut_av = 0.
tt_bl_av = 0.
tt_lt_av = 0.
tt_ut_av = 0.
colwat_av = 0.
dmdry_av = 0.
dmdry_bl_tot = 0.
dmdry_lt_tot = 0.
dmdry_ut_tot = 0.
mass_cg_av = 0.
mass_dp_av = 0.
mass_dn_av = 0.
mass_ut_lt_av = 0.
mass_lt_bl_av = 0.
sub_ut_lt_av = 0.
sub_lt_bl_av = 0.
mextra_bl_av = 0.
f_dp_av = 0.
f_cg_av = 0.
f_caput_av = 0.
f_caplt_av = 0.
f_colrh_dp_av = 0.
f_colrh_cg_av = 0.
rv_bl_av = 0.
rv_ut_av = 0.
rv_lt_av = 0.
!------------------------------------------------------------------------
!  Loop  (In timestep loop)
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  colwat = 0.
  dmdry = 0.
  do k = 1,nlev
    xxx = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k))) + mextra(i,j,k)
    colwat = colwat + rv(i,j,k)*xxx
    dmdry = dmdry + xxx
    if (k == 1) dmdry_bl_tot = dmdry_bl_tot + xxx
    if (k == 2) dmdry_lt_tot = dmdry_lt_tot + xxx
    if (k == 3) dmdry_ut_tot = dmdry_ut_tot + xxx
    if (k == 1) mextra_bl_av = mextra_bl_av + mextra(i,j,1)
  end do
  colwat_av = colwat_av + colwat
  dmdry_av = dmdry_av + dmdry
  mass_cg_av = mass_cg_av + mass_cg(i,j)
  mass_dp_av = mass_dp_av + mass_dp(i,j)
  mass_dn_av = mass_dn_av + mass_dn(i,j)
  mass_ut_lt_av = mass_ut_lt_av + mass_ut_lt(i,j)
  mass_lt_bl_av = mass_lt_bl_av + mass_lt_bl(i,j)
  sub_ut_lt_av = sub_ut_lt_av + sub_ut_lt(i,j)
  sub_lt_bl_av = sub_lt_bl_av + sub_lt_bl(i,j)
  f_dp_av = f_dp_av + f_dp(i,j)
  f_cg_av = f_cg_av + f_cg(i,j)
  f_caput_av = f_caput_av + f_caput(i,j)
  f_caplt_av = f_caplt_av + f_caplt(i,j)
  f_colrh_dp_av = f_colrh_dp_av + f_colrh_dp(i,j)
  f_colrh_cg_av = f_colrh_cg_av + f_colrh_cg(i,j)
  rain_an_av = rain_an_av + rain_an(i,j)
  rain_dp_av = rain_dp_av + rain_dp(i,j)
  rain_dr_av = rain_dr_av + rain_dr(i,j)
  rain_cg_av = rain_cg_av + rain_cg(i,j)
  sst_av = sst_av + sst(i,j)
  evap_dp_av = evap_dp_av + evap_dp(i,j)
  down_dp_av = down_dp_av + down_dp(i,j)
  evap_cg_av = evap_cg_av + evap_cg(i,j)
  rh_bl_av = rh_bl_av + rh(i,j,1)
  rh_lt_av = rh_lt_av + rh(i,j,2)
  rh_ut_av = rh_ut_av + rh(i,j,3)
  tt_bl_av = tt_bl_av + tt(i,j,1)
  tt_lt_av = tt_lt_av + tt(i,j,2)
  tt_ut_av = tt_ut_av + tt(i,j,3)
  rv_bl_av = rv_bl_av + rv(i,j,1)
  rv_lt_av = rv_lt_av + rv(i,j,2)
  rv_ut_av = rv_ut_av + rv(i,j,3)
end do
end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
!print*,'After loop mass_ut_lt_av = ',mass_ut_lt_av
rh_bl_av = rh_bl_av/float(nlon*nlat)
rh_lt_av = rh_lt_av/float(nlon*nlat)
rh_ut_av = rh_ut_av/float(nlon*nlat)
rain_an_av = rain_an_av/(float(nlon*nlat))
rain_dp_av = rain_dp_av/(float(nlon*nlat))
rain_dr_av = rain_dr_av/(float(nlon*nlat))
sst_av = sst_av/(float(nlon*nlat))
rain_cg_av = rain_cg_av/(float(nlon*nlat))
evap_cg_av = evap_cg_av/(float(nlon*nlat))
evap_dp_av = evap_dp_av/(float(nlon*nlat))
down_dp_av = down_dp_av/(float(nlon*nlat))
colwat_av = colwat_av/(float(nlat)*float(nlon))
dmdry_av = dmdry_av/(float(nlat)*float(nlon))
colwat_end = colwat_av
dmdry_end = dmdry_av
dmdry_bl_end = dmdry_bl_tot/(float(nlat)*float(nlon))
dmdry_lt_end = dmdry_lt_tot/(float(nlat)*float(nlon))
dmdry_ut_end = dmdry_ut_tot/(float(nlat)*float(nlon))
mass_cg_av = mass_cg_av/(float(nlat)*float(nlon))
mass_dp_av = mass_dp_av/(float(nlat)*float(nlon))
mass_dn_av = mass_dn_av/(float(nlat)*float(nlon))
mass_ut_lt_av = mass_ut_lt_av/(float(nlat)*float(nlon))
mass_lt_bl_av = mass_lt_bl_av/(float(nlat)*float(nlon))
sub_ut_lt_av = sub_ut_lt_av/(float(nlat)*float(nlon))
sub_lt_bl_av = sub_lt_bl_av/(float(nlat)*float(nlon))
mextra_bl_av = mextra_bl_av/(float(nlat)*float(nlon))
f_dp_av = f_dp_av/float(nlon*nlat)
f_cg_av = f_cg_av/float(nlon*nlat)
tt_bl_av = tt_bl_av/float(nlon*nlat)
tt_lt_av = tt_lt_av/float(nlon*nlat)
tt_ut_av = tt_ut_av/float(nlon*nlat)
f_caput_av = f_caput_av/float(nlon*nlat)
f_caplt_av = f_caplt_av/float(nlon*nlat)
f_colrh_dp_av = f_colrh_dp_av/float(nlon*nlat)
f_colrh_cg_av = f_colrh_cg_av/float(nlon*nlat)
rv_bl_av = rv_bl_av/float(nlon*nlat)
rv_lt_av = rv_lt_av/float(nlon*nlat)
rv_ut_av = rv_ut_av/float(nlon*nlat)
!------------------------------------------------------------------------
!
!  DG: colwat_error   (In timestep loop)
!  colwat_end = colwat_start + watflux_av - rain_dp_av + evap_dp_av + down_dp_av
!
!------------------------------------------------------------------------
colwat_error = colwat_start - colwat_end + watflux_av - rain_an_av - rain_dp_av -  &
               rain_dr_av - rain_cg_av + evap_dp_av + evap_cg_av + down_dp_av
if (abs(colwat_error) > colwat_tolerance) then
  print*,'-----------  colwat error at end of timestep ------'
  print*,'colwat_error = ',colwat_error
  print*,'colwat_start = ',colwat_start
  print*,'colwat_end = ',colwat_end
  print*,'watflux_av = ',watflux_av
  print*,'rain_an_av = ',rain_an_av
  print*,'rain_dp_av = ',rain_dp_av
  print*,'rain_dr_av = ',rain_dr_av
  print*,'rain_cg_av = ',rain_cg_av
  print*,'evap_cg_av = ',evap_cg_av
  print*,'evap_dp_av = ',evap_dp_av
  print*,'down_dp_av = ',down_dp_av
  stop '-------  colwat error at end of timestep ---------------'
endif
!------------------------------------------------------------------------
!
!   Check for error in total dry mass  (In timestep loop)
!
!------------------------------------------------------------------------
drymass_error = dmdry_end - dmdry_start
if (abs(drymass_error) > dmdry_tolerance) then
  print*,'-----------  drymass error at end of timestep ------'
  print*,'drymass_error = ',drymass_error
  print*,'dmdry_start = ',dmdry_start
  print*,'dmdry_end = ',dmdry_end
  stop 'drymass error at end of timestep'
endif
!stop 'before error check' ! problem after here
!------------------------------------------------------------------------
!
!   Check for conservation of dry mass at every level
!   dmdry_bl_end = dmdry_bl_start - mass_cg_av - mass_dp_av + mass_dn_av - mass_lt_bl_av
!   dmdry_lt_end = dmdry_lt_start + mass_cg_av - mass_dn_av + mass_ut_lt_av - mass_lt_bl_av
!   dmdry_ut_end = dmdry_ut_start + mass_dp_av - mass_ut_lt_av
!
!------------------------------------------------------------------------
xxx = - mass_cg_av - mass_dp_av + mass_dn_av + mass_lt_bl_av + sub_lt_bl_av 
drymass_bl_error = dmdry_bl_start - dmdry_bl_end + xxx
!
xxx = mass_cg_av - mass_dn_av + mass_ut_lt_av - sub_lt_bl_av + sub_ut_lt_av - mass_lt_bl_av 
drymass_lt_error = dmdry_lt_start - dmdry_lt_end + xxx
!
xxx = mass_dp_av - mass_ut_lt_av - sub_ut_lt_av 
drymass_ut_error = dmdry_ut_start - dmdry_ut_end + xxx
!------------------------------------------------------------------------
!  Check errors  (In timestep loop)
!------------------------------------------------------------------------
if (abs(drymass_bl_error) > dmdry_tolerance) then
  print*,'--------- drymass_bl_error  at end of timestep ---------------'
  print*,'drymass_bl_error = ',drymass_bl_error 
  print*,'drymass_lt_error = ',drymass_lt_error 
  print*,'drymass_ut_error = ',drymass_ut_error 
  print*,'dmdry_bl_start = ',dmdry_bl_start
  print*,'dmdry_bl_end = ',dmdry_bl_end 
  print*,'dmdry_bl change = ',dmdry_bl_end-dmdry_bl_start
  print*,'mass_cg_av = ',mass_cg_av
  print*,'mass_dp_av = ',mass_dp_av
  print*,'mass_dn_av = ',mass_dn_av
  print*,'mass_lt_bl_av = ',mass_lt_bl_av
  print*,'sub_lt_bl_av = ',sub_lt_bl_av
  stop '------ drymass_bl_error  at end of timestep -----------'
endif
if (abs(drymass_lt_error) > dmdry_tolerance) then
  print*,'--------- drymass_lt_error  at end of timestep ---------------'
  print*,'drymass_lt_error = ',drymass_lt_error 
  print*,'dmdry_lt_start = ',dmdry_lt_start
  print*,'mass_cg_av = ',mass_cg_av
  print*,'mass_dp_av = ',mass_dp_av
  print*,'mass_dn_av = ',mass_dn_av
  print*,'mass_lt_bl_av = ',mass_lt_bl_av
  print*,'mass_ut_lt_av = ',mass_ut_lt_av
  print*,'dmdry_lt_end = ',dmdry_lt_end 
  stop 'drymass_lt_error  at end of timestep'
endif
if (abs(drymass_ut_error) > dmdry_tolerance) then
  print*,'--------- drymass_ut_error  at end of timestep ---------------'
  print*,'drymass_ut_error = ',drymass_ut_error 
  print*,'dmdry_ut_start = ',dmdry_ut_start
  print*,'mass_cg_av = ',mass_cg_av
  print*,'mass_dp_av = ',mass_dp_av
  print*,'mass_dn_av = ',mass_dn_av
  print*,'mass_lt_bl_av = ',mass_lt_bl_av
  print*,'mass_ut_lt_av = ',mass_ut_lt_av
  print*,'dmdry_ut_end = ',dmdry_ut_end 
  stop 'drymass_ut_error  at end of timestep'
endif
!stop 'before WATER'  ! after here
!------------------------------------------------------------------------
!
!  In timestep loop
!  DG: Write out files:
!
!  WATER   (In timestep loop)
!
!------------------------------------------------------------------------
xxx = sec_day/tstep
if (nstep == 1) open(unit = 99, file = dir_temp//'WATER')
days = float(nstep)*(tstep/sec_day)
write(99,"(8F13.6,E12.3,11F13.6)")       &
days,                       & ! 1.day
colwat_av,                  & ! 2.
watflux_av,                 & ! 3. change in column water/one timestep
colwat_end-colwat_start,    & ! 4.
rain_dp_av*xxx,             & ! 5.
evap_dp_av*xxx,             & ! 6.
rain_cg_av*xxx,             & ! 7.
evap_cg_av*xxx,             & ! 8.
colwat_error,               & ! 9.
rain_mmday,                 & ! 10.
colrh_av,                   & ! 11.
rh_bl_av,                   & ! 12.
rh_lt_av,                   & ! 13.
rh_ut_av,                   & ! 14.
tt_bl_av,                   & ! 15.
tt_lt_av,                   & ! 16.
tt_ut_av,                   & ! 17.
down_dp_av*xxx,             & ! 18.
rain_an_av*xxx,             & ! 19.
rain_dr_av*xxx                ! 20.
if (nstep == nsteps) close(99)
!------------------------------------------------------------------------
!
!  WATER.ONE   (In timestep loop)
!
!------------------------------------------------------------------------
xxx = sec_day/tstep
if (nstep == 1) open(unit = 91, file = dir_temp//'WATER.ONE')
write(91,"(33E12.4)")       &
days,                       & ! 1. day
senflux(i_one,j_one),       & ! 2  J/m2
watflux(i_one,j_one),       & ! 3. J/m2
rain_dp(i_one,j_one)*xxx,   & ! 4.
evap_dp(i_one,j_one)*xxx,   & ! 5.
down_dp(i_one,j_one)*xxx,   & ! 6.
rain_cg(i_one,j_one)*xxx,   & ! 7.
evap_cg(i_one,j_one)*xxx,   & ! 8.
rain_an(i_one,j_one)*xxx,   & ! 9.
rain_dr(i_one,j_one)*xxx,   & ! 10.
colrh(i_one,j_one),         & ! 11.
rh(i_one,j_one,1),          & ! 12.
rh(i_one,j_one,2),          & ! 13.
rh(i_one,j_one,3),          & ! 14.
tt(i_one,j_one,1),          & ! 15.
tt(i_one,j_one,2),          & ! 16.
tt(i_one,j_one,3),          & ! 17.
mass_cg(i_one,j_one),       & ! 18.
mass_dp(i_one,j_one),       & ! 19.
mass_dn(i_one,j_one),       & ! 20.
sub_ut_lt(i_one,j_one),     & ! 21.
sub_lt_bl(i_one,j_one),     & ! 22.
mass_ut_lt(i_one,j_one),    & ! 23.
mass_lt_bl(i_one,j_one),    & ! 24.
mextra(i_one,j_one,1),      & ! 25.
mextra(i_one,j_one,2),      & ! 26.
mextra(i_one,j_one,3),      & ! 27.
rv(i_one,j_one,1),          & ! 28.
rv(i_one,j_one,2),          & ! 29.
rv(i_one,j_one,3),          & ! 30.
zz(i_one,j_one,1),          & ! 31.
zz(i_one,j_one,2),          & ! 32.
zz(i_one,j_one,3)             ! 33.
if (nstep == nsteps) close(91)
!------------------------------------------------------------------------
!
!  MASS (In timestep loop)
!
!------------------------------------------------------------------------
if (nstep == 1) open(unit = 98, file = dir_temp//'MASS')
write(98,"(17F15.4)")  &
days,              &  ! 1.day
mass_cg_av,        &  ! 2.
mass_dp_av,        &  ! 3.
mass_dn_av,        &  ! 4.
sub_lt_bl_av,      &  ! 5.
sub_ut_lt_av,      &  ! 6.
mass_lt_bl_av,     &  ! 7.
mass_ut_lt_av,     &  ! 8.
mextra_bl_av,      &  ! 9.
mextra_lt_av,      &  ! 10.
mextra_ut_av,      &  ! 11.
dmdry_bl_start,    &  ! 12.
dmdry_lt_start,    &  ! 13.
dmdry_ut_start,    &  ! 14.
rv_bl_av,          &  ! 15.
rv_lt_av,          &  ! 16.
rv_ut_av              ! 17.
if (nstep == nsteps) close(98)
! stop 'before CAPE'
!------------------------------------------------------------------------
!
!  CAPE  (In timestep loop)
!
!------------------------------------------------------------------------
print*,'writing on CAPE'
if (nstep == 1) open(unit = 55, file = dir_temp//'CAPE')
write(55,"(15F15.4)")  &
days,          &  ! 1.day
bowen_av,      &  ! 2.bowen_av
!ps_av,        &  ! 3.ps_av ! something wrong with ps_av
xxx,           &  ! 3.ps_av 
f_dp_av,       &  ! 4.f_dp
f_cg_av,       &  ! 5.f_dp
f_caput_av,    &  ! 6.f_caput
f_caplt_av,    &  ! 7.f_caplt
f_colrh_dp_av, &  ! 8.f_colrh
f_colrh_cg_av, &  ! 9.f_colrh
cape_lt_av,    &  ! 10.cape_lt
cape_ut_av,    &  ! 11.cape_ut
zz_bl_av,      &  ! 12.zz_bl_av
zz_lt_av,      &  ! 13.zz_lt_av
zz_ut_av,      &  ! 14.zz_ut_av
dmdry_av          ! 15.dmdry_av
if (nstep == nsteps) close(55)
!stop 'before FLUX'
!------------------------------------------------------------------------
!
!  FLUX  (In timestep loop)
!
!------------------------------------------------------------------------
if (nstep == 1) open(unit = 96, file = dir_temp//'FLUX')
write(96,"(5F15.4)")  &
days,                     &  ! 1.day
bowen_av,                 &  ! 2.bowen_av
drv_av*sec_day,           &  ! 3.drv_av
dtt_av*(sec_day/tstep),   &  ! 4.dtt_av
sst_av-tkelvin               ! 5.
if (nstep == nsteps) close(96)

!stop 'after Flux'
!------------------------------------------------------------------------
!
!  WKK
!  In timestep loop
!
!  Calculate prect_wk:
!  This is tricky
!  want to average rain over 6 hour intervals
!  must truncate series so have complete number of days
!  Use day(n): number of days since start of simulation
!  NRECS_WK = nsteps - ndays_spinup*nsteps_day  ! exclude first spinup days
!  nsteps = 800
!  ndays_spinup = 5
!  nsteps_day = 24
!  NRECS_WK = 680
!  
!  
!------------------------------------------------------------------------
if (nstep == 1) num_wk = 0
if (nstep > (ndays_spinup*nsteps_day)) then
!stop 'entering prect_wk definition'
!------------------------------------------------------------------------
!  Update num_wk and check for error
!------------------------------------------------------------------------
  num_wk = num_wk + 1
  if (num_wk > NRECS_WK) then
    print*,'NRECS_WK incorrect'
    print*,'NRECS_WK = ',NRECS_WK
    print*,'num_wk = ',num_wk
    print*,'nstep = ',nstep
    stop 'NRECS_WK too small'
  endif
!------------------------------------------------------------------------
!  tim_wk is the time array in days
!  NRECS_WK is hard-wired to be correct (including complete days only).
!  Drop the last part of the record if involves incomplete days
!------------------------------------------------------------------------
  tim_wk(num_wk) = float(num_wk - 1)/float(nsteps_day)
  print*,'prect_wk calculation num_wk NRECS_WK tim_wk(num_wk) = ',num_wk,NRECS_WK,tim_wk(num_wk)
!------------------------------------------------------------------------
!  Loop over lats/lons
!------------------------------------------------------------------------
  do i = 1,nlon
  do j = 1,NLATS_WK
!------------------------------------------------------------------------
!  - What is this 39 number? The -16 latitude
!    starts at index 40.
!  - Defining lats_wk and lons_wk in this way is overkill, but no harm
!  - NLATS_WK = 21, nlat = 41; j = 21 is the equator, want to go from j = 11 to 31
!------------------------------------------------------------------------
    lons_wk(i) = lon(i)
    jj = j+10 ! From 11 to 31
    lats_wk(j) = lat(jj) 
    prect_wk(i,j,num_wk) = rain_tt(i,jj)*(sec_day/tstep)
!------------------------------------------------------------------------
! end loops over lat/lon
!------------------------------------------------------------------------
  end do
  end do
!------------------------------------------------------------------------
! endif for NRECS_WK
!------------------------------------------------------------------------
!stop 'after prect_wk defined'
endif  ! endif for nstep big enough
!------------------------------------------------------------------------
!
!  Archive hov (In timestep loop)
!
!------------------------------------------------------------------------
do i = 1,nlon
  xxx = 0.
  do j = 1,nlat
    if (abs(lat(j)) < lat_hov) then
    hov_15(i,nstep) = hov_15(i,nstep) + rain_tt(i,j)*(sec_day/tstep)
    hov_colrh(i,nstep) = hov_colrh(i,nstep) + colrh(i,j)
    xxx = xxx + 1.
    endif
  end do
  hov_15(i,nstep) = hov_15(i,nstep)/xxx
  hov_colrh(i,nstep) = hov_colrh(i,nstep)/xxx
end do
!------------------------------------------------------------------------
!
!  Mean quantities  (In timestep loop)
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  wdelta_bl_mn(i,j) = wdelta_bl_mn(i,j) + wdelta_bl(i,j)
  wconv_bl_mn(i,j) = wconv_bl_mn(i,j) + wconv_bl(i,j)
  rain_mn(i,j) = rain_mn(i,j) + rain_tt(i,j)*(sec_day/tstep)
  rain_an_mn(i,j) = rain_an_mn(i,j) + rain_an(i,j)*(sec_day/tstep)
  rain_dp_mn(i,j) = rain_dp_mn(i,j) + rain_dp(i,j)*(sec_day/tstep)
  rain_cg_mn(i,j) = rain_cg_mn(i,j) + rain_cg(i,j)*(sec_day/tstep)
  evap_cg_mn(i,j) = evap_cg_mn(i,j) + evap_cg(i,j)*(sec_day/tstep)
  cape_lt_mn(i,j) = cape_lt_mn(i,j) + cape_lt(i,j)
  cape_ut_mn(i,j) = cape_ut_mn(i,j) + cape_ut(i,j)
  colrh_mn(i,j) = colrh_mn(i,j) + colrh(i,j)
  t_bl_mn(i,j) = t_bl_mn(i,j) + tt(i,j,1)
  t_lt_mn(i,j) = t_lt_mn(i,j) + tt(i,j,2)
  t_ut_mn(i,j) = t_ut_mn(i,j) + tt(i,j,3)
  rh_bl_mn(i,j) = rh_bl_mn(i,j) + rh(i,j,1)
  rh_lt_mn(i,j) = rh_lt_mn(i,j) + rh(i,j,2)
  rh_ut_mn(i,j) = rh_ut_mn(i,j) + rh(i,j,3)
  sub_ut_lt_mn(i,j) = sub_ut_lt_mn(i,j) + sub_ut_lt(i,j)
  sub_lt_bl_mn(i,j) = sub_lt_bl_mn(i,j) + sub_lt_bl(i,j)
  mass_ut_lt_mn(i,j) = mass_ut_lt_mn(i,j) + mass_ut_lt(i,j)
  mass_lt_bl_mn(i,j) = mass_lt_bl_mn(i,j) + mass_lt_bl(i,j)
  mass_dp_mn(i,j) = mass_dp_mn(i,j) + mass_dp(i,j)
  mass_cg_mn(i,j) = mass_cg_mn(i,j) + mass_cg(i,j)
  mass_dn_mn(i,j) = mass_dn_mn(i,j) + mass_dn(i,j)
  mextra_ut_mn(i,j) = mextra_ut_mn(i,j) + mextra(i,j,3)
  mextra_lt_mn(i,j) = mextra_lt_mn(i,j) + mextra(i,j,2)
  mextra_bl_mn(i,j) = mextra_bl_mn(i,j) + mextra(i,j,1)
  f_dp_mn(i,j) = f_dp_mn(i,j) + f_dp(i,j)
  f_cg_mn(i,j) = f_cg_mn(i,j) + f_cg(i,j)
  f_caput_mn(i,j) = f_caput_mn(i,j) + f_caput(i,j)
  f_caplt_mn(i,j) = f_caplt_mn(i,j) + f_caplt(i,j)
  f_colrh_dp_mn(i,j) = f_colrh_dp_mn(i,j) + f_colrh_dp(i,j)
  f_colrh_cg_mn(i,j) = f_colrh_cg_mn(i,j) + f_colrh_cg(i,j)
  senflux_mn(i,j) = senflux_mn(i,j) + senflux(i,j)
  watflux_mn(i,j) = watflux_mn(i,j) + watflux(i,j)
  sst_mn(i,j) = sst_mn(i,j) + sst(i,j)
  sigma1_mn(i,j) = sigma1_mn(i,j) + sigma1(i,j)
  sigma2_mn(i,j) = sigma2_mn(i,j) + sigma2(i,j)
  do k = 1,nlev
    mconv_mn(i,j,k) = mconv_mn(i,j,k) + mconv(i,j,k)
    zz_mn(i,j,k) = zz_mn(i,j,k) + zz(i,j,k)
  end do
end do
end do
!------------------------------------------------------------------------
!
!  Add current timestep to rain_pdf
!
!------------------------------------------------------------------------
  do i = 1,nlon
  do j = 1,nlat
!-----------------------------------------------------------
!  Determine nr_want
!-----------------------------------------------------------
    nr_want = 0
    do nr = 1,num_rain_int
!------------------------------------------------------------------------------
!  Define rain1 and rain2
!------------------------------------------------------------------------------
    if (nr == 1) then
      rain1 = 0.
    else
      rain1 = 0.5*(grid_pdf(nr-1)+grid_pdf(nr))
    endif
    if (nr == num_rain_int) then
      rain2 = grid_pdf(nr)
    else
      rain2 = 0.5*(grid_pdf(nr)+grid_pdf(nr+1))
    endif
!------------------------------------------------------------------------------
!   Find nr_want
!------------------------------------------------------------------------------
    rr = rain_tt(i,j)*(sec_day/tstep)
    if ((rr >= rain1).and.(rr <= rain2)) then
      nr_want = nr
    endif
!------------------------------------------------------------------------------
!  end loop over nr
!------------------------------------------------------------------------------
    end do
!------------------------------------------------------------------------------
!  Add to rain_pdf
!------------------------------------------------------------------------------
  if (nr_want /= 0) rain_pdf(nr_want) = rain_pdf(nr_want) + 1.
!------------------------------------------------------------------------------
! end do over i,j
!------------------------------------------------------------------------------
end do
end do
!------------------------------------------------------------------------
!
!  EVV
!  - look for events
!  each event
!   - sum over ii,jj
!   - if abs(ii-i) + abs(jj-j) < 15
!
!------------------------------------------------------------------------
convert = sec_day/tstep
do i = 15,nlon-15 !  avoid points near the boundaries for simplicity
do j = 1,nlat
  if (abs(lat(j)) < 15.) then
  if (rain_tt(i,j)*convert > rain_event_min) then
!    print*,'EVV rain_tt(i,j)*yyy = ',rain_tt(i,j)*yyy 
    do ii = 1,nlon
    do jj = 1,nlat
       xxx = abs(float(i-ii))
       yyy = abs(float(j-jj))
       if ((xxx + yyy) < 8.) then
         zzz = sqrt(xxx*xxx + yyy*yyy)
         zzz = zzz*1.5*110. ! assume one grid spacing equals 1.5*110 km
!         print*,'i ii xxx = ',i,ii,xxx
!         print*,'j jj yyy = ',j,jj,yyy
!         print*,'zzz = ',zzz
         do nr = 1,num_radial
           xxx = abs(zzz - rdist(nr))
           if (xxx <= 0.5*dist_int) then
!             print*,'Adding to nr = ',nr
             num_rad(nr) = num_rad(nr) + 1.
             var_rad(1,nr) = var_rad(1,nr) + zzz
             var_rad(2,nr) = var_rad(2,nr) + rain_tt(ii,jj)*convert
             var_rad(3,nr) = var_rad(3,nr) - mass_lt_bl(ii,jj)
             var_rad(4,nr) = var_rad(4,nr) - sub_lt_bl(ii,jj) 
             var_rad(5,nr) = var_rad(5,nr) - mass_lt_bl(ii,jj) - sub_lt_bl(ii,jj) 
           endif
         end do
       endif
    end do
    end do
  endif
  endif
end do
end do
!------------------------------------------------------------------------
!
!  Move data in stuff file back one timestep
!
!------------------------------------------------------------------------
if (nstep > nsteps_stuff) then
  do i = 1,nlon
  do j = 1,nlat
  do ns = 1,nsteps_stuff-1
  do nv = 1,nvar_stuff
    stuff(i,j,ns,nv) = stuff(i,j,ns+1,nv)
  end do
  end do
  end do
  end do
endif
!------------------------------------------------------------------------
!
!  Add current timestep to last timestep of stuff
!
!------------------------------------------------------------------------
print*,'Defining stuff'
yyy = sec_day/tstep
ns = nsteps_stuff
if (nstep < nsteps_stuff) ns = nstep
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
  stuff(i,j,ns,i_prect) = rain_tt(i,j)*yyy
  stuff(i,j,ns,i_t_bl) = tt(i,j,1)
  stuff(i,j,ns,i_t_lt) = tt(i,j,2)
  stuff(i,j,ns,i_t_ut) = tt(i,j,3)
  stuff(i,j,ns,i_cape_ut) = cape_ut(i,j)
  stuff(i,j,ns,i_cape_lt) = cape_lt(i,j)
  stuff(i,j,ns,i_mass_dp) = mass_dp(i,j)
  stuff(i,j,ns,i_mass_dn) = mass_dn(i,j)
  stuff(i,j,ns,i_mass_cg) = mass_cg(i,j)
  stuff(i,j,ns,i_f_dp) = f_dp(i,j)
  stuff(i,j,ns,i_f_cg) = f_cg(i,j)
  stuff(i,j,ns,i_f_caput) = f_caput(i,j)
  stuff(i,j,ns,i_f_caplt) = f_caplt(i,j)
  stuff(i,j,ns,i_f_colrh_dp) = f_colrh_dp(i,j)
  stuff(i,j,ns,i_f_colrh_cg) = f_colrh_cg(i,j)
  stuff(i,j,ns,i_colrh) = colrh(i,j)
  stuff(i,j,ns,i_rh_bl) = rh(i,j,1)
  stuff(i,j,ns,i_rh_lt) = rh(i,j,2)
  stuff(i,j,ns,i_rh_ut) = rh(i,j,3)
  stuff(i,j,ns,i_mconv1) = mconv(i,j,1)
  stuff(i,j,ns,i_mconv2) = mconv(i,j,2)
  stuff(i,j,ns,i_mconv3) = mconv(i,j,3)
  stuff(i,j,ns,i_mass_lt_bl) = mass_lt_bl(i,j)
  stuff(i,j,ns,i_mass_ut_lt) = mass_ut_lt(i,j)
  stuff(i,j,ns,i_sub_lt_bl) = sub_lt_bl(i,j)
  stuff(i,j,ns,i_sub_ut_lt) = sub_ut_lt(i,j)
  stuff(i,j,ns,i_sum_lt_bl) = sub_lt_bl(i,j) + mass_lt_bl(i,j)
  stuff(i,j,ns,i_sum_ut_lt) = sub_ut_lt(i,j) + mass_ut_lt(i,j)
  stuff(i,j,ns,i_rain_cg) = rain_cg(i,j)*yyy
  stuff(i,j,ns,i_rain_an) = rain_an(i,j)*yyy
  stuff(i,j,ns,i_zz_bl) = zz(i,j,1)
  stuff(i,j,ns,i_zz_lt) = zz(i,j,2)
  stuff(i,j,ns,i_zz_ut) = zz(i,j,3)
  stuff(i,j,ns,i_down_dp) = down_dp(i,j)*yyy
  stuff(i,j,ns,i_evap_dp) = evap_dp(i,j)*yyy
  stuff(i,j,ns,i_mconv_diff) = mconv(i,j,2) - mconv(i,j,1)
  stuff(i,j,ns,i_mextra_ut) = mextra(i,j,3)
  stuff(i,j,ns,i_mextra_lt) = mextra(i,j,2)
  stuff(i,j,ns,i_mextra_bl) = mextra(i,j,1)
  stuff(i,j,ns,i_rain_dr) = rain_dr(i,j)*yyy
  stuff(i,j,ns,i_evap_cg) = evap_cg(i,j)*yyy
  stuff(i,j,ns,i_colrh_diff) = colrh(i,j) - colrh_prev(i,j)
  stuff(i,j,ns,i_rain_tt_diff) = (rain_tt(i,j) - rain_tt_prev(i,j))*yyy
  stuff(i,j,ns,i_mass_ut_lt_diff) = mass_ut_lt(i,j) - mass_ut_lt_prev(i,j)
  stuff(i,j,ns,i_mass_lt_bl_diff) = mass_lt_bl(i,j) - mass_lt_bl_prev(i,j)
  stuff(i,j,ns,i_sub_ut_lt_diff) = sub_ut_lt(i,j) - sub_ut_lt_prev(i,j)
  stuff(i,j,ns,i_sub_lt_bl_diff) = sub_lt_bl(i,j) - sub_lt_bl_prev(i,j)
  stuff(i,j,ns,i_rv_bl) = rv(i,j,1)
  stuff(i,j,ns,i_dp_prod) = dp_prod(i,j)
  stuff(i,j,ns,i_mextra_diff) = mextra(i,j,2) - mextra(i,j,1)
  stuff(i,j,ns,i_mextra_all) = mextra(i,j,3) + mextra(i,j,2) + mextra(i,j,1)
  stuff(i,j,ns,i_mconv_all) = mconv(i,j,3) + mconv(i,j,2) + mconv(i,j,1)
!------------------------------------------------------------------------------
! end do over i,j
!------------------------------------------------------------------------------
end do
end do
!------------------------------------------------------------------------
!
!    Look for events at time midpoint of stuff, and add to events (in timestep loop)
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  if (nstep > nsteps_events) then
  if (stuff(i,j,ns_mid,i_prect) > rain_event_min) then
    num_events = num_events + 1
    do ns = 1,nsteps_stuff
    do nv = 1,nvar_stuff
      events(ns,nv) = events(ns,nv) + stuff(i,j,ns,nv)  
    end do
    end do
  endif
  endif
end do
end do
!------------------------------------------------------------------------
!
!  MMM
!  - nlon_mjo = 121
!  - nlon = 240
!  - add current timestep of stuff to mjo_mean
!  - i_mjo can vary between 1 and 120 (ii = 61 to 180)
!  If i_mjo = 1, then actual center of mjo is ii = 61
!  want ii to vary from 1 to 121
!  i_mjo   actual ii index    start ii    end ii
!    2          61              1         121
!   121        180            120         240
!------------------------------------------------------------------------
if (i_mjo /= 0) then
  do i = 1,nlon_mjo
  do j = 1,nlat
  do nv = 1,nvar_stuff
     ii = i + i_mjo - 2
     mjo_mean(i,j,nv) = mjo_mean(i,j,nv) + stuff(ii,j,nsteps_stuff,nv)
  end do
  end do
  end do
endif
!------------------------------------------------------------------------
!
!  **********  End loop over time steps  *******************8
!
!------------------------------------------------------------------------
end do
!stop 'end loop over timesteps'  ! problem after here
print*,'end loop over timesteps'  
!------------------------------------------------------------------------
!
!  EVV
!  Normalize and write out
!
!------------------------------------------------------------------------
open(unit = 10, file = dir_temp//'RADIAL')
  do nr = 1,num_radial
    do nv = 1,num_var
      var_rad(nv,nr) = var_rad(nv,nr)/num_rad(nr)
    end do
    write(10,"(6F14.8)") rdist(nr),(var_rad(nv,nr),nv=1,num_var)
  end do
close(10)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
do nnn = 1,num_speed_mjo
do num = 1,num_pdf
  speed1 = speed_pdf_grid(num) - 0.5*pdf_int
  speed2 = speed_pdf_grid(num) + 0.5*pdf_int
  if ((speed_mjo(nnn) > speed1).and.(speed_mjo(nnn) < speed2)) then
    speed_pdf(num) = speed_pdf(num) + 1.
    speed_pdf_tot = speed_pdf_tot + 1.
  endif
end do
end do
print*,'=========================================================='
if (num_speed_mjo_av > 0) then
  print*,'AVERAGE MJO speed = ',speed_mjo_av/float(num_speed_mjo_av)
else
  print*,'AVERAGE MJO speed undefined'
endif
!------------------------------------------------------------------------
!
!  Normalize speed_pdf and write it out
!
!------------------------------------------------------------------------
do num = 1,num_pdf
  if (speed_pdf_tot > 0.5) then
    speed_pdf(num) = speed_pdf(num)/speed_pdf_tot
  else
    speed_pdf(num) = bad
  endif
end do
open(unit = 10, file = dir_temp//'SPEED_PDF')
do num = 1,num_pdf
  write(10,"(2F14.8)") speed_pdf_grid(num),speed_pdf(num)
end do
close(10)
!------------------------------------------------------------------------
!
!  Normalize mjo_mean
!
!------------------------------------------------------------------------
if (num_mjo > 0) then
  do i = 1,nlon_mjo
  do j = 1,nlat
  do nv = 1,nvar_stuff
     mjo_mean(i,j,nv) = mjo_mean(i,j,nv)/float(num_mjo)
  end do
  end do
  end do
else
  print*,'Can Not define mjo_mean num_mjo = 0'
endif
!------------------------------------------------------------------------
!
!  Use mjo_men and friction to get balanced wind.
!
!------------------------------------------------------------------------
if (num_mjo > 0) then
  print*,'Calling friction'
  call friction
else
  print*,'Can not call friction since num_mjo = 0'
endif
!------------------------------------------------------------------------
!
!  Determine Center of Mass from mjo_mean
!  - For variables that change sign, and sum to zero over larger areas,
!    the longitude expectation value will be weird, since xxx small.
!
!------------------------------------------------------------------------
print*,'Start cofm num_mjo = ',num_mjo
if (num_mjo > 0) then
print*,'Defining cofm'
do nv = 1,nvar_stuff
  cofm(nv) = 0.
  xxx = 0.
!  print*,'starting nv = ',nv
  do i = 1,nlon_mjo
  do j = 1,nlat
!  mconv3, mass_lt_bl, mass_ut_lt, sum_lt_bl, sum_ut_lt
    if (abs(lat(j)) < cofm_lat) then
    if ((nv == 22).or.(nv == 23).or.(nv == 24).or.(nv == 34).or.(nv == 35)) then
      if (mjo_mean(i,j,nv) > 0.) then
        yyy = 0.
      else
        yyy = -mjo_mean(i,j,nv)
      endif
      cofm(nv) = cofm(nv) + float(i)*yyy
      xxx = xxx + yyy
    elseif (nv == 27) then  ! mass_dp growth term
      zzz = -mjo_mean(i,j,34)*mjo_mean(i,j,10)
      if (zzz < 0.) zzz = 0.
      cofm(nv) = cofm(nv) + float(i)*zzz
      xxx = xxx + zzz
    else
      cofm(nv) = cofm(nv) + float(i)*mjo_mean(i,j,nv)
      xxx = xxx + mjo_mean(i,j,nv)
    endif
    endif
!    if (nv == 22) print*,'nv i j mjo_mean(i,j,nv) = ',nv,i,j,mjo_mean(i,j,nv)
!    if (nv == 22) print*,'xxx cofm(nv) = ',xxx,cofm(nv)
  end do
  end do
  if (abs(xxx) > 0.001) then
    cofm(nv) = cofm(nv)/xxx
  else
    cofm(nv) = bad
  endif
end do
else
  print*,'Can not define cofm since mjo_mean undefined'
endif
print*,'finished cofm'
!------------------------------------------------------------------------
!
!  Write it out
!
!------------------------------------------------------------------------
print*,'opening COFM'
open(unit = 10, file = dir_temp//'COFM')
do nv = 1,nvar_stuff
!  write(10,"(3F12.2,I4,A12)") 1.5*cofm(nv),lat(21),1.5*(cofm(nv)-cofm(1)),nv,name_var(nv)
  write(10,"(3F12.2,I4,A12)") 1.5*cofm(nv),lat(21),1.5*(cofm(nv)-cofm(7)),nv,name_var(nv) ! wrt mass_dp
end do
close(10)
print*,'WROTE OUT COFM'
!------------------------------------------------------------------------
!
!  Determine divergent wind
!
!------------------------------------------------------------------------
!call massflux
!------------------------------------------------------------------------
!
!  Define mjo_mean_lon
!
!------------------------------------------------------------------------
do i = 1,nlon_mjo
do nv = 1,nvar_stuff
  mjo_mean_lon_15S_15N(i,nv) = 0.
  mjo_mean_lon_10S_10N(i,nv) = 0.
  mjo_mean_lon_5S_5N(i,nv) = 0.
  num_15 = 0
  num_10 = 0
  num_5 = 0
  do j = 1,nlat
  if (abs(lat(j)) < 15.) then
     mjo_mean_lon_15S_15N(i,nv) = mjo_mean_lon_15S_15N(i,nv) + mjo_mean(i,j,nv)
     num_15 = num_15 + 1
  endif
  if (abs(lat(j)) < 10.) then
     mjo_mean_lon_10S_10N(i,nv) = mjo_mean_lon_10S_10N(i,nv) + mjo_mean(i,j,nv)
     num_10 = num_10 + 1
  endif
  if (abs(lat(j)) < 5.) then
     mjo_mean_lon_5S_5N(i,nv) = mjo_mean_lon_5S_5N(i,nv) + mjo_mean(i,j,nv)
     num_5 = num_5 + 1
  endif
  end do
  mjo_mean_lon_15S_15N(i,nv) = mjo_mean_lon_15S_15N(i,nv)/float(num_15)
  mjo_mean_lon_10S_10N(i,nv) = mjo_mean_lon_10S_10N(i,nv)/float(num_10)
  mjo_mean_lon_5S_5N(i,nv) = mjo_mean_lon_5S_5N(i,nv)/float(num_5)
end do
end do
!------------------------------------------------------------------------
!
!  i_lead for 15S-15N
!  - assume longitude index with largest mass_dp is leading edge
!  - v = integral(dp_prod*dx)/mass_dp(max)
!
!------------------------------------------------------------------------
mass_dp_max = 0.
i_lead = 0
do i = 1,nlon_mjo
  if (mjo_mean_lon_15S_15N(i,i_mass_dp) > mass_dp_max) then
    mass_dp_max = mjo_mean_lon_15S_15N(i,i_mass_dp)
    i_lead = i
  endif
end do
mjo_sp = 0.
!do i = i_lead+1,nlon_mjo
do i = i_lead,nlon_mjo
mjo_sp = mjo_sp + mjo_mean_lon_15S_15N(i,i_dp_prod)*dx ! don't use value at i_lead-1 since negative
end do
mjo_sp = mjo_sp/(mass_dp_max*3600.)
print*,'--------  ESTIMATING mjo_sp 15S-15N ---------------'
print*,'i_lead = ',i_lead
print*,'lon(i_lead) = ',lon(i_lead)
print*,'mass_dp_max = ',mass_dp_max
print*,'mjo_sp = ',mjo_sp
print*,'-------------------------------------------'
!------------------------------------------------------------------------
!
!  i_lead for 15S-15N
!  - assume longitude index with largest mass_dp is leading edge
!  - v = integral(dp_prod*dx)/mass_dp(max)
!
!------------------------------------------------------------------------
mass_dp_max = 0.
i_lead = 0
do i = 1,nlon_mjo
  if (mjo_mean_lon_10S_10N(i,i_mass_dp) > mass_dp_max) then
    mass_dp_max = mjo_mean_lon_10S_10N(i,i_mass_dp)
    i_lead = i
  endif
end do
mjo_sp = 0.
!do i = i_lead+1,nlon_mjo
do i = i_lead,nlon_mjo
mjo_sp = mjo_sp + mjo_mean_lon_10S_10N(i,i_dp_prod)*dx ! don't use value at i_lead-1 since negative
end do
mjo_sp = mjo_sp/(mass_dp_max*3600.)
print*,'--------  ESTIMATING mjo_sp 10S-10N ---------------'
print*,'i_lead = ',i_lead
print*,'lon(i_lead) = ',lon(i_lead)
print*,'mass_dp_max = ',mass_dp_max
print*,'mjo_sp = ',mjo_sp
print*,'-------------------------------------------'
!------------------------------------------------------------------------
!
!  i_lead for 5S-5N
!  - assume longitude index with largest mass_dp is leading edge
!  - v = integral(dp_prod*dx)/mass_dp(max)
!
!------------------------------------------------------------------------
mass_dp_max = 0.
i_lead = 0
do i = 1,nlon_mjo
  if (mjo_mean_lon_5S_5N(i,i_mass_dp) > mass_dp_max) then
    mass_dp_max = mjo_mean_lon_5S_5N(i,i_mass_dp)
    i_lead = i
  endif
end do
mjo_sp = 0.
do i = i_lead+1,nlon_mjo
  mjo_sp = mjo_sp + mjo_mean_lon_5S_5N(i,i_dp_prod)*dx
end do
mjo_sp = mjo_sp/(mass_dp_max*3600.)
print*,'--------  ESTIMATING mjo_sp 5S-5N ------------'
print*,'i_lead = ',i_lead
print*,'lon(i_lead) = ',lon(i_lead)
print*,'mass_dp_max = ',mass_dp_max
print*,'mjo_sp = ',mjo_sp
print*,'-------------------------------------------'
!------------------------------------------------------------------------
!
!   Normalize rain pdf (after timesteps)
!
!------------------------------------------------------------------------
ntot = 0
do nr = 1,num_rain_int
  ntot = ntot + nint(rain_pdf(nr))
end do
do nr = 1,num_rain_int
  if (ntot == 0) stop 'ntot is zero'
  rain_pdf(nr) = rain_pdf(nr)/float(ntot)
  if (rain_pdf(nr) < 0.00001) rain_pdf(nr) = 0.000001
end do
open(unit = 10, file = dir_temp//'RAIN_PDF')
do nr = 1,num_rain_int
  write(10,"(2F14.8)") grid_pdf(nr),rain_pdf(nr)
end do
close(10)
!------------------------------------------------------------------------
!
!  Write out HOV file
!
!------------------------------------------------------------------------
day_max = float(nsteps)*(tstep/sec_day)
open(unit = 10, file = dir_temp//'HOV.RAIN.15')
write(10,*) nsteps
write(10,*) nlon
write(10,*) day_max
do nstep = 1,nsteps
  write(10,"(240F8.2)") (hov_15(i,nstep),i=1,nlon)
end do
close(10)
!
open(unit = 10, file = dir_temp//'HOV.COLRH.15')
write(10,*) nsteps
write(10,*) nlon
write(10,*) day_max
do nstep = 1,nsteps
  write(10,"(240F8.2)") (hov_colrh(i,nstep),i=1,nlon)
end do
close(10)
!------------------------------------------------------------------------
!
!  Normalize Mean quantities
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  wdelta_bl_mn(i,j) = wdelta_bl_mn(i,j)/float(nsteps)
  wconv_bl_mn(i,j) = wconv_bl_mn(i,j)/float(nsteps)
  rain_mn(i,j) = rain_mn(i,j)/float(nsteps)
  if (rain_mn(i,j) < 0.000001) rain_mn(i,j) = 0.000001
  rain_an_mn(i,j) = rain_an_mn(i,j)/float(nsteps)
  rain_dp_mn(i,j) = rain_dp_mn(i,j)/float(nsteps)
  rain_cg_mn(i,j) = rain_cg_mn(i,j)/float(nsteps)
  evap_cg_mn(i,j) = evap_cg_mn(i,j)/float(nsteps)
  cape_lt_mn(i,j) = cape_lt_mn(i,j)/float(nsteps)
  cape_ut_mn(i,j) = cape_ut_mn(i,j)/float(nsteps)
  colrh_mn(i,j) = colrh_mn(i,j)/float(nsteps)
  t_bl_mn(i,j) = t_bl_mn(i,j)/float(nsteps) 
  t_lt_mn(i,j) = t_lt_mn(i,j)/float(nsteps)
  t_ut_mn(i,j) = t_ut_mn(i,j)/float(nsteps)
  rh_bl_mn(i,j) = rh_bl_mn(i,j)/float(nsteps)
  rh_lt_mn(i,j) = rh_lt_mn(i,j)/float(nsteps)
  rh_ut_mn(i,j) = rh_ut_mn(i,j)/float(nsteps)
  f_colrh_dp_mn(i,j) = f_colrh_dp_mn(i,j)/float(nsteps)
  f_colrh_cg_mn(i,j) = f_colrh_cg_mn(i,j)/float(nsteps)
  f_caput_mn(i,j) = f_caput_mn(i,j)/float(nsteps)
  f_caplt_mn(i,j) = f_caplt_mn(i,j)/float(nsteps)
  senflux_mn(i,j) = senflux_mn(i,j)/float(nsteps)
  watflux_mn(i,j) = watflux_mn(i,j)/float(nsteps)
  sst_mn(i,j) = sst_mn(i,j)/float(nsteps)
  sigma1_mn(i,j) = sigma1_mn(i,j)/float(nsteps)
  sigma2_mn(i,j) = sigma2_mn(i,j)/float(nsteps)
  f_dp_mn(i,j) = f_dp_mn(i,j)/float(nsteps)
  f_cg_mn(i,j) = f_cg_mn(i,j)/float(nsteps)
  sub_ut_lt_mn(i,j) = sub_ut_lt_mn(i,j)/float(nsteps)
  sub_lt_bl_mn(i,j) = sub_lt_bl_mn(i,j)/float(nsteps)
  mass_ut_lt_mn(i,j) = mass_ut_lt_mn(i,j)/float(nsteps)
  mass_lt_bl_mn(i,j) = mass_lt_bl_mn(i,j)/float(nsteps)
  mass_dp_mn(i,j) = mass_dp_mn(i,j)/float(nsteps)
  mass_cg_mn(i,j) = mass_cg_mn(i,j)/float(nsteps)
  mass_dn_mn(i,j) = mass_dn_mn(i,j)/float(nsteps)
  mextra_ut_mn(i,j) = mextra_ut_mn(i,j)/float(nsteps)
  mextra_lt_mn(i,j) = mextra_lt_mn(i,j)/float(nsteps)
  mextra_bl_mn(i,j) = mextra_bl_mn(i,j)/float(nsteps)
  do k = 1,nlev
    mconv_mn(i,j,k) = mconv_mn(i,j,k)/float(nsteps)
    zz_mn(i,j,k) = zz_mn(i,j,k)/float(nsteps)
  end do
end do
end do
!stop 'before writing out MEANS' ! problem after here
!------------------------------------------------------------------------
!
!  MMM
!  Write out MJO MEANS
!
!------------------------------------------------------------------------
print*,'Writing out MJO Mean quantities'
open(unit = 10, file = dir_temp//'MJO')
write(10,*) nlon_mjo
write(10,*) nlat
do i = 1,nlon_mjo
  write(10,*) lon_mjo(i)
end do
write(10,*)
do j = nlat,1,-1
  write(10,*) lat(j)
end do
write(10,*)
do nv = 1,nvar_stuff
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,nv),i=1,nlon_mjo)
end do
end do
!
!  BL Winds
!
do j = nlat,1,-1
  write(10,"(121E12.4)") (uu(i,j,1),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (vv(i,j,1),i=1,nlon_mjo)
end do
!
!  LT Winds
!
do j = nlat,1,-1
  write(10,"(121E12.4)") (uu(i,j,2),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (vv(i,j,2),i=1,nlon_mjo)
end do
!
!  BL Winds
!
do j = nlat,1,-1
  write(10,"(121E12.4)") (uu(i,j,3),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (vv(i,j,3),i=1,nlon_mjo)
end do
close(10)
!------------------------------------------------------------------------
!
!  WINDS
!
!------------------------------------------------------------------------
print*,'writing out WINDS'
open(unit = 10, file = dir_temp//'WINDS')
!do i = 1,nlon_mjo
do i = 1,nlon_mjo,3  ! output only every third longitude value
do j = nlat,1,-1
  write(10,"(8F12.4)")    &
    lon(i),lat(j),        &
    uu(i,j,1),vv(i,j,1),  &
    uu(i,j,2),vv(i,j,2),  &
    uu(i,j,3),vv(i,j,3)
end do
end do
close(10)
!------------------------------------------------------------------------
!
! ROSSBY
!
!------------------------------------------------------------------------
print*,'writing out ROSSBY'
open(unit = 10, file = dir_temp//'ROSSBY')
do j = nlat,1,-1
  write(10,"(4F12.4)")        &
     lat(j),                  &
     0.001*ros_radius(j,1),   &
     0.001*ros_radius(j,2),   &
     0.001*ros_radius(j,3)
end do
close(10)
!------------------------------------------------------------------------
!
!  MMM
!  Write out 4 select variables for a 4 panel figure MJO_WANT
!
!integer, parameter :: i_prect = 1
!integer, parameter :: i_rain_cg = 26
!integer, parameter :: i_mconv_diff = 31
!integer, parameter :: i_sum_lt_bl = 34
!
!------------------------------------------------------------------------
print*,'Writing out MJO Mean quantities'
open(unit = 10, file = dir_temp//'MJO_WANT')
write(10,*) nlon_mjo
write(10,*) nlat
do i = 1,nlon_mjo
  write(10,*) lon_mjo(i)
end do
write(10,*)
do j = nlat,1,-1
  write(10,*) lat(j)
end do
write(10,*)
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,i_prect),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,i_mass_dp),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,i_mconv2),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (-mjo_mean(i,j,i_mconv_diff),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (-mjo_mean(i,j,i_mass_lt_bl),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (-mjo_mean(i,j,i_sum_lt_bl),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,i_rain_cg),i=1,nlon_mjo)
end do
close(10)
!------------------------------------------------------------------------
!
!  MMM
!  Write out 4 select variables for a 4 panel figure MJO_PS
!
!------------------------------------------------------------------------
print*,'Writing out MJO Mean quantities'
open(unit = 10, file = dir_temp//'MJO_PS')
write(10,*) nlon_mjo
write(10,*) nlat
do i = 1,nlon_mjo
  write(10,*) lon_mjo(i)
end do
write(10,*)
do j = nlat,1,-1
  write(10,*) lat(j)
end do
write(10,*)
do j = nlat,1,-1
  write(10,"(121E12.4)") (mjo_mean(i,j,i_mextra_all),i=1,nlon_mjo)
end do
do j = nlat,1,-1
  write(10,"(121E12.4)") (uu(i,j,1),i=1,nlon_mjo)
end do
close(10)
!------------------------------------------------------------------------
!
!  MMM
!  Write out MJO_LON_15S_15N MEANS
!
!------------------------------------------------------------------------
print*,'Writing out MJO_LON Mean quantities'
open(unit = 10, file = dir_temp//'MJO_LON_15S_15N_20')
do i = 1,nlon_mjo
  write(10,"(21E12.4)") lon_mjo(i),(mjo_mean_lon_15S_15N(i,nv),nv=1,20)
end do
close(10)
!
open(unit = 10, file = dir_temp//'MJO_LON_15S_15N_40')
do i = 1,nlon_mjo
  write(10,"(21E12.4)") lon_mjo(i),(mjo_mean_lon_15S_15N(i,nv),nv=21,40)
end do
close(10)
!
open(unit = 10, file = dir_temp//'MJO_LON_15S_15N_60')
do i = 1,nlon_mjo
  write(10,"(13E12.4)") lon_mjo(i),(mjo_mean_lon_15S_15N(i,nv),nv=41,nvar_stuff)
end do
close(10)
!------------------------------------------------------------------------
!
!  MMM
!  Write out MJO_LON_5S_5N MEANS
!
!------------------------------------------------------------------------
print*,'Writing out MJO_LON Mean quantities'
open(unit = 10, file = dir_temp//'MJO_LON_5S_5N_20')
do i = 1,nlon_mjo
  write(10,"(21E12.4)") lon_mjo(i),(mjo_mean_lon_5S_5N(i,nv),nv=1,20)
end do
close(10)
!
open(unit = 10, file = dir_temp//'MJO_LON_5S_5N_40')
do i = 1,nlon_mjo
  write(10,"(21E12.4)") lon_mjo(i),(mjo_mean_lon_5S_5N(i,nv),nv=21,40)
end do
close(10)
!
open(unit = 10, file = dir_temp//'MJO_LON_5S_5N_60')
do i = 1,nlon_mjo
  write(10,"(13E12.4)") lon_mjo(i),(mjo_mean_lon_5S_5N(i,nv),nv=41,nvar_stuff)
end do
close(10)
!------------------------------------------------------------------------
!
!  Write out 2 Weights
!
!------------------------------------------------------------------------
open(unit = 10, file = dir_temp//'WEIGHTS')
write(10,*) nlon
write(10,*) nlat
do i = 1,nlon
  write(10,*) lon(i)
end do
write(10,*)
do j = nlat,1,-1
  write(10,*) lat(j)
end do
write(10,*)
do j = nlat,1,-1
  write(10,"(240E12.4)") (weight_save(i,j),i=1,nlon)
end do
close(10)
print*,'wrote out WEIGHTS'
!------------------------------------------------------------------------
!
!  Write out Mean quantities
!
!------------------------------------------------------------------------
print*,'Writing out Mean quantities'
open(unit = 10, file = dir_temp//'MEANS')
write(10,*) nlon
write(10,*) nlat
do i = 1,nlon
  write(10,*) lon(i)
end do
write(10,*)
do j = nlat,1,-1
  write(10,*) lat(j)
end do
write(10,*)
do j = nlat,1,-1
  write(10,"(240E12.4)") (rain_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (dlog10(rain_mn(i,j)),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rain_dp_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rain_an_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rain_cg_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (evap_cg_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (t_bl_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (t_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (t_ut_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rh_bl_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rh_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (rh_ut_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_dp_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_colrh_dp_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (colrh_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_caput_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (cape_ut_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_cg_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_colrh_cg_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (f_caplt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (cape_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (sub_ut_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (sub_lt_bl_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mass_ut_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mass_lt_bl_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mass_dp_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mass_dn_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mass_cg_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mextra_ut_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mextra_lt_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mextra_bl_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mconv_mn(i,j,1),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mconv_mn(i,j,2),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mconv_mn(i,j,3),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (mconv_mn(i,j,1)-mconv_mn(i,j,2),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (wconv_bl_mn(i,j)*xxx,i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (wdelta_bl_mn(i,j)*xxx,i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (dlog10(sigma1_mn(i,j)),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (dlog10(sigma2_mn(i,j)),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (zz_mn(i,j,1),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (zz_mn(i,j,2),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (zz_mn(i,j,3),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (watflux_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (senflux_mn(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (sst_mn(i,j)-tkelvin,i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (weight_save_upper(i,j),i=1,nlon)
end do
do j = nlat,1,-1
  write(10,"(240E12.4)") (weight_save_lower(i,j),i=1,nlon)
end do
close(10)
print*,'wrote out MEANS'
!stop 'wrote out MEANS' ! before div by zero without -O3
!------------------------------------------------------------------------
!
!  Normalize events
!
!------------------------------------------------------------------------
!stop 'before defining events'  ! before div by zero without -O3
print*,'num_events = ',num_events
do ns = 1,nsteps_stuff
do nv = 1,nvar_stuff
   if (num_events > 0) then
     events(ns,nv) = events(ns,nv)/float(num_events)
   else
     events(ns,nv) = 0.
   endif
end do
end do
!stop 'defined events' ! before div by zero
!------------------------------------------------------------------------
!
!  Write out events
!
!------------------------------------------------------------------------
open(unit = 10, file = dir_temp//'EVENTS1')
do ns = 1,nsteps_stuff
  write(10,"(20E16.6)") float(ns-1)*1.0-12.0,(events(ns,nv),nv=1,19)
end do
close(10)
!
open(unit = 10, file = dir_temp//'EVENTS2')
do ns = 1,nsteps_stuff
  write(10,"(20E16.6)") float(ns-1)*1.0-12.0,(events(ns,nv),nv=20,38)
end do
close(10)
!
open(unit = 10, file = dir_temp//'EVENTS3')
do ns = 1,nsteps_stuff
  write(10,"(15E16.6)") float(ns-1)*1.0-12.0,(events(ns,nv),nv=39,nvar_stuff)
end do
close(10)

print*,'wrote out EVENTS'
print*,'NRECS_WK = ',NRECS_WK

!------------------------------------------------------------------------
!
!  Define div
!
!------------------------------------------------------------------------
do nd = 1,num_div
  pp_div(nd) = 1000. - float(nd-1)*50.  ! need num_div = 18
end do
do ns = 1,nsteps_stuff
  mm_ut = ((ph(3) - ph(4))/g) + events(ns,i_mextra_ut)   ! kg/m2
  mm_lt = ((ph(2) - ph(3))/g) + events(ns,i_mextra_lt)   ! kg/m2
  mm_bl = ((ph(1) - ph(2))/g) + events(ns,i_mextra_bl)   ! kg/m2
  div_ut = -24.*events(ns,i_mconv3)/mm_ut
  div_lt = -24.*events(ns,i_mconv2)/mm_lt
  div_bl = -24.*events(ns,i_mconv1)/mm_bl
  xxx = events(ns,i_mextra_ut) + events(ns,i_mextra_lt) + events(ns,i_mextra_bl)
  pph_1 = ph(1) + g*xxx   ! surface pressure
  xxx = events(ns,i_mextra_ut) + events(ns,i_mextra_lt) 
  pph_2 = ph(2) + g*xxx  ! pressure at top of BL
  xxx = events(ns,i_mextra_ut) 
  pph_3 = ph(3) + g*xxx  ! pressure at top of LT
  pph_4 = ph(4)
  pph_1 = 0.01*pph_1
  pph_2 = 0.01*pph_2
  pph_3 = 0.01*pph_3
  pph_4 = 0.01*pph_4
  do nd = 1,num_div
    if (pp_div(nd) > pph_2) then
      div(ns,nd) = div_bl
    elseif (pp_div(nd) > pph_3) then
      div(ns,nd) = div_lt
    else
      div(ns,nd) = div_ut
    endif
  end do
end do
!------------------------------------------------------------------------
!
!  Write out div
!
!------------------------------------------------------------------------
open(unit = 10, file = dir_temp//'EVENTS_DIV')
write(10,*) nsteps_stuff
write(10,*) num_div
do ns = 1,nsteps_stuff
  write(10,*) float(ns-1)*1.0-12.0
end do
write(10,*)
do nd = num_div,1,-1
  write(10,*) pp_div(nd)
end do
write(10,*)
do nd = num_div,1,-1
  write(10,"(25E12.4)") (div(ns,nd),ns=1,nsteps_stuff)
end do
!------------------------------------------------------------------------
!
!   *********************   AFTER LOOPS  *******************************
!
!   WKK
!   Writing out the PRECT WK File
!
!   First step: Create the file.
!
!------------------------------------------------------------------------
print*,'----------  Starting section to write out WK File ------------'
!stop 'before first nf90_create' ! before div by zero with -O3

status = nf90_create(FILE_NAME_PRECT_WK, nf90_clobber, ncid_wk)
if (status /= nf90_noerr) then
  stop 'create file incorrect'
endif

print*,'Called nf90_create for WK file ncid_wk = ',ncid_wk,nf90_clobber
print*,'FILE_NAME_PRECT_WK = ',FILE_NAME_PRECT_WK

!stop 'after first nf90_create' !  before div by zero with -O3
!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Define the dimensions. The record dimension is defined to have
! unlimited length - it can grow as needed. In this example it is
! the time dimension.
!
!  The function NF90_DEF_DIM adds a new dimension to an open netCDF dataset in 
!  define mode. It returns (as an argument) a dimension ID, given the netCDF ID, 
!  the dimension name, and the dimension length. At most one unlimited length 
!  dimension, called the record dimension, may be defined for each netCDF dataset.
!
!  ncid_wk should be defined
!  LON_NAME should be defined
!  NLONS_WK should be defined
!  lon_dimid is the output
!
!------------------------------------------------------------------------
print*,'Starting section defining the dimensions'
print*,'LON_NAME = ',LON_NAME
print*,'NLONS_WK = ',NLONS_WK
print*,'lon_dimid = ',lon_dimid
!stop 'before first nf90_def_dim'  ! before zero problem

status = nf90_def_dim(ncid_wk, LON_NAME, NLONS_WK, lon_dimid)
if (status /= nf90_noerr) then
  stop 'def_dim for lon incorrect'
endif

status = nf90_def_dim(ncid_wk, LAT_NAME, NLATS_WK, lat_dimid)
if (status /= nf90_noerr) then
  stop 'def_dim for lat incorrect'
endif

status = nf90_def_dim(ncid_wk, REC_NAME, NRECS_WK, rec_dimid)
if (status /= nf90_noerr) then
  stop 'def_dim for tim incorrect'
endif

print*,'lon_dimid NLONS_WK = ',lon_dimid,NLONS_WK
print*,'lat_dimid NLATS_WK = ',lat_dimid,NLATS_WK
print*,'rec_dimid NRECS_WK = ',rec_dimid,NRECS_WK

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Define the coordinate variables. We will only define coordinate
! variables for lat and lon.  Ordinarily we would need to provide
! an array of dimension IDs for each variable's dimensions, but
! since coordinate variables only have one dimension, we can
! simply provide the address of that dimension ID (lat_dimid) and
! similarly for (lon_dimid).
!
! The function NF90_DEF_VAR adds a new variable to an open netCDF dataset 
! in define mode. It returns (as an argument) a variable ID, given the netCDF 
! ID, the variable name, the variable type, the number of dimensions, and a 
! list of the dimension IDs.
!
! NF90_EBADNAME The specified variable name is the name of another existing variable.
!
! I am getting: Note: The following floating-point exceptions are signalling: 
!  IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
!
!------------------------------------------------------------------------
print*,'--- Starting section defining the coordinate variables'
print*,'should be defined ncid_wk = ',ncid_wk  ! input
print*,'should be defined LON_NAME ',LON_NAME  ! input
print*,'variable type NF90_REAL = ',NF90_REAL ! variable type; input predefined try NF90_FLOAT?
print*,'variable type NF90_FLOAT = ',NF90_FLOAT ! variable type; input predefined try NF90_FLOAT?
print*,'should be defined lon_dimid = ',lon_dimid ! dimension ID
print*,'should be undefined lon_varid = ',lon_varid  ! Should be undefined; Output

print*,'Starting def_var for longitude ncid_wk = ',ncid_wk

status = nf90_def_var(ncid_wk, LON_NAME, NF90_REAL, lon_dimid, lon_varid)
if (status /= nf90_noerr) then
  print*,'status = ',status
  stop 'def_var for lon incorrect'
endif

print*,'Starting def_var for latitude'
status = nf90_def_var(ncid_wk, LAT_NAME, NF90_REAL, lat_dimid, lat_varid)
if (status /= nf90_noerr) then
  stop 'def_var for lat incorrect'
endif

print*,'Starting def_var for time'
status = nf90_def_var(ncid_wk, REC_NAME, NF90_REAL, rec_dimid, rec_varid)
if (status /= nf90_noerr) then
  stop 'def_var for rec incorrect'
endif
print*,'lon_varid = ',lon_varid
print*,'lat_varid = ',lat_varid
print*,'rec_varid = ',rec_varid

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! The dimids array is used to pass the dimids of the dimensions of
! the netCDF variables. Both of the netCDF variables we are creating
! share the same four dimensions. In Fortran, the unlimited
! dimension must come last on the list of dimids.
!
!------------------------------------------------------------------------
print*,'lon_dimid = ',lon_dimid
print*,'lat_dimid = ',lat_dimid
print*,'rec_dimid = ',rec_dimid
dimids_wk = (/ lon_dimid, lat_dimid, rec_dimid /)

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Define the netCDF variables for the data array
!
!------------------------------------------------------------------------
print*,'Defining rain array'
status = nf90_def_var(ncid_wk, RAIN_NAME, NF90_REAL, dimids_wk, rain_varid)
print*,'rain_varid = ',rain_varid
if (status /= nf90_noerr) then
  stop 'def_var not written correctly'
endif

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! End define mode.
!
!------------------------------------------------------------------------
status =  nf90_enddef(ncid_wk)
if (status /= nf90_noerr) then
  stop 'enddef not written correctly'
endif

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Write longitude coordinate variable data to netCDF file.
!
!------------------------------------------------------------------------
status = nf90_put_var(ncid_wk, lon_varid, lons_wk)
if (status /= nf90_noerr) then
  stop 'longitude coordinate data not written correctly'
endif
print*,'longitude coordinate data written correctly'

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Write latitude coordinate variable data to netCDF file.
!
!------------------------------------------------------------------------
status = nf90_put_var(ncid_wk, lat_varid, lats_wk)
if (status /= nf90_noerr) then
  stop 'latitude coordinate data not written correctly'
endif
print*,'latitude coordinate data written correctly'


!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
! Write time coordinate variable data to netCDF file.
!
!------------------------------------------------------------------------
status = nf90_put_var(ncid_wk, rec_varid, tim_wk)

if (status /= nf90_noerr) then
  print*,'time coordinate data not written correctly'
  stop 'time coordinate data not written correctly'
endif
print*,'time coordinate data written correctly'

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
!  Not sure what this does; may not be needed
!
!------------------------------------------------------------------------
  count = (/ NLONS_WK, NLATS_WK, NRECS_WK /)
  start = (/ 1, 1, 1 /)

!------------------------------------------------------------------------
!
!   Writing out the PRECT WK File
!
!------------------------------------------------------------------------
print*,'Writing the data'
status = nf90_put_var(ncid_wk, rain_varid, prect_wk, start = start, count = count)

if (status /= nf90_noerr) then
  print*,'prect_wk not written correctly'
  stop 'prect_wk not written correctly'
endif
print*,'prect_wk written correctly'

!------------------------------------------------------------------------
!
! Close the file. This causes netCDF to flush all buffers and make
! sure your data are really written to disk.
!
!------------------------------------------------------------------------
print*,'Closing the WK netcdf file'
status = nf90_close(ncid_wk)

if (status /= nf90_noerr) then
  print*,'file not closed correctly'
  stop 'file not closed correctly'
endif
print *,"*** SUCCESS writing netcdf file"


!------------------------------------------------------------------------
!
!  end
!
!------------------------------------------------------------------------
end program toy
!------------------------------------------------------------------------
!
!**********************************************************************
!***************  MV: move_vert
!**********************************************************************
!
!------------------------------------------------------------------------
subroutine move_vert(i,j,k1,k2,dmdry_move,nmove)
!------------------------------------------------------------------------
!
!   Move mass in the vertical direction for LT and UT levels
!   - checks for colwat and dmdry errors
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
integer, intent(in) :: i
integer, intent(in) :: j
integer, intent(in) :: k1  ! originating level
integer, intent(in) :: k2  ! receiving level
double precision, intent(inout) :: dmdry_move ! mass moved
integer, intent(in) :: nmove ! which process is calling move_vert
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: num_neg,k
!------------------------------------------------------------------------
!
!  MV: Local real variables
!
!------------------------------------------------------------------------
double precision :: rv1_initial,rv2_initial,hm1_initial,hm2_initial,rv2_new
double precision :: p2_bot,p2_top,rv2mass_new,p1_top,p1_bot,rv2mass_initial
double precision :: hm2mass_new,hm2_new,dm2_want,dm2_extra,dm2_have,dm1_have
double precision :: dmdry1_initial,dmdry1_new,dm1_extra,hmmass_move
double precision :: dmdry2_initial,dmdry2_new
double precision :: rvmass_move,dmdry_want
double precision :: dm1_want,mextra_initial
double precision :: colwat_initial,dmdry_initial
double precision :: colwat_final,dmdry_final,xxx,colwat,dmdry_error
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (dmdry_move /= dmdry_move) then
  print*,'dmdry_move = ',dmdry_move
  stop '------- dmdry_move NaN on entry to move_vert'
endif
!------------------------------------------------------------------------
!  MV: Define colwat_initial,dmdry_initial
!------------------------------------------------------------------------
colwat = 0.
dmdry = 0.
do k = 1,nlev
  xxx = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k))) + mextra(i,j,k)
  dmdry = dmdry + xxx
  colwat = colwat + rv(i,j,k)*xxx
  if (k == k1) dmdry1_initial = xxx
  if (k == k2) dmdry2_initial = xxx
end do
colwat_initial = colwat
dmdry_initial = dmdry
mextra_initial = mextra(i,j,k1)
!------------------------------------------------------------------------
!  Check for max value
!------------------------------------------------------------------------
if (dmdry_move > 0.5*dmdry1_initial) dmdry_move = 0.5*dmdry1_initial
!------------------------------------------------------------------------
!  Check that dmdry_move not too big.
!  Not really a problem is dmdry_move > 0.5*dmdry2_initial, since might mean
!   dmdry2_initial is too small and trying to make it bigger.
!------------------------------------------------------------------------
!if ((dmdry_move > 0.5*dmdry1_initial).or.(dmdry_move > 0.5*dmdry2_initial).or.  &
if ((dmdry_move > 0.5*dmdry1_initial).or.(dmdry1_initial < 0.).or.(dmdry2_initial < 0.)) then
  print*,'-------- WARNING: dmdry_move in move_vert quite big -----'
  print*,'i = ',i
  print*,'j = ',j
  print*,'originating level k1 = ',k1
  print*,'receivinglevel k2 = ',k2
  print*,'process = ',nmove
  print*,'process 1 = Subsidence from UT -> LT'
  print*,'process 2 = Subsidence from LT -> BL'
  print*,'process 3 = Dyn Descent from UT -> LT'
  print*,'process 4 = Dyn Ascent from LT -> UL'
  print*,'process 5 = Dyn Descent from LT -> BL'
  print*,'process 6 = Dyn Ascent from BL -> LT'
  print*,'mextra(i,j,k1) = ',mextra(i,j,k1)
  print*,'mextra(i,j,k2) = ',mextra(i,j,k2)
  print*,'dmdry_move = ',dmdry_move
  print*,'dmdry1_initial = ',dmdry1_initial
  print*,'dmdry2_initial = ',dmdry2_initial
  print*,'rain_dp (mm/day)= ',rain_dp(i,j)*(sec_day/tstep)
  print*,'mass_dp = ',mass_dp(i,j)
  print*,'mass_dn = ',mass_dn(i,j)
  print*,'omega1(i,j) = ',omega1(i,j)
  print*,'omega2(i,j) = ',omega2(i,j)
  print*,'sigma1(i,j) = ',sigma1(i,j)
  print*,'sigma2(i,j) = ',sigma2(i,j)
  print*,'sigma_scale/sigma1(i,j) = ',sigma_scale/sigma1(i,j)
  print*,'sigma_scale/sigma2(i,j) = ',sigma_scale/sigma2(i,j)
  stop '--------- dmdry_move too big at start of move_vert ----'
  if (dmdry1_initial > dmdry2_initial) then
    dmdry_move = 0.5*dmdry2_initial
  else
    dmdry_move = 0.5*dmdry1_initial
  endif
endif
!------------------------------------------------------------------------
!  MV: Define originating(k1) and receiving(k2) initial properties
!------------------------------------------------------------------------
rv1_initial = rv(i,j,k1)
rv2_initial = rv(i,j,k2)
hm1_initial = hm(i,j,k1)
hm2_initial = hm(i,j,k2)
p1_bot = ph(k1)
p1_top = ph(k1+1)
p2_bot = ph(k2)
p2_top = ph(k2+1)
!------------------------------------------------------------------------
!  MV: Define total rv and hm being moved.
!------------------------------------------------------------------------
rvmass_move = dmdry_move*rv(i,j,k1)
hmmass_move = dmdry_move*hm(i,j,k1)
dmdry1_new = dmdry1_initial  - dmdry_move ! 
!------------------------------------------------------------------------
!  Diagnostics
!------------------------------------------------------------------------
if ((i == i_want).and.(j == j_want).and.(nprint_vert == 1)) then
  print*,'rvmass_move = ',rvmass_move
  print*,'rv1_initial = ',rv1_initial
  print*,'rv2_initial = ',rv2_initial
  print*,'Moving from k1 to k2'
  print*,'k1 k2 = ',k1,k2
  print*,'p1_bot p1_top = ',p1_bot,p1_top
  print*,'p2_bot p2_top = ',p2_bot,p2_top
  print*,'dmdry_move = ',dmdry_move
  print*,'dmdry1_initial = ',dmdry1_initial
  print*,'dmdry2_initial = ',dmdry2_initial
  print*,'before adjust mextra(i,j,k1) = ',mextra(i,j,k1)
  print*,'before adjust mextra(i,j,k2) = ',mextra(i,j,k2)
endif
!------------------------------------------------------------------------
!  MV: Define new total drymass/rv/hm of the receiving layer (layer 2)
!  - the total initial plus the amount moving.
!------------------------------------------------------------------------
dmdry2_initial = (p2_bot - p2_top)/(g*(1.0+rv2_initial)) + mextra(i,j,k2)
dmdry2_new = dmdry2_initial + dmdry_move
rv2mass_initial = dmdry2_initial*rv2_initial
rv2mass_new = rv2mass_initial + rvmass_move
hm2mass_new = hm2_initial*dmdry2_initial + hmmass_move
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if ((dmdry1_new < 0.).or.(dmdry2_new < 0.)) then
  print*,'--------  dmdry1_new negative ---'
  print*,'i = ',i
  print*,'j = ',j
  print*,'originating level k1 = ',k1
  print*,'receivinglevel k2 = ',k2
  print*,'process = ',nmove
  print*,'process 1 = Subsidence from UT -> LT'
  print*,'process 2 = Subsidence from LT -> BL'
  print*,'process 5 = Dyn Descent from LT -> BL'
  print*,'process 6 = Dyn Ascent from LT -> BL'
  print*,'p1_bot p1_top = ',p1_bot,p1_top
  print*,'p2_bot p2_top = ',p2_bot,p2_top
  print*,'mextra(i,j,3) - mextra(i,j,2) = ',mextra(i,j,3)-mextra(i,j,2)
  print*,'f_sub = ',f_sub
  print*,'dmdry_move = ',dmdry_move
  print*,'dmdry1_initial = ',dmdry1_initial
  print*,'dmdry2_initial = ',dmdry2_initial
  print*,'dmdry1_new = ',dmdry1_new
  print*,'dmdry2_new = ',dmdry2_new
  print*,'before adjust mextra(i,j,k1) = ',mextra(i,j,k1)
  print*,'before adjust mextra(i,j,k2) = ',mextra(i,j,k2)
  stop '------ dmdry1_new negative in move_vert'
endif
!------------------------------------------------------------------------
!  Normalize rv2/hm2 by the total new mass
!------------------------------------------------------------------------
rv2_new = rv2mass_new/dmdry2_new  ! new rv of receiving grid cell
hm2_new = hm2mass_new/dmdry2_new  ! New hm of receiving grid cell
rv(i,j,k2) = rv2_new
hm(i,j,k2) = hm2_new
if (rv(i,j,k2) /= rv(i,j,k2)) then
  print*,'------- NaN problem ---'
  print*,'rv(i,j,k2) = ',rv(i,j,k2) 
  print*,'rv(i,j,k1) = ',rv(i,j,k1) 
  print*,'i j k2 = ',i,j,k2
  print*,'dmdry2_new = ',dmdry2_new
  print*,'rv2_new = ',rv2_new
  print*,'dmdry2_initial = ',dmdry2_initial
  print*,'dmdry_move = ',dmdry_move
  print*,'dmdry1_initial = ',dmdry1_initial
  stop '-------- rv(i,j,k2) NaN in move_vert ---------'
endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (rv2_new < 0.) then
  print*,'---------- rv2_new < 0 in move_vert ----------'
  print*,'rv2_new = ',rv2_new
  print*,'rv(i,j,k2) = ',rv(i,j,k2) 
  print*,'rv(i,j,k1) = ',rv(i,j,k1) 
  print*,'i j k2 = ',i,j,k2
  print*,'dmdry2_new = ',dmdry2_new
  print*,'dmdry2_initial = ',dmdry2_initial
  print*,'dmdry_move = ',dmdry_move
  print*,'dmdry1_initial = ',dmdry1_initial
  stop '---------- rv2_new < 0 in move_vert ----------'
endif
!------------------------------------------------------------------------
!  Want to determine the new mextra of receiving layer.
!  - should be zero in the case of uplift from the BL->LT
!------------------------------------------------------------------------
dm1_want = (p1_bot - p1_top)/g
dm2_want = (p2_bot - p2_top)/g
dm2_have = (1.0 + rv2_new)*dmdry2_new
dm1_have = (1.0 + rv1_initial)*dmdry1_new
!------------------------------------------------------------------------
!  MV: define new mextra
!------------------------------------------------------------------------
dm1_extra = dm1_have - dm1_want
dm2_extra = dm2_have - dm2_want
mextra(i,j,k1) = dm1_extra/(1.0 + rv1_initial)
mextra(i,j,k2) = dm2_extra/(1.0 + rv2_new)
!------------------------------------------------------------------------
!  Diagnostics
!------------------------------------------------------------------------
if ((i == i_want).and.(j == j_want).and.(nprint_vert == 1)) then
  print*,'Moving from UT to LT'
  print*,'dm1_want = ',dm1_want
  print*,'dm1_have = ',dm1_have
endif
!------------------------------------------------------------------------
! MV: Diagnostics
!------------------------------------------------------------------------
if ((i == i_want).and.(j == j_want).and.(nprint_vert == 1)) then
   print*,'Moving up from BL to LT'
   print*,'dm1_want = ',dm1_want
   print*,'dm1_have = ',dm1_have
   print*,'Pressure adjustment dm1_have*g = ',dm1_have*g
endif
!------------------------------------------------------------------------
!  MV: Define colwat_final,dmdry_final
!------------------------------------------------------------------------
colwat = 0.
dmdry = 0.
do k = 1,nlev
  xxx = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k))) + mextra(i,j,k)
  colwat = colwat + rv(i,j,k)*xxx
  dmdry = dmdry + xxx
end do
colwat_final = colwat
dmdry_final = dmdry
!------------------------------------------------------------------------
!   MV: Diagnostics
!------------------------------------------------------------------------
if ((i == i_want).and.(j == j_want).and.(nprint_vert == 1)) then
  print*,'After move:'
  print*,'mextra(i,j,1) = ',mextra(i,j,1)
  print*,'mextra(i,j,2) = ',mextra(i,j,2)
  print*,'mextra(i,j,3) = ',mextra(i,j,3)
endif
!------------------------------------------------------------------------
!
!  MV: Check for dmdry colwat conservation
!
!------------------------------------------------------------------------
colwat_error = colwat_final - colwat_initial
dmdry_error = dmdry_final - dmdry_initial
if ((abs(colwat_error) > colwat_tolerance).or.(abs(dmdry_error) > dmdry_tolerance)) then
   print*,'------- Water or drymass Conservation error in move_vert ---------------------'
   print*,'colwat_error = ',colwat_error
   print*,'dmdry_error = ',dmdry_error
   print*,'k1 k2 = ',k1,k2
   print*,'mextra_initial = ',mextra_initial
   print*,'colwat_initial = ',colwat_initial
   print*,'dmdry_initial = ',dmdry_initial
   print*,'colwat_final = ',colwat_final
   print*,'dmdry_final = ',dmdry_final
   print*,'dmdry_move = ',dmdry_move
   stop 'ERROR in move_vert exceeds tolerance'
endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine move_vert
!------------------------------------------------------------------------
!
!  **********************************************************************
!  ****************  TW: trade_wind   ***********************************
!  **********************************************************************
!
!------------------------------------------------------------------------
subroutine trade_wind
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!   Local integer variables in trade_wind
!
!------------------------------------------------------------------------
integer :: i,j
!------------------------------------------------------------------------
!
!   Local real variables in trade_wind
!
!------------------------------------------------------------------------
double precision :: u_tw,dx,dy,xxx,rv_want,hm_want,fx
!------------------------------------------------------------------------
!
!   Advect rv,hm
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  dx = (2.0*pi*re*cosd(lat(j))/float(nlon))  ! distance of 1 degree along x direction
  dy = (2.0*pi*re/float(nlat-1))*(2.0*dlat/360.0)
  xxx = abs((lat(j)*dy)/dist_tw)
  u_tw = u_trade_wind*exp(-xxx*xxx)
  fx = (u_tw*tstep)/dx  ! fraction of a grid cell in zonal direction
  if (i == nlon) then
    rv_want = (1.0-fx)*rv(i,j,1) + fx*rv(1,j,1)
    hm_want = (1.0-fx)*hm(i,j,1) + fx*hm(1,j,1)
  else
    rv_want = (1.0-fx)*rv(i,j,1) + fx*rv(i+1,j,1)
    hm_want = (1.0-fx)*hm(i,j,1) + fx*hm(i+1,j,1)
  endif
  rv(i,j,1) = rv_want
  hm(i,j,1) = hm_want
end do
end do
!------------------------------------------------------------------------
!
!  End
!
!------------------------------------------------------------------------
end subroutine trade_wind
!------------------------------------------------------------------------
!
!  **********************************************************************
!  ****************  SP: spread_mass  ***********************************
!  **********************************************************************
!
!  Gaussian spreading:
!  ------------------
!  - from one grid point, loop over all other grid points
!  - Calculate the distance
!  - Calculate the Gaussian factor
!  - Calculate the Gaussian weighted mean mextra about a grid point.
!  - If local mextra is larger than this Gaussian mean, redistribute some fraction of this
!    excess according to the deficit from the mean, and the Gaussian weighting.
!  - Similarly for a deficit.
!
!  Steps:
!  - calculate mean pressure pav of adjacent grid cells
!  - From zzz at ppp, and hyd relation calculate heights at pav
!
!------------------------------------------------------------------------
subroutine spread_mass
!------------------------------------------------------------------------
!
!   Move mass in the zonal (ndir = 1) or meridional (ndir = 2) direction
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!   Local integer variables in spread_mass
!
!------------------------------------------------------------------------
integer :: i,j,k
integer :: ii,jj,kk
!------------------------------------------------------------------------
!
!   Local real variables in spread_mass
!
!------------------------------------------------------------------------
double precision :: p1_bot,p1_top,dx,dy,mextra_move,zzz,zzz_lt,zzz_bl,zzz_ut,f_move
double precision :: rv1,hm1,dmdry_new,dmm
double precision :: weight_x_east(0:nlon-1,nlat,nlev),weight_y(nlat,nlat,nlev)
double precision :: weight_x_west(0:nlon-1,nlat,nlev)
double precision :: dm1,dm2,dmdry1,dmdry2,dmadd,rvadd,hmadd
double precision :: dmdry1_want,dmdry2_want,dm_move
double precision :: colwat_bl,dmdry_bl
double precision :: colwat_lt,dmdry_lt
double precision :: colwat_ut,dmdry_ut
double precision :: colwat_final,dmdry_final,yyy,xxx,colwat,dmdry_error
double precision :: dmdry_error_bl,colwat_error_bl
double precision :: dmdry_error_lt,colwat_error_lt
double precision :: dmdry_error_ut,colwat_error_ut
double precision :: pav,pp1,pp2,tt1,tt2,dpp1,dpp2,tv1,tv2,rho1,rho2,dz1,dz2,zz1,zz2
double precision :: mextra_mean,weight_total,mextra_diff,weight_total_eq
double precision :: dmmass_initial,rvmass_initial,hmmass_initial
double precision :: dmmass_final,rvmass_final,hmmass_final
double precision :: dist(nlon,nlat),weight(nlon,nlat)
double precision :: dm_add(nlon,nlat,nlev)
double precision :: rv_add(nlon,nlat,nlev)
double precision :: hm_add(nlon,nlat,nlev)
integer :: used(nlon,nlat),jjj,iii,di,i_eq,j_eq,i1,i2,j1,j2
!------------------------------------------------------------------------
!
!  Update pph - pressure at boundaries
!
!------------------------------------------------------------------------
!stop 'starting spread_mass'  ! no problem to here
do i = 1,nlon
do j = 1,nlat
  xxx = mextra(i,j,1) + mextra(i,j,2) + mextra(i,j,3)
  pph(i,j,1) = ph(1) + g*xxx
  xxx = mextra(i,j,2) + mextra(i,j,3)
  pph(i,j,2) = ph(2) + g*xxx
  xxx = mextra(i,j,3)
  pph(i,j,3) = ph(3) + g*xxx
  pph(i,j,4) = ph(4)
end do
end do
!------------------------------------------------------------------------
!
!  Update ppp - pressure at midpoints
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
do k = 1,nlev
  ppp(i,j,k) = 0.5*(pph(i,j,k) + pph(i,j,k+1))
end do
end do
end do
!------------------------------------------------------------------------
!
!  Calculate Rossby radius.
!  - zzz factor: take into account the reduction in the Rossby radius with
!    distance from the equator.
!------------------------------------------------------------------------
do j = 1,nlat
! bl
  if (abs(lat(j)) < lat_rossby_bl) then
    zzz_bl = 1.0  ! close to equator ros_radius has its max value
    xxx = abs(lat(j))/lat_rossby_bl
    zzz_bl = 1.0 + ros_radius_add*(1.- xxx*xxx) ! close to equator ros_radius has its max value
  else
    zzz_bl = lat_rossby_bl/abs(lat(j))
  endif
! lt
  if (abs(lat(j)) < lat_rossby_lt) then
    zzz_lt = 1.0  ! close to equator ros_radius has its max value
    xxx = abs(lat(j))/lat_rossby_lt
    zzz_lt = 1.0 + ros_radius_add*(1.- xxx*xxx) ! close to equator ros_radius has its max value
  else
    zzz_lt = lat_rossby_lt/abs(lat(j))
  endif
!  ut
  if (abs(lat(j)) < lat_rossby_ut) then
!    zzz_ut = 1.0  ! close to equator ros_radius has its max value
    xxx = abs(lat(j))/lat_rossby_ut
    zzz_ut = 1.0 + ros_radius_add*(1.- xxx*xxx) ! close to equator ros_radius has its max value
!    print*,'j lat(j) zzz_ut = ',j,lat(j),zzz_ut
  else
    zzz_ut = lat_rossby_ut/abs(lat(j))
  endif
!
  ros_radius(j,1) = dist_bl*zzz_bl
  ros_radius(j,2) = dist_lt*zzz_lt
  ros_radius(j,3) = dist_ut*zzz_ut
end do
!------------------------------------------------------------------------
!
!  SP: Define colwat,dmdry of each layer
!
!------------------------------------------------------------------------
do k = 1,nlev
  colwat = 0.
  dmdry = 0.
  do i = 1,nlon
  do j = 1,nlat
    xxx = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k))) + mextra(i,j,k)
    dmdry = dmdry + xxx
    colwat = colwat + rv(i,j,k)*xxx
  end do
  end do
  colwat = colwat/(float(nlat)*float(nlon))
  dmdry = dmdry/(float(nlat)*float(nlon))
  if (k == 1) colwat_bl = colwat
  if (k == 1) dmdry_bl = dmdry
  if (k == 2) colwat_lt = colwat
  if (k == 2) dmdry_lt = dmdry
  if (k == 3) colwat_ut = colwat
  if (k == 3) dmdry_ut = dmdry
end do
!------------------------------------------------------------------------
!
!  Initialize
!
!------------------------------------------------------------------------
do k = 1,nlev
do i = 1,nlon
do j = 1,nlat
   dm_add(i,j,k) = 0.
   rv_add(i,j,k) = 0.
   hm_add(i,j,k) = 0.
end do
end do
end do
!------------------------------------------------------------------------
!
!  Calculate weight_y(j,jj,k)
!  - for jj >= j
!
!------------------------------------------------------------------------
dy = ((pi*re)/180.)   ! distance of 1 degree along the y direction in m
do k = 1,nlev
do j = 1,nlat
do jj = j,nlat  ! Note starting from j here 
  weight_y(j,jj,k) = 0.
  if (jj == j) then
    weight_y(j,jj,k) = 1.
  else
    if (jj-1 == 0) stop 'jj-1 problem'
    xxx = dy/(0.5*(ros_radius(jj,k) + ros_radius(jj-1,k)))
    weight_y(j,jj,k) = weight_y(j,jj-1,k)*exp(-xxx*xxx)
  endif
end do
end do
end do
!------------------------------------------------------------------------
!
!  Assume weight_y(j,jj,k) is symmetrical in j,jj
!
!------------------------------------------------------------------------
do k = 1,nlev
do j = 1,nlat
do jj = 1,nlat
  if (jj < j) weight_y(j,jj,k) = weight_y(jj,j,k)
!  if ((k == 1).and.(j == 30)) print*,'jj lat(jj) weight_y(j,jj,k) = ',jj,lat(jj),weight_y(j,jj,k)
end do
end do
end do
!stop  'got weight_y'  ! No problem to here
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
do j = 1,nlat
do k = 1,nlev
do i = 0,nlon-1
  weight_x_east(i,j,k) = 0.
  weight_x_west(i,j,k) = 0.
end do
end do
end do
!------------------------------------------------------------------------
!
!  Calculate weight_x(0:nlon-1,nlat,nlev)
!  - the first index refers to the grid point separation along the x direction
!    of the two points. 
!
!------------------------------------------------------------------------
!stop 'starting weight_x calculation'  ! No problem to here
do j = 1,nlat
!  print*,'j = ',j
  dx = (2.0*pi*re*cosd(lat(j))/float(nlon))  ! distance of 1 degree along x direction
!  print*,'dx = ',dx
  do k = 1,nlev
    xxx = dx/(ros_radius(j,k)*zonal_amp*float(apply_zonal_amp(k))) ! for zonal_amp = 0.9, xxx bigger
    yyy = (dx*zonal_amp*float(apply_zonal_amp(k)))/ros_radius(j,k) ! for zonal_amp = 0.9, yyy smaller
    do i = 0,nlon-1
!      print*,'i = ',i
      if (i == 0) then
        weight_x_east(i,j,k) = 1.
        weight_x_west(i,j,k) = 1.
      else
        if (weight_x_east(i-1,j,k) > 1.0E-10) then
          weight_x_east(i,j,k) = weight_x_east(i-1,j,k)*exp(-xxx*xxx)
        else
          weight_x_east(i,j,k) = 0.0
        endif
        if (weight_x_west(i-1,j,k) > 1.0E-10) then
          weight_x_west(i,j,k) = weight_x_west(i-1,j,k)*exp(-yyy*yyy)
        else
          weight_x_west(i,j,k) = 0.0
        endif
      endif
!      print*,'weight_x_east(i,j,k) = ',weight_x_east(i,j,k)
!      print*,'weight_x_west(i,j,k) = ',weight_x_west(i,j,k)
    end do
  end do
end do
!stop 'got weight_x'
!------------------------------------------------------------------------
!
!  weight_y and weight_x now determined.
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
! Loop over heights
!
!------------------------------------------------------------------------
do k = 1,nlev
!------------------------------------------------------------------------
!
!  Get weight_total_eq: for a location on the equator
!
!------------------------------------------------------------------------
  i_eq = int(nlon/2)
  j_eq = int(nlat/2)
!  stop 'calling get_weight'
  call get_weight(i_eq,j_eq,k,used,weight,mextra_mean,weight_total_eq,  &
           weight_x_east,weight_x_west,weight_y)
!------------------------------------------------------------------------
!
!  Loop over i,j
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!
! Get weight function for i,j,k
!
!------------------------------------------------------------------------
!  stop 'calling get_weight'
  call get_weight(i,j,k,used,weight,mextra_mean,weight_total,  &
           weight_x_east,weight_x_west,weight_y)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
!  if ((i == 5).and.(j == 5)) then
!    print*,'weight_total mextra_mean = ',weight_total,mextra_mean
!    do ii = 1,nlon
!    do jj = 1,nlat
!      print*,'ii jj weight(ii,jj) = ',ii,jj,weight(ii,jj)
!    end do
!    end do
!  endif
!------------------------------------------------------------------------
!
! Save upper tropospheric weight at one location
!
!------------------------------------------------------------------------
  if ((j == j_weight_save_upper).and.(i == i_weight_save_upper).and.(nstep == 1).and.(k == 3)) then
  do ii = 1,nlon
  do jj = 1,nlat
      weight_save_upper(ii,jj) = weight(ii,jj)
  end do
  end do
  endif
!------------------------------------------------------------------------
!
! Save lower tropospheric weight at one location
!
!------------------------------------------------------------------------
  if ((j == j_weight_save_lower).and.(i == i_weight_save_lower).and.(nstep == 1).and.(k == 2)) then
  do ii = 1,nlon
  do jj = 1,nlat
      weight_save_lower(ii,jj) = weight(ii,jj)
  end do
  end do
  endif
  do ii = 1,nlon
  do jj = 1,nlat
      weight_save(ii,jj) = weight_save_lower(ii,jj) + weight_save_upper(ii,jj)
  end do
  end do
!------------------------------------------------------------------------
!
!  Inside loop i,j,k
!
!  Determine mextra_diff for i,j,k
!  - This is the mean difference with respect to the weighted surrounding area
!  - Inside loop over i,j,k
!
!------------------------------------------------------------------------
  mextra_mean = mextra_mean/weight_total
  mextra_diff = mextra(i,j,k) - mextra_mean
!------------------------------------------------------------------------
!
!  Inside loop i,j,k
!
!  Calculate f_move
!  - decrease depending on local Rossby radius.
!
!------------------------------------------------------------------------
   f_move = f_move_eq*(weight_total/weight_total_eq)
!   if (i == 1) print*,'f_move lat(j) = ',f_move,lat(j)
!------------------------------------------------------------------------
!
!  Inside loop i,j,k
!
!   Move mass,rv,hm for all ii,jj for which used(ii,jj) = 1 
!  - Inside loop over i,j,k
!  - Changing the order of the ii,jj looping here DOES affect rain variance. 
!  YYYY1 = YYYY3 but YYYY2 is different
!
!------------------------------------------------------------------------
! YYYY1
  do ii = 1,nlon
  do jj = 1,nlat
! YYYY2 
!  do ii = nlon,1,-1
!  do jj = nlat,1,-1
! YYYY3
!  do jj = 1,nlat
!  do ii = 1,nlon
  if (used(ii,jj) == 1) then
!------------------------------------------------------------------------
!
!  Inside loop i,j,k
!
!  - local values of rv/hm do not change
!  - mextra_diff and dmadd should be the same sign.
!
!------------------------------------------------------------------------
    dmadd = f_move*mextra_diff*weight(ii,jj)/weight_total
    if (weight(ii,jj) < 0.) stop 'weight(ii,jj) negative'
!------------------------------------------------------------------------
!   mextra_diff > 0: Move mass from i,j to ii,jj
!------------------------------------------------------------------------
    if (mextra_diff > 0.) then
      rvadd =  dmadd*rv(i,j,k)
      hmadd = dmadd*hm(i,j,k)
      if (dmadd < 0.) stop 'dmadd less than zero'
      if (rvadd < 0.) stop 'rvadd less than zero'
      if (hmadd < 0.) stop 'hmadd less than zero'
      dm_add(i,j,k) = dm_add(i,j,k) - dmadd
      rv_add(i,j,k) = rv_add(i,j,k) - rvadd
      hm_add(i,j,k) = hm_add(i,j,k) - hmadd
      dm_add(ii,jj,k) = dm_add(ii,jj,k) + dmadd
      rv_add(ii,jj,k) = rv_add(ii,jj,k) + rvadd
      hm_add(ii,jj,k) = hm_add(ii,jj,k) + hmadd
!------------------------------------------------------------------------
!   mextra_diff < 0: Move mass from ii,jj to i,j
!------------------------------------------------------------------------
    else 
      rvadd =  dmadd*rv(ii,jj,k)
      hmadd = dmadd*hm(ii,jj,k)
      if (dmadd > 0.) stop 'dmadd larger than zero'
      if (rvadd > 0.) stop 'rvadd larger than zero'
      if (hmadd > 0.) stop 'hmadd larger than zero'
      dm_add(i,j,k) = dm_add(i,j,k) - dmadd  ! adding to i,j
      rv_add(i,j,k) = rv_add(i,j,k) - rvadd  ! adding to i,j
      hm_add(i,j,k) = hm_add(i,j,k) - hmadd  ! adding to i,j
      dm_add(ii,jj,k) = dm_add(ii,jj,k) + dmadd  ! removing from ii,jj
      rv_add(ii,jj,k) = rv_add(ii,jj,k) + rvadd  ! removing from ii,jj 
      hm_add(ii,jj,k) = hm_add(ii,jj,k) + hmadd  ! removing from ii,jj
    endif
!------------------------------------------------------------------------
!
!  end loop over ii,jj
!
!------------------------------------------------------------------------
  endif   ! for used(ii,jj) = 1
  end do
  end do
!------------------------------------------------------------------------
!
!  end loop over k,i,j
!
!------------------------------------------------------------------------
end do
end do
end do
!------------------------------------------------------------------------
!
!
!  Make adjustments: NEW loop over i,j,k
!
!------------------------------------------------------------------------
do k = 1,nlev
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!  Find new dry mass
!------------------------------------------------------------------------
  p1_bot = ph(k)
  p1_top = ph(k+1)
  rv1 = rv(i,j,k)
  hm1 = hm(i,j,k)
  dm1 = ((p1_bot - p1_top)/g) 
  dmmass_initial = (dm1/(1.0+rv1)) + mextra(i,j,k)
  rvmass_initial = rv1*dmmass_initial
  hmmass_initial = hm1*dmmass_initial
  dmmass_final = dmmass_initial + dm_add(i,j,k)
  rvmass_final = rvmass_initial + rv_add(i,j,k) 
  hmmass_final = hmmass_initial + hm_add(i,j,k)
  rv(i,j,k) = rvmass_final/dmmass_final
  hm(i,j,k) = hmmass_final/dmmass_final
  dmm = ((p1_bot - p1_top)/g) 
  dmm = dmm/(1.0+rv(i,j,k))
  mextra(i,j,k) = dmmass_final - dmm
  mconv(i,j,k) = dm_add(i,j,k)
!------------------------------------------------------------------------
!  end loop over i,j,k
!------------------------------------------------------------------------
end do
end do
end do
!------------------------------------------------------------------------
!
!   *****************   ERROR CHECK  **************************
!
!  SP: Define colwat_final,dmdry_final
!
!------------------------------------------------------------------------
do kk = 1,nlev
  colwat = 0.
  dmdry = 0.
  do ii = 1,nlon
  do jj = 1,nlat
    xxx = (ph(kk) - ph(kk+1))/(g*(1.0+rv(ii,jj,kk))) + mextra(ii,jj,kk)
    dmdry = dmdry + xxx
    colwat = colwat + rv(ii,jj,kk)*xxx
  end do
  end do
  colwat = colwat/(float(nlat)*float(nlon))
  dmdry = dmdry/(float(nlat)*float(nlon))
  if (kk == 1) colwat_error_bl = colwat - colwat_bl
  if (kk == 1) dmdry_error_bl = dmdry - dmdry_bl
  if (kk == 2) colwat_error_lt = colwat - colwat_lt
  if (kk == 2) dmdry_error_lt = dmdry - dmdry_lt
  if (kk == 3) colwat_error_ut = colwat - colwat_ut
  if (kk == 3) dmdry_error_ut = dmdry - dmdry_ut
end do
!------------------------------------------------------------------------
!
!  Check for dmdry colwat conservation
!  - There are problems with mass/colwat conservation here, I think mainly having
!    to do with the extremely large number of calculations being made, and the fact
!    that small numbers are always being added to big numbers. For some reason the
!    errors aren't random but accumulate.
!
!------------------------------------------------------------------------
!print*,'spread_mass: dmdry_error_bl dmdry_error_lt dmdry_error_ut = ',dmdry_error_bl,dmdry_error_lt,dmdry_error_ut
!print*,'spread_mass: colwat_error_bl colwat_error_lt colwat_error_ut = ',colwat_error_bl,colwat_error_lt,colwat_error_ut
if ((abs(colwat_error_bl) > colwat_tolerance).or.(abs(dmdry_error_bl) > dmdry_tolerance)) then
   print*,'-------------- Conservation error in spread_mass -------------------'
   print*,'colwat_error_bl = ',colwat_error_bl
   print*,'dmdry_error_bl = ',dmdry_error_bl
   stop 'ERROR in spread_mass exceeds tolerance'
endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine spread_mass
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  get_weight  ******************************************
! *********************************************************************
!
!------------------------------------------------------------------------
subroutine get_weight(i,j,k,used,weight,mextra_mean,weight_total,  &
           weight_x_east,weight_x_west,weight_y)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
integer, intent(in) :: i
integer, intent(in) :: j
integer, intent(in) :: k
integer, intent(out) :: used(nlon,nlat)
double precision, intent(out) :: weight(nlon,nlat)
double precision, intent(out) :: mextra_mean
double precision, intent(out) :: weight_total
double precision, intent(in) :: weight_x_east(0:nlon-1,nlat,nlev)
double precision, intent(in) :: weight_x_west(0:nlon-1,nlat,nlev)
double precision, intent(in) :: weight_y(nlat,nlat,nlev)
!------------------------------------------------------------------------
!
!  Local variables
!
!------------------------------------------------------------------------
integer :: ii,jj,di
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
mextra_mean = 0. 
weight_total = 0.0d0
!------------------------------------------------------------------------
!
!  Loop over ii,jj
!  - A change in the oder of this loop does not affect the weights
!
!------------------------------------------------------------------------
do ii = 1,nlon
do jj = 1,nlat
!------------------------------------------------------------------------
!
!  Initialize used and weight
!
!------------------------------------------------------------------------
  used(ii,jj) = 0
  weight(ii,jj) = 0.0d0
!------------------------------------------------------------------------
!
!  Determine di
!  - i is the "origin" longitude index. If i < ii, multiply weight by some factor.
!  Suppose i = 360 and ii = 1: want to use zonal_amp
!  Suppose i = 181 and ii = 1: want to use zonal_amp
!  Suppose i = 180 and ii = 1: do not want to use zonal_amp
!  Suppose i = 360 and ii = 1: do want to use zonal_amp
!  Suppose i = 360 and ii = 179: do want to use zonal_amp
!------------------------------------------------------------------------
    di = abs(i-ii)
    if (abs(i+nlon-ii) < di) di = abs(i+nlon-ii)
    if (abs(i-nlon-ii) < di) di = abs(i-nlon-ii)
    if (di > (nlon-1)) stop 'di out of range'
    if ((ii > i).and.(ii < (i + 180))) then 
       weight(ii,jj) = weight_x_east(di,jj,k)*weight_y(j,jj,k)
    elseif ((i > ii).and.((ii + 360 - i) <= 180)) then
       weight(ii,jj) = weight_x_east(di,jj,k)*weight_y(j,jj,k)
    else
       weight(ii,jj) = weight_x_west(di,jj,k)*weight_y(j,jj,k)
    endif
!------------------------------------------------------------------------
!
!  end loop over ii,jj
!
!------------------------------------------------------------------------
end do
end do
!------------------------------------------------------------------------
!
!  Add them up
!  - The order of the loops here does affect weight_total. Apparently matters
!    if start with the bigger numbers first or the smaller numbers. This does
!    then affect the rainfall variance after a few weeks. Apparently is a result
!    of adding numbers of very different magnitudes.
!
!------------------------------------------------------------------------
!  BBBB1
do ii = 1,nlon
do jj = 1,nlat
!  BBBB2
!  do jj = 1,nlat
!  do ii = 1,nlon
  if (weight(ii,jj) > weight_min) then
    used(ii,jj) = 1
    weight_total = weight_total + weight(ii,jj)
    mextra_mean = mextra_mean + mextra(ii,jj,k)*weight(ii,jj)
  else
    used(ii,jj) = 0
  endif
end do
end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine get_weight
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
! *********************************************************************
! *************  Evaporation  *****************************************
! *********************************************************************
!
!  Evaporation proportional to rvmass_avail = (rv_evap_max - rv_lt)*dmdry
!
!------------------------------------------------------------------------
subroutine evap(kk)
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
integer, intent(in) :: kk  ! height level
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: i,j,k,nprint
!------------------------------------------------------------------------
!
!   Local real variables in evap
!
!------------------------------------------------------------------------
double precision :: t_initial,wat_error,rvmass_max,rvmass_initial,rvmass_evap,rh_evap_max
double precision :: rv_initial,rv_final,dmdry_error,dmdry_final,dmdry_initial,hm_initial
double precision :: km_final,km_initial,pppp,ri_final,rl_final,t_final,rt_final,f_evap
double precision :: rain_evap_avail,rhh,rvmass_evap_max
!------------------------------------------------------------------------
!
!  To avoid creating instability, only do evaporation if level 2 PT
!  is some threshold larger than level 1 PT.
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
if (kk == 2) then
  rh_evap_max = rh_cg_evap_max
endif
if (kk == 3) then
  rh_evap_max = rh_dp_evap_max
endif
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
if (pt(i,j,2) > (pt(i,j,1) + pt_add)) then
  rvmass_evap = 0.
  rv_initial = rv(i,j,2)
  hm_initial = hm(i,j,2)
  t_initial = tt(i,j,2)  ! though t may not be consistent with hm/rv
!------------------------------------------------------------------------
!  Find f_evap and rain_evap_avail
!------------------------------------------------------------------------
  rhh = rv_initial/rsat(t_initial,pp(2))
  if (kk == 2) then 
    f_evap = sigmoidal(rhh,f_evap_cg_half,f_evap_cg_scale,f_evap_cg_low,f_evap_cg_add)
    rain_evap_avail = f_evap*(rain_cg(i,j) - evap_cg(i,j))
  else
    f_evap = sigmoidal(rhh,f_evap_dp_half,f_evap_dp_scale,f_evap_dp_low,f_evap_dp_add)
    rain_evap_avail = f_evap*(rain_dp(i,j) + rain_an(i,j) - evap_dp(i,j) - down_dp(i,j))
  endif
  if (rain_evap_avail < 0.) stop 'rain_evap_avail is negative'
!------------------------------------------------------------------------
!  Find rvmass_evap_max: max allowed evaporation
!------------------------------------------------------------------------
  dmdry_initial = (ph(2) - ph(3))/(g*(1.00000+rv_initial)) + mextra(i,j,2)
  rvmass_max = rh_evap_max*rsat(t_initial,pp(2))*dmdry_initial
  rvmass_initial = rv_initial*dmdry_initial
  rvmass_evap_max = (rvmass_max - rvmass_initial)
  if (rvmass_evap_max < 0.) rvmass_evap_max = 0.
  rvmass_evap = rain_evap_avail
  if (rvmass_evap > rvmass_evap_max) rvmass_evap = rvmass_evap_max
! print*,'rvmass_max = ',rvmass_max
! print*,'rvmass_initial = ',rvmass_initial
! print*,'rvmass_evap = ',rvmass_evap
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'============  In evaporation loop before evap nstep = ',nstep
    print*,'initial tt(i,j,2) = ',tt(i,j,2)
    print*,'rv_initial = ',rv_initial
    print*,'hm_initial = ',hm_initial
    print*,'rvmass_evap = ',rvmass_evap
    print*,'rain_dp(i,j) = ',rain_dp(i,j)
    print*,'evap_dp(i,j) = ',evap_dp(i,j)
    print*,'down_dp(i,j) = ',down_dp(i,j)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (rvmass_evap < 0.) then 
    print*,'kk = ',kk
    print*,'rvmass_evap = ',rvmass_evap
    stop 'rvmass_evap negative'
  endif
!------------------------------------------------------------------------
!  Find new T,rv if rvmass_evap > 0.
!------------------------------------------------------------------------
  if (rvmass_evap > 0.) then
!   print*,'entering evap'
    if (kk == 2) then
      evap_cg(i,j) = evap_cg(i,j) + rvmass_evap
    else
      evap_dp(i,j) = evap_dp(i,j) + rvmass_evap
      if (evap_dp(i,j) > (rain_dp(i,j)+rain_an(i,j)-down_dp(i,j))) stop 'evap_dp larger than rain_dp'
    endif
    rv_final = (rvmass_initial + rvmass_evap)/dmdry_initial
    dm_want = (ph(2) - ph(3))/g  
    dm_final = dmdry_initial*(1.00000 + rv_initial) + rvmass_evap
    dm_extra = dm_final - dm_want
    mextra(i,j,2) = dm_extra/(1.00000 + rv_final)
    km_initial = hm_initial - (1.000000 + rv_initial)*g*zz(i,j,2)
    km_final = km_initial + cl*tt(i,j,2)*(rv_final - rv_initial)
    t_final = tt(i,j,2) !  Need T initial guess
    rt_final = rv_final
    rl_final = 0.
    ri_final = 0.
    call t_from_km_allow_sat(km_final,rv_final,ppp(i,j,2),t_final)
    rv(i,j,2) = rv_final
    if (rv(i,j,2) /= rv(i,j,2)) stop 'rv(i,j,2)) NaN'
    tt(i,j,2) = t_final
    hm(i,j,2) = hmoist(t_final,rt_final,rv_final,ri_final,zz(i,j,2))
!------------------------------------------------------------------------
!  
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'============  In evaporation loop after evap nstep = ',nstep
    print*,'New tt(i,j,2) = ',tt(i,j,2)
    print*,'rv_final = ',rv_final
    print*,'Change in rv rv_final-rv_initial = ',rv_final-rv_initial
    print*,'km_initial = ',km_initial
    print*,'km_final = ',km_final
    print*,'----------'
  endif
!------------------------------------------------------------------------
!  water drymass error checks
!------------------------------------------------------------------------
    dmdry_final = (ph(2) - ph(3))/(g*(1.00000 + rv_final)) + mextra(i,j,2)
    dmdry_error = dmdry_final - dmdry_initial
    wat_error = rv_final*dmdry_final - rv_initial*dmdry_initial - rvmass_evap
    if ((abs(dmdry_error) > dmdry_tolerance).or.(abs(wat_error) > 1.0E-09)) then
      print*,'---------- dmdry or water error from evaporation -------'
      print*,'wat_error = ',wat_error
      print*,'rain_cg(i,j,) = ',rain_cg(i,j)
      print*,'evap_cg(i,j,) = ',evap_cg(i,j)
      print*,'rvmass_evap = ',rvmass_evap
      print*,'rv_initial = ',rv_initial
      print*,'rv_final = ',rv_final
      print*,'dmdry_error = ',dmdry_error
      print*,'dmdry_final = ',dmdry_final
      print*,'dmdry_initial = ',dmdry_initial
      stop '---------- dmdry error from evaporation -------'
    endif
!------------------------------------------------------------------------
!  endif for rvmass_evap > 0.
!------------------------------------------------------------------------
  endif
!------------------------------------------------------------------------
!  end loop over i,j
!------------------------------------------------------------------------
endif
end do
end do 
!------------------------------------------------------------------------
!  
!  Diagnostics  
!  
!------------------------------------------------------------------------
  if (nprint_dp == 1) then
    print*,'============  After evaporation timestep nstep = ',nstep
    print*,'tt(i_want,j_want,1) = ',tt(i_want,j_want,1)
    print*,'tt(i_want,j_want,2) = ',tt(i_want,j_want,2)
    print*,'tt(i_want,j_want,3) = ',tt(i_want,j_want,3)
    print*,'pt(i_want,j_want,1) = ',pt(i_want,j_want,1)
    print*,'pt(i_want,j_want,2) = ',pt(i_want,j_want,2)
    print*,'pt(i_want,j_want,3) = ',pt(i_want,j_want,3)
    print*,'hm(i_want,j_want,2) = ',hm(i_want,j_want,2)
    print*,'----------'
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine evap
!------------------------------------------------------------------------
!
!  **********************************************************************
!  ****************  DR: drizzle  ***************************************
!  **********************************************************************
!
!------------------------------------------------------------------------
subroutine drizzle
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: i,j,nprint
!------------------------------------------------------------------------
!
!   Local real variables in convect
!
!------------------------------------------------------------------------
double precision :: rvsp,rv_avail,rv_initial,rv_final,ri_final,km_initial
double precision :: km_wat,km_final,t_final
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (do_dr == 1) then
do i = 1,nlon
do j = 1,nlat
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  rvsp = rsat(tt(i,j,1),ppp(i,j,1))
  rv_avail = rv(i,j,1) - rh_dr_max*rvsp
  if (rv_avail > 0.) then
    dmdry = (ph(1) - ph(2))/(g*(1.0 + rv(i,j,1))) + mextra(i,j,1)
    rv_initial = rv(i,j,1)
    rv_final = rv(i,j,1) - (tstep/tau_dr)*rv_avail
    ri_final = 0.
    km_initial = hm(i,j,1) - (1. + rv_initial)*g*zz(i,j,1)  
    km_wat = cl*tt(i,j,1)
    km_final = km_initial - km_wat*(rv_initial - rv_final)
    call t_from_km_allow_sat(km_final,rv_final,ppp(i,j,1),t_final)
    rv(i,j,1) = rv_final
    tt(i,j,1) = t_final
    rain_dr(i,j) = (rv_initial - rv_final)*dmdry
    hm(i,j,1) = hmoist(t_final,rv_final,rv_final,ri_final,zz(i,j,1))
    dm_want = (ph(1) - ph(2))/g  
    dm_final = (1.0 + rv_final)*dmdry  ! total mass after rain removal
    dm_extra = dm_final - dm_want
    mextra(i,j,1) = dm_extra/(1.0 + rv_final)
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end do
end do
endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine drizzle
!------------------------------------------------------------------------
!
!  **********************************************************************
!  ****************  CN: convect  ***************************************
!  **********************************************************************
!
!------------------------------------------------------------------------
subroutine convect(kk)
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
integer, intent(in) :: kk  ! height level
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: i,j,k,nprint,ip1,im1
!------------------------------------------------------------------------
!
!   Local real variables in convect
!
!------------------------------------------------------------------------
double precision :: colwat_bf_cn,dmdry_bf_cn,colwat_af_cn,dmdry_af_cn ,mdp_initial,mcg_initial
double precision :: rain_cn,dmdry_cn,dmdry_error,dmdry_initial,dmdry_final,rv_max
double precision :: dpp,hm_final,hmp,hm_initial,km_final,kmp,rh_max,km_wat
double precision :: ri_final,rip,rl_final,rlp,rt_final,rtp,rv_bl,rv_final,rv_initial
double precision :: rvp,rvsp,tp,tp_final,zzz,xxx,dmdry_bl,dmdry_bl_af_cn,dmdry_bl_bf_cn
double precision :: dmdry_bl_error,xxx_ip1,xxx_def,xxx_im1,yyy
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
print*,'Starting Convection'
do i = 1,nlon
do j = 1,nlat
!if (nstep == 24) print*,'i j = ',i,j
!------------------------------------------------------------------------
!  Initialize 
!------------------------------------------------------------------------
   if (kk == 3) mass_dp(i,j) = 0.
   if (kk == 3) dp_prod(i,j) = 0.
   if (kk == 2) mass_cg(i,j) = 0.
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  nprint = 1
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'---entering convect --- nstep = ',nstep
    print*,'kk = ',kk
    print*,'entering mass_dp(i,j) = ',mass_dp(i,j)
    print*,'entering mass_cg(i,j) = ',mass_cg(i,j)
    print*,'initial T tt(i,j,kk) = ',tt(i,j,kk)
  endif
!------------------------------------------------------------------------
!  Calculate colwat before convection: colwat_bf_cn
!------------------------------------------------------------------------
    colwat_bf_cn = 0.
    dmdry_bf_cn = 0.
    do k = 1,nlev
      dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k)))
      dmdry = dmdry + mextra(i,j,k)
      if (k == 1) dmdry_bl = dmdry
      colwat_bf_cn = colwat_bf_cn + rv(i,j,k)*dmdry
      dmdry_bf_cn = dmdry_bf_cn + dmdry
      if (k == 1) dmdry_bl_bf_cn = dmdry
    end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'f_dp = ',f_dp(i,j)
    print*,'f_cg = ',f_cg(i,j)
  endif
!------------------------------------------------------------------------
!  Define BL parcel properties 
!------------------------------------------------------------------------
  dmdry_cn = 0.
  rv_bl = rv(i,j,1)
  rvp = rv(i,j,1)
  rtp = rvp
  rip = 0.
  rlp = 0.
  hmp = hm(i,j,1)
!------------------------------------------------------------------------
!  From hmp conservation, find tp,rlp,rip in UT/LT
!  - parcel assumed saturated at detrainment
!  - In this call to t_from_km, the tp guess is used to calculate
!    rvsp, so the max RH is with respect to the background atmosphere.
!------------------------------------------------------------------------
  rv_max = rh_max_detrain*rsat(tt(i,j,kk),ppp(i,j,kk))
  if (rvp > rv_max) rvp = rv_max
  rip = 0.
  rlp = rtp - rvp
  if (rlp < 0.) rlp = 0.
  zzz = zz(i,j,kk)
  kmp = hmp - (1. + rvp)*g*zzz  ! kmp when moved to LT/UT
  call t_from_km_allow_sat(kmp,rvp,ppp(i,j,kk),tp)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (kk == 2) then
!    print*,'----- rtp rv_max rlp = ',rtp,rv_max,rlp
!    print*,'tt(i,j,kk) = ',tt(i,j,kk)
!    print*,'tt(i,j,1) = ',tt(i,j,1)
!    print*,'ppp(i,j,1) = ',ppp(i,j,1)
!    print*,'ppp(i,j,kk) = ',ppp(i,j,kk)
!    print*,'rsat(tt(i,j,1),ppp(i,j,1)) = ',rsat(tt(i,j,1),ppp(i,j,1))
!    print*,'rsat(tt(i,j,2),ppp(i,j,2)) = ',rsat(tt(i,j,2),ppp(i,j,2))
!    stop
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  rvsp = rsat(tt(i,j,kk),ppp(i,j,kk))
  if (kk == 3) then 
!    print*,'--------------'
!    print*,'tt(i,j,kk) = ',tt(i,j,kk)
!    print*,'zz(i,j,kk) = ',zz(i,j,kk)
!    print*,'ppp(i,j,kk) = ',ppp(i,j,kk)
!    print*,'kmp = ',kmp
!    print*,'hmp = ',hmp
!    print*,'rvp = ',rvp
!    print*,'rlp = ',rlp
!    print*,'rip = ',rip
!    print*,'tp RH at detrainment = ',tp,rvp/rvsp
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'parcel T at new zz t_p = ',tp
    print*,'parcel height = ',zzz
    print*,'Background T tt(i,j,kk) = ',tt(i,j,kk)
  endif
!------------------------------------------------------------------------
!
! Find BL mass
!
!------------------------------------------------------------------------
  dpp = (ph(1) - ph(2)) + g*mextra(i,j,1)
  dmdry_bl = (ph(1) - ph(2))/(g*(1.0+rv_bl)) + mextra(i,j,1)
!------------------------------------------------------------------------
!  Need f_dp for CG -> DP conversion
!------------------------------------------------------------------------
  f_colrh_dp(i,j) = sigmoidal(colrh(i,j),colrh_dp_half,colrh_dp_scale,colrh_dp_low,colrh_dp_add)
  f_caput(i,j) = sigmoidal(cape_ut(i,j),caput_half,caput_scale,caput_low,caput_add)
  f_dp(i,j) = f_colrh_dp(i,j)*f_caput(i,j)
!------------------------------------------------------------------------
!
!  DP
!
!  wconv_bl ~ 3/50 ~ 0.06  (3 mm/day divided by timesteps/day)
!  maxx f_dp ~ 0.03  (since multiplied by tstep/tau)
!
!------------------------------------------------------------------------
  if (kk == 3) then 
!------------------------------------------------------------------------
! Calculate mass_dp(i,j)  DPPPP
! Should be a lagged relation between DP mass and vertical motion, in phase with colrh
!------------------------------------------------------------------------
  ip1 = i + 1
  if (ip1 == nlon+1) ip1 = 1
  im1 = i - 1
  if (im1 == 0) im1 = nlon
  xxx_def = -mass_lt_bl(i,j) - sub_lt_bl(i,j)  ! DEF
!  xxx_ip1 = -mass_lt_bl(ip1,j) - sub_lt_bl(ip1,j)  ! destroys MJO
!  xxx_im1 = -mass_lt_bl(im1,j) - sub_lt_bl(im1,j) ! hmm error
!  xxx = 0.5*(xxx_im1 + xxx_def)
!  xxx = 0.65*xxx_ip1 + 0.35*xxx_def
  xxx = xxx_def
  mass_dp(i,j) = mass_dp_old(i,j) 
  if (xxx > 0.) then 
    mass_dp(i,j) = mass_dp_old(i,j) + dp_inc*xxx*f_dp(i,j)
    dp_prod(i,j) = dp_inc*xxx*f_dp(i,j)  ! diagnostics
  else
    yyy =  dp_dec*xxx*(1.0-f_dp(i,j))
    if (abs(yyy) > mass_dp_old(i,j)) yyy = -mass_dp_old(i,j)
    mass_dp(i,j) = mass_dp_old(i,j) + yyy
    dp_prod(i,j) = yyy
  endif
  if (mass_dp(i,j) < 0.) mass_dp(i,j) = 0.
  dmdry_cn = mass_dp(i,j)
!------------------------------------------------------------------------
!  endif for dp
!------------------------------------------------------------------------
  endif
!------------------------------------------------------------------------
!
!  CG Congestus
!
!------------------------------------------------------------------------
  if (kk == 2) then 
!------------------------------------------------------------------------
! Calculate CG properties (for diagnostics even if kk = 2)
!------------------------------------------------------------------------
  f_colrh_cg(i,j) = sigmoidal(colrh(i,j),colrh_cg_half,colrh_cg_scale,colrh_cg_low,colrh_cg_add)
  f_caplt(i,j) = sigmoidal(cape_lt(i,j),caplt_half,caplt_scale,caplt_low,caplt_add)
  f_cg(i,j) = f_colrh_cg(i,j)*f_caplt(i,j)  ! Don't need tstep for CG since prop to mass_lt_bl
!------------------------------------------------------------------------
!  Dyn production: 
!------------------------------------------------------------------------
  ip1 = i + 1
  if (ip1 == nlon+1) ip1 = 1
  im1 = i - 1
  if (im1 == 0) im1 = nlon
  xxx_def = -mass_lt_bl(i,j) - sub_lt_bl(i,j)  ! DEF
!  xxx_ip1 = -mass_lt_bl(ip1,j) - sub_lt_bl(ip1,j)  ! destroys MJO
!  xxx_im1 = -mass_lt_bl(im1,j) - sub_lt_bl(im1,j) ! hmm error
!  xxx = 0.5*(xxx_im1 + xxx_def)
!  xxx = 0.65*xxx_ip1 + 0.35*xxx_def
  xxx = xxx_def
  if (xxx > 0.) then 
    mass_cg(i,j) = cg_amp_d*xxx*f_cg(i,j)
  else
    mass_cg(i,j) = 0.
  endif
  mass_cg(i,j) = mass_cg(i,j) + cg_amp_c*f_cg(i,j)
  dmdry_cn = mass_cg(i,j)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  endif
!------------------------------------------------------------------------
!  diagnostics
!------------------------------------------------------------------------
    nprint = 1
    if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'-------- calculated dmdry_cn for cg -----------------'
    print*,'initial mass_cg(i,j) = ',mass_cg(i,j)
    print*,'colrh(i,j) = ',colrh(i,j)
    print*,'f_cg = ',f_cg(i,j)
    print*,'dmdry_bl = ',dmdry_bl
    print*,'dmdry_cn = ',dmdry_cn
    Print*,'mass_lt_bl(i,j) = ',mass_lt_bl(i,j)
    print*,'f_cg = ',f_cg(i,j)
    print*,'f_caplt = ',f_caplt(i,j)
    print*,'f_caput = ',f_caput(i,j)
    print*,'f_colrh_cg = ',f_colrh_cg(i,j)
    print*,'f_colrh_dp = ',f_colrh_dp(i,j)
    print*,'----------------'
    endif
!------------------------------------------------------------------------
!
!  END FOR DP OR CG
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!  Check limits
!------------------------------------------------------------------------
!if (nstep == 24) print*,'checking limits'
  if (dmdry_cn < 0.) dmdry_cn = 0.
  if (dmdry_cn > dmdry_bl_fract*dmdry_bl) then
!    print*,'------ WARNING: dmdry_cn too large a fraction of BL dmdry -----'
    print*,'W:dmdry_cn big'
!    print*,'kk = ',kk
!    print*,'i j = ',i,j
!    print*,'dmdry_cn = ',dmdry_cn
!    print*,'dmdry_bl = ',dmdry_bl
!    print*,'f_dp = ',f_dp(i,j)
!    print*,'f_cg = ',f_cg(i,j)
!    Print*,'mass_lt_bl(i,j) = ',mass_lt_bl(i,j)
!    print*,'cape_ut(i,j) = ',cape_ut(i,j)
!    print*,'cape_lt(i,j) = ',cape_lt(i,j)
!    print*,'f_caput = ',f_caput(i,j)
!    print*,'f_caplt = ',f_caplt(i,j)
!    print*,'f_colrh_dp = ',f_colrh_dp(i,j)
!    print*,'f_colrh_cg = ',f_colrh_cg(i,j)
!    print*,'mass_dp(i,j) = ',mass_dp(i,j)
!    print*,'mass_cg(i,j) = ',mass_cg(i,j)
    dmdry_cn = dmdry_bl_fract*dmdry_bl
!    print*,'--------- dmdry_cn reduced ------------'
!    stop '--------- dmdry_cn too big ------------'
  endif
!------------------------------------------------------------------------
! Remove air from BL
!------------------------------------------------------------------------
  mextra(i,j,1) = mextra(i,j,1) - dmdry_cn
!------------------------------------------------------------------------
!  Check for negative dmdry_cn
!------------------------------------------------------------------------
  if ((dmdry_cn < 0.).or.(ph(1) < ph(2))) then
    print*,'*********  start of deep convect ********'
    print*,'dmdry_cn = ',dmdry_cn
    print*,'ph(1) ',ph(1)
    print*,'ph(2) ',ph(2)
    print*,'p1 = ',p1
    print*,'rv_bl = ',rv_bl
    print*,'kmp = ',kmp
    print*,'tp = ',tp
    print*,'zzz = ',zzz
    stop 'dmdry_cn negative or ps too low at start of deep convection'
  endif
!------------------------------------------------------------------------
!  DC: Produce rain from condensate
!  - This is the rain from assuming updraft parcel is saturated at detrainment.
!------------------------------------------------------------------------
  if ((rip + rlp) < 0.) stop 'negative rain_cn'
  rain_cn = dmdry_cn*(rip+rlp)
!------------------------------------------------------------------------
!  Check for NaN
!------------------------------------------------------------------------
  if (rain_cn /= rain_cn) then
    print*,'rain_cn = ',rain_cn 
    print*,'dmdry_cn = ',dmdry_cn
    print*,'rip = ',rip
    print*,'rlp = ',rlp
    print*,'xxx = ',xxx
    print*,'Initial UT T tt(i,j,3) = ',tt(i,j,3) 
    print*,'T of detrained parcel tp = ',tp
    print*,'rvp = ',rvp
    print*,'rv(i,j,1) = ',rv(i,j,1)
    print*,'rtp = ',rtp
    print*,'rlp = ',rlp
    print*,'rip = ',rip
    print*,'rvsp = ',rvsp
    print*,'before rain hm(i,j,1) = ',hm(i,j,1)
    print*,'after rain hmp = ',hmp
    print*,'hmp = ',hmp
    print*,'kmp = ',kmp
    print*,'dmdry_cn = ',dmdry_cn
    print*,'zzz = ',zzz
    stop 'rain_cn not a number'
  endif
!  if (nstep == 24) print*,'Initial rain_cn = ',rain_cn
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!if (nstep == 24) print*,'rain mass'
  if (kk == 3) then
    rain_dp(i,j) = rain_cn
    mass_dp(i,j) = dmdry_cn
!    print*,'Initial rain_dp(i,j) = ',rain_dp(i,j)
  else
    rain_cg(i,j) = rain_cn
    mass_cg(i,j) = dmdry_cn
  endif
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
  if (rain_dp(i,j) > 10.) then
    print*,'------  Large rain_dp  in convect ---'
    print*,'rain_dp(i,j) in mm/day = ',rain_dp(i,j)*(sec_day/tstep)
!    stop '--- really big rain_dp -----'
  endif
!------------------------------------------------------------------------
!  DC: Define new updraft parcel rtp/hmp after condensate removal 
!  Assume all rain is water.
!  hmp is the original hm in the BL
!  can't use tp here since not update/accurate.
!------------------------------------------------------------------------
  km_wat = cl*tt(i,j,kk)
  hmp = hmp - (rip+rlp)*(km_wat + g*zz(i,j,kk))
  rtp = rvp
  rlp = 0.
  rip = 0.
!------------------------------------------------------------------------
!  Check tp after condensate removal
!------------------------------------------------------------------------
  if (tp > 350.) then
    stop 'tp too hot after condensate removal'
  endif
!------------------------------------------------------------------------
!  DC: Detrainment Diagnostics
!------------------------------------------------------------------------
!if (nstep == 24) print*,'Calling rsat tp = ',tp
  rvsp = rsat(tp,pp(kk))
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
  print*,'------  Detrainment Diagnostics -----'
  print*,'Initial grid point T tt(i,j,kk) = ',tt(i,j,kk) 
  print*,'T of detrained parcel tp = ',tp
  print*,'rvp = ',rvp
  print*,'rtp = ',rtp
  print*,'rlp = ',rlp
  print*,'rip = ',rip
  print*,'rvsp = ',rvsp
  print*,'before rain hm(i,j,1) = ',hm(i,j,1)
  print*,'after rain hmp = ',hmp
  print*,'kmp = ',kmp
  print*,'dmdry_cn = ',dmdry_cn
  print*,'rain_dp(i,j) in mm/day = ',rain_dp(i,j)*sec_day
  endif
!------------------------------------------------------------------------
!  DC: Mix detrained parcel with detrainment grid cell and determine new UT parcel properties.
!------------------------------------------------------------------------
  rv_initial = rv(i,j,kk)
  hm_initial = hm(i,j,kk)
  dmdry_initial = (ph(kk) - ph(kk+1))/(g*(1.0+rv_initial)) + mextra(i,j,kk)
  dmdry_final = dmdry_cn + dmdry_initial
  rv_final = (rvp*dmdry_cn + rv_initial*dmdry_initial)/dmdry_final
  rt_final = rv_final
  rl_final = 0.
  ri_final = 0.
  hm_final = (hmp*dmdry_cn + hm_initial*dmdry_initial)/dmdry_final
  km_final = hm_final - (1. + rv_final)*g*zzz
  tp_final = tt(i,j,kk) !  T have initial guess
!  if (nstep == 24) print*,'rv_final = ',rv_final
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'------- After detrainment in convect --'
    print*,'initial grid point T = ',tt(i,j,kk)
    print*,'grid point tp_final = ',tp_final
  endif
!------------------------------------------------------------------------
!   Determine new final water partitioning: allow supersaturation for lag in anvil prod
!   - Do not produce any rain here.
!------------------------------------------------------------------------
!  call t_from_km(km_final,rt_final,pp(kk),rv_final,rl_final,ri_final,tp_final)
  call t_from_km_allow_sat(km_final,rv_final,ppp(i,j,kk),tp_final)
!  rain_cn = dmdry_final*(ri_final+rl_final)  ! turn off rain production
  rt_final = rv_final
  rl_final = 0.
  ri_final = 0.
  hm_final = hmoist(tp_final,rt_final,rv_final,ri_final,zzz)
!  if (nstep == 24) print*,'hm_final = ',hm_final
!------------------------------------------------------------------------
!  Assign final grid cell values at detrainment
!------------------------------------------------------------------------
  rv(i,j,kk) = rv_final
  tt(i,j,kk) = tp_final
  hm(i,j,kk) = hm_final
  if (rv(i,j,kk) /= rv(i,j,kk)) stop 'rv(i,j,kk) NaN'
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if (rv_final < 0.) then
    print*,'rv_final = ',rv_final
    stop '------rv_final negative in convect --------'
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'------- After supersaturation check in convect --'
    print*,'Additional anvil rain from supersaturation = ',rain_cn
    print*,'Final T tp_final = ',tp_final
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
   if (tp_final > 350.) then
     print*,'tp_final = ',tp_final
     stop 'tp_final coming out of t_from_km in deep convection too hot'
   endif
!------------------------------------------------------------------------
!  DC: Adjust mextra(i,j,kk) in UT or LT
!------------------------------------------------------------------------
   dm_want = (ph(kk) - ph(kk+1))/g  
   dm_final = (1.0 + rv_final)*dmdry_final  ! total mass after mixing
   dm_extra = dm_final - dm_want
   mextra(i,j,kk) = dm_extra/(1.0 + rv_final)
!------------------------------------------------------------------------
!   Diagnostics
!------------------------------------------------------------------------
   if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
     print*,'mextra(i,j,kk) = ',mextra(i,j,kk)
   endif
!------------------------------------------------------------------------
!  DC: Check for negative mextra at UT level (should not happen)
!------------------------------------------------------------------------
!   if (mextra(i,j,3) < 0.) then
   if (mextra(i,j,3) < -1000000.) then
     print*,'mextra(i,j,kk) = ',mextra(i,j,kk)
     print*,'dmdry_cn = ',dmdry_cn
     print*,'kk = ',kk
     stop 'mextra in UT negative after convection'
   endif
!------------------------------------------------------------------------
!  DC: Calculate colwat after convection: colwat_af_dp
!  - If water is conserved here, then deep convection should be conserving water
!------------------------------------------------------------------------
  if ((i == i_want).and.(j == j_want).and.(nprint_dp == 1)) then
    print*,'------------------------'
    print*,'colwat_af_cn = ',colwat_af_cn
    print*,'Change in colwat available for precip = ',colwat_bf_cn - colwat_af_cn
    print*,'rain_dp(i,j) = ',rain_dp(i,j)
    print*,'water dp error = ',colwat_bf_cn - colwat_af_cn - rain_cn
  endif
!------------------------------------------------------------------------
!  endif for mass_cg or mass_dp larger than threshold 
!------------------------------------------------------------------------
! endif
!------------------------------------------------------------------------
!  DC: Colwat after convection
!------------------------------------------------------------------------
  colwat_af_cn = 0.
  dmdry_af_cn = 0.
  do k = 1,nlev
    dmdry = (ph(k) - ph(k+1))/(g*(1.0+rv(i,j,k))) + mextra(i,j,k)
    colwat_af_cn = colwat_af_cn + rv(i,j,k)*dmdry
    dmdry_af_cn = dmdry_af_cn + dmdry
    if (k == 1) dmdry_bl_af_cn = dmdry
  end do
!------------------------------------------------------------------------
!  convect: Stop if conservation error
!------------------------------------------------------------------------
  if (kk == 3) colwat_error = colwat_af_cn - colwat_bf_cn + rain_dp(i,j)
  if (kk == 2) colwat_error = colwat_af_cn - colwat_bf_cn + rain_cg(i,j)
  dmdry_error = dmdry_af_cn - dmdry_bf_cn
!  if (nstep == 24) print*,'dmdry_error = ',dmdry_error
  if ((abs(colwat_error) > colwat_tolerance).or.(abs(dmdry_error) > dmdry_tolerance)) then
     print*,'-----------  conservation problem in deep convection ---'
     print*,'kk = ',kk
     print*,'colwat_error = ',colwat_error
     print*,'dmdry_error = ',dmdry_error
     print*,'colwat_bf_cn = ',colwat_bf_cn
     print*,'colwat_af_cn = ',colwat_af_cn
     print*,'dmdry_af_cn = ',dmdry_af_cn
     print*,'dmdry_bf_cn = ',dmdry_bf_cn
     print*,'rain_cg(i,j) = ',rain_cg(i,j)
     print*,'rain_dp(i,j) = ',rain_dp(i,j)
     print*,'dmdry_cn = ',dmdry_cn
     stop 'conservation error after deep convection'
  endif
!------------------------------------------------------------------------
!
!  convect: Check for BL mass error
!
!------------------------------------------------------------------------
  if (kk == 3) then
     dmdry_bl_error = dmdry_bl_af_cn - dmdry_bl_bf_cn + mass_dp(i,j)
  else
     dmdry_bl_error = dmdry_bl_af_cn - dmdry_bl_bf_cn + mass_cg(i,j)
  endif
  if (abs(dmdry_bl_error) > 1.0E-08) then
    print*,'----------  dmdry_bl_error in convect ---'
    print*,'dmdry_cn = ',dmdry_cn
    print*,'kk = ',kk
    print*,'dmdry_bl_bf_cn = ',dmdry_bl_bf_cn
    print*,'dmdry_bl_af_cn = ',dmdry_bl_af_cn
    print*,'dmdry_bl_error  = ',dmdry_bl_error
    print*,'mass_dp(i,j) = ',mass_dp(i,j)
    print*,'mass_cg(i,j) = ',mass_cg(i,j)
    print*,'dmdry_cn = ',dmdry_cn
    print*,'f_dp = ',f_dp(i,j)
    print*,'f_cg = ',f_cg(i,j)
    stop '------------  dmdry_bl_error in convect  -----------------'
  endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
!  mdp_error = mdp_initial = mass_dp(i,j)
!  mcg_error = mcg_initial = mass_cg(i,j)
!------------------------------------------------------------------------
!  end loop over nlon,nlat
!------------------------------------------------------------------------
end do
end do
print*,'finished convect'
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine convect
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!  **********************************************************************
!  ****************  HT: height  ***************************************
!  **********************************************************************
!
!  Calculated zz: mean height of each layer. Needed for hmoist
!  - also ppp,pt
!
!------------------------------------------------------------------------
subroutine height
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!   Local variables
!
!------------------------------------------------------------------------
integer :: i,j,k
double precision :: ztop,zbot,pbot,ptop,rv_in,hm_in,tt_out,xxx,zmid,pmid
!------------------------------------------------------------------------
!
!  Define pph - pressure at boundaries
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
  xxx = mextra(i,j,1) + mextra(i,j,2) + mextra(i,j,3)
  pph(i,j,1) = ph(1) + g*xxx
  xxx = mextra(i,j,2) + mextra(i,j,3)
  pph(i,j,2) = ph(2) + g*xxx
  xxx = mextra(i,j,3)
  pph(i,j,3) = ph(3) + g*xxx
  pph(i,j,4) = ph(4)
end do
end do
!------------------------------------------------------------------------
!
!  Define ppp - pressure at midpoints
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
do k = 1,nlev
  ppp(i,j,k) = 0.5*(pph(i,j,k) + pph(i,j,k+1))
end do
end do
end do
!------------------------------------------------------------------------
!
!  Get zz and tt at midpoints consistent with hm,rv, and hyd balance.
!
!------------------------------------------------------------------------
do i = 1,nlon
do j = 1,nlat
do k = 1,nlev
  if (k == 1) zh(i,j,k) = 0.
  pbot = pph(i,j,k)
  ptop = pph(i,j,k+1)
  pmid = ppp(i,j,k)
  zbot = zh(i,j,k)
  hm_in = hm(i,j,k)
  rv_in = rv(i,j,k)
  tt_out = tt(i,j,k)
!  print*,'i j k = ',i,j,k
  call t_from_hm(pbot,pmid,ptop,zbot,zmid,ztop,hm_in,rv_in,tt_out,i,j)
  tt(i,j,k) = tt_out
  zz(i,j,k) = zmid
  zh(i,j,k+1) = ztop
  pt(i,j,k) = tt(i,j,k)*(100000./ppp(i,j,k))**(rdd/cpd)
end do
end do
end do
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine height
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  Given the hm/rv of a pressure layer, find the mean T and width of the
!    layer using hydrostatic balance and the definition of hm.
!
!  Procedure:
!  - make initial guess of T based on hm/rv/zbot (using km and t_from_km)
!  - use virtual T to get density
!  - use density,hyd relation, and (pbot-ptop) to get dz
!  - use dz to get z_av
!  - use z_av to get hm value
!  - determine error in hm.
!  - iterate
!  - use hydrostatic relation
!  - Define range of allowed T and iterate to solution
!  
!------------------------------------------------------------------------
subroutine t_from_hm(pbot,pmid,ptop,zbot,zmid,ztop,hm_in,rv_in,tt_out,i,j)
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
double precision, intent(in) :: pbot
double precision, intent(in) :: pmid
double precision, intent(in) :: ptop
double precision, intent(in) :: zbot
double precision, intent(out) :: zmid
double precision, intent(out) :: ztop
double precision, intent(in) :: hm_in
double precision, intent(in) :: rv_in
double precision, intent(inout) :: tt_out
integer, intent(in) :: i
integer, intent(in) :: j
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: niter,iter
double precision :: dtt,rt_in,ri_in,dpp,kmm,t1,t2,tmid,hmerr,tv,rho,dz
double precision :: hmerr1,hmerr2,zzz,hmm,xxx,hmerr_new
!------------------------------------------------------------------------
!  Diagnostics
!------------------------------------------------------------------------
!print*,'--------  incoming values ---'
!print*,'incoming pbot = ',pbot
!print*,'incoming ptop = ',ptop
!print*,'incoming zbot = ',zbot
!print*,'incoming hm_in = ',hm_in
!print*,'incoming rv_in = ',rv_in
!print*,'incoming tt_out = ',tt_out
if (zbot /= zbot) stop 'zbot NaN problem at start of t_from_hm'
if (rv_in /= rv_in) stop 'rv_in NaN problem at start of t_from_hm'
if (hm_in /= hm_in) stop 'hm_in NaN problem at start of t_from_hm'
!------------------------------------------------------------------------
!
! Set tmid equal to old tt value
!
!------------------------------------------------------------------------
  hmerr_new = 100.  ! arbitrary initial large value so enters loop
  dtt = dtt_start
  dpp = pbot - ptop
  tmid = tt_out
  rt_in = rv_in
  ri_in = 0.
!------------------------------------------------------------------------
!
!   Start iteration
!
!------------------------------------------------------------------------
do iter = 1,niter_hm
if (hmerr_new > hmerr_want) then
!------------------------------------------------------------------------
! Define bracketing T
!------------------------------------------------------------------------
!  print*,'iter = ',iter
!  print*,'hmerr = ',hmerr
!  print*,'tmid = ',tmid
  t1 = tmid - 0.5*dtt
  t2 = tmid + 0.5*dtt
!------------------------------------------------------------------------
!
!   Assume t1: get hmerr1
!
!------------------------------------------------------------------------
  tv = t1*(1. + epsi*rv_in)/(1. + rv_in)
  rho = pmid/(rdd*tv)
  dz = dpp/(g*rho)
  zzz = zbot + 0.5*dz
  hmm = hmoist(t1,rt_in,rv_in,ri_in,zzz)
!  hm(i,j,k) = hmoist(tt(i,j,k),rtp,rv(i,j,k),rip,zz(i,j,k))
  hmerr1 = hmm-hm_in
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
!  print*,'t1 = ',t1
!  print*,'dz = ',dz
!  print*,'hmm = ',hmm
!  print*,'hmerr1 = ',hmerr1
!  print*,'zzz = ',zzz
!  print*,'tv = ',tv
!  print*,'dpp = ',dpp
!------------------------------------------------------------------------
!
!   Assume t2: get hmerr2
!
!------------------------------------------------------------------------
  tv = t2*(1. + epsi*rv_in)/(1. + rv_in)
  rho = pmid/(rdd*tv)
  dz = dpp/(g*rho)
  zzz = zbot + 0.5*dz
  hmm = hmoist(t2,rt_in,rv_in,ri_in,zzz)
  hmerr2 = hmm-hm_in
!------------------------------------------------------------------------
!  Check that hmerr of different sign
!  - should have hmerr1 < 0. (since t1 < real T)
!  - should have hmerr2 > 0. (since t2 > real T)
!------------------------------------------------------------------------
  if (((hmerr1*hmerr2) > 0.).and.(nprint_hmerr == 1)) then 
    print*,'----------  hmerr problem ----'
    print*,'pmid = ',pmid
    print*,'hmm for t2 = ',hmm
    print*,'hm_in = ',hm_in
    print*,'iter = ',iter
    print*,'hmerr1 for t1 = ',hmerr1
    print*,'hmerr2 for t2 = ',hmerr2
    print*,'t1 = ',t1
    print*,'t2 = ',t2
    print*,'tmid = ',tmid
    print*,'dz = ',dz
    print*,'incoming pbot = ',pbot
    print*,'incoming ptop = ',ptop
    print*,'incoming zbot = ',zbot
    print*,'incoming hm_in = ',hm_in
    print*,'incoming rv_in = ',rv_in
    print*,'incoming tt_out = ',tt_out
    stop 'hmerr of the same sign'
  endif
!------------------------------------------------------------------------
!
!  Use hmerr to get new tmid
!  - xxx should be > 0.
!------------------------------------------------------------------------
  if ((hmerr1*hmerr2) < 0.) then 
    xxx = -hmerr1/(hmerr2-hmerr1)
    tmid = t1 + xxx*dtt
    dtt = dtt_reduce*dtt
  else
    if (abs(hmerr1) < abs(hmerr2)) then
      tmid = t1
    else
      tmid = t2
    endif
  endif
!------------------------------------------------------------------------
!  Find new hmerr/dz/zz from best tmid
!------------------------------------------------------------------------
  tv = tmid*(1. + epsi*rv_in)/(1. + rv_in)
  rho = pmid/(rdd*tv)
  dz = dpp/(g*rho)
  zzz = zbot + 0.5*dz
  hmm = hmoist(tmid,rt_in,rv_in,ri_in,zzz)
  hmerr_new = abs(hmm-hm_in)
  if (hmerr_new < hmerr_want) goto 60
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
if (hmerr_new > hmerr_want) then
!  print*,'-------  in loop --------'
!  print*,'xxx = ',xxx
!  print*,'hmerr1 = ',hmerr1
!  print*,'hmerr2 = ',hmerr2
!  print*,'New hmerr = ',hmerr_new
endif
!------------------------------------------------------------------------
!
!  end iterative loop
!
!------------------------------------------------------------------------
endif
end do
!------------------------------------------------------------------------
!
!  Check for convergence
!
!------------------------------------------------------------------------
if (hmerr > hmerr_want) then
  stop 'convergence in t_from_hm not achieved' 
endif
!------------------------------------------------------------------------
!
!  Assign output values
!
!------------------------------------------------------------------------
60 continue
ztop = zbot + dz
zmid = zzz
tt_out = tmid
!------------------------------------------------------------------------
!  Check for unrealistic values of hmm
!------------------------------------------------------------------------
if (hmm > 410000.) then
  print*,'--------- unrealistic hmm in t_from_hm ----'
  print*,'hmm = ',hmm
  print*,'tmid = ',tmid
  print*,'rv_in = ',rv_in
  print*,'rt_in = ',rt_in
  print*,'zzz = ',zzz
  print*,'i = ',i
  print*,'j = ',j
  stop '--- bad hmm ----'
endif
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
!print*,'Number of iterations = ',iter
!print*,'------  outgoing values ---'
!print*,'outgoing ztop = ',ztop
!print*,'outgoing zmid = ',zmid
!print*,'outgoing tt_out = ',tt_out
!print*,'hmerr_new = ',hmerr_new
!stop
!------------------------------------------------------------------------
!
!  End subroutine t_from_hm
!
!------------------------------------------------------------------------
end subroutine t_from_hm
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  Balance condition  (see Pu and Dickinson)
!
!  du/dt = f*v - (1/rho)*dp/dx - eps*u 
!  dv/dt  = -f*u - (1/rho)*dp/dy - eps*v  
!
! f = 2*omega(sin(phi)
! Assume du/dt,dv/dt = 0
! Know rho,dp/dx,eps,f
!
! (1)   f*v - (1/rho)*dp/dx - eps*u = 0
! (2)  -f*u - (1/rho)*dp/dy - eps*v = 0
!
!  From (1) eps*u = f*v - (1/rho)*dp/dx
!           u = (f/eps)*v - (1/eps)*(1/rho)*dp/dx
!  Plug into 2:
!     -f*[(f/eps)*v - (1/eps)*(1/rho)*dp/dx] - (1/rho)*dp/dy - eps*v = 0
!     -(f**2/eps)*v + (f/eps)*(1/rho)*dp/dx - (1/rho)*dp/dy - eps*v = 0
!     -(f**2/eps)*v + (f/eps)*(1/rho)*dp/dx - (1/rho)*dp/dy - eps*v = 0
!     v*[eps + (f**2/eps)] = (f/eps)*(1/rho)*dp/dx - (1/rho)*dp/dy
!     v*{rho*[eps + (f**2/eps)]} = (f/eps)*dp/dx - dp/dy
!     Define AA = rho*[eps + (f**2/eps)]
!     v*AA = (f/eps)*dp/dx - dp/dy
!
!  Solutions:
!     v = (1/AA)*{(f/eps)*dp/dx - dp/dy}
!     u = (f/eps)*v - (1/eps)*(1/rho)*dp/dx
!
!------------------------------------------------------------------------
subroutine friction
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
use model
implicit none
!------------------------------------------------------------------------
!
!   Local integer variables
!
!------------------------------------------------------------------------
integer :: i,j,k,jp1,jm1,im1,ip1
double precision :: ps(nlon_mjo,nlat,nlev),ff,pmid,tv,rho,dpdx,dpdy,AA,dx,dy
!------------------------------------------------------------------------
!
!  Define ps at surface (in Pa)
!  convert kg/m2 to N/m2
!
!------------------------------------------------------------------------
do i = 1,nlon_mjo
do j = 1,nlat
    ps(i,j,1) = g*mjo_mean(i,j,i_mextra_all) 
    ps(i,j,2) = g*(mjo_mean(i,j,i_mextra_ut) + mjo_mean(i,j,i_mextra_lt))
    ps(i,j,3) = g*mjo_mean(i,j,i_mextra_ut)
end do
end do
!------------------------------------------------------------------------
!
!     Define AA = rho*[frict + (f**2/frict)]
!     v = (1/AA)*{(f/frict)*dp/dx - dp/dy}
!     u = (f/frict)*v - (1/frict)*(1/rho)*dp/dx
!
!
!------------------------------------------------------------------------
do k = 1,nlev
do i = 1,nlon_mjo
do j = 1,nlat
  ff = 2.0*omega*sind(lat(j))
  pmid = ppp(i,j,1)
  tv = tt(i,j,1)*(1. + epsi*rv(i,j,1))/(1. + rv(i,j,1))
  rho = pmid/(rdd*tv)
  AA = rho*(frict + ((ff*ff)/frict))
  dx = (2.0*pi*re/float(nlon))
  dy = (2.0*pi*re/float(nlat-1))*(2.0*dlat/360.0)
!  dpdx
  if (i == 1) then
    dpdx = (ps(i+1,j,k) - ps(i,j,k))/dx
  elseif (i == nlon_mjo) then
    dpdx = (ps(i,j,k) - ps(i-1,j,k))/dx
  else
    dpdx = (ps(i+1,j,k) - ps(i-1,j,k))/(2.0*dx)
  endif
!  dpdy
  if (j == 1) then
    dpdy = (ps(i,j+1,k) - ps(i,j,k))/dy
  elseif (j == nlat) then
    dpdy = (ps(i,j,k) - ps(i,j-1,k))/dy
  else
    dpdy = (ps(i,j+1,k) - ps(i,j-1,k))/(2.0*dy)
  endif
  vv(i,j,k) = (1./AA)*((ff/frict)*dpdx - dpdy)
  uu(i,j,k) = (ff/frict)*vv(i,j,k) - (1.0/frict)*(1/rho)*dpdx
end do
end do
end do
!------------------------------------------------------------------------
!
!  Do 3x3 averaging
!
!------------------------------------------------------------------------
print*,'start 3 by 3 averaging'
do k = 1,nlev
!------------------------------------------------------------------------
! uu
!------------------------------------------------------------------------
  do i = 1,nlon_mjo
  do j = 1,nlat
    im1 = i-1
    if (im1 == 0) im1 = nlon_mjo
    ip1 = i+1
    if (ip1 == nlon_mjo+1) ip1 = 1
    jm1 = j-1
    if (jm1 == 0) jm1 = nlat
    jp1 = j+1
    if (jp1 == nlat+1) jp1 = 1
    temp(i,j) = uu(i,jp1,k) + uu(im1,jp1,k) + uu(ip1,jp1,k) + &
                uu(i,j,k) + uu(im1,j,k) + uu(ip1,j,k) +       &
                uu(i,jm1,k) + uu(im1,jm1,k) + uu(ip1,jm1,k) 
    temp(i,j) = temp(i,j)/9.0
  end do
  end do
  do i = 1,nlon_mjo
  do j = 1,nlat
    uu(i,j,k) = temp(i,j)
  end do
  end do
!------------------------------------------------------------------------
! vv
!------------------------------------------------------------------------
  do i = 1,nlon_mjo
  do j = 1,nlat
    im1 = i-1
    if (im1 == 0) im1 = nlon_mjo
    ip1 = i+1
    if (ip1 == nlon_mjo+1) ip1 = 1
    jm1 = j-1
    if (jm1 == 0) jm1 = nlat
    jp1 = j+1
    if (jp1 == nlat+1) jp1 = 1
    temp(i,j) = vv(i,jp1,k) + vv(im1,jp1,k) + vv(ip1,jp1,k) + &
                vv(i,j,k) + vv(im1,j,k) + vv(ip1,j,k) +       &
                vv(i,jm1,k) + vv(im1,jm1,k) + vv(ip1,jm1,k) 
    temp(i,j) = temp(i,j)/9.0
  end do
  end do
  do i = 1,nlon_mjo
  do j = 1,nlat
    vv(i,j,k) = temp(i,j)
  end do
  end do
!------------------------------------------------------------------------
!  end loop over k
!------------------------------------------------------------------------
end do
!------------------------------------------------------------------------
!
!  End subroutine friction
!
!------------------------------------------------------------------------
end subroutine friction
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
