C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.16 2005/05/24 21:03:08 molod Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
      subroutine fizhi_step_diag(myid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
     .  radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
     .  turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
     .  lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
     .  im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
C***********************************************************************
      implicit none

      integer myid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
      _RL p(im2,jm2,Nbi,Nbj)
      _RL uphy(im2,jm2,Nrphys,Nbi,Nbj)
      _RL vphy(im2,jm2,Nrphys,Nbi,Nbj)
      _RL thphy(im2,jm2,Nrphys,Nbi,Nbj)
      _RL sphy(im2,jm2,Nrphys,Nbi,Nbj)
      _RL qq(im2,jm2,Nrphys,Nbi,Nbj),pk(im2,jm2,Nrphys,Nbi,Nbj)
      _RL dp(im2,jm2,Nrphys,Nbi,Nbj)
      _RL radswt(im2,jm2,Nbi,Nbj),radswg(im2,jm2,Nbi,Nbj)
      _RL swgclr(im2,jm2,Nbi,Nbj),osr(im2,jm2,Nbi,Nbj)
      _RL osrclr(im2,jm2,Nbi,Nbj),st4(im2,jm2,Nbi,Nbj)
      _RL dst4(im2,jm2,Nbi,Nbj),tgz(im2,jm2,Nbi,Nbj)
      _RL tg0(im2,jm2,Nbi,Nbj),radlwg(im2,jm2,Nbi,Nbj)
      _RL lwgclr(im2,jm2,Nbi,Nbj)
      _RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
      _RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
      _RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
      _RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
      _RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
      _RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
      _RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
      _RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
      _RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
      _RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
      _RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
      _RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
      _RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)

      integer  i,j,L
      _RL pinv(im2,jm2), qbar(im2,jm2),tmpdiag(im2,jm2)
#ifdef ALLOW_DIAGNOSTICS
      logical  diagnostics_is_on
      external diagnostics_is_on
#endif

C **********************************************************************        

#ifdef ALLOW_DIAGNOSTICS
      do j=jm1,jm2
      do i=im1,im2
      pinv(i,j) = 1.0 / p(i,j,bi,bj)
      enddo
      enddo

c Surface Pressure (mb)
c ---------------------------------
      if(diagnostics_is_on('PS      ',myid) ) then
       call diagnostics_fill(p(1,1,bi,bj),'PS      ',0,1,3,bi,bj,myid)
      endif

c Incident Solar Radiation (W/m**2)
c ---------------------------------
      if(diagnostics_is_on('RADSWT  ',myid) ) then
       call diagnostics_fill(radswt,'RADSWT  ',0,1,3,bi,bj,myid)
      endif

