!```````````````````````````````````````````````````````````````````
module postp_comp_mod

!+

implicit none


contains


!```````````````````````````````````````````````````````````````````
subroutine comp_tempk(tempK, theta, PI)

!+ Computes temperature [K], given theta and normalized Exner function

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cp

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: theta, PI

real(f4), dimension(:,:,:), intent(OUT) :: tempK

!== Local declarations:

! NONE

!=== Executable statements:

tempK = theta*PI/cp 

end subroutine comp_tempk


!```````````````````````````````````````````````````````````````````
subroutine comp_SBPStempk(Q_tempK, density, cp)

!+ Computes SBPS temperature [K], given SBPS Q, density, and cp

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:,:), intent(IN) :: density, cp 

real(f4), dimension(:,:,:,:), intent(INOUT) :: Q_tempK

!== Local declarations:

! NONE

!=== Executable statements:

! Do not calculate nonphysical cases:
where (Q_tempK > 1.e-6_f4) Q_tempK = Q_tempK/(density*cp)

end subroutine comp_SBPStempk


!```````````````````````````````````````````````````````````````````
subroutine comp_press(pressure, PI)

!+ Computes pressure [Pa], given normalized Exner function

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cpor, p00, cp

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: PI

real(f4), dimension(:,:,:), intent(OUT) :: pressure

!== Local declarations:

! NONE

!=== Executable statements:

pressure  = p00*((pi/cp)**cpor)

end subroutine comp_press


!```````````````````````````````````````````````````````````````````
subroutine comp_nh_pressp(NHprp, theta, orig_PI, orig_press, rtgt, tgt_igd,  &
     ifn, i_AZ)

!+ Computes non-hydrostatic pressure perturbation [Pa]

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: g_areoid, p00, cp, cpor 
use postp_main_memory_mod, only: dzm, zt

implicit none

!=== Argument declarations:

integer, intent(IN) :: tgt_igd, ifn
real(f4), dimension(:,:,:), intent(IN) :: orig_press, orig_PI, theta
real(f4), dimension(:,:), intent(IN) :: rtgt
integer, dimension(:) :: i_AZ

real(f4), dimension(:,:,:), intent(OUT) :: NHprp

!== Local declarations:

real(f4), dimension(size(theta,1)) :: new_PI

real(f4) :: new_press
integer :: i, j, k, n1, kk

!=== Executable statements:

n1 = size(theta, 1)

do j=1,size(NHprp, 3)
   do i=1,size(NHprp, 2)

      ! Hydrostatic integration to obtain hydrostatic PI values:
      new_PI(n1) = orig_PI(n1,i,j)
      do k=n1,2,-1
         new_PI(k-1) = new_PI(k)+g_areoid*rtgt(i,j)/(dzm(k-1,tgt_igd,ifn)*  &
              0.5_f4*(theta(k,i,j)+theta(k-1,i,j)))
      end do

      ! Finish up:
      do k=1,size(NHprp, 1)
         kk = i_AZ(k)
         new_press = p00*((new_PI(kk)/cp)**cpor)   ! Convert new PI to pressure
         NHprp(k,i,j) = orig_press(kk,i,j)-new_press  ! Non-hydrostatic pressure
                                                      !  perturbation [Pa]
      end do

   end do
end do

end subroutine comp_nh_pressp


!```````````````````````````````````````````````````````````````````
subroutine comp_density(density, theta, PI)

!+ Computes "real" density [kg/m^3], given theta and normalized Exner function

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use body_specific_constants_mod, only: cp, p00, cpor, R_air

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: theta, PI

real(f4), dimension(:,:,:), intent(OUT) :: density

!== Local declarations:

real(f4) :: c1, c2

!=== Executable statements:

c1 = cpor-ONE_f4
c2 = p00*(cp**(-c1))/R_air

density = c2*(PI**c1)/theta

end subroutine comp_density


!```````````````````````````````````````````````````````````````````
subroutine comp_brunt_vaisala(en, theta, PI, rtg, tgt_igd, ifn, i_AZ)

!+ Computes the (dry) Brunt-Vaisala frequency

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4, ONE_f4
use body_specific_constants_mod, only: g_areoid
use postp_main_memory_mod, only: zt

implicit none

!=== Argument declarations:

integer, intent(IN) :: tgt_igd, ifn
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: theta, PI
real(f4), dimension(:,:), intent(IN) :: rtg

real(f4), dimension(:,:,:), intent(OUT) :: en

!== Local declarations:

integer :: i, j, k, kk
real(f4) :: en2
real(f4), dimension(size(theta,1)) :: vctr11, vctr12, vctr32, vctr1, vctr2

!=== Executable statements:

! Calculate potential temperature profile:

do j=1,size(theta, 3)
   do i=1,size(theta, 2)

      vctr11 = theta(:,i,j); vctr12 = theta(:,i,j); vctr32 = ZERO_f4

      do k=2,size(theta, 1)-1   ! g/dz_2kt
         vctr1(k) = g_areoid/((zt(k+1,tgt_igd,ifn)-zt(k-1,tgt_igd,ifn))*  &
              rtg(i,j))
      end do

      do k=2,size(theta, 1)-1
         !   g    [ dtheta_2kt    dvapterm_2kt ]
         ! -----  [  -------   -   ----------  ]
         ! dz_2kt [ theta_kt                   ]
         en2 = vctr1(k)*((vctr12(k+1)-vctr12(k-1))/vctr12(k)-  &
              (vctr32(k+1)-vctr32(k-1)))
         vctr2(k) = sign(ONE_f4, en2)*sqrt(abs(en2))
      end do
      vctr2(1) = vctr2(2)
      vctr2(size(theta, 1)) = vctr2(size(theta, 1)-1)

      do k=1,size(en, 1)
         kk = i_AZ(k)
         en(k,i,j) = vctr2(kk)
      end do

   end do
end do

end subroutine comp_brunt_vaisala


!```````````````````````````````````````````````````````````````````
subroutine comp_sfc_air_density(sfc_air_density, sfc_air_tempk, PI)

!+ Computes surface air density [kg/m^3], given surface air temperature and
!   normalized Exner function. Intended for surface interaction use (especially
!   aeolian), and will tend to slightly over(under)-estimate the "real" air
!   density value during the night(day). 

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use body_specific_constants_mod, only: cpi, p00, cpor, R_air

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: PI
real(f4), dimension(:,:,:), intent(IN) :: sfc_air_tempk

real(f4), dimension(:,:,:), intent(OUT) :: sfc_air_density

!== Local declarations:

integer :: i, j, ip
real(f4) :: c1, sfc_pi

!=== Executable statements:

c1 = p00/R_air

do ip=1,size(sfc_air_density, 3)
   do j=1,size(PI, 3)
      do i=1,size(PI, 2)
         sfc_pi = 0.5_f4*(PI(1,i,j)+PI(2,i,j))
         sfc_air_density(i,j,ip) = c1*((sfc_pi*cpi)**cpor)/sfc_air_tempk(i,j,ip)
      end do
   end do
end do

end subroutine comp_sfc_air_density


!```````````````````````````````````````````````````````````````````
subroutine comp_density_lin(density, dn0, pip_over_pi0, thp_over_th0)

!+ Computes "linearized" density [kg/m^3]

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cv, R_air

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: dn0, pip_over_pi0, thp_over_th0

real(f4), dimension(:,:,:), intent(OUT) :: density

!== Local declarations:

real(f4) :: c1

!=== Executable statements:

c1 = cv/R_air

density = dn0+(c1*pip_over_pi0-thp_over_th0)*dn0

end subroutine comp_density_lin


!```````````````````````````````````````````````````````````````````
subroutine comp_pert_over_base(pob, pert, base, pert_type)

!+ Computes quotient of var perturbation and var base state

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4

implicit none

!=== Argument declarations:

integer, intent(IN) :: pert_type
real(f4), dimension(:,:,:), intent(IN) :: pert, base

real(f4), dimension(:,:,:), intent(OUT) :: pob

!== Local declarations:

! NONE

!=== Executable statements:

select case (pert_type)
case(0)    ! pert is full variable
   pob = pert/base-ONE_f4
case(1)    ! pert is only perturbation of variable
   pob = pert/base
end select

end subroutine comp_pert_over_base


!```````````````````````````````````````````````````````````````````
subroutine comp_vertsum(var_vsum, var, rtgt, kmin, kmax, tgt_igd, ifn)

!+ Computes vertical sum (weighted by layer thickness) of a variable's values

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use postp_main_memory_mod, only: dzt

implicit none

!=== Argument declarations:

integer, intent(IN) :: kmin, kmax, tgt_igd, ifn
real(f4), dimension(:,:,:), intent(IN) :: var
real(f4), dimension(:,:), intent(IN) :: rtgt

real(f4), dimension(:,:), intent(OUT) :: var_vsum

!== Local declarations:

integer :: i, j, k
real(f4) :: vsum

!=== Executable statements:

do j=1,size(var, 3)
   do i=1,size(var, 2)
      vsum = ZERO_f4
      do k=kmin,kmax
         vsum = vsum+var(k,i,j)/dzt(k,tgt_igd,ifn)*rtgt(i,j) 
      end do
      var_vsum(i,j) = vsum
   end do
end do

end subroutine comp_vertsum


!```````````````````````````````````````````````````````````````````
subroutine comp_rtgX(topXin_rtgXout)

!+ Computes rtgX vertical stretch factors

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use postp_main_memory_mod, only: ztop

implicit none

!=== Argument declarations:

real(f4), dimension(:,:), intent(INOUT) :: topXin_rtgXout

!== Local declarations:

! NONE

!=== Executable statements:

topXin_rtgXout = ONE_f4-topXin_rtgXout/ztop 

end subroutine comp_rtgX


!```````````````````````````````````````````````````````````````````
subroutine comp_speed(ax, ay, speed)

!+ Computes vector magnitude, given x- and y-components

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: ax, ay

real(f4), dimension(:,:,:), intent(OUT) :: speed

!== Local declarations:

! NONE

!=== Executable statements:

speed = sqrt(ax**2+ay**2) 

end subroutine comp_speed


!```````````````````````````````````````````````````````````````````
subroutine comp_dir(ax, ay, dir)

!+ Computes vector direction [deg_from_N], given x- and y-components

use univ_kind_defs_mod, only: f4
use univ_convert_units_mod, only: rad_to_deg_f4
use univ_math_constants_mod, only: ZERO_f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: ax, ay

real(f4), dimension(:,:,:), intent(OUT) :: dir

!== Local declarations:

integer :: i, j, k
real(f4) :: u, v, ftmp

!=== Executable statements:

do j=1,size(ax, 3)
   do i=1,size(ax, 2)
      do k=1,size(ax, 1)
         ! Direction [deg_from_N]: 
         u = ax(k,i,j); v = ay(k,i,j)
         if (abs(u) > 1.e-5_f4) then
            ftmp = atan2(v, u)*rad_to_deg_f4
         else
            if (abs(v) > 1.e-5_f4) then
               ftmp = 180._f4-90._f4*(v/abs(v))
            else
               ftmp = 270._f4
            end if
         end if
         ftmp = 270._f4-ftmp

         if (ftmp < ZERO_f4) ftmp = ftmp+360._f4
         if (ftmp >= 360._f4) ftmp = ftmp-360._f4

         dir(k,i,j) = ftmp
      end do
   end do
end do

end subroutine comp_dir


!```````````````````````````````````````````````````````````````````
subroutine comp_avgu(u_raw, u_avg, i_AZ, i_X, i_Y)

!+ Averages u wind component values to "T" points

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_HALF_f4
use postp_main_memory_mod, only: undef_value_f4

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ, i_X, i_Y
real(f4), dimension(:,:,:), intent(IN) :: u_raw
real(f4), dimension(:,:,:), intent(OUT) :: u_avg

!== Local declarations:

integer :: i, j, k, ii, jj, kk

!=== Executable statements:

do j=1,size(u_avg, 3)
   jj = i_Y(j)
   do i=1,size(u_avg, 2)
      ii = i_X(i)
      if (ii > 1) then
         do k=1,size(u_avg, 1)
            kk = i_AZ(k)
            u_avg(k,i,j) = ONE_HALF_f4*(u_raw(kk,ii,jj)+u_raw(kk,ii-1,jj))
         end do
      else
         u_avg(:,i,j) = undef_value_f4
      end if
   end do
end do

end subroutine comp_avgu


!```````````````````````````````````````````````````````````````````
subroutine comp_avgv(v_raw, v_avg, i_AZ, i_X, i_Y)

!+ Averages v wind component values to "T" points

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_HALF_f4
use postp_main_memory_mod, only: undef_value_f4

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ, i_X, i_Y
real(f4), dimension(:,:,:), intent(IN) :: v_raw
real(f4), dimension(:,:,:), intent(OUT) :: v_avg

!== Local declarations:

integer :: i, j, k, ii, jj, kk

!=== Executable statements:

do j=1,size(v_avg, 3)
   jj = i_Y(j)
   if (jj > 1) then
      do i=1,size(v_avg, 2)
         ii = i_X(i)
         do k=1,size(v_avg, 1)
            kk = i_AZ(k)
            v_avg(k,i,j) = ONE_HALF_f4*(v_raw(kk,ii,jj)+v_raw(kk,ii,jj-1))
         end do
      end do
   else
      v_avg(:,:,j) = undef_value_f4
   end if
end do

end subroutine comp_avgv


!```````````````````````````````````````````````````````````````````
subroutine comp_avgw(w_raw, w_avg, i_AZ)

!+ Averages (vertically) staggered w wind component values to "T" points

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_HALF_f4
use univ_exception_mod, only: exception, FATAL
use univ_maxlen_string_mod, only: MAX_PROC_NAME_LEN
use postp_main_memory_mod, only: undef_value_f4

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: w_raw

real(f4), dimension(:,:,:), intent(OUT) :: w_avg

!== Local declarations:

character(len=MAX_PROC_NAME_LEN), parameter ::  &
     procedure_name = 'comp_avgw'

integer :: i, j, k, kk

!=== Executable statements:

do j=1,size(w_avg, 3)
   do i=1,size(w_avg, 2)
      if (i_AZ(1) == 1) w_avg(1,i,j) = undef_value_f4
      do k=1,size(w_avg, 1)
         kk = i_AZ(k)
         if (kk >= 2) w_avg(k,i,j) = ONE_HALF_f4*  &
              (w_raw(kk,i,j)+w_raw(kk-1,i,j)) 
      end do
   end do
end do

end subroutine comp_avgw


!```````````````````````````````````````````````````````````````````
subroutine comp_areo_rotate(ca, cb, xt, yt, i_X, i_Y, plat, plon)

!+ Rotates OPS vector components to true areographical orientation

use univ_kind_defs_mod, only: f4
use univ_OPS_proj_mod, only: xy_ll_OPS
use body_specific_constants_mod, only: radius
use polarst_mod, only: uvtoueve

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_X, i_Y
real(f4), dimension(:), intent(IN) :: xt
real(f4), dimension(:), intent(IN) :: yt
real(f4), intent(IN) :: plat,plon

real(f4), dimension(:,:,:), intent(INOUT) :: ca, cb

!== Local declarations:

integer :: i, j, k
real(f4) :: qlat, qlon, u, v

!=== Executable statements:

do j=1,size(ca, 3)
   do i=1,size(ca, 2)
      do k=1,size(ca, 1)
         call xy_ll_OPS(qlat, qlon, plat, plon, xt(i_X(i)), yt(i_Y(j)),  &
              radius)
         u = ca(k,i,j); v = cb(k,i,j)
         call uvtoueve(u, v, ca(k,i,j), cb(k,i,j), qlat, qlon, plat, plon)
      end do
   end do
end do

end subroutine comp_areo_rotate


!```````````````````````````````````````````````````````````````````
subroutine comp_sfc_T2th(theta_sfc, sfc_T, PI)

!+ Computes sfc potential temperature, given sfc PI and temperature

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cp

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: PI
real(f4), dimension(:,:,:), intent(IN) :: sfc_T

real(f4), dimension(:,:,:), intent(OUT) :: theta_sfc

!== Local declarations:

integer :: i, j, ip
real(f4) :: PI_sfc

!=== Executable statements:

do ip=1,size(sfc_T, 3)
   do j=1,size(sfc_T, 2)
      do i=1,size(sfc_T, 1)
         PI_sfc = (PI(1,i,j)+PI(2,i,j))*0.5_f4
         theta_sfc(i,j,ip) = sfc_T(i,j,ip)*cp/PI_sfc   ! Surface theta value
      end do
   end do
end do

end subroutine comp_sfc_T2th


!```````````````````````````````````````````````````````````````````
subroutine comp_TItvar_sfc(TItvar_sfc, sbps_density, sbps_thrmcap,  &
     sbps_thrmcnd, sbps_lyrthkn)

!+ Computes *surface* time-varying thermal inertia

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:,:), intent(IN) :: sbps_density, sbps_thrmcap,  &
     sbps_thrmcnd, sbps_lyrthkn

real(f4), dimension(:,:,:), intent(OUT) :: TItvar_sfc

!== Local declarations:

integer :: i, j, ip, ilev
real(f4) :: vol_thrmcap

!=== Executable statements:

do ip=1,size(sbps_lyrthkn, 4)
   do j=1,size(sbps_lyrthkn, 3)
      do i=1,size(sbps_lyrthkn, 2)

         ! Find current surface (layer index); sbps_lyrthkn > 0. if layer at or
         !  beneath current surface:
         ilev = size(sbps_lyrthkn, 1)  ! Try top SBPS layer first 
         do while (sbps_lyrthkn(ilev,i,j,ip) < tiny(vol_thrmcap))
            ilev = ilev-1                                       
         end do

         vol_thrmcap = sbps_density(ilev,i,j,ip)*sbps_thrmcap(ilev,i,j,ip)

         ! Calculate surface layer thermal inertia:
         TItvar_sfc(i,j,ip) = sqrt(vol_thrmcap*sbps_thrmcnd(ilev,i,j,ip))
      end do
   end do
end do

end subroutine comp_TItvar_sfc


!```````````````````````````````````````````````````````````````````
subroutine comp_CO2Tsat_sfc(Tsat, press)

!+ Computes sfc CO2 saturation temperatures

use univ_kind_defs_mod, only: f4
use PRAMS_SBPS_mod, only: Tsat_CO2ice

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: press

real(f4), dimension(:,:), intent(OUT) :: Tsat

!== Local declarations:

integer :: i, j
real(f4) :: sfc_press

!=== Executable statements:

do j=1,size(Tsat, 2)
   do i=1,size(Tsat, 1)
      sfc_press = (press(1,i,j)+press(2,i,j))*0.5_f4   ! Calc sfc pressure
      Tsat(i,j) = Tsat_CO2ice(sfc_press)
   end do
end do

end subroutine comp_CO2Tsat_sfc


!```````````````````````````````````````````````````````````````````
subroutine comp_CO2Tsat_lyr(Tsat, press)

!+ Computes 3D lyr CO2 saturation temperatures

use univ_kind_defs_mod, only: f4
use PRAMS_SBPS_mod, only: Tsat_CO2ice

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: press

real(f4), dimension(:,:,:), intent(OUT) :: Tsat

!== Local declarations:

integer :: i, j, k

!=== Executable statements:

do j=1,size(Tsat, 3)
   do i=1,size(Tsat, 2)
      do k=1,size(Tsat, 1)
         Tsat(k,i,j) = Tsat_CO2ice(press(k,i,j))
      end do
   end do
end do

end subroutine comp_CO2Tsat_lyr


!```````````````````````````````````````````````````````````````````
subroutine comp_modlyr_Z(output_Z, topo, rtg, zt, action, i_AZ)

!esta es la subrutina llamada por postp_varlib
!+ Computes model layer height info

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:), intent(IN) :: topo, rtg
real(f4), dimension(:), intent(IN) :: zt
character(len=*), intent(IN) :: action

real(f4), dimension(:,:,:), intent(OUT) :: output_Z

!== Local declarations:

integer :: i, j, k, kk

!=== Executable statements:

if (action == 'AGL') then    ! Output height above ground level

   do j=1,size(output_Z, 3)
      do i=1,size(output_Z, 2)
         do k=1,size(output_Z, 1)
            kk = i_AZ(k)
            output_Z(k,i,j) = zt(kk)*rtg(i,j)
         end do
      end do
   end do

else if (action == 'AREO') then    ! Output areopotential height
   
   do j=1,size(output_Z, 3)
      do i=1,size(output_Z, 2)
         do k=1,size(output_Z, 1)
            kk = i_AZ(k)
            output_Z(k,i,j) = zt(kk)*rtg(i,j)+topo(i,j)
         end do
      end do
   end do

end if

end subroutine comp_modlyr_Z

subroutine comp_zm(output_Z, topo, rtg, zm, action, i_AZ)

!esta es la subrutina llamada por postp_varlib
!+ Computes model layer height info

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:), intent(IN) :: topo, rtg
real(f4), dimension(:), intent(IN) :: zm
character(len=*), intent(IN) :: action

real(f4), dimension(:,:,:), intent(OUT) :: output_Z

!== Local declarations:

integer :: i, j, k, kk

!=== Executable statements:

if (action == 'AGL') then    ! Output height above ground level

   do j=1,size(output_Z, 3)
      do i=1,size(output_Z, 2)
         do k=1,size(output_Z, 1)
            kk = i_AZ(k)
            output_Z(k,i,j) = zm(kk)*rtg(i,j)
         end do
      end do
   end do

else if (action == 'AREO') then    ! Output areopotential height

   do j=1,size(output_Z, 3)
      do i=1,size(output_Z, 2)
         do k=1,size(output_Z, 1)
            kk = i_AZ(k)
            output_Z(k,i,j) = zm(kk)*rtg(i,j)+topo(i,j)
         end do
      end do
   end do

end if

end subroutine comp_zm


!```````````````````````````````````````````````````````````````````
subroutine comp_richnum(i_AZ, u_avg, v_avg, theta, nsfc_theta, topo,  &
     zt, ztop, richnum, richnum_tmp)

