#include <params.h>

      module mo_physlic

      contains

      subroutine physlic( ip, lat, dtime, calday, ncdate, &
                          ncsec, nstep, ps, oro, phis, &
                          pt0, sh0, t1, sh1, u, &
                          v, ts, taux, tauy, hflx, &
                          shflx, qrad, etadot, omga, cgs, &
                          kvh, zmu, zmd, zeu, hketa, &
                          hkbeta, cldtop, cldbot, zm, zi, &
                          cmfdqr, cldfr, concld, nrain, nevapr, &
                          cwat, pdmfuen, pdmfude, pdmfden, pdmfdde, &
                          as, ptadv, shadv, vdiffbj, vdiffaj, &
                          convj, cloudj, snow, fsds, soilw, &
                          precip, ts_avg, fsds_avg, plonl, platl, elapsed )
!-----------------------------------------------------------------------
! 	... physics parameterizations and output fields
!-----------------------------------------------------------------------

      use mo_mpi,       only : thisnode, base_lat
      use mo_grid,      only : plev, plevp, pcnst, pnats
      use plevs,        only : plevs0
      use mo_constants, only : rgrav, latwts
      use mo_control,   only : diconvccm, arconvccm, dicldphys, chemistry, &
                               tropchem, rdomega, advqth, addqrad, divdiff, arvdiff, &
                               xactive_drydep, xactive_h2o
      use mo_histout,   only : outfld, hfile, moz_file_cnt, sim_file_cnt
      use mo_surface,   only : sflxdr
      use mo_chemdr,    only : chemdr
      use mo_vdiff,     only : vdiffdr
      use mo_cloud,     only : clouddr
      use mo_aphys,     only : aphys, aphysar
      use chem_mods,    only : adv_mass
      use mass_diags,   only : fdiags, advn_flux, hsa_fac
      use mo_virtem,    only : virtem
      use mo_z0,        only : zmid, zint, zr2a, z2phi
      use mo_pottemp,   only : t2pt, pt2t
      use mo_mass,      only : gamdry2
      use mo_lghtning,  only : do_lghtning, lght_no_prod
      use eslookup,     only : aqsat
      use mo_local_time,only : local_time, local_time_char
      use m_tracname,   only : tracnam

      implicit none

