#include <misc.h>
#include <params.h>
      subroutine tphysbc(lat     ,clat    ,clon    ,coslat  ,
     $                   oro     ,pmidm1  ,pintm1  ,pdelm1  ,omga    ,
     $                   um1     ,vm1     ,tm1     ,qm1     ,psm1    ,
     $                   pmlnm1  ,pilnm1  ,zm      ,asdir   ,asdif   ,
     $                   aldir   ,aldif   ,snowh   ,cnt     ,cnb     ,
     $                   qrs     ,qrl     ,flwds   ,fsns    ,tvm     ,
     $                   thm     ,cmfmc   ,lwup    ,precst  ,srfrad  ,
     $                   sols    ,soll    ,solsd   ,solld   )
C-----------------------------------------------------------------------
C
C Tendency physics BEFORE coupling to land, sea, and ice models.
C Call physics subroutines and compute the following:
C     o cloud calculations (cloud fraction, emissivity, etc.)
C     o radiation calculations
C Pass surface fields for separate surface flux calculations
C Dump appropriate fields to history file.
C
C----------------------------Code History-------------------------------
C
C Original version:  CCM1
C Standardized:      J. Rosinski, June 1992
C                    T. Acker, March 1996
C Reviewed:          J. Kiehl, J. Hack, August 1992
C Include Ice:       J. Kiehl, January 1993
C Reviewed:          J. Kiehl, April 1996
C
C-----------------------------------------------------------------------
c
c $Id: tphysbc.F,v 1.18.2.1 2000/09/03 22:19:07 rosinski Exp $
c
C-----------------------------------------------------------------------
      use precision
      use pmgrid
      use history, only: nhtfrq
C-----------------------------------------------------------------------
#include <implicit.h>
C------------------------------Parameters-------------------------------
#include <comcon.h>
C-----------------------------------------------------------------------
#include <comctl.h>
C-----------------------------------------------------------------------
#include <comtim.h>
C------------------------------Arguments--------------------------------
C
C Input arguments
C
      integer lat               ! Model latitude index
      real(r8) clat                 ! Current latitude(radians)
      real(r8) clon(plon)           ! Centered longitude(radians)
      real(r8) coslat               ! Cosine latitude
      real(r8) oro(plond)           ! Sfc type flag (ocean,land,sea ice)
      real(r8) pmidm1(plond,plev)   ! Midpoint pressures
      real(r8) pintm1(plond,plevp)  ! Interface pressures
      real(r8) pdelm1(plond,plev)   ! Layer pressure depth
      real(r8) omga(plond,plev)     ! Vertical pressure velocity
      real(r8) um1(plond,plev)      ! U-wind (time n-1)
      real(r8) vm1(plond,plev)      ! V-wind (time n-1)
      real(r8) tm1(plond,plev)      ! Temperature
      real(r8) qm1(plond,plev,pcnst+pnats) ! Specific humidity
      real(r8) psm1(plond)          ! Surface pressure
      real(r8) pmlnm1(plond,plev)   ! ln(pmidm1)
      real(r8) pilnm1(plond,plevp)  ! ln(pintm1)
      real(r8) zm(plond,plev)       ! Height of midpoints (above surface)
      real(r8) asdir(plond)         ! Albedo: shortwave, direct
      real(r8) asdif(plond)         ! Albedo: shortwave, diffuse
      real(r8) aldir(plond)         ! Albedo: longwave, direct
      real(r8) aldif(plond)         ! Albedo: longwave, diffuse
      real(r8) cmfmc(plond,plev)    ! Convective mass flux--m sub c
      real(r8) lwup(plond)          ! Surface longwave up flux
C
C Output arguments
C
      real(r8) snowh(plond)         ! Snow depth (liquid water equivalent)
      real(r8) cnt(plond)           ! Top level of convective activity
      real(r8) cnb(plond)           ! Lowest level of convective activity
      real(r8) qrs(plond,plev)      ! Shortwave heating rate
      real(r8) qrl(plond,plev)      ! Longwave  heating rate
      real(r8) flwds(plond)         ! Surface longwave down flux
      real(r8) fsns(plond)          ! Surface solar absorbed flux
      real(r8) tvm(plond,plev)      ! Virtual temperature
      real(r8) thm(plond,plev)      ! Potential temperature
      real(r8) precst(plond)        ! Total snow precip. rate
      real(r8) srfrad(plond)        ! Net surface radiative flux (watts/m**2)
      real(r8) sols(plond)          ! Direct beam solar rad. onto srf (sw)
      real(r8) soll(plond)          ! Direct beam solar rad. onto srf (lw)
      real(r8) solsd(plond)         ! Diffuse solar radiation onto srf (sw)
      real(r8) solld(plond)         ! Diffuse solar radiation onto srf (lw)