!+ Calculates bulk Richardson number 
!
!  Valid for bottom boundary of each layer {thus the value for (absolute) k=1
!   is undefined}

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4, ONE_HALF_f4
use univ_exception_mod, only: exception, FATAL
use univ_maxlen_string_mod, only: MAX_PROC_NAME_LEN
use postp_main_memory_mod, only: undef_value_f4
use body_specific_constants_mod, only: g_areoid

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: theta, u_avg, v_avg
real(f4), dimension(:,:,:), intent(IN) :: nsfc_theta
real(f4), dimension(:,:), intent(IN) :: topo
real(f4), dimension(:), intent(IN) :: zt
real(f4), intent(IN) :: ztop
real(f4), dimension(:,:,:), intent(INOUT) :: richnum_tmp

real(f4), dimension(:,:,:), intent(OUT) :: richnum

!== Local declarations:

character(len=MAX_PROC_NAME_LEN), parameter ::  &
     procedure_name = 'comp_richnum'

integer :: i, j, k, ip
real(f4) :: rtg, dz, dtheta, theta_avg, du, dv

!=== Executable statements:

ip = 1

do j=1,size(theta, 3)
   do i=1,size(theta, 2)

      rtg = ONE_f4-topo(i,j)/ztop

      k = 1; richnum_tmp(k,i,j) = undef_value_f4  ! Not enough info to calculate

      do k=1,size(richnum_tmp, 1)-1
         dz = (zt(k+1)-zt(k))*rtg    ! dz of flux layer
         dtheta = theta(k+1,i,j)-theta(k,i,j)
         theta_avg = ONE_HALF_f4*(theta(k+1,i,j)+theta(k,i,j))
         du = u_avg(k+1,i,j)-u_avg(k,i,j)
         dv = v_avg(k+1,i,j)-v_avg(k,i,j)
         richnum_tmp(k+1,i,j) = g_areoid*dz*dtheta/(theta_avg*(du*du+dv*dv))
      end do

      ! Return a (sub)set of values:
      do k=1,size(richnum, 1)
         richnum(k,i,j) = richnum_tmp(i_AZ(k),i,j)
      end do

   end do
end do

end subroutine comp_richnum


!```````````````````````````````````````````````````````````````````
subroutine comp_MOreduct_X(speed, u, v, theta, nsfc_theta, ustar, tstar,  &
     z0net, PI, topo, zlev2, ztop, zAGL_tgt, tempk_rdc, u_rdc, v_rdc, vtype)

!+ Extrapolates the k=2 wind components (u and v) [m/s] *OR* temperature
!   towards the ground, using MO similarity theory reduction
!
!  NOTE: *NOT* needed for (vtype = 'T'):  u, v, u_rdc, v_rdc, ustar
!        *NOT* needed for (vtype = 'U'):  tstar, PI, tempk_rdc, v, v_rdc
!        *NOT* needed for (vtype = 'V'):  tstar, PI, tempk_rdc, u, u_rdc

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4, ZERO_f4
use body_specific_constants_mod, only: g_areoid, vonk, cp

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: speed, u, v, theta, PI
real(f4), dimension(:,:,:), intent(IN) :: ustar, tstar, z0net, nsfc_theta
real(f4), dimension(:,:), intent(IN) :: topo
real(f4), intent(IN) :: zlev2, ztop, zAGL_tgt
character(len=1), intent(IN) :: vtype

real(f4), dimension(:,:), intent(OUT) :: tempk_rdc,u_rdc,v_rdc

!== Local declarations:

! Coefficients for PRAMS (`starsl`; Louis 1981) atm structure functions:
real(f4), parameter :: b = 5._f4, d = 5._f4, csm = 7.5_f4, csh = 5._f4 

integer :: i, j, ip
real(f4) :: rtg, zAGL_lv2, z0, windspd, richnum, c1, rtmp, rtmpv, rtmpu,  &
     ustar_x, ustar_y, rtmpt, PI_tgtlev, fh, fm, a2, cm, ch, cosine1, sine1
logical :: var_is_U, var_is_T, var_is_V

!=== Executable statements:
 
var_is_T = ( vtype == 'T' )
var_is_U = ( vtype == 'U' )
var_is_V = ( vtype == 'V' )

! Monin-Obukhov similarity theory-based reduction scheme from Arya (p166,
!  Eq 11.12):

do j=1,size(speed, 3)
   do i=1,size(speed, 2)

      rtg = ONE_f4-topo(i,j)/ztop
      zAGL_lv2 = zlev2*rtg    ! Height AGL of k=2 sigmaz level

      a2 = (vonk/log(zAGL_lv2/z0))**2
      
      if (var_is_T) then       !
         rtmpt = ZERO_f4           !
      else if (var_is_U) then   ! Initialize patch-weighted reduced vars
         rtmpu = ZERO_f4           !
      else if (var_is_V) then   !
         rtmpv = ZERO_f4           !   
      end if

!%      do ip=1,n5
      ip = 1
      
      z0 = z0net(i,j,ip)   ! Surface aerodynamic roughness length
 
      ! Drag coefficient for neutral stability (here, same for both h & m):
      a2 = (vonk/log(zAGL_lv2/z0))**2

      ! Adjust differently based on static stability:

      rtmp = theta(2,i,j)-nsfc_theta(i,j,ip)

      if (rtmp > ZERO_f4) then     ! Stable case

         windspd = max(speed(2,i,j), 0.1_f4)   ! [m/s]
         
         richnum = g_areoid*zAGL_lv2*rtmp/(0.5_f4*  &
              (theta(2,i,j)+nsfc_theta(i,j,ip))*  &
              windspd*windspd)    ! Richardson # (used to approximate z/L)

         fm = ONE_f4+(2._f4*b*richnum/sqrt(ONE_f4+d*richnum))
         fm = ONE_f4/fm    ! Momentum structure function

         fh = ONE_f4+(3._f4*b*richnum*sqrt(ONE_f4+d*richnum))
         fh = ONE_f4/fh    ! Heat (thermal) structure function

      else        ! Unstable case

         windspd = max(speed(2,i,j), ONE_f4)   ! [m/s]

         richnum = g_areoid*zAGL_lv2*rtmp/(0.5_f4*  &
              (theta(2,i,j)+nsfc_theta(i,j,ip))*  &
              windspd*windspd)    ! Richardson # (used to approximate z/L)

         c1 = b*a2*sqrt(zAGL_lv2/z0*abs(richnum))

         cm = csm*c1
         fm = ONE_f4-(2._f4*b*richnum/(ONE_f4+2._f4*cm))   ! Momentum structure function
         ch = csh*c1
         fh = ONE_f4-(3._f4*b*richnum/(ONE_f4+3._f4*ch))  ! Heat (thermal) structure function

      end if
         
      if (var_is_T) then

         c1 = (log(zAGL_tgt/z0)-fh)/vonk
         rtmp = nsfc_theta(i,j,ip)+tstar(i,j,ip)*c1
         !%       rtmpt = rtmpt+rtmp*patfrac(i,j,ip)
         rtmpt = rtmp

      else      ! Velocity of some type

         c1 = (log(zAGL_tgt/z0)-fm)/vonk

         if (var_is_U) then

            cosine1 = u(2,i,j)/windspd    ! Determine partitioning factor
                                          !  for wind component

            ustar_x = cosine1*ustar(i,j,ip)   ! Component of u* in x-dir

            rtmp = ustar_x*c1
            !%      rtmpu = rtmpu+rtmp*patfrac(i,j,ip)
            rtmpu = rtmp

         else if (var_is_V) then

            sine1 = v(2,i,j)/windspd         ! Determine partitioning factor
                                             !  for wind component

            ustar_y = sine1*ustar(i,j,ip)   ! Component of u* in y-dir

            rtmp = ustar_y*c1
            !%      rtmpv = rtmpv+rtmp*patfrac(i,j,ip)
            rtmpv = rtmp

         end if

      end if
      
!%      end do

      if (var_is_T) then

         c1 = (zAGL_lv2-zAGL_tgt)/zAGL_lv2
         PI_tgtlev = ((PI(1,i,j)+PI(2,i,j))*0.5_f4)*c1 + PI(2,i,j)*(ONE_f4-c1)

         tempk_rdc(i,j) = rtmpt*PI_tgtlev/cp   ! MO-reduced temperature
                                               !  (convert theta to temp)

      else if (var_is_U) then
         u_rdc(i,j) = rtmpu    ! MO-reduced wind component
      else if (var_is_V) then
         v_rdc(i,j) = rtmpv    ! MO-reduced wind component
      end if

   end do
end do

end subroutine comp_MOreduct_X


!```````````````````````````````````````````````````````````````````
subroutine comp_ustar_neutral(speed, z0net, topo, zlev2, ztop, ustar_neutral)

!+ Uses the k=2 wind speed [m/s] to estimate ustar if
!   the atmospheric stability were neutral, using M-O similarity theory

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use univ_maxlen_string_mod, only: MAX_PROC_NAME_LEN
use univ_exception_mod, only: exception, FATAL
use body_specific_constants_mod, only: vonk

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: speed
real(f4), dimension(:,:,:), intent(IN) :: z0net
real(f4), dimension(:,:), intent(IN) :: topo
real(f4), intent(IN) :: zlev2, ztop

real(f4), dimension(:,:,:), intent(OUT) :: ustar_neutral

!== Local declarations:

character(len=MAX_PROC_NAME_LEN), parameter ::  &
     procedure_name = 'comp_ustar_neutral'

integer :: i, j, ip
real(f4) :: rtg, zAGL_lv2, z0

!=== Executable statements:

! Monin-Obukhov similarity theory-based reduction scheme from Arya (p166,
!  Eq 11.12):

do j=1,size(ustar_neutral, 2)
   do i=1,size(ustar_neutral, 1)
      
      rtg = ONE_f4-topo(i,j)/ztop
      zAGL_lv2 = zlev2*rtg    ! Height AGL of k=2 sigmaz level

      do ip=1,size(ustar_neutral, 3)
         z0 = z0net(i,j,ip)   ! Surface aerodynamic roughness length
 
         ! Neutral stability ustar [m/s]:
         ustar_neutral(i,j,ip) = vonk*speed(2,i,j)/log(zAGL_lv2/z0)

      end do

   end do
end do

end subroutine comp_ustar_neutral