!-----------------------------------------------------------------------
! 	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ip                        ! longitude index for this slice
      integer, intent(in) :: lat                       ! latitude index for this slice
      integer, intent(in) :: platl                     ! latitude tile dimension
      integer, intent(in) :: plonl                     ! longitude tile dimension
      integer, intent(in) :: ncdate                    ! date at beginning of current timestep in yymmdd format
      integer, intent(in) :: ncsec                     ! seconds relative to ncdate
      integer, intent(in) :: nstep                     ! timestep counter

      real, intent(in)    :: dtime                     ! timestep (s)
      real, intent(in)    :: calday                    ! calendar day in range [1., 366.)
      real, intent(in)    :: ps(plonl)                 ! surface pressure (Pa)
      real, intent(in)    :: oro(plonl)                ! orography
      real, intent(in)    :: phis(plonl)               ! surface geopotential (m^2/s^2)
      real, intent(in)    :: ts(plonl)                 ! surface temperature (K)
      real, intent(in)    :: ts_avg(plonl)             ! average surface temperature (K)
      real, intent(in)    :: taux(plonl)               ! x surface stress (n)
      real, intent(in)    :: tauy(plonl)               ! y surface stress (n)
      real, intent(in)    :: hflx(plonl)               ! surface sensible heat flux (w/m2)
      real, intent(in)    :: shflx(plonl)              ! surface water vapor flux (kg/m2/s)
      real, intent(in)    :: snow(plonl)               ! snow height (m)
      real, intent(in)    :: fsds(plonl)               ! srf direct radiation (W/m^2)
      real, intent(in)    :: fsds_avg(plonl)           ! average srf direct radiation (W/m^2)
      real, intent(inout) :: soilw(plonl)              ! soil moisture fraction
      real, intent(in)    :: u(plonl,plev)             ! u wind (m/s)
      real, intent(in)    :: v(plonl,plev)             ! v wind (m/s)
      real, intent(in)    :: etadot(plonl,plevp)       ! vertical velocity (etadot)
      real, intent(in)    :: pt0(plonl,plev)           ! potential temperature from input data
                                                       ! interpolated to beginning of current timestep, (k)
      real, intent(in)    :: sh0(plonl,plev)           ! specific humidity from input data interpolated
                                                       ! to beginning of current timestep, (kg/kg)
      real, intent(in)    :: t1(plonl,plev)            ! temperature from input data interpolated to
                                                       ! end of current timestep, (k)
      real, intent(in)    :: sh1(plonl,plev)           ! specific humidity from input data interpolated
                                                       ! to end of current timestep, (kg/kg)
      real, intent(in)    :: zmu(plonl,plev)           ! mu2 from conv_ccm, kg/m2/s
      real, intent(in)    :: zmd(plonl,plev)           ! md2 from conv_ccm, kg/m2/s
      real, intent(in)    :: zeu(plonl,plev)           ! eu2 from conv_ccm, 1/s
      real, intent(in)    :: hketa(plonl,plev)         ! convective mass flux, hack convection (kg/m^2/s)
      real, intent(in)    :: hkbeta(plonl,plev)        ! overshoot parameter, hack convection (fraction)
      real, intent(in)    :: cwat(plonl,plev)          ! cloud water (kg/kg)
      real, intent(in)    :: pdmfuen(plonl,plev)       ! entrainment into updraft  
      real, intent(in)    :: pdmfude(plonl,plev)       ! detrainment into updraft  
      real, intent(in)    :: pdmfden(plonl,plev)       ! entrainment into downdraft
      real, intent(in)    :: pdmfdde(plonl,plev)       ! detrainment into downdraft
      real, intent(in)    :: qrad(plonl,plev)          ! radiative heating tendency (k/s)
      real, intent(inout) :: cgs(plonl,plevp)          ! counter-gradient coefficient
      real, intent(inout) :: kvh(plonl,plevp)          ! vertical diffusion coefficient
      real, intent(inout) :: cldfr(plonl,plevp)        ! cloud fraction
      real, intent(inout) :: zi(plonl,plevp)           ! potential height above surface at interfaces (m)

      real, intent(inout) :: as(plonl,plev,pcnst)      ! advected species concentrations (kg/kg)
      real, intent(inout) :: cldtop(plonl)             ! cloud top level index
      real, intent(inout) :: cldbot(plonl)             ! cloud bottom level index
      real, intent(inout) :: precip(plonl)             ! srf precipitation
      real, intent(inout) :: elapsed(10)               ! timing
      real, intent(inout) :: ptadv(plonl,plev)         ! potential temperature from input data t0
                                                       ! advected one timestep on input.  vdiff, qrad,
                                                       ! convection and cloud tendencies added on output.
      real, intent(inout) :: shadv(plonl,plev)         ! specific humidity from input data sh0 ad
                                                       ! advected one timestep on input.  vdiff,
                                                       ! convection and cloud tendencies added on output.
      real, intent(inout) :: zm(plonl,plev)            ! potential height above surface at midpoints (m)
      real, intent(inout) :: omga(plonl,plev)          ! vertical velocity (dp/dt)
      real, intent(inout) :: cmfdqr(plonl,plev)        ! dq/dt due to convective rainout 
      real, intent(inout) :: concld(plonl,plev)        ! convective cloud fraction
      real, intent(inout) :: nevapr(plonl,plev)        ! rate of evaporation of precipitation (1/s)
      real, intent(inout) :: nrain(plonl,plev)         ! rate of release of stratiform precip (1/s)

      real, intent(out)   :: vdiffaj(plonl,pcnst)
      real, intent(out)   :: vdiffbj(plonl,pcnst)
      real, intent(out)   :: convj(plonl,pcnst)
      real, intent(out)   :: cloudj(plonl,pcnst)

