source: codes/icosagcm/trunk/src/dcmip/dcmip2016_simple_physics.f90

MODULE dcmip2016_simple_physics_mod




CONTAINS
!-----------------------------------------------------------------------
!
!  Date:  April 26, 2016 (Version 6)
!
!  Simple Physics Package
!
!  SIMPLE_PHYSICS includes large-scale precipitation, surface fluxes and
!  boundary-leyer mixing. The processes are time-split in that order.
!  A partially implicit formulation is used to foster numerical
!  stability. The routine assumes that the model levels are ordered
!  in a top-down approach, e.g. level 1 denotes the uppermost full model
!  level.
!
!  This routine is based on an implementation which was developed for
!  the NCAR Community Atmosphere Model (CAM). Adjustments for other
!  models may be necessary.
!
!  The routine provides both updates of the state variables u, v, T, q
!  (these are local copies of u,v,T,q within this physics routine) and
!  also collects their time tendencies. The latter might be used to
!  couple the physics and dynamics in a process-split way. For a
!  time-split coupling, the final state should be given to the
!  dynamical core for the next time step.
!
! Test:      0 = Reed and Jablonowski (2011) tropical cyclone test
!            1 = Moist baroclinic instability test
! RJ2012_precip:
!         true  = Turn on Reed and Jablonowski (2012) precip scheme 
!         false = Turn off Reed and Jablonowski (2012) precip scheme
! TC_PBL_mod:
!         true  = Turn on George Bryan PBL mod for tropical cyclone test
!         false = Turn off George Bryan PBL mod (i.e., run as in Reed and Jablonowski (2012))
!
!  SUBROUTINE SIMPLE_PHYSICS(pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test, RJ2012_precip, TC_PBL_mod)
!
!  Input variables:
!     pcols  - number of atmospheric columns (#)
!     pver   - number of model levels (#)
!     dtime  - time step (s)
!     lat    - latitude (radians)
!     t      - temperature at model levels (K)
!     q      - specific humidity at model levels (gm/gm)
!     u      - zonal wind at model levels (m/s)
!     v      - meridional wind at model levels (m/s)
!     pmid   - pressure at model levels (Pa)
!     pint   - pressure at interfaces (Pa)
!     pdel   - layer thickness (Pa)
!     rpdel  - reciprocal of layer thickness (1/Pa)
!     ps     - surface pressure (Pa)
!     test   - test case to use for sea-surface temperatures
!     RJ2012_precip - RJ2012 precip flag
!     TC_PBL_mod    - PCL modification for TC test 
!
!  Output variables:
!     Increments are added into t, q, u, v, pmid, pint, pdel, rpdel and ps
!     which are returned to the routine from which SIMPLE_PHYSICS was
!     called.  Precpitation is returned via precl.
!
!  Change log:
!  v2: removal of some NCAR CAM-specific 'use' associations
!  v3: corrected precl(i) computation, the precipitation rate is now
!      computed via a vertical integral, the previous single-level
!      computation in v2 was a bug
!  v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing
!      the temperature variable: '-t(i,1)'
!  v4: modified and enhanced parameter list to make the routine truly
!      standalone, the number of columns and vertical levels have been
!      added: pcols, pver
!  v4: 'ncol' has been removed, 'pcols' is used instead
!  v5: the sea surface temperature (SST) field Tsurf is now an array,
!      the SST now depends on the latitude
!  v5: addition of the latitude array 'lat' and the flag 'test' in the
!      parameter list
!      if test = 0: constant SST is used, correct setting for the
!                   tropical cyclone test
!      if test = 1: newly added latitude-dependent SST is used,
!                   correct setting for the moist baroclinic wave test
!                   with simple-physics
!  v6: addition of flags for a modified PBL for the TC test and
!      to turn off large-scale condensation scheme when using Kessler physics
!      Included virtual temperature in density calculation in PBL scheme
!      Also, included the virtual temperature, instead of temperature, for
!      the calculation of rho in the PBL scheme
!      (v6_1) Minor specification and generalization fixes.
!      
! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone
!            simulations of intermediate complexity: A test case for AGCMs,
!            J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099
!-----------------------------------------------------------------------

SUBROUTINE SIMPLE_PHYSICS(pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test, RJ2012_precip, TC_PBL_mod)

  ! use physics_types     , only: physics_dme_adjust   ! This is for CESM/CAM
  ! use cam_diagnostics,    only: diag_phys_writeout   ! This is for CESM/CAM

   implicit none

   integer, parameter :: r8 = selected_real_kind(12)

!
! Input arguments - MODEL DEPENDENT
!
   integer, intent(in)  :: pcols        ! Set number of atmospheric columns
   integer, intent(in)  :: pver         ! Set number of model levels
   real(r8), intent(in) :: dtime        ! Set model physics timestep
   real(r8), intent(in) :: lat(pcols)   ! Latitude
   integer, intent(in)  :: test         ! Test number
   logical, intent(in)  :: RJ2012_precip
   logical, intent(in)  :: TC_PBL_mod

!
! Input/Output arguments
!
!  pcols is the maximum number of vertical columns per 'chunk' of atmosphere
!
   real(r8), intent(inout) :: t(pcols,pver)      ! Temperature at full-model level (K)
   real(r8), intent(inout) :: q(pcols,pver)      ! Specific Humidity at full-model level (kg/kg)
   real(r8), intent(inout) :: u(pcols,pver)      ! Zonal wind at full-model level (m/s)
   real(r8), intent(inout) :: v(pcols,pver)      ! Meridional wind at full-model level (m/s)
   real(r8), intent(inout) :: pmid(pcols,pver)   ! Pressure is full-model level (Pa)
   real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa)
   real(r8), intent(inout) :: pdel(pcols,pver)   ! Layer thickness (Pa)
   real(r8), intent(in   ) :: rpdel(pcols,pver)  ! Reciprocal of layer thickness (1/Pa)
   real(r8), intent(inout) :: ps(pcols)          ! Surface Pressue (Pa)

!
! Output arguments
!
   real(r8), intent(out) :: precl(pcols)         ! Precipitation rate (m_water / s)

!
!---------------------------Local workspace-----------------------------
!

! Integers for loops

   integer  i,k                         ! Longitude, level indices

! Physical Constants - Many of these may be model dependent

   real(r8) gravit                      ! Gravity
   real(r8) rair                        ! Gas constant for dry air
   real(r8) cpair                       ! Specific heat of dry air
   real(r8) latvap                      ! Latent heat of vaporization
   real(r8) rh2o                        ! Gas constant for water vapor
   real(r8) epsilo                      ! Ratio of gas constant for dry air to that for vapor
   real(r8) zvir                        ! Constant for virtual temp. calc. =(rh2o/rair) - 1
   real(r8) a                           ! Reference Earth's Radius (m)
   real(r8) omega                       ! Reference rotation rate of the Earth (s^-1)
   real(r8) pi                          ! pi

! Simple Physics Specific Constants 

!++++++++
   real(r8) Tsurf(pcols)                ! Sea Surface Temperature (constant for tropical cyclone)
!++++++++                                 Tsurf needs to be dependent on latitude for the
                                        ! moist baroclinic wave test, adjust

   real(r8) SST_TC                      ! Sea Surface Temperature for tropical cyclone test
   real(r8) T0                          ! Control temp for calculation of qsat
   real(r8) e0                          ! Saturation vapor pressure at T0 for calculation of qsat
   real(r8) rhow                        ! Density of Liquid Water
   real(r8) p0                          ! Constant for calculation of potential temperature
   real(r8) Cd0                         ! Constant for calculating Cd from Smith and Vogl 2008
   real(r8) Cd1                         ! Constant for calculating Cd from Smith and Vogl 2008
   real(r8) Cm                          ! Constant for calculating Cd from Smith and Vogl 2008
   real(r8) v20                         ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
   real(r8) C                           ! Drag coefficient for sensible heat and evaporation
   real(r8) T00                         ! Horizontal mean T at surface for moist baro test
   real(r8) u0                          ! Zonal wind constant for moist baro test
   real(r8) latw                        ! halfwidth for  for baro test
   real(r8) eta0                        ! Center of jets (hybrid) for baro test
   real(r8) etav                        ! Auxiliary variable for baro test
   real(r8) q0                          ! Maximum specific humidity for baro test
   real(r8) kappa                       ! von Karman constant

! Physics Tendency Arrays
  real(r8) dtdt(pcols,pver)             ! Temperature tendency
  real(r8) dqdt(pcols,pver)             ! Specific humidity tendency
  real(r8) dudt(pcols,pver)             ! Zonal wind tendency
  real(r8) dvdt(pcols,pver)             ! Meridional wind tendency

! Temporary variables for tendency calculations

   real(r8) tmp                         ! Temporary
   real(r8) qsat                        ! Saturation vapor pressure
   real(r8) qsats                       ! Saturation vapor pressure of SST

! Variables for Boundary Layer Calculation

   real(r8) wind(pcols)                 ! Magnitude of Wind
   real(r8) Cd(pcols)                   ! Drag coefficient for momentum
   real(r8) Km(pcols,pver+1)            ! Eddy diffusivity for boundary layer calculations 
   real(r8) Ke(pcols,pver+1)            ! Eddy diffusivity for boundary layer calculations
   real(r8) rho                         ! Density at lower/upper interface
   real(r8) za(pcols)                   ! Heights at midpoints of first model level
   real(r8) zi(pcols,pver+1)            ! Heights at model interfaces
   real(r8) dlnpint                     ! Used for calculation of heights
   real(r8) pbltop                      ! Top of boundary layer
   real(r8) zpbltop                     ! Top of boundary layer for George Bryan Modifcation
   real(r8) pblconst                    ! Constant for the calculation of the decay of diffusivity 
   real(r8) CA(pcols,pver)              ! Matrix Coefficents for PBL Scheme 
   real(r8) CC(pcols,pver)              ! Matrix Coefficents for PBL Scheme 
   real(r8) CE(pcols,pver+1)            ! Matrix Coefficents for PBL Scheme
   real(r8) CAm(pcols,pver)             ! Matrix Coefficents for PBL Scheme
   real(r8) CCm(pcols,pver)             ! Matrix Coefficents for PBL Scheme
   real(r8) CEm(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
   real(r8) CFu(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
   real(r8) CFv(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
   real(r8) CFt(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
   real(r8) CFq(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme


! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to
! conserve the mass of the dry air

   real(r8) qini(pcols,pver)            ! Initial specific humidity

!===============================================================================
!
! Physical Constants - MAY BE MODEL DEPENDENT
!
!===============================================================================
   gravit = 9.80616_r8                   ! Gravity (9.80616 m/s^2)
   rair   = 287.0_r8                     ! Gas constant for dry air: 287 J/(kg K)
   cpair  = 1.0045e3_r8                  ! Specific heat of dry air: here we use 1004.5 J/(kg K)
   latvap = 2.5e6_r8                     ! Latent heat of vaporization (J/kg)
   rh2o   = 461.5_r8                     ! Gas constant for water vapor: 461.5 J/(kg K)
   epsilo = rair/rh2o                    ! Ratio of gas constant for dry air to that for vapor
   zvir   = (rh2o/rair) - 1._r8          ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608
   a      = 6371220.0_r8                 ! Reference Earth's Radius (m)
   omega  = 7.29212d-5                   ! Reference rotation rate of the Earth (s^-1)
   pi     = 4._r8*atan(1._r8)            ! pi

!===============================================================================
!
! Local Constants for Simple Physics
!
!===============================================================================
      C        = 0.0011_r8      ! From Simth and Vogl 2008
      SST_TC   = 302.15_r8      ! Constant Value for SST
      T0       = 273.16_r8      ! control temp for calculation of qsat
      e0       = 610.78_r8      ! saturation vapor pressure at T0 for calculation of qsat
      rhow     = 1000.0_r8      ! Density of Liquid Water 
      Cd0      = 0.0007_r8      ! Constant for Cd calc. Simth and Vogl 2008
      Cd1      = 0.000065_r8    ! Constant for Cd calc. Simth and Vogl 2008
      Cm       = 0.002_r8       ! Constant for Cd calc. Simth and Vogl 2008
      v20      = 20.0_r8        ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
      p0       = 100000.0_r8    ! Constant for potential temp calculation
      pbltop   = 85000._r8      ! Top of boundary layer in p
      zpbltop  = 1000._r8       ! Top of boundary layer in z
      pblconst = 10000._r8      ! Constant for the calculation of the decay of diffusivity
      T00      = 288.0_r8         ! Horizontal mean T at surface for moist baro test
      u0       = 35.0_r8          ! Zonal wind constant for moist baro test
      latw     = 2.0_r8*pi/9.0_r8 ! Halfwidth for  for baro test
      eta0     = 0.252_r8         ! Center of jets (hybrid) for baro test
      etav     = (1._r8-eta0)*0.5_r8*pi ! Auxiliary variable for baro test
      q0       = 0.021_r8         ! Maximum specific humidity for baro test
      kappa    = 0.4_r8         ! von Karman constant

!===============================================================================
!
! Definition of local arrays
!
!===============================================================================
!
! Calculate hydrostatic height
!
     do i=1,pcols
        dlnpint = log(ps(i)) - log(pint(i,pver))  ! ps(i) is identical to pint(i,pver+1), note: this is the correct sign (corrects typo in JAMES paper) 
        za(i) = rair/gravit*t(i,pver)*(1._r8+zvir*q(i,pver))*0.5_r8*dlnpint
        zi(i,pver+1) = 0.0_r8
     end do
!
! Set Initial Specific Humidity
!
     qini(:pcols,:pver) = q(:pcols,:pver)
!
! Set Sea Surface Temperature (constant for tropical cyclone)
! Tsurf needs to be dependent on latitude for moist baroclinic wave test
! Tsurf needs to be constant for tropical cyclone test
!
     if (test .eq. 1) then ! Moist Baroclinic Wave Test
        do i=1,pcols
           Tsurf(i) = (T00 + pi*u0/rair * 1.5_r8 * sin(etav) * (cos(etav))**0.5_r8 *                 &
                     ((-2._r8*(sin(lat(i)))**6 * ((cos(lat(i)))**2 + 1._r8/3._r8) + 10._r8/63._r8)* &
                     u0 * (cos(etav))**1.5_r8  +                                                    &
                     (8._r8/5._r8*(cos(lat(i)))**3 * ((sin(lat(i)))**2 + 2._r8/3._r8) - pi/4._r8)*a*omega*0.5_r8 ))/ &
                     (1._r8+zvir*q0*exp(-(lat(i)/latw)**4))

        end do
     else if (test .eq. 0) then ! Tropical Cyclone Test
        do i=1,pcols
           Tsurf(i) = SST_TC
        end do
     end if

!===============================================================================
!
! Set initial physics time tendencies and precipitation field to zero
!
!===============================================================================
     dtdt(:pcols,:pver)  = 0._r8            ! initialize temperature tendency with zero
     dqdt(:pcols,:pver)  = 0._r8            ! initialize specific humidity tendency with zero
     dudt(:pcols,:pver)  = 0._r8            ! initialize zonal wind tendency with zero
     dvdt(:pcols,:pver)  = 0._r8            ! initialize meridional wind tendency with zero
     precl(:pcols) = 0._r8                  ! initialize precipitation rate with zero

!===============================================================================
!
! Large-Scale Condensation and Precipitation
!
!===============================================================================

      if (RJ2012_precip) then
!
! Calculate Tendencies
!
      do k=1,pver
         do i=1,pcols
            qsat = epsilo*e0/pmid(i,k)*exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0))
            if (q(i,k) > qsat) then
               tmp  = 1._r8/dtime*(q(i,k)-qsat)/(1._r8+(latvap/cpair)*(epsilo*latvap*qsat/(rair*t(i,k)**2)))
               dtdt(i,k) = dtdt(i,k)+latvap/cpair*tmp
               dqdt(i,k) = dqdt(i,k)-tmp
               precl(i)  = precl(i)+tmp*pdel(i,k)/(gravit*rhow)
            end if
         end do
      end do
!
! Update moisture and temperature fields from Larger-Scale Precipitation Scheme
!
      do k=1,pver
         do i=1,pcols
            t(i,k) =  t(i,k) + dtdt(i,k)*dtime
            q(i,k) =  q(i,k) + dqdt(i,k)*dtime
         end do
      end do

!===============================================================================
! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC
!
! note: The variables, as done in many GCMs, are written to the history file
!       after the moist physics process.  This ensures that the moisture fields
!       are somewhat in equilibrium.
!===============================================================================
  !  call diag_phys_writeout(state)   ! This is for CESM/CAM
    
      end if
     
!===============================================================================
!
! Turbulent mixing coefficients for the PBL mixing of horizontal momentum,
! sensible heat and latent heat
!
! We are using Simplified Ekman theory to compute the diffusion coefficients
! Kx for the boundary-layer mixing. The Kx values are calculated at each time step
! and in each column.
!
!===============================================================================!
! Compute magnitude of the wind and drag coeffcients for turbulence scheme
!
     do i=1,pcols
        wind(i) = sqrt(u(i,pver)**2+v(i,pver)**2)
     end do
     do i=1,pcols
        if( wind(i) .lt. v20) then
           Cd(i) = Cd0+Cd1*wind(i) 
        else
           Cd(i) = Cm
        endif
     end do

     if (TC_PBL_mod) then !Bryan TC PBL Modification 
     do k=pver,1,-1
        do i=1,pcols
           dlnpint = log(pint(i,k+1)) - log(pint(i,k))
           zi(i,k) = zi(i,k+1)+rair/gravit*t(i,k)*(1._r8+zvir*q(i,k))*dlnpint
           if( zi(i,k) .le. zpbltop) then
              Km(i,k) = kappa*sqrt(Cd(i))*wind(i)*zi(i,k)*(1._r8-zi(i,k)/zpbltop)*(1._r8-zi(i,k)/zpbltop)
              Ke(i,k) = kappa*sqrt(C)*wind(i)*zi(i,k)*(1._r8-zi(i,k)/zpbltop)*(1._r8-zi(i,k)/zpbltop) 
           else
              Km(i,k) = 0.0_r8
              Ke(i,k) = 0.0_r8
           end if 
        end do
     end do     
     else ! Reed and Jablonowski (2012) Configuration
     do k=1,pver
        do i=1,pcols
           if( pint(i,k) .ge. pbltop) then
              Km(i,k) = Cd(i)*wind(i)*za(i) 
              Ke(i,k) = C*wind(i)*za(i)
           else
              Km(i,k) = Cd(i)*wind(i)*za(i)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2)
              Ke(i,k) = C*wind(i)*za(i)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2)
           end if
        end do
     end do
     end if
!===============================================================================
! Update the state variables u, v, t, q with the surface fluxes at the
! lowest model level, this is done with an implicit approach
! see Reed and Jablonowski (JAMES, 2012)
!
! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1
! Tsurf needs to be dependent on latitude for the
! moist baroclinic wave test 
!===============================================================================
     do i=1,pcols
        qsats = epsilo*e0/ps(i)*exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0))
        dudt(i,pver) = dudt(i,pver) + (u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))-u(i,pver))/dtime
        dvdt(i,pver) = dvdt(i,pver) + (v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))-v(i,pver))/dtime
        u(i,pver) = u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
        v(i,pver) = v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
        dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))-t(i,pver))/dtime 
        t(i,pver) = (t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))  
        dqdt(i,pver) = dqdt(i,pver) +((q(i,pver)+C*wind(i)*qsats*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))-q(i,pver))/dtime
        q(i,pver) = (q(i,pver)+C*wind(i)*qsats*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))
     end do