!```````````````````````````````````````````````````````````````````
subroutine comp_3Daer_optdepth(aer_numconc, rtgt, zt, zm, aer_OD, aer_type,  &
     aer_bin, wvlen_type, calc_mode, aer_OD_tmp, i_AZ)

!+ Calculate aerosol optical depth in 3D space (supports multiple radius bins)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use univ_convert_units_mod, only: cm3_to_m3_f4
use postp_main_memory_mod, only: dust_Qe_vis, dust_xsarea_fact, H2oi_Qe_vis,  &
     H2Oi_xsarea_fact, CO2i_Qe_vis, CO2i_xsarea_fact, dust_qe_IR, H2Oi_qe_IR,  &
     CO2i_qe_IR

implicit none

!=== Argument declarations:

integer, dimension(:), intent(IN) :: i_AZ
integer, intent(IN) :: aer_type, aer_bin, wvlen_type
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc  ! [#/cm^3]
real(f4), dimension(:,:), intent(IN) :: rtgt  ! [none]
real(f4), dimension(:), intent(IN) :: zt, zm  ! [m]
character(len=3), intent(IN) :: calc_mode   ! 'cum' or 'lyr'

real(f4), dimension(:), intent(INOUT) :: aer_OD_tmp
real(f4), dimension(:,:,:), intent(INOUT) :: aer_OD   ! [none]

!== Local declarations:

integer :: k, i, j, kk
real(f4) :: tau_bin, Qe_b, xsarea_b

integer, save :: ncall = 0

!=== Executable statements:

if (ncall == 0) then   ! Only do once
   call postp_init_optprop1    ! Obtain needed aerosol optical properties
   ncall = 1
end if

if (aer_type == 1) then    ! Dust
   if (wvlen_type == 1) then
      Qe_b = dust_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then
      Qe_b = dust_Qe_IR(aer_bin)
   end if
   xsarea_b = dust_xsarea_fact(aer_bin)
else if (aer_type == 2) then    ! H2O ice
   if (wvlen_type == 1) then   
      Qe_b = H2Oi_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then   
      Qe_b = H2Oi_Qe_IR(aer_bin)
   end if
   xsarea_b = H2Oi_xsarea_fact(aer_bin)
else if (aer_type == 3) then   ! CO2 ice
   if (wvlen_type == 1) then   
      Qe_b = CO2i_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then   
      Qe_b = CO2i_Qe_IR(aer_bin)
   end if
   xsarea_b = CO2i_xsarea_fact(aer_bin)
end if

! Add this bin's aerosol optical depth values (aer_OD MUST be initialized to
!  zero outside this routine):

if (calc_mode == 'cum') then

   do j=1,size(aer_OD, 3)
      do i=1,size(aer_OD, 2)
         
         tau_bin = ZERO_f4   ! Set optical depth at top of model to zero
         do k=size(aer_OD_tmp, 1)-1,2,-1
            tau_bin = tau_bin+Qe_b*xsarea_b*  &
                 aer_numconc(k,i,j)/cm3_to_m3_f4*rtgt(i,j)*(zt(k+1)-zt(k))
            tau_bin = max(1.e-20_f4, tau_bin)    ! Optical depth for this z,bin
            aer_OD_tmp(k) = aer_OD_tmp(k)+tau_bin   ! Total optical depth
         end do
         aer_OD_tmp(size(aer_OD_tmp, 1)) = ZERO_f4
         aer_OD_tmp(1) = aer_OD_tmp(2) ! Avoid vis/interp problems

         do k=1,size(aer_OD, 1)
            kk = i_AZ(k)
            aer_OD(k,i,j) = aer_OD(k,i,j)+aer_OD_tmp(kk)
         end do

      end do
   end do

else if (calc_mode == 'lyr') then

   do j=1,size(aer_OD, 3)
      do i=1,size(aer_OD, 2)
         
         do k=2,size(aer_OD_tmp, 1)-1
            tau_bin = Qe_b*xsarea_b*  &
                 aer_numconc(k,i,j)/cm3_to_m3_f4   ! Leave as "per unit volume"
            tau_bin = max(1.e-20_f4, tau_bin)    ! Optical depth for this z,bin
            aer_OD_tmp(k) = aer_OD_tmp(k)+tau_bin   ! Total optical depth
         end do
         aer_OD_tmp(size(aer_OD_tmp, 1)) = ZERO_f4
         aer_OD_tmp(1) = aer_OD_tmp(2) ! Avoid vis/interp problems      

         do k=1,size(aer_OD, 1)
            kk = i_AZ(k)
            aer_OD(k,i,j) = aer_OD(k,i,j)+aer_OD_tmp(kk)
         end do

      end do
   end do

end if

end subroutine comp_3Daer_optdepth


!```````````````````````````````````````````````````````````````````
subroutine comp_2Daer_optdepth(aer_numconc, rtgt, zt, zm, aer_OD, aer_type,  &
     aer_bin, wvlen_type)

!+ Calculate aerosol optical depth in 2D space (column total optical depth;
!   supports multiple radius bins)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use univ_convert_units_mod, only: cm3_to_m3_f4
use postp_main_memory_mod, only: dust_Qe_vis, dust_xsarea_fact, H2oi_Qe_vis,  &
     H2Oi_xsarea_fact, CO2i_Qe_vis, CO2i_xsarea_fact, dust_qe_IR, H2Oi_qe_IR,  &
     CO2i_qe_IR

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, aer_bin, wvlen_type
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc  ! [#/cm^3]
real(f4), dimension(:,:), intent(IN) :: rtgt  ! [none]
real(f4), dimension(:), intent(IN) :: zt, zm   ! [m]

real(f4), dimension(:,:), intent(INOUT) :: aer_OD   ! [none]

!== Local declarations:

integer :: k, i, j
real(f4) :: tau_bin, Qe_b, xsarea_b

integer, save :: ncall = 0

!=== Executable statements:

if (ncall == 0) then   ! Only do once
   call postp_init_optprop1    ! Obtain needed aerosol optical properties
   ncall = 1
end if

if (aer_type == 1) then    ! Dust
   if (wvlen_type == 1) then
      Qe_b = dust_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then
      Qe_b = dust_Qe_IR(aer_bin)
   end if
   xsarea_b = dust_xsarea_fact(aer_bin)
else if (aer_type == 2) then    ! H2O ice
   if (wvlen_type == 1) then   
      Qe_b = H2Oi_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then   
      Qe_b = H2Oi_Qe_IR(aer_bin)
   end if
   xsarea_b = H2Oi_xsarea_fact(aer_bin)
else if (aer_type == 3) then   ! CO2 ice
   if (wvlen_type == 1) then   
      Qe_b = CO2i_Qe_vis(aer_bin)
   else if (wvlen_type == 2) then   
      Qe_b = CO2i_Qe_IR(aer_bin)
   end if
   xsarea_b = CO2i_xsarea_fact(aer_bin)
end if

! Add this bin's aerosol optical depth values (aer_OD MUST be initialized to
!  zero outside this routine):

do j=1,size(aer_numconc, 3)
   do i=1,size(aer_numconc, 2)
         
      tau_bin = ZERO_f4   ! Set optical depth at top of model to zero
      do k=size(aer_numconc, 1)-1,2,-1
         tau_bin = tau_bin+Qe_b*xsarea_b*  &
              aer_numconc(k,i,j)/cm3_to_m3_f4*rtgt(i,j)*(zt(k+1)-zt(k))
         tau_bin = max(1.e-20_f4, tau_bin)      ! Optical depth for this z,bin
      end do
      aer_OD(i,j) = aer_OD(i,j)+tau_bin   ! Total optical depth
      
   end do
end do

end subroutine comp_2Daer_optdepth


!```````````````````````````````````````````````````````````````````
subroutine comp_supsat_H2Ov(tempK, rv, pvapX, supsatX)

!+ Computes supersaturation of water vapor (over ice)

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cp, molmass_H2O
use univ_phys_constants_mod, only: R_univ_gas_f4
use univ_convert_units_mod, only: g_cm3_to_kg_m3_f8

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: tempK,  &   ! [K]
                                          pvapX,  &     ! [Pa]
                                          rv          ! [g/cm^3]

real(f4), dimension(:,:,:), intent(OUT) :: supsatX

!== Local declarations:

real(f4) :: rvap

!=== Executable statements:

rvap = R_univ_gas_f4/molmass_H2O

supsatX = (rv*g_cm3_to_kg_m3_f8*rvap*tempK-pvapX)/pvapX

! Adjust value to be more easily understood
!  (i.e., supsat=0.01 --> 1%, supsat=-0.3 --> -30%
supsatX = supsatX*100._f4

end subroutine comp_supsat_H2Ov


!```````````````````````````````````````````````````````````````````
subroutine comp_supsat_CO2i(air_density, CO2_gas_frac, tempK, pvapi, supsati)

!+ Computes supersaturation of CO2 gas (over ice)

use univ_kind_defs_mod, only: f4
use body_specific_constants_mod, only: cp, molmass_CO2
use univ_phys_constants_mod, only: R_univ_gas_f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: tempK,  &   ! [K]
                                          pvapi,  &     ! [Pa]
                                          air_density,  &  ! [kg/m^3]
                                          CO2_gas_frac   ! [none]

real(f4), dimension(:,:,:), intent(OUT) :: supsati   ! [% supersaturation]

!== Local declarations:

real(f4) :: rvap

!=== Executable statements:

rvap = R_univ_gas_f4/molmass_CO2

supsati = (air_density*CO2_gas_frac*rvap*tempK-pvapi)/pvapi 