!-----------------------------------------------------------------------
! 	... local variables
!-----------------------------------------------------------------------
      real, parameter :: trop_min = 5.e3                ! tropopause min height (m)
      real, parameter :: trop_max = 20.e3               ! tropopause max height (m)
      real, parameter :: max_dtdz = 2.e-3               ! max value of dt/dz for tropopause location (K/m)
      real, parameter :: m2km     = 1.e-3               ! meters to kilometers

      integer :: &
        i, k, m, n, &  ! indices
        istat, &       ! return status
        astat, &       ! allocation status
        file, &        ! history file index
        hndx           ! work index
      integer :: &
        ioro(plonl)    ! nint(oro)

      real :: vfac, dtdz
      real, allocatable :: nas(:,:,:)

      real :: &
        coro(plonl)             ! orography for cloud routine
      real, dimension(plonl,plev) :: &
        pmid, &                ! pressure at layer midpoints
        pdel, &                ! pressure difference across layers
        pmln, &                ! log of midpoint pressure
        rpdel, &               ! 1/pdel
        rpdeli, &              ! 1./(pmid(k+1)-pmid(k))
        t0, &                  ! temperature corresponding to pt0
        tadv, &                ! temperature corresponding to ptadv
        tv, &                  ! virtual temperature
        rh, &                  ! relative humidity
        zma, &                 ! absolute potential height at midpoints
        phima, &               ! absolute geopotential at midpoints
        xtmp, &                ! temp storage for output tracers
        cme                    ! rate of cond-evap within the cloud
      real, dimension(plonl,plevp) :: &
        pint, &                ! pressure at layer interfaces
        piln, &                ! log of interface pressure
        zia                    ! absolute potential height at interfaces
      real, dimension(plonl) :: &
        cldtot, &              ! diagnosed tot random overlap cld cover
        wind_speed             ! wind speed at srf

!-----------------------------------------------------------------------
!     	... vertical diffusion
!-----------------------------------------------------------------------
      real :: &
        sflx(plonl,pcnst)      ! surface flux for advected species
      real, dimension(plonl) :: &
        pblh, &                ! planetary boundary layer height
        tpert, &               ! convective temperature excess
        qpert                  ! convective humidity excess

!-----------------------------------------------------------------------
!     	... convective mass flux
!-----------------------------------------------------------------------
      real, dimension(plonl,plev) :: &
        cmfmc, &               ! moist convection cloud mass flux
        zdu, &                 ! du2 from conv_ccm, 1/s
        dlf                    ! detraining cloud water from convection
      real, dimension(plonl) :: &
        precc, &               ! convective precip (m/s)
        precl                  ! large scale precip (m/s)

      real :: tmp(plonl,plev)
      real :: wrk(plonl)
      real :: wrk2d(plonl,1)
      real :: etimer

      character(len=32) :: fldname
      character(len=8)  :: cdate(2)
      character(len=10) :: ctime(2)

!-----------------------------------------------------------------------
!  	... initialize some variables
!-----------------------------------------------------------------------
      ioro(:plonl) = nint( oro(:plonl) )
      coro(:plonl) = oro(:plonl)

      call plevs0( ps, pint, pmid, pdel, plonl )

      do k = 1,plev
         pmln(:plonl,k)  = log( pmid(:plonl,k) )
         rpdel(:plonl,k) = 1. / pdel(:plonl,k)
         piln(:plonl,k)  = log( pint(:plonl,k) )
      end do
      piln(:plonl,plevp) = log( pint(:plonl,plevp) )
      do k = 1,plev-1
         rpdeli(:plonl,k) = 1. / ( pmid(:plonl,k+1) - pmid(:plonl,k) )
      end do
      rpdeli(:plonl,plev) = 0.

      if( advqth ) then
         call pt2t( pt0, pmid, t0, plonl )
         call pt2t( ptadv, pmid, tadv, plonl )
         call virtem( tadv, shadv, tv, plonl )
      else
         call virtem( t1, sh1, tv, plonl )
      end if