!===============================================================================
! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012)
!===============================================================================
! Calculate Diagonal Variables for Implicit PBL Scheme
!
      do k=1,pver-1
         do i=1,pcols
            rho = (pint(i,k+1)/(rair*(t(i,k+1)*(1._r8+zvir*q(i,k+1))+t(i,k)*(1._r8+zvir*q(i,k)))/2.0_r8)) 
            CAm(i,k) = rpdel(i,k)*dtime*gravit*gravit*Km(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))    
            CCm(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Km(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
            CA(i,k) = rpdel(i,k)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
            CC(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho/(pmid(i,k+1)-pmid(i,k))
         end do
      end do
      do i=1,pcols
         CAm(i,pver) = 0._r8
         CCm(i,1) = 0._r8
         CEm(i,pver+1) = 0._r8
         CA(i,pver) = 0._r8
         CC(i,1) = 0._r8
         CE(i,pver+1) = 0._r8
         CFu(i,pver+1) = 0._r8
         CFv(i,pver+1) = 0._r8
         CFt(i,pver+1) = 0._r8
         CFq(i,pver+1) = 0._r8 
      end do
      do i=1,pcols
         do k=pver,1,-1
            CE(i,k) = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) 
            CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
            CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1))/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
            CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1))/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
            CFt(i,k) = ((p0/pmid(i,k))**(rair/cpair)*t(i,k)+CA(i,k)*CFt(i,k+1))/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) 
            CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1))/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
       end do
      end do