! Adjust value to be more easily understood
!  (i.e., supsat=0.01 --> 1%, supsat=-0.3 --> -30%
supsati = supsati*100._f4

end subroutine comp_supsat_CO2i


!```````````````````````````````````````````````````````````````````
subroutine comp_gas_vapor_press(tempK, pvapi, pvapl, gasname)

!+ Calculates the vapor pressure (over liquid,ice) for specified gas

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use univ_convert_units_mod, only: dyne_o_cm2_to_Pa_f4, mb_to_Pa_f4, atm_to_Pa_f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: tempK
character(len=1), intent(IN) :: gasname

real(f4), dimension(:,:,:), intent(OUT) :: pvapi, pvapl

!== Local declarations:

! Define coefficients in Buck's formulation for saturation vapor pressures:

real(f4), parameter :: BAI = 6.1115e3_f4, BBI = 23.036_f4,  &
                             BCI = 279.82_f4, BDI = 333.7_f4,  &
                             BAL = 6.1121e3_f4, BBL = 18.729_f4,  &
                             BCL = 257.87_f4, BDL = 227.3_f4

!=== Executable statements:

if (gasname == 'W') then     ! H2O vapor
   
   ! Saturation vapor pressure over liquid water and water ice
   !  (from Buck [J. Atmos. Sci., 20, 1527, 1981]):
   
   pvapl = tempK-273.16_f4        ! Temperature in deg_Celsius
   pvapl = BAL*exp((BBL-pvapl/BDL)*pvapl/(pvapl+BCL))
   pvapl = pvapl*dyne_o_cm2_to_Pa_f4   ! [dynes/cm^2 -> Pa] conversion
            
   !orig%         pvapi(ik,igas) = BAI*exp((BBI-tt/BDI)*tt/(tt+BCI))
   pvapi = 6.11_f4*exp(22.5_f4*(ONE_f4-(273.16_f4/tempK))) !Tony alt??
   pvapi = pvapi*mb_to_Pa_f4         ! " [mb -> Pa] conversion

   ! Check to see whether temperature is ouside range of validity
   !  for parameterizations (???); threshold is 1.e-10 in cgs units:
   where (pvapl <= 1.e-11_f4) pvapl = 1.e-11_f4

else if (gasname == 'C') then     ! carbon dioxide gas

   ! Vapor pressures for CO2 taken from [Kasting, Icarus 94, 1991 (Atm)]:

   where (tempK >= 216.56_f4)
      pvapi = 10._f4**(3.128082_f4-(867.2124_f4/tempK)+  &
           18.65612e-3_f4*tempK-72.48820e-6_f4*tempK**2+  &
           93.e-9_f4*tempK**3)
      pvapi = pvapi*atm_to_Pa_f4
      pvapl = pvapi
   elsewhere
      pvapi = 10._f4**(6.760956_f4-(1284.07_f4/(tempK-4.718_f4))+  &
           1.256e-4_f4*(tempK-143.15_f4))
      pvapi = pvapi*atm_to_Pa_f4  !??? include conversion???correct???
      pvapl = pvapi
   end where

end if

end subroutine comp_gas_vapor_press


!```````````````````````````````````````````````````````````````````
subroutine comp_sfc_stress(ustar, sfc_air_density, sfc_stress)

!+ Computes aerodynamic surface shear stress [Pa], given ustar and surface air
!   density

use univ_kind_defs_mod, only: f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: ustar, sfc_air_density

real(f4), dimension(:,:,:), intent(OUT) :: sfc_stress

!== Local declarations:

integer :: ip

!=== Executable statements:

do ip=1,size(sfc_stress, 3)
   sfc_stress(:,:,ip) = ustar(:,:,ip)*ustar(:,:,ip)*sfc_air_density(:,:,ip)
end do

end subroutine comp_sfc_stress


!```````````````````````````````````````````````````````````````````
subroutine comp_ustar_t(tempk_air_sfc, massdensity_air_sfc,  &
     given_ustar_t_table, ustar_t)

!+ Computes minimum (saltation) threshold friction velocity [m/s] (from a
!   pre-filled table).

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4
use postp_main_memory_mod, only: ustar_t_table

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: tempk_air_sfc,  &  ! [K]
                                          massdensity_air_sfc  ! [kg/m^3]
type(ustar_t_table), intent(IN) :: given_ustar_t_table

real(f4), dimension(:,:,:), intent(OUT) :: ustar_t

!== Local declarations:

integer :: i, j, i_itp1, i_itp2, ip
real(f4) :: ftmp, ftmp1, f_itp1, f_itp2, ftmp2

!=== Executable statements:

! Determine current ustar_t:

!%do ip=1,n5
ip = 1
do j=1,size(tempk_air_sfc, 2)
   do i=1,size(tempk_air_sfc, 1)

      ftmp = min(max(massdensity_air_sfc(i,j,ip),  &
           given_ustar_t_table%ust_tab_rhof_min),  &
           given_ustar_t_table%ust_tab_rhof_max-  &
           given_ustar_t_table%ust_tab_rhof_inc*1.001_f4)
      ftmp1 = (ftmp-given_ustar_t_table%ust_tab_rhof_min)/  &
           given_ustar_t_table%ust_tab_rhof_inc
      i_itp1 = int(ftmp1)+1
      f_itp1 = ftmp1-real(i_itp1-1, kind=f4)

      ftmp = min(max(tempk_air_sfc(i,j,ip),  &
           given_ustar_t_table%ust_tab_T_min),  &
           given_ustar_t_table%ust_tab_T_max-  &
           given_ustar_t_table%ust_tab_T_inc*1.001_f4)
      ftmp2 = (ftmp-given_ustar_t_table%ust_tab_T_min)/  &
           given_ustar_t_table%ust_tab_T_inc
      i_itp2 = int(ftmp2)+1
      f_itp2 = ftmp2-real(i_itp2-1, kind=f4)
      
      ftmp1 = given_ustar_t_table%ustar_t_min(i_itp1,i_itp2)*(ONE_f4-f_itp1)+  &
           given_ustar_t_table%ustar_t_min(i_itp1+1,i_itp2)*f_itp1
      ftmp2 =  &
           given_ustar_t_table%ustar_t_min(i_itp1,i_itp2+1)*(ONE_f4-f_itp1)+  &
           given_ustar_t_table%ustar_t_min(i_itp1+1,i_itp2+1)*f_itp1

      ustar_t(i,j,ip) = ftmp1*(ONE_f4-f_itp2)+ftmp2*f_itp2
   end do
end do
!%end do

end subroutine comp_ustar_t


!```````````````````````````````````````````````````````````````````
subroutine postp_init_optprop1

!+

use univ_kind_defs_mod, only: f4
use postp_main_memory_mod, only: r, nbins, dust_Qe_vis, dust_xsarea_fact,  &
     H2Oi_Qe_vis, H2Oi_xsarea_fact, CO2i_Qe_vis, CO2i_xsarea_fact,  &
     dust_Qe_IR, H2Oi_Qe_IR, CO2i_Qe_IR, dust_op_case
use aertxrad_mem_mod, only: naerspc_rad, spc_xsarea_fact, sw_abp_optprop,  &
     lw_abp_optprop, spc_ig, spc_ctype, aeradID_dust, aeradID_H2Oi, aeradID_CO2i
use aerad_init_mod, only: init_optprop

implicit none

!=== Argument declarations:

! NONE

!== Local declarations:

integer, parameter :: vis_wvband = 6, IR_wvband = 4

integer :: ip, ib, ig, ic_fvws_r, ic_fvws_b, iw, ia

!=== Executable statements:

! Ingest optical property data:
call init_optprop(r, dust_op_case)

! Copy particle cross-sectional area values:
do ip=1,naerspc_rad ! Loop over all possible RT-active airborne particle species
   ig = spc_ig(ip)  ! Relative group ID#
   ia = spc_ctype(ip)   ! airborne particle type
   
   do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
      select case (ia)
      case (aeradID_dust)
         dust_xsarea_fact(ib) = spc_xsarea_fact(ib,ip)  ! [m^2]
      case (aeradID_H2Oi)
         H2Oi_xsarea_fact(ib) = spc_xsarea_fact(ib,ip)  ! [m^2]
      case (aeradID_CO2i)
         CO2i_xsarea_fact(ib) = spc_xsarea_fact(ib,ip)  ! [m^2]
      end select
   end do
end do

!********* "Shortwave" spectral intervals **************************************

! Set wavelength-, radius-, and specie-dependent airborne/sfc particle optical
!  properties:

iw = vis_wvband
do ip=1,naerspc_rad ! Loop over all possible RT-active airborne particle species

   ic_fvws_r = sw_abp_optprop(ip)%search_beg_idx_fvws_r
   ic_fvws_b = sw_abp_optprop(ip)%search_beg_idx_fvws_b
   ig = spc_ig(ip)  ! Relative group ID#
   ia = spc_ctype(ip)   ! airborne particle type

   select case (ia)
   case (aeradID_dust)   
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         dust_Qe_vis(ib) = sw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   case (aeradID_H2Oi)
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         H2Oi_Qe_vis(ib) = sw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   case (aeradID_CO2i)
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         CO2i_Qe_vis(ib) = sw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   end select

end do

!********* "Longwave" spectral intervals **************************************

! Set wavelength-, radius-, and specie-dependent airborne/sfc particle optical
!  properties:

iw = IR_wvband
do ip=1,naerspc_rad ! Loop over all possible RT-active airborne particle species

   ic_fvws_r = lw_abp_optprop(ip)%search_beg_idx_fvws_r
   ic_fvws_b = lw_abp_optprop(ip)%search_beg_idx_fvws_b
   ig = spc_ig(ip)  ! Relative group ID#
   ia = spc_ctype(ip)   ! airborne particle type

   select case (ia)
   case (aeradID_dust)   
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         dust_Qe_IR(ib) = lw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   case (aeradID_H2Oi)
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         H2Oi_Qe_IR(ib) = lw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   case (aeradID_CO2i)
      do ib=1,nbins(ig)  ! Loop over # of airborne particle radius bins
         CO2i_Qe_IR(ib) = lw_abp_optprop(ip)%Q_ext(ic_fvws_r,ic_fvws_b,ib,iw)
      end do
   end select

end do

end subroutine postp_init_optprop1


!```````````````````````````````````````````````````````````````````
subroutine comp_3Daer_mnrad(aer_numconc, aer_mnrad, aer_totnum, aer_type,  &
     aer_bin, phase, i_AZ, n_AZ_abs)

!+ Calculate mean particle radius (number-wise) (supports multiple radius bins)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4, ZERO_f4
use postp_main_memory_mod, only: r, dust_ig, H2Oi_ig, CO2i_ig

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, aer_bin, phase, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc   ! [#/cm^3]

real(f4), dimension(:,:,:), intent(INOUT) :: aer_mnrad,  &  ! [cm]
                                             aer_totnum   ! [#/cm^3]

!== Local declarations:

integer :: k, i, j, kk, n1
real(f4) :: r_bin

!=== Executable statements:

n1 = size(aer_mnrad, 1)

if (phase == 0) then          ! Initialize summing arrays

   aer_mnrad = ZERO_f4    ! Fill entire arrays at once
   aer_totnum = ZERO_f4   !

else if (phase == 1) then    ! Sum over aerosol bins

   select case (aer_type)
   case(1)    ! Dust
      r_bin = r(aer_bin,dust_ig)   ! Avg radius in this bin [cm]
   case(2)    ! H2O ice
      r_bin = r(aer_bin,H2Oi_ig)   ! Avg radius in this bin [cm]
   case(3)    ! CO2 ice
      r_bin = r(aer_bin,CO2i_ig)   ! Avg radius in this bin [cm]
   end select

   ! Add this bin's aerosol contribution ('aer_totnum' and 'aer_mnrad' 
   !  MUST be initialized to zero with 'phase 0'):

   do j=1,size(aer_mnrad, 3)
      do i=1,size(aer_mnrad, 2)
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               aer_mnrad(k,i,j) = aer_mnrad(k,i,j)+aer_numconc(k,i,j)*r_bin
               aer_totnum(k,i,j) = aer_totnum(k,i,j)+aer_numconc(k,i,j)
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) then
            aer_mnrad(n1,i,j) = ZERO_f4
            aer_totnum(n1,i,j) = ONE_f4
         end if
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_mnrad(1,i,j) = aer_mnrad(2,i,j)
               aer_totnum(1,i,j) = aer_totnum(2,i,j)
            else
               aer_mnrad(1,i,j) = ZERO_f4
               aer_totnum(1,i,j) = ZERO_f4
            end if
         end if
      end do
   end do

else if (phase == 2) then    ! Finish computation

   do j=1,size(aer_mnrad, 3)
      do i=1,size(aer_mnrad, 2)
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               aer_mnrad(k,i,j) = aer_mnrad(k,i,j)/  &
                    max(1.e-12_f4, aer_totnum(k,i,j))
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) then
            if (n1 >= 2) then
               aer_mnrad(n1,i,j) = aer_mnrad(n1-1,i,j)
            else
               aer_mnrad(n1,i,j) = ZERO_f4
            end if
         end if
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_mnrad(1,i,j) = aer_mnrad(2,i,j)
            else
               aer_mnrad(1,i,j) = ZERO_f4
            end if
         end if
      end do
   end do

end if

end subroutine comp_3Daer_mnrad


!```````````````````````````````````````````````````````````````````
subroutine comp_3Daer_variance(aer_numconc, aer_mnrad, aer_totnum,  &
     aer_variance, aer_type, aer_bin, phase, i_AZ, n_AZ_abs)

!+ Calculate variance of particle radii (number-wise) (supports multiple radius
!   bins)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use postp_main_memory_mod, only: r, dust_ig, H2Oi_ig, CO2i_ig

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, aer_bin, phase, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc,  &   ! [#/cm^3]
                                          aer_mnrad,  &  ! [cm]
                                          aer_totnum   ! [#/cm^3]

real(f4), dimension(:,:,:), intent(INOUT) :: aer_variance   ! [cm^2]

!== Local declarations:

integer :: k, i, j, kk, n1
real(f4) :: r_bin

!=== Executable statements:

n1 = size(aer_variance, 1)

if (phase == 0) then          ! Initialize summing array

   aer_variance = ZERO_f4    ! Fill entire array at once

else if (phase == 1) then    ! Sum over aerosol bins

   select case (aer_type)
   case(1)    ! Dust
      r_bin = r(aer_bin,dust_ig)   ! Avg radius in this bin [cm]
   case(2)    ! H2O ice
      r_bin = r(aer_bin,H2Oi_ig)   ! Avg radius in this bin [cm]
   case(3)    ! CO2 ice
      r_bin = r(aer_bin,CO2i_ig)   ! Avg radius in this bin [cm]
   end select

   ! Add this bin's aerosol contribution ('aer_variance' 
   !  MUST be initialized to zero with 'phase 0'):

   do j=1,size(aer_variance, 3)
      do i=1,size(aer_variance, 2)
         
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               aer_variance(k,i,j) = aer_variance(k,i,j)+  &
                    aer_numconc(k,i,j)*(r_bin-aer_mnrad(k,i,j))**2
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) aer_variance(n1,i,j) = ZERO_f4
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_variance(1,i,j) = aer_variance(2,i,j)
            else
               aer_variance(1,i,j) = ZERO_f4
            end if
         end if

      end do
   end do

else if (phase == 2) then    ! Finish computation

   do j=1,size(aer_variance, 3)
      do i=1,size(aer_variance, 2)
         
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               aer_variance(k,i,j) = aer_variance(k,i,j)/  &
                    max(1.e-12_f4, aer_totnum(k,i,j))
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) then
            if (n1 >= 2) then
               aer_variance(n1,i,j) = aer_variance(n1-1,i,j)
            else
               aer_variance(n1,i,j) = ZERO_f4
            end if
         end if
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_variance(1,i,j) = aer_variance(2,i,j)
            else
               aer_variance(1,i,j) = ZERO_f4
            end if
         end if

      end do
   end do