!----------------------------------------------------------------
!	... form geopotential height above surface
!           at midpoints and interfaces
!----------------------------------------------------------------
      call zmid( piln, pmln, tv, zm, plonl )
      call zint( piln, pmln, tv, zm, zi, plonl )
!----------------------------------------------------------------
!	... form absolute geopotential height
!           at midpoints and interfaces
!----------------------------------------------------------------
      call zr2a( zm, plev, phis, zma, plonl )
      call zr2a( zi, plevp, phis, zia, plonl )
      call z2phi( zma, plev, phima, plonl )

      if( .not. rdomega ) then
         do k = 1,plev
!----------------------------------------------------------------
!	... temporary change; using etadot from ffsl
!----------------------------------------------------------------
            omga(:plonl,k) = etadot(:plonl,k)
         end do
      end if

!----------------------------------------------------------------
! 	... fluxes and delta to history file
!----------------------------------------------------------------
      do file = 1,moz_file_cnt
         if( fdiags(file) > 0 ) then
	    vfac = hsa_fac*latwts(base_lat+lat)*rgrav
	    do n = 1,fdiags(file)
!----------------------------------------------------------------
! 	... mass fluxes to history file
!----------------------------------------------------------------
	       do m = 6,8
	          fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1)+m-1)
	          select case( m )
	             case( 6 )
			tmp(:,:) = advn_flux(n,file)%hist_xflx(:,:,lat,ip)
	             case( 7 )
			tmp(:,:) = advn_flux(n,file)%hist_yflx(:,:,lat,ip)
	             case( 8 )
			tmp(:,:) = advn_flux(n,file)%hist_vflx(:,:,lat,ip)
	          end select
                  call outfld( fldname, tmp, plonl, ip, lat, file )
	       end do
!----------------------------------------------------------------
! 	... advection delta to history file
!----------------------------------------------------------------
	       fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1))
               tmp(:,:) = advn_flux(n,file)%mass_delta(:,:,lat,ip)
               call outfld( fldname, tmp, plonl, ip, lat, file )
!----------------------------------------------------------------
! 	... surface pressure delta to history file
!----------------------------------------------------------------
	       m = advn_flux(n,file)%index
	       advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = &
                 (as(:plonl,:plev,m)*pdel(:plonl,:plev) &
                  - advn_flux(n,file)%adv_mass(:plonl,:plev,lat,ip))*vfac
	       fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1)+1)
               tmp(:,:) = advn_flux(n,file)%mass_delta(:,:,lat,ip)
               call outfld( fldname, tmp, plonl, ip, lat, file )
            end do
         end if
      end do

!----------------------------------------------------------------
! 	... surface fluxes / dry deposition
!----------------------------------------------------------------
      if( xactive_drydep ) then
         wind_speed(:) = sqrt( u(:,plev)*u(:,plev) + v(:,plev)*v(:,plev) )
      end if
      call sflxdr( lat, ip, dtime, calday, ncdate, &
                   ncsec, ioro, as, pmid, rpdel, &
		   tv, zi, ts, sflx, snow, &
                   fsds, pint(1,plevp), wind_speed, shadv(1,plev), precip, &
                   t1(1,plev), soilw, shflx, ts_avg, fsds_avg, plonl )

!================================================================
!  	... vertical diffusion
!================================================================
      call gamdry2( lat, pcnst, as, pdel, vdiffbj, plonl )
      do file = 1,moz_file_cnt
         if( fdiags(file) > 0 ) then
	    do n = 1,fdiags(file)
	       m = advn_flux(n,file)%index
	       advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = as(:plonl,:plev,m)
	    end do
         end if
      end do
      call vdiffdr( divdiff, arvdiff, lat, ip, u, &
		    v, tadv, t1, pmid, pint, &
		    rpdel, rpdeli, dtime, zm, taux, &
		    tauy, hflx, sflx, shflx, cgs, &
		    kvh, ptadv, as, shadv, pblh, &
		    tpert, qpert, plonl )
      call gamdry2( lat, pcnst, as, pdel, vdiffaj, plonl )
      do file = 1,moz_file_cnt
         if( fdiags(file) > 0 ) then
	    do n = 1,fdiags(file)
	       m = advn_flux(n,file)%index
	       advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = &
                       (as(:plonl,:plev,m) - advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip)) &
                       *pdel(:plonl,:plev)*vfac
	       fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1)+3)
               tmp(:,:) = advn_flux(n,file)%mass_delta(:,:,lat,ip)
               call outfld( fldname, tmp, plonl, ip, lat, file )
	    end do
         end if
      end do