c Net Solar Radiation at the Ground (W/m**2)
c ------------------------------------------
      if(diagnostics_is_on('RADSWG  ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'RADSWG  ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Net Clear Sky Solar Radiation at the Ground (W/m**2)
c ----------------------------------------------------
      if(diagnostics_is_on('SWGCLR  ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'SWGCLR  ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Outgoing Solar Radiation at top (W/m**2)
c -----------------------------------------
      if(diagnostics_is_on('OSR     ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'OSR     ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Outgoing Clear Sky Solar Radiation at top (W/m**2)
c ---------------------------------------------------
      if(diagnostics_is_on('OSRCLR  ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'OSRCLR  ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Upward Longwave Flux at the Ground (W/m**2)
c -------------------------------------------
      if(diagnostics_is_on('LWGUP   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = st4(i,j,bi,bj) 
     .                 + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'LWGUP   ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Net Longwave Flux at the Ground (W/m**2)
c ----------------------------------------
      if(diagnostics_is_on('RADLWG  ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = radlwg(i,j,bi,bj) + 
     .                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'RADLWG  ',0,1,3,bi,bj,myid)
      endif
                                                                                
c Net Longwave Flux at the Ground Clear Sky (W/m**2)
c --------------------------------------------------
      if(diagnostics_is_on('LWGCLR  ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = lwgclr(i,j,bi,bj) + 
     .                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'LWGCLR  ',0,1,3,bi,bj,myid)
      endif
                                                                                
C **********************************************************************        
      do L=1,Nrphys

c Total Diabatic U-Tendency (m/sec/day)
c -------------------------------------
      if(diagnostics_is_on('DIABU   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = (moistu (i,j,L,bi,bj)+turbu(i,j,L,bi,bj) )*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'DIABU   ',L,1,3,bi,bj,myid)
      endif

c Total Diabatic V-Tendency (m/sec/day)
c -------------------------------------
      if(diagnostics_is_on('DIABV   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = (moistv (i,j,L,bi,bj)+turbv(i,j,L,bi,bj) )*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'DIABV   ',L,1,3,bi,bj,myid)
      endif

c Total Diabatic T-Tendency (deg/day)
c -----------------------------------
      if(diagnostics_is_on('DIABT   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     .   ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +
     .      lwdt(i,j,L,bi,bj) +
     .      dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +
     .      swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )
     .      * pk(i,j,L,bi,bj)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'DIABT   ',L,1,3,bi,bj,myid)
      endif

c Total Diabatic Q-Tendency (g/kg/day)
c ------------------------------------
      if(diagnostics_is_on('DIABQ   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     . ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) * 
     .                                      pinv(i,j)*86400*1000
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'DIABQ   ',L,1,3,bi,bj,myid)
      endif
    
c Longwave Heating (deg/day)
c --------------------------
      if(diagnostics_is_on('RADLW   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     . ( lwdt(i,j,l,bi,bj) + 
     .            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
     .                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'RADLW   ',L,1,3,bi,bj,myid)
      endif

c Longwave Heating Clear-Sky (deg/day)
c ------------------------------------
      if(diagnostics_is_on('LWCLR   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     . ( lwdtclr(i,j,l,bi,bj) + 
     .            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
     .                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'LWCLR   ',L,1,3,bi,bj,myid)
      endif
                                                                                
c Solar Radiative Heating (deg/day)
c ---------------------------------
      if(diagnostics_is_on('RADSW   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     .  + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
     .                   pk(i,j,l,bi,bj)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'RADSW   ',L,1,3,bi,bj,myid)
      endif
                                                                                
c Clear Sky Solar Radiative Heating (deg/day)
c -------------------------------------------
      if(diagnostics_is_on('SWCLR   ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = 
     .  + swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
     .                   pk(i,j,l,bi,bj)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'SWCLR   ',L,1,3,bi,bj,myid)
      endif
                                                                                
c Averaged U-Field (m/sec)
c ------------------------
      if(diagnostics_is_on('UWND    ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = uphy(i,j,L,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'UWND    ',L,1,3,bi,bj,myid)
      endif

c Averaged V-Field (m/sec)
c ------------------------
      if(diagnostics_is_on('VWND    ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = vphy(i,j,L,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'VWND    ',L,1,3,bi,bj,myid)
      endif

c Averaged T-Field (deg)
c ----------------------
      if(diagnostics_is_on('TMPU    ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'TMPU    ',L,1,3,bi,bj,myid)
      endif
                                                                                
c Averaged QQ-Field (m/sec)**2
c ----------------------------
      if(diagnostics_is_on('TKE     ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = qq(i,j,L,bi,bj)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'TKE     ',L,1,3,bi,bj,myid)
      endif
                                                                                
c Averaged Q-Field (g/kg)
c -----------------------
      if(diagnostics_is_on('SPHU    ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
        tmpdiag(i,j) = sphy(i,j,L,bi,bj) * 1000.
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'SPHU    ',L,1,3,bi,bj,myid)
      endif
                                                                                
      enddo

C **********************************************************************        

c Vertically Averaged Moist-T Increment (K/day)
c ---------------------------------------------
      if(diagnostics_is_on('VDTMOIST',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = 0.0
       enddo
       enddo
       do L=1,Nrphys
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = qbar(i,j) + 
     .             moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
       enddo
       enddo
       enddo
       do j=jm1,jm2
       do i=im1,im2
       tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'VDTMOIST',0,1,3,bi,bj,myid)
      endif

c Vertically Averaged Turb-T Increment (K/day)
c --------------------------------------------
      if(diagnostics_is_on('VDTTURB ',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = 0.0
       enddo
       enddo
       do L=1,Nrphys
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = qbar(i,j) + 
     .             turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
       enddo
       enddo
       enddo
       do j=jm1,jm2
       do i=im1,im2
       tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'VDTTURB ',0,1,3,bi,bj,myid)
      endif

c Vertically Averaged RADLW Temperature Increment (K/day)
c -------------------------------------------------------
      if(diagnostics_is_on('VDTRADLW',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = 0.0
       enddo
       enddo
       do L=1,Nrphys
       do j=jm1,jm2
       do i=im1,im2
        qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
     .  dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
     .             *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
       enddo
       enddo
       enddo
       do j=jm1,jm2
       do i=im1,im2
       tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'VDTRADLW',0,1,3,bi,bj,myid)
      endif

c Vertically Averaged RADSW Temperature Increment (K/day)
c -------------------------------------------------------
      if(diagnostics_is_on('VDTRADSW',myid) ) then
       do j=jm1,jm2
       do i=im1,im2
       qbar(i,j) = 0.0
       enddo
       enddo
       do L=1,Nrphys
       do j=jm1,jm2
       do i=im1,im2
        qbar(i,j) = qbar(i,j) + 
     .             swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
       enddo
       enddo
       enddo
       do j=jm1,jm2
       do i=im1,im2
       tmpdiag(i,j) = qbar(i,j) *
     .             radswt(i,j,bi,bj) * pinv(i,j) * pinv(i,j) * 86400
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'VDTRADSW',0,1,3,bi,bj,myid)
      endif

#endif
      return
      end