end if

end subroutine comp_3Daer_variance


!```````````````````````````````````````````````````````````````````
subroutine comp_H2Ov_cond_coldepth(H2Ov_mc, H2Ov_pd, rtgt, zm)

!+ Compute depth [cm] of precipitable water in vapor column  

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use univ_convert_units_mod, only: m_to_cm_f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: H2Ov_mc
real(f4), dimension(:), intent(IN) :: zm
real(f4),dimension(:,:), intent(IN) :: rtgt

real(f4), dimension(:,:), intent(OUT) :: H2Ov_pd

!== Local declarations:

real(f4), parameter :: H2O_rho = 1._f4   ! Density of water [g/cm^3]

integer :: k,i,j
real(f4) :: H2Ov_lyrmass, tmpsum

!=== Executable statements:

do j=1,size(H2Ov_mc, 3)
   do i=1,size(H2Ov_mc, 2)
      tmpsum = ZERO_f4
      do k=2,size(H2Ov_mc, 1)
         H2Ov_lyrmass = H2Ov_mc(k,i,j)*rtgt(i,j)*(zm(k)-zm(k-1))*m_to_cm_f4
         tmpsum = tmpsum+H2Ov_lyrmass/H2O_rho
      end do
      H2Ov_pd(i,j) = tmpsum    ! [cm]
   end do
end do

end subroutine comp_H2Ov_cond_coldepth


!```````````````````````````````````````````````````````````````````
subroutine comp_mass_mr(mixing_ratio, X_density, air_density, X_is_in_cgs)

!+ Computes mass mixing ratio (g_of_X / g_of_air) in a unit volume

use univ_kind_defs_mod, only: f4
use univ_convert_units_mod, only: kg_to_g_f4, m3_to_cm3_f4

implicit none

!=== Argument declarations:

real(f4), dimension(:,:,:), intent(IN) :: air_density,  &   ! [kg/m^3]
     X_density  ! In units of either [g/cm^3] or [kg/m^3], but set
                !  'X_is_in_cgs' appropriately
logical, intent(IN) :: X_is_in_cgs

real(f4), dimension(:,:,:), intent(OUT) ::  &
     mixing_ratio ! [mass_of_X / mass_of_air]

!== Local declarations:

! NONE

!=== Executable statements:

if (X_is_in_cgs) then    ! X is in cgs units; convert air density to cgs

   mixing_ratio = X_density/(air_density*kg_to_g_f4/m3_to_cm3_f4)

else         ! X is in MKS units; no need to convert air density

   mixing_ratio = X_density/air_density

end if

end subroutine comp_mass_mr


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_numc2massc(aer_massconc, aer_numconc, aer_type, aer_bin,  &
     phase, i_AZ, n_AZ_abs)

!+ Computes aerosol mass concentration [g/cm^3] from [#/cm^3]; can sum bins

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use postp_main_memory_mod, only: dust_ig, H2Oi_ig, CO2i_ig, rmass

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, phase, aer_bin, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc

real(f4), dimension(:,:,:), intent(INOUT) :: aer_massconc

!== Local declarations:

integer :: i, j, k, n1, kk
real(f4) :: rmass_bin

!=== Executable statements:

if (phase == 0) then          ! Initialize summing array

   aer_massconc = ZERO_f4    ! Fill entire array at once

else if (phase == 1) then

   select case (aer_type)
   case(1)    ! Dust
      rmass_bin = rmass(aer_bin,dust_ig)   ! Avg particle mass in this bin [g]
   case(2)    ! H2O ice
      rmass_bin = rmass(aer_bin,H2Oi_ig)   ! Avg particle mass in this bin [g]
   case(3)    ! CO2 ice
      rmass_bin = rmass(aer_bin,CO2i_ig)   ! Avg particle mass in this bin [g]
   end select

   ! Add this bin's aerosol contribution ('aer_massconc' MUST be initialized
   !  to zero with 'phase 0' *before* starting this phase):

   n1 = size(aer_massconc, 1)
   do j=1,size(aer_massconc, 3)
      do i=1,size(aer_massconc, 2)
         
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               aer_massconc(k,i,j) = aer_massconc(k,i,j)+  &
                    aer_numconc(k,i,j)*rmass_bin
            end if
         end do

         if (i_AZ(n1) == n_AZ_abs) aer_massconc(n1,i,j) = ZERO_f4
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_massconc(1,i,j) = aer_massconc(2,i,j)
            else
               aer_massconc(1,i,j) = ZERO_f4
            end if
         end if

      end do
   end do

end if

end subroutine comp_aer_numc2massc


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_numc2massc_2D(aer_massconc, aer_numconc, aer_type,  &
     aer_bin, phase)

!+ Computes aerosol mass concentration [g/cm^3] from [#/cm^3]; can sum bins

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use postp_main_memory_mod, only: dust_ig, H2Oi_ig, CO2i_ig, rmass

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, phase, aer_bin
real(f4), dimension(:,:), intent(IN) :: aer_numconc

real(f4), dimension(:,:), intent(INOUT) :: aer_massconc

!== Local declarations:

real(f4) :: rmass_bin

!=== Executable statements:

if (phase == 0) then          ! Initialize summing array

   aer_massconc = ZERO_f4    ! Fill entire array at once

else if (phase == 1) then

   select case (aer_type)
   case(1)    ! Dust
      rmass_bin = rmass(aer_bin,dust_ig)   ! Avg particle mass in this bin [g]
   case(2)    ! H2O ice
      rmass_bin = rmass(aer_bin,H2Oi_ig)   ! Avg particle mass in this bin [g]
   case(3)    ! CO2 ice
      rmass_bin = rmass(aer_bin,CO2i_ig)   ! Avg particle mass in this bin [g]
   end select

   ! Add this bin's aerosol contribution ('aer_massconc' MUST be initialized
   !  to zero with 'phase 0' *before* starting this phase):

   aer_massconc = aer_massconc+aer_numconc*rmass_bin

end if

end subroutine comp_aer_numc2massc_2D


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_pshell_massc(ncore_shells, aer_massconc_csh1,  &
     aer_massconc_csh2, aer_numconc, aer_massconc_psh, aer_type, aer_bin,  &
     phase, i_AZ, n_AZ_abs)

!+ Computes aerosol mass concentration of the "primary" shell (e.g., of water
!   ice for a "water ice particle") [g/cm^3]; can sum bins

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4
use postp_main_memory_mod, only: dust_ig, H2Oi_ig, CO2i_ig, rmass

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, phase, aer_bin, ncore_shells, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: aer_numconc
real(f4), dimension(:,:,:), intent(IN) :: aer_massconc_csh1, aer_massconc_csh2

real(f4), dimension(:,:,:), intent(INOUT) :: aer_massconc_psh

!== Local declarations:

integer :: i, j, k, ic, kk, n1
real(f4) :: rmass_bin, core_mc, particle_mc, pshell_mc

!=== Executable statements:

if (phase == 0) then          ! Initialize summing array

   aer_massconc_psh = ZERO_f4    ! Fill entire array at once