!----------------------------------------------------------------
! 	... update the temperature field to be consistent with the potential
!           temperature field which has just been modified by vertical diffusion.
!----------------------------------------------------------------
      if( advqth ) then
         call pt2t( ptadv, pmid, tadv, plonl )
         if( addqrad ) then
!----------------------------------------------------------------
!  	... add radiative heating tendencies to advected/diffused t
!----------------------------------------------------------------
            do k = 1,plev
               do i = 1,plonl
                  tadv(i,k) = tadv(i,k) + qrad(i,k)*dtime
               end do
            end do
            call t2pt( tadv, pmid, ptadv, plonl )
         end if
      end if

      do file = 1,moz_file_cnt
         if( fdiags(file) > 0 ) then
	    do n = 1,fdiags(file)
	       m = advn_flux(n,file)%index
	       advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = as(:plonl,:plev,m)
	    end do
         end if
      end do
      if( arconvccm ) then
!----------------------------------------------------------------
! 	... archived ccm scheme (zhang/hack)
!----------------------------------------------------------------
         call aphysar( lat, nstep, dtime, pdel, rpdel, &
                       zmu, zmd, zeu, hketa, hkbeta, &
		       as, plonl )
      else if( diconvccm ) then
!----------------------------------------------------------------
! 	... diagnosed ccm scheme (zhang/hack)
!----------------------------------------------------------------
         call aphys( lat, ip, nstep, dtime, pmid, &
	             pint, pdel, rpdel, zm, phima, &
		     zi, tpert, qpert, phis, pblh, &
		     tadv, shadv, as, cmfmc, zdu, &
		     cmfdqr, precc, cldtop, cldbot, dlf, plonl )
!----------------------------------------------------------------
!  	... update theta (advected q is updated in aphys)
!----------------------------------------------------------------
         call t2pt( tadv, pmid, ptadv, plonl )
      end if
      precip(:) = precc(:)
      call gamdry2( lat, pcnst, as, pdel, convj, plonl )
      do file = 1,moz_file_cnt
         if( fdiags(file) > 0 ) then
	    do n = 1,fdiags(file)
	       m = advn_flux(n,file)%index
	       advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = &
                  (as(:plonl,:plev,m) - advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip)) &
                   *pdel(:plonl,:plev)*vfac
	       fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1)+2)
               tmp(:,:) = advn_flux(n,file)%mass_delta(:,:,lat,ip)
               call outfld( fldname, tmp, plonl, ip, lat, file )
	    end do
         end if
      end do

!================================================================
!   	... cloud physics
!================================================================
      if( dicldphys ) then
         call clouddr( ip, lat, dtime, oro, pint, &
		       pmid, pdel, rpdeli, ts, tadv, &
		       shadv, omga, cmfmc, zdu, cldtop, &
		       cldbot, dlf, t0, sh0, cmfdqr, &
                       as(1,1,pcnst), cldfr, concld, cldtot, cme, &
		       nevapr, nrain, precl, plonl, platl )
!----------------------------------------------------------------
!   	... update theta
!----------------------------------------------------------------
         call t2pt( tadv, pmid, ptadv, plonl )
         precip(:) = precip(:) + precl(:)
      else