!
! Calculate the updated temperaure and specific humidity and wind tendencies
!
! First we need to calculate the tendencies at the top model level
!
      do i=1,pcols
            dudt(i,1)  = dudt(i,1)+(CFu(i,1)-u(i,1))/dtime
            dvdt(i,1)  = dvdt(i,1)+(CFv(i,1)-v(i,1))/dtime
            u(i,1)    = CFu(i,1)
            v(i,1)    = CFv(i,1)
            dtdt(i,1)  = dtdt(i,1)+(CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)-t(i,1))/dtime
            t(i,1)    = CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)
            dqdt(i,1)  = dqdt(i,1)+(CFq(i,1)-q(i,1))/dtime
            q(i,1)  = CFq(i,1)
      end do

      do i=1,pcols
         do k=2,pver
            dudt(i,k)  = dudt(i,k)+(CEm(i,k)*u(i,k-1)+CFu(i,k)-u(i,k))/dtime
            dvdt(i,k)  = dvdt(i,k)+(CEm(i,k)*v(i,k-1)+CFv(i,k)-v(i,k))/dtime
            u(i,k)    = CEm(i,k)*u(i,k-1)+CFu(i,k) 
            v(i,k)    = CEm(i,k)*v(i,k-1)+CFv(i,k)
            dtdt(i,k)  = dtdt(i,k)+  &
                        ((CE(i,k)*t(i,k-1)*(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k))*(pmid(i,k)/p0)**(rair/cpair)-t(i,k))/dtime 
            t(i,k)    = (CE(i,k)*t(i,k-1)*(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k))*(pmid(i,k)/p0)**(rair/cpair)
            dqdt(i,k)  = dqdt(i,k)+(CE(i,k)*q(i,k-1)+CFq(i,k)-q(i,k))/dtime
            q(i,k)  = CE(i,k)*q(i,k-1)+CFq(i,k)
         end do
      end do

!===============================================================================
!
! Dry Mass Adjustment - THIS PROCESS WILL BE MODEL SPECIFIC
!
! note: Care needs to be taken to ensure that the model conserves the total
!       dry air mass. Add your own routine here.
!===============================================================================
  !  call physics_dme_adjust(state, tend, qini, dtime)   ! This is for CESM/CAM

   return
end subroutine SIMPLE_PHYSICS 


END MODULE dcmip2016_simple_physics_mod