else if (phase == 1) then

   select case (aer_type)
   case(1)    ! Dust
      rmass_bin = rmass(aer_bin,dust_ig)   ! Avg particle mass in this bin [g]
   case(2)    ! H2O ice
      rmass_bin = rmass(aer_bin,H2Oi_ig)   ! Avg particle mass in this bin [g]
   case(3)    ! CO2 ice
      rmass_bin = rmass(aer_bin,CO2i_ig)   ! Avg particle mass in this bin [g]
   end select

   ! Add this bin's aerosol contribution ('aer_massconc' MUST be initialized
   !  to zero with 'phase 0' *before* starting this phase):

   n1 = size(aer_massconc_psh, 1)
   do j=1,size(aer_massconc_psh, 3)
      do i=1,size(aer_massconc_psh, 2)
         
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               particle_mc = aer_numconc(k,i,j)*rmass_bin
               pshell_mc = particle_mc
               do ic=1,ncore_shells
                  if (ic == 1) then
                     core_mc = aer_massconc_csh1(k,i,j)
                  else if (ic == 2) then
                     core_mc = aer_massconc_csh2(k,i,j)
                  end if
                  pshell_mc = pshell_mc-core_mc
               end do
               pshell_mc = max(ZERO_f4, pshell_mc)   ! Clip if value is negative
               aer_massconc_psh(k,i,j) = aer_massconc_psh(k,i,j)+pshell_mc
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) aer_massconc_psh(n1,i,j) = ZERO_f4
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               aer_massconc_psh(1,i,j) = aer_massconc_psh(2,i,j)
            else
               aer_massconc_psh(1,i,j) = ZERO_f4
            end if
         end if

      end do
   end do

end if

end subroutine comp_aer_pshell_massc


!```````````````````````````````````````````````````````````````````
subroutine comp_3D_esum(sumA, X, phase, i_AZ, n_AZ_abs)

!+ Tool for summing successive 3D array element values

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ZERO_f4

implicit none

!=== Argument declarations:

integer, intent(IN) :: phase, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), dimension(:,:,:), intent(IN) :: X

real(f4), dimension(:,:,:), intent(INOUT) :: sumA

!== Local declarations:

integer :: i, j, k, n1, kk

!=== Executable statements:

if (phase == 0) then          ! Initialize summing array

   sumA = ZERO_f4    ! Fill entire array at once

else if (phase == 1) then

   ! Add this array's contribution ('sumA' MUST be initialized
   !  to zero with 'phase 0' *before* starting this phase):

   n1 = size(sumA, 1)
   do j=1,size(sumA, 3)
      do i=1,size(sumA, 2)
         
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               sumA(k,i,j) = sumA(k,i,j)+X(k,i,j)
            end if
         end do
         if (i_AZ(n1) == n_AZ_abs) sumA(n1,i,j) = ZERO_f4
         if (i_AZ(1) == 1) then
            ! k=1 val (avoid display intrp problems):
            if (n1 >= 2) then
               sumA(1,i,j) = sumA(2,i,j)
            else
               sumA(1,i,j) = ZERO_f4
            end if
         end if

      end do
   end do

end if

end subroutine comp_3D_esum


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_vfall(rhoa, t, ig, ib, vf)

!+ Evaluates particle fall velocities

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: pi_f4, ONE_f4, ZERO_f4, ONE_NINTH_f4,  &
     ONE_SIXTH_f4, ONE_THIRD_f4, TWO_THIRDS_f4
use univ_convert_units_mod, only: J_to_erg_f4, kg_to_g_f4, m_to_cm_f4,  &
     kg_m3_to_g_cm3_f4
use postp_main_memory_mod, only: ngroups, nbins, r, max_nbins, maxnzp, rhop, vol
use body_specific_constants_mod, only: R_air, g_areoid

implicit none

!=== Argument declarations:

integer, intent(IN) :: ib, ig
real(f4), dimension(:,:,:), intent(IN) :: rhoa,  &  ! Air density [kg/m^3]
                                          t    ! Air temperature [K]

real(f4), dimension(:,:,:), intent(OUT) :: vf    ! Fall speed [cm/s]

!== Local declarations:

! This routine evaluates particle fall velocities, vf(k) [cm s^-1]
!  and Reynolds' numbers based on fall velocities, re(j,i,k) [dimensionless].
!  indices correspond to vertical level <k>, bin index <i>, and aerosol
!  group <j>.
!
!  Non-spherical particles are treated through shape factors <ishape>
!  and <eshape>.
!  
!  General method is to first use Stokes' flow to estimate fall
!  velocity, then calculate reynolds' number, then use "y function" 
!  (defined in Pruppacher and Klett) to reevaluate reynolds' number,
!  from which the fall velocity is finally obtained.
!
!  This routine requires that vertical profiles of temperature <t>,
!  air density <rhoa>, and viscosity <rmu> are defined (i.e., initatm.f
!  must be called before this).  The vertical profile with ix = iy = 1
!  is used.
!
!  ******* MODIFIED FOR POSTPROCESSOR *******

real(f4) :: f1, f2, f3, exx, exy, xcc, xa, vg, rhoa_cgs, rmfp, r_shape, rkn,  &
     expon, rfix, x, y, ex, b0, bb1, bb2, bb3, bxx, z
integer :: k, i, j
real(f4), parameter :: POWMAX = 85._f4,  &
     R_air_cgs = R_air*J_to_erg_f4/kg_to_g_f4, g_cgs = g_areoid*m_to_cm_f4
real(f4) :: bpm, rmu, re, eshape
integer :: ishape

!=== Executable statements:

!###
ishape = 1 
eshape = ONE_f4
!###

! First evaluate factors that depend upon particle shape (used in correction
!  factor <bpm> below):  

select case (ishape)
case(1)     ! Spherical
   f1 = ONE_f4
   f2 = ONE_f4
case(2)     ! Hexagonal
   ! Taken from Turco et al. (Planet. Space Sci. Rev. 30, 1147-1181, 1982)
   !  with diffuse reflection of air molecules assumed:
   f2 = (pi_f4*ONE_NINTH_f4/tan(pi_f4*ONE_SIXTH_f4))**ONE_THIRD_f4*  &
        eshape**ONE_SIXTH_f4
case(3)     ! Spheroidal
   ! Taken FROM WHERE?:
   f2 = TWO_THIRDS_f4**ONE_THIRD_f4*eshape**ONE_SIXTH_f4
end select

! Yields <f3> = 1. for <eshape> = 1 :
f3 = 1.39_f4/sqrt((1.14_f4+0.25_f4/eshape)*(0.89_f4+eshape*0.5_f4))
f2 = f2*f3

if (eshape > ONE_f4) then

   ! For the Stokes regime there is no separate data for hexagonal plates or
   ! columns, so we use prolate spheroids.  This is from Fuchs' book:
   exx = eshape**2-ONE_f4
   exy = sqrt(exx)
   xcc = 1.333_f4*exx/((2._f4*eshape**2-ONE_f4)*log(eshape+exy)/exy  &
        -eshape)
   xa = 2.666_f4*exx/((2._f4*eshape**2-3._f4)*log(eshape+exy)/exy  &
        +eshape)
   f1 = eshape**(-ONE_THIRD_f4)*(xcc+2._f4*xa)*ONE_THIRD_f4
      
else if (eshape < ONE_f4) then

   ! Use oblate spheroids for disks (eshape < 1.); Also from Fuchs' book:
   bxx = ONE_f4/eshape
   exx = bxx**2-ONE_f4
   exy = sqrt(exx)
   xcc = 1.333_f4*exx/(bxx*(bxx**2-2._f4)*atan(exy)/exy+bxx)
   xa = 2.666_f4*exx/(bxx*(3._f4*bxx**2-2._f4)*atan(exy)/exy-bxx)
   f1 = bxx**ONE_THIRD_f4*(xcc+2._f4*xa)*ONE_THIRD_f4

end if

do j=1,size(t, 3)
   do i=1,size(t, 2)
      do k=2,size(t, 1)      !  Loop over column

         ! This is just <rhoa> in Cartesian coordinates:       
         rhoa_cgs = rhoa(k,i,j)*kg_m3_to_g_cm3_f4

         ! <vg> is the mean thermal velocity of air molecules [cm/s]:
         vg = sqrt(8._f4/pi_f4*R_air_cgs*t(k,i,j))

         ! <rmfp> is the mean free path of air molecules [cm]:
         ! Viscosity (CO2 atm) [g/(cm.s)]:
         !@alt_formulation   rmu(k) = 5.65e-8*((t(k)**2.5)/t(k)+269.)
         rmu = 1.59e-6_f4*t(k,i,j)**1.5_f4/(t(k,i,j)+244.4_f4)*10._f4
         rmfp = 2._f4*rmu/(rhoa_cgs*vg)

         ! <r_shape> is the radius of particle used to calculate <re>:
         select case (ishape)
         case(1)    ! Spherical
            r_shape = r(ib,ig)
         case(2)    ! Hexagonal
            r_shape = r(ib,ig)*0.85_f4*eshape**(-ONE_THIRD_f4)
         case(3)    ! Spheroidal
            r_shape = r(ib,ig)*eshape**(-ONE_THIRD_f4)
         end select

         ! <rkn> is the Knudsen number:
         rkn = rmfp/r(ib,ig)

         ! <bpm> is correction term for particle shape and non-continuum
         !  effects. It is also used to calculate coagulation kernels and
         !  diffusion coefficients:        
         expon = -0.87_f4/rkn
         expon = max(-POWMAX, expon)
         bpm = ONE_f4+f1*f2*(1.246_f4*rkn+0.42_f4*rkn*exp(expon))
         
         ! These are first guesses for fall velocity and Reynolds' #; 
         !  valid for Reynolds' number < 0.01
         
         vf(k,i,j) = 2._f4*ONE_NINTH_f4*rhop(k,ib,ig)*(r(ib,ig)**2)  &
              *g_cgs*bpm/(f1*rmu)
         
         re = 2._f4*rhoa_cgs*r_shape*vf(k,i,j)/rmu
         
         ! <rfix> is used in drag coefficient:
         rfix = vol(ib,ig)*rhop(k,ib,ig)*g_cgs*rhoa_cgs/rmu**2
         
         if ((re >= 0.01_f4) .and. (re <= 300._f4)) then
            
            ! This is "regime 2" in Pruppacher and Klett (chap. 10):

            if (ishape == 1) then     ! Spherical

               x = log(24._f4*re/bpm)
               y = -3.18657_f4+x*(0.992696_f4-x*(0.00153193_f4  &
                    + x*(0.000987059_f4+x*(0.000578878_f4  &
                    - x*(8.55176e-05_f4-x*3.27815e-06_f4)))))
               if (y < -675._f4) then
                  y = -675._f4
               else if (y >= 741._f4) then
                  y = 741._f4
               end if

               re = exp(y)*bpm

            else if (eshape <= ONE_f4) then
               
               if (ishape == 2) then
                  x = log10(16._f4*rfix/(3._f4*sqrt(3._f4*ONE_f4)))
               else if (ishape == 3) then
                  x = log10(8._f4*rfix/pi_f4)
               end if

               if (eshape <= 0.2_f4) then
                  b0 = -1.33_f4
                  bb1 = 1.0217_f4
                  bb2 = -0.049018_f4
                  bb3 = ZERO_f4
               else if (eshape <= 0.5_f4) then
                  ex = (eshape-0.2_f4)/0.3_f4
                  b0 = -1.33_f4+ex*(-1.3247_f4+1.33_f4*ONE_f4)
                  bb1 = 1.0217_f4+ex*(1.0396_f4-1.0217_f4*ONE_f4)
                  bb2 = -0.049018_f4+ex*(-0.047556_f4+0.049018_f4*ONE_f4)
                  bb3 = ZERO_f4+ex*(-0.002327_f4*ONE_f4)
               else
                  ex = (eshape-0.5_f4)/0.5_f4
                  b0 = -1.3247_f4+ex*(-1.310_f4+1.3247_f4*ONE_f4)
                  bb1 = 1.0396_f4+ex*(0.98968_f4-1.0396_f4*ONE_f4)
                  bb2 = -0.047556_f4+ex*(-0.042379_f4+0.047556_f4*ONE_f4)
                  bb3 = -0.002327_f4+ex*(0.002327_f4*ONE_f4)
               end if
               y = b0+x*(bb1+x*(bb2+x*bb3))
               
               re = 10._f4**y*bpm
               
            else if (eshape > ONE_f4) then
               
               x = log10(2._f4*rfix/eshape)
               
               if (eshape <= 2._f4) then
                  ex = eshape-ONE_f4
                  b0 = -1.310_f4+ex*(-1.11812_f4+1.310_f4*ONE_f4)
                  bb1 = 0.98968_f4+ex*(0.97084_f4-0.98968_f4*ONE_f4)
                  bb2 = -0.042379_f4+ex*(-0.058810_f4+0.042379_f4*ONE_f4)
                  bb3 = ZERO_f4+ex*(0.002159_f4*ONE_f4)
               else if (eshape <= 10._f4) then
                  ex = (eshape-2._f4)*0.125_f4
                  b0 = -1.11812_f4+ex*(-0.90629_f4+1.11812_f4*ONE_f4)
                  bb1 = 0.97084_f4+ex*(0.90412_f4-0.97084_f4*ONE_f4)
                  bb2 = -0.058810_f4+ex*(-0.059312_f4+0.058810_f4*ONE_f4)
                  bb3 = 0.002159_f4+ex*(0.0029941_f4-0.002159_f4*ONE_f4)
               else
                  ex = 10._f4/eshape
                  b0 = -0.79888_f4+ex*(-0.90629_f4+0.79888_f4*ONE_f4)
                  bb1 = 0.80817_f4+ex*(0.90412_f4-0.80817_f4*ONE_f4)
                  bb2 = -0.030528_f4+ex*(-0.059312_f4+0.030528_f4*ONE_f4)
                  bb3 = ZERO_f4+ex*(0.0029941_f4*ONE_f4)
               end if
               y = b0+x*(bb1+x*(bb2+x*bb3))
               
               re = 10._f4**y*bpm
               
            end if

            ! Adjust <vf> for non-sphericicity:

            vf(k,i,j) = re*rmu/(2._f4*r_shape*rhoa_cgs)

         end if

         if (re > 300._f4) then

            z = ((1.e6_f4*rhoa_cgs**2)/  &
                 (g_cgs*rhop(k,ib,ig)*rmu**4))**ONE_SIXTH_f4
            b0 = (24._f4*vf(k,i,j)*rmu)*0.01_f4
            x  = log(z*b0)
            y  = -5.00015_f4+x*(5.23778_f4-x*(2.04914_f4-x*(0.475294_f4  &
                 - x*(0.0542819_f4-x*0.00238449_f4))))
            if (y < -675._f4) then
               y = -675._f4
            else if(y >= 741._f4) then
               y = 741._f4
            end if

            re = z*exp(y)*bpm
            vf(k,i,j) = re*rmu/(2._f4*r(ib,ig)*rhoa_cgs)

         end if

      end do
   end do
end do

end subroutine comp_aer_vfall


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_coremc_kludge(max_log10_sumbpmc, sumbpmc, binpmc,  &
     binpfracmc, aer_type, i_AZ, n_AZ_abs)

!+ Computes aerosol mass concentration kludge (clipping-type)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4, ZERO_f4

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type, n_AZ_abs
integer, dimension(:), intent(IN) :: i_AZ
real(f4), intent(IN) :: max_log10_sumbpmc
real(f4), dimension(:,:,:), intent(IN) :: sumbpmc, binpmc

real(f4), dimension(:,:,:), intent(INOUT) :: binpfracmc  ! IN: Dust core mass
                                                   ! OUT: "New" dust core mass
                                                   !   OR "New" H2O_ice mass

!== Local declarations:

integer :: i, j, k, n1, kk
real(f4) :: core_mass_frac

!=== Executable statements:

n1 = size(binpfracmc, 1)

if (aer_type == 1) then    ! binpfracmc (OUT) is dust core mass

   do j=1,size(binpfracmc, 3)
      do i=1,size(binpfracmc, 2)
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               if (binpmc(k,i,j)-binpfracmc(k,i,j) < 1.e-3_f4) then ! Too much core
                  if (max_log10_sumbpmc-log10(max(tiny(core_mass_frac),  &
                       sumbpmc(k,i,j))) <= 2._f4) then
                     core_mass_frac = 0.03_f4 ! Massive cloud (thick ice shell?)
                  else
                     core_mass_frac = 0.9_f4
                  end if
                  binpfracmc(k,i,j) = binpmc(k,i,j)*core_mass_frac
               else   ! No change
                  binpfracmc(k,i,j) = binpfracmc(k,i,j)
               end if
            end if
         end do
      end do
   end do

else if (aer_type == 2) then    ! binpfracmc (OUT) is H2O_ice mass

   do j=1,size(binpfracmc, 3)
      do i=1,size(binpfracmc, 2)
         do k=1,n1
            kk = i_AZ(k)
            if (kk >= 2 .and. kk <= n_AZ_abs-1) then
               if (binpmc(k,i,j)-binpfracmc(k,i,j) < 1.e-3_f4) then ! Too much core
                  if (max_log10_sumbpmc-log10(max(tiny(core_mass_frac),  &
                       sumbpmc(k,i,j))) <= 2) then
                     core_mass_frac = 0.03_f4 ! Massive cloud (thick ice shell?)
                  else
                     core_mass_frac = 0.9_f4
                  end if
                  binpfracmc(k,i,j) = binpmc(k,i,j)*(ONE_f4-core_mass_frac)
               else   ! No kludge "needed"
                  binpfracmc(k,i,j) = binpmc(k,i,j)-binpfracmc(k,i,j)
               end if
            end if
         end do
      end do
   end do

end if

if (i_AZ(n1) == n_AZ_abs) binpfracmc(n1,:,:) = ZERO_f4
if (i_AZ(1) == 1) then
   ! k=1 val (avoid display intrp problems):
   if (n1 >= 2) then
      binpfracmc(1,:,:) = binpfracmc(2,:,:)
   else
      binpfracmc(1,:,:) = ZERO_f4
   end if
end if

end subroutine comp_aer_coremc_kludge


!```````````````````````````````````````````````````````````````````
subroutine comp_aer_coremc_kludge_2D(max_log10_sumbpmc, sumbpmc, binpmc,  &
     binpfracmc, aer_type)