!----------------------------------------------------------------
!   	... special code to set the cloud top
!----------------------------------------------------------------
	 do i = 1,plonl
	    cldtop(i) = real( plev )
            cldbot(i) = real( 1 )
	    do k = 1,plev
	       if( hketa(i,k) /= 0. .and. zma(i,k) > 0. ) then
		  cldtop(i) = real( k )
		  exit
	       end if
	    end do
            do k = plev,1,-1
               if( hketa(i,k) /= 0. .and. zma(i,k) > 0. ) then
                  cldbot(i) = real( k )
                  exit
               end if
            end do
	 end do
      end if
!----------------------------------------------------------------
!   	... cloud lightning no production
!----------------------------------------------------------------
      if( do_lghtning ) then
         call lght_no_prod( ncdate, ncsec, calday, cldtop, &
			    cldbot, zm, zi, t1, plonl, &
                            lat, ip )
      end if
      call gamdry2( lat, pcnst, as, pdel, cloudj, plonl )
!----------------------------------------------------------------
!   	... chemistry
!----------------------------------------------------------------
      if( tropchem ) then
	 do file = 1,moz_file_cnt
            if( fdiags(file) > 0 ) then
	       do n = 1,fdiags(file)
	          m = advn_flux(n,file)%index
	          advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = as(:plonl,:plev,m)
	       end do
            end if
         end do
	 if( dicldphys ) then
	    if( .not. xactive_h2o ) then
               call chemdr( as, pint, nstep, calday, ncdate, &
                            ncsec, lat, ip, dtime, ps, &
                            pmid, pdel, coro, ts, zma, &
                            zi, zia, phis, cldfr, cldtop, cmfdqr, &
                            nrain, nevapr, as(1,1,pcnst), t1, u, &
                            v, sh1, snow, zm, plonl, platl )
	    else
               call chemdr( as, pint, nstep, calday, ncdate, &
			    ncsec, lat, ip, dtime, ps, &
			    pmid, pdel, coro, ts, zma, &
			    zi, zia, phis, cldfr, cldtop, cmfdqr, &
			    nrain, nevapr, as(1,1,pcnst), t1, u, &
                            v, shadv, snow, zm, plonl, platl )
	    end if
	 else
	    if( .not. xactive_h2o ) then
               call chemdr( as, pint, nstep, calday, ncdate, &
			    ncsec, lat, ip, dtime, ps, &
			    pmid, pdel, coro, ts, zma, &
			    zi, zia, phis, cldfr, cldtop, cmfdqr, &
			    nrain, nevapr, cwat, t1, u, &
                            v, sh1, snow, zm, plonl, platl )
	    else
               call chemdr( as, pint, nstep, calday, ncdate, &
			    ncsec, lat, ip, dtime, ps, &
			    pmid, pdel, coro, ts, zma, &
			    zi, zia, phis, cldfr, cldtop, cmfdqr, &
			    nrain, nevapr, cwat, t1, u, &
                            v, shadv, snow, zm, plonl, platl )
	    end if
	 end if
	 do file = 1,moz_file_cnt
            if( fdiags(file) > 0 ) then
	       do n = 1,fdiags(file)
	          m = advn_flux(n,file)%index
	          advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = &
                       (as(:plonl,:plev,m) - advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip)) &
                        *pdel(:plonl,:plev)*vfac
	          fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(16,2)+8*(n-1)+4)
                  tmp(:,:) = advn_flux(n,file)%mass_delta(:,:,lat,ip)
                  call outfld( fldname, tmp, plonl, ip, lat, file )
	       end do
            end if
         end do
      else
	 do file = 1,moz_file_cnt
            if( fdiags(file) > 0 ) then
	       do n = 1,fdiags(file)
	          advn_flux(n,file)%mass_delta(:plonl,:plev,lat,ip) = 0.
	       end do
            end if
	 end do
      end if