C
C---------------------------Local workspace-----------------------------
C
      integer  i,k              ! Longitude, level indices
      real(r8) clwp(plond,plev)     ! Presribed cloud liq. h2o path 
      real(r8) rel(plond,plev)      ! Liquid cloud particle effective radius
      real(r8) rei(plond,plev)      ! Ice effective drop size (microns)
      real(r8) fice(plond,plev)     ! Fractional ice content within cloud
      real(r8) cld(plond,plevp)     ! Fraction cloud cover
      real(r8) effcld(plond,plevp)  ! Effective cloud=cld*emis
      real(r8) emis(plond,plev)     ! Cloud longwave emissivity
      real(r8) clc(plond)           ! Total convective cloud (cloud scheme)
      real(r8) rpdeli(plond,plev)   ! 1./(pmid(k+1)-pmid(k))
      real(r8) rpdel(plond,plev)    ! 1./(pint(k+1)-pint(k))
      real(r8) cltot(plond)         ! Diagnostic total cloud cover
      real(r8) cllow(plond)         !       "     low  cloud cover
      real(r8) clmed(plond)         !       "     mid  cloud cover
      real(r8) clhgh(plond)         !       "     hgh  cloud cover
      real(r8) coszrs(plond)        ! Cosine solar zenith angle
C        
      real(r8) actlwp(plond,plev)   ! Actual cloud liquid h2o path 
c      real(r8) totlwp(plond)       ! Actual column cloud liquid h2o path 
c      real(r8) sum                 ! Sum of actlwp
      real(r8) cloudup(plond,plev)  ! Cloud fraction during omega < 0.
C                                 
      real(r8) tpw(plond)           ! Total precipitable water
      real(r8) hl (plond)           ! Liquid water scale height
C
C-----------------------------------------------------------------------
C
C Cosine solar zenith angle for current time step
C
      call zenith (calday  ,clat    , clon, coszrs  )
C
C Precompute frequently used variables
C
      do k=1,plev
         do i=1,plon
            rpdel(i,k) = 1./pdelm1(i,k)
         end do
      end do
C
C 1./delta-p needed by cloud routines
C
      if (dosw .or. dolw) then
         do k=1,plev-1
            do i=1,plon
               rpdeli(i,k) = 1./(pmidm1(i,k+1) - pmidm1(i,k))
            end do
         end do
C
C Compute cloud properties for input to radiation
C
         call cldint(pmidm1  ,rpdeli  ,tm1     ,qm1     ,omga    ,
     $               cnt     ,cnb     ,pilnm1  ,pmlnm1  ,tvm     ,
     $               zm      ,cld     ,clc     ,clwp    ,emis    ,
     $               effcld  ,cltot   ,cllow   ,clmed   ,clhgh   ,
     $               oro     ,rel     ,rei     ,fice    ,pdelm1  ,
     $               tpw     ,hl      ,cmfmc   ,psm1    ,snowh   )
C
C Dump cloud field information to history tape buffer (diagnostics)
C
         call outfld('CLOUD   ',cld    ,plond,lat)
         call outfld('EFFCLD  ',effcld ,plond,lat)
         call outfld('CNVCLD  ',clc    ,plond,lat)
         call outfld('CLDTOT  ',cltot  ,plond,lat)
         call outfld('CLDLOW  ',cllow  ,plond,lat)
         call outfld('CLDMED  ',clmed  ,plond,lat)
         call outfld('CLDHGH  ',clhgh  ,plond,lat)
         call outfld('LWSH    ',hl     ,plond,lat)
C
C Special diagnostic cloud water fields:
C
         call outfld('SETLWP  ',clwp   ,plond,lat)
C
C Outfield special monthly cloud diagnostics
C      
         if (nhtfrq(1).eq.0) then
            do k=1,plev
               do i=1,plon
                  cloudup(i,k) = 0.
                  if (omga(i,k).lt.0.) then
                     cloudup(i,k) = cld(i,k)
                  end if
               end do
            end do
            call outfld('CLOUDUP ',cloudup ,plond,lat)
         end if
C
C Compute actual liquid water path
C
         do k=1,plev
           do i=1,plon
             actlwp(i,k) = cld(i,k)*clwp(i,k)
           end do
         end do
         call outfld('CLDLWP  ',actlwp ,plond,lat)
c         do i=1,plon
c           sum = 0.
c           do k=1,plev
c             sum = sum + actlwp(i,k)
c           end do
c           totlwp(i) = sum
c         end do
c         call outfld('TOTLWP  ',totlwp ,plond,lat)
C
C Complete radiation calculations
C
         call radctl(         clat    ,coslat  ,lat     ,lwup    ,
     $               pmidm1  ,pintm1  ,pmlnm1  ,pilnm1  ,tm1     ,
     $               qm1     ,cld     ,effcld  ,clwp    ,coszrs  ,
     $               asdir   ,asdif   ,aldir   ,aldif   ,fsns    ,
     $               qrs     ,qrl     ,flwds   ,rel     ,rei     ,
     $               fice    ,sols    ,soll    ,solsd   ,solld   )
      end if
C
C Compute net surface radiative flux for use by surface temperature code.
C Note that units have already been converted to mks in RADCTL.  Since
C fsns and flwds are in the buffer, array values will be carried across
C timesteps when the radiation code is not invoked.
C
      do i=1,plon
         srfrad(i) = fsns(i) + flwds(i)
      end do
      call outfld('SRFRAD  ',srfrad,plond,lat)
C
C Save atmospheric fields to force surface models
C
      call srfxfer(lat           ,psm1        ,um1(1,plev) ,
     $                  vm1(1,plev)   ,tm1(1,plev),
     $             qm1(1,plev,1) ,thm(1,plev) ,zm(1,plev)  ,
     $                  pmidm1        ,rhoh2o     ,
     $             rpdel(1,plev) ,fsns        ,precst      )
C
      return
C
      end
 