!+ Computes aerosol mass concentration kludge (clipping-type)

use univ_kind_defs_mod, only: f4
use univ_math_constants_mod, only: ONE_f4, ZERO_f4

implicit none

!=== Argument declarations:

integer, intent(IN) :: aer_type
real(f4), intent(IN) :: max_log10_sumbpmc
real(f4), dimension(:,:), intent(IN) :: sumbpmc, binpmc

real(f4), dimension(:,:), intent(INOUT) :: binpfracmc  ! IN: Dust core mass
                                                   ! OUT: "New" dust core mass
                                                   !   OR "New" H2O_ice mass

!== Local declarations:

integer :: i, j
real(f4) :: core_mass_frac

!=== Executable statements:

if (aer_type == 1) then    ! binpfracmc (OUT) is dust core mass

   do j=1,size(binpfracmc, 2)
      do i=1,size(binpfracmc, 1)
         if (binpmc(i,j)-binpfracmc(i,j) < 1.e-3_f4) then ! Too much core
            if (max_log10_sumbpmc-log10(max(tiny(core_mass_frac),  &
                 sumbpmc(i,j))) <= 2._f4) then
               core_mass_frac = 0.03_f4   ! Massive cloud (thick ice shells?)
            else
               core_mass_frac = 0.9_f4
            end if
            binpfracmc(i,j) = binpmc(i,j)*core_mass_frac
         else   ! No change
            binpfracmc(i,j) = binpfracmc(i,j)
         end if
      end do
   end do

else if (aer_type == 2) then    ! binpfracmc (OUT) is H2O_ice mass

   do j=1,size(binpfracmc, 2)
      do i=1,size(binpfracmc, 1)
         if (binpmc(i,j)-binpfracmc(i,j) < 1.e-3_f4) then ! Too much core
            if (max_log10_sumbpmc-log10(max(tiny(core_mass_frac),  &
                 sumbpmc(i,j))) <= 2) then
               core_mass_frac = 0.03_f4   ! Massive cloud (thick ice shells?)
            else
               core_mass_frac = 0.9_f4
            end if
            binpfracmc(i,j) = binpmc(i,j)*(ONE_f4-core_mass_frac)
         else   ! No kludge "needed"
            binpfracmc(i,j) = binpmc(i,j)-binpfracmc(i,j)
         end if
      end do
   end do

end if

end subroutine comp_aer_coremc_kludge_2D


!```````````````````````````````````````````````````````````````````
subroutine comp_sfc_value_extrp(z_lyr_mid_AGL, a3D_array, s2D_array, mode)

!+ Estimate surface value of 'a3D_array' through extrapolation.

use univ_kind_defs_mod, only: f4
use univ_exception_mod, only: exception, FATAL
use univ_maxlen_string_mod, only: MAX_PROC_NAME_LEN

implicit none

!=== Argument declarations:

integer, intent(IN) :: mode
real(f4), dimension(:,:,:), intent(IN) :: z_lyr_mid_AGL, a3D_array

real(f4), dimension(:,:,:), intent(OUT) :: s2D_array

!== Local declarations:

character(len=MAX_PROC_NAME_LEN), parameter ::  &
     procedure_name = 'comp_sfc_value_extrp'

integer :: i, j, ip
real(f4) :: slope

!=== Executable statements:

select case (mode)
case (1)  ! Requires a3D_array and z_lyr_mid_AGL of the "_faz_" genre
   do ip=1,size(s2D_array, 3)
      do j=1,size(a3D_array, 3)
         do i=1,size(a3D_array, 2)
            slope = (a3D_array(3,i,j)-a3D_array(2,i,j))/  &
                 (z_lyr_mid_AGL(3,i,j)-z_lyr_mid_AGL(2,i,j))         
            s2D_array(i,j,ip) = a3D_array(2,i,j)-z_lyr_mid_AGL(2,i,j)*slope
         end do
      end do
   end do
end select

end subroutine comp_sfc_value_extrp


end module postp_comp_mod