!----------------------------------------------------------------
!     	... output fields corresponding to time index nstep
!----------------------------------------------------------------
!----------------------------------------------------------------
!     	... cloud height
!----------------------------------------------------------------
      do i = 1,plonl
	 k = max( 1,nint( cldtop(i) ) )
	 wrk(i) = m2km * zm(i,k)
      end do
      do file = 1,sim_file_cnt
         call outfld( 'PCTOP', wrk,  plonl,  ip, lat, file )
      end do
      do i = 1,plonl
         k = min( max( 1,nint( cldbot(i) ) ), plev )
         wrk(i) = m2km * zm(i,k)
      end do
      do file = 1,sim_file_cnt
         call outfld( 'PCBOT', wrk,  plonl,  ip, lat, file )
         call outfld( 'CLOUD', cldfr(1,2),  plonl,  ip, lat, file )
         call outfld( 'CLDTOT', cldtot,  plonl,  ip, lat, file )
         call outfld( 'PS', ps, plonl, ip, lat, file )
         call outfld( 'PBLH', pblh, plonl, ip, lat, file )
         call outfld( 'PRECL', precl, plonl, ip, lat, file )
         call outfld( 'PRECC', precc, plonl, ip, lat, file )
         call outfld( 'PRECT', precip, plonl, ip, lat, file )
         call outfld( 'U', u, plonl, ip, lat, file )
         call outfld( 'V', v, plonl, ip, lat, file )
         if( xactive_h2o ) then
            call outfld( 'Q', shadv,  plonl, ip, lat, file )
         else
            call outfld( 'Q', sh1,  plonl, ip, lat, file )
         end if
         call outfld( 'T', t1,   plonl, ip, lat, file )
         call outfld( 'ORO', oro,  plonl, ip, lat, file )
         call outfld( 'SNOWH', snow,  plonl, ip, lat, file )
         call outfld( 'ETADOT', etadot(1,2), plonl, ip, lat, file )
         call outfld( 'OMEGA', omga, plonl, ip, lat, file )
         call outfld( 'CMFDQR', cmfdqr, plonl, ip, lat, file )
         if( local_time(file) > 0 ) then
            fldname = 'PS_' // local_time_char(file) // '_LT'
            call outfld( fldname, ps, plonl, ip, lat, file )
            fldname = 'T_' // local_time_char(file) // '_LT'
            call outfld( fldname, t1, plonl, ip, lat, file )
         end if
      end do

!--------------------------------------------------------
!	... the transported species are now output in
!           the chemmain subroutine except for cloud water
!--------------------------------------------------------
      if( .not. chemistry .and. pcnst > 1 ) then
	 do file = 1,moz_file_cnt
	    if( hfile(file)%wrhstts ) then
               do m = 1,hfile(file)%histout_cnt(1,1)
	          fldname = hfile(file)%hist_inst(hfile(file)%histout_ind(1,1)+m-1)
	          hndx    = hfile(file)%inst_map(hfile(file)%histout_ind(1,1)+m-1)
                  call outfld( fldname, as(1,1,hndx), plonl, ip, lat, file )
               end do
	    end if
            do m = 1,hfile(file)%histout_cnt(1,2)
	       fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(1,2)+m-1)
	       hndx    = hfile(file)%timav_map(hfile(file)%histout_ind(1,2)+m-1)
               call outfld( fldname, as(1,1,hndx), plonl, ip, lat, file )
            end do
         end do
      end if
      do file = 1,sim_file_cnt
#ifdef DI_CLOUD_PHYS
         call outfld( tracnam(pcnst), as(1,1,pcnst), plonl, ip, lat, file )
#else
         call outfld( tracnam(pcnst), cwat, plonl, ip, lat, file )
#endif
      end do
      if( advqth ) then
         do file = 1,sim_file_cnt
            call outfld( 'XQ', shadv, plonl, ip, lat, file )
            call outfld( 'XTHETA', ptadv, plonl, ip, lat, file )
         end do
      end if

!--------------------------------------------------------
! 	... free dynamic memory
!--------------------------------------------------------
!     if( chemistry .and. pnats /= 0 ) then
! deallocate( nas )
!     end if
#ifdef DEBUG
      write(*,'(''physlic: leaving for j,ip,node = '',3i4)') lat,ip,thisnode
#endif

      end subroutine physlic

      end module mo_physlic