pkg/fizhi/fizhi_step_diag.F

C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.6 2004/07/26 18:45:17 molod Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
      subroutine fizhi_step_diag(myThid,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

#ifdef ALLOW_DIAGNOSTICS
#include "SIZE.h"
#include "diagnostics_SIZE.h"
#include "diagnostics.h"
#endif

      integer myThid,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),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)

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

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

c Incident Solar Radiation (W/m**2)
c ---------------------------------
      if (iradswt.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iradswt,bi,bj)= qdiag(i,j,iradswt,bi,bj) +
     .                                                 radswt(i,j,bi,bj)
      enddo
      enddo
      endif
                                                                                
c Net Solar Radiation at the Ground (W/m**2)
c ------------------------------------------
      if (iradswg.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iradswg,bi,bj) = qdiag(i,j,iradswg,bi,bj) + 
     .                               radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
      enddo
      enddo
      endif
                                                                                
c Net Clear Sky Solar Radiation at the Ground (W/m**2)
c ----------------------------------------------------
      if (iswgclr.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iswgclr,bi,bj) = qdiag(i,j,iswgclr,bi,bj) + 
     .                               swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
      enddo
      enddo
      endif
                                                                                
c Outgoing Solar Radiation at top (W/m**2)
c -----------------------------------------
      if (iosr.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iosr,bi,bj) = qdiag(i,j,iosr,bi,bj) + 
     .                            (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
      enddo
      enddo
      endif
                                                                                
c Outgoing Clear Sky Solar Radiation at top (W/m**2)
c ---------------------------------------------------
      if (iosrclr.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iosrclr,bi,bj) = qdiag(i,j,iosrclr,bi,bj) + 
     .                         (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)         
      enddo
      enddo
      endif
                                                                                
c Upward Longwave Flux at the Ground (W/m**2)
c -------------------------------------------
      if (ilwgup.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,ilwgup,bi,bj) = qdiag(i,j,ilwgup,bi,bj) + st4(i,j,bi,bj) 
     .                 + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
      enddo
      enddo
      endif
                                                                                
c Net Longwave Flux at the Ground (W/m**2)
c ----------------------------------------
      if (iradlwg.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iradlwg,bi,bj) =  qdiag(i,j,iradlwg,bi,bj) + 
     .                                              radlwg(i,j,bi,bj) + 
     .                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
      enddo
      enddo
      endif
                                                                                
c Net Longwave Flux at the Ground Clear Sky (W/m**2)
c --------------------------------------------------
      if (ilwgclr.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,ilwgclr,bi,bj) = qdiag(i,j,ilwgclr,bi,bj) + 
     .                                               lwgclr(i,j,bi,bj) + 
     .                   dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
      enddo
      enddo
      endif
                                                                                
      nradswt = nradswt + 1
      nradswg = nradswg + 1
      nswgclr = nswgclr + 1
      nosr    = nosr    + 1
      nosrclr = nosrclr + 1
      nradlwg = nradlwg + 1
      nlwgclr = nlwgclr + 1
      nlwgup  = nlwgup  + 1
                                                                                
C **********************************************************************        
      do L=1,Nrphys

c Total Diabatic U-Tendency (m/sec/day)
c -------------------------------------
      if( idiabu.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,idiabu+L-1,bi,bj) = qdiag(i,j,idiabu+L-1,bi,bj)
     .            + ( moistu (i,j,L,bi,bj) + turbu(i,j,L,bi,bj) )*86400
      enddo
      enddo
      endif
                                                                  
c Total Diabatic V-Tendency (m/sec/day)
c -------------------------------------
      if( idiabv.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,idiabv+L-1,bi,bj) = qdiag(i,j,idiabv+L-1,bi,bj)
     .            + ( moistv (i,j,L,bi,bj) + turbv(i,j,L,bi,bj) )*86400
      enddo
      enddo
      endif

c Total Diabatic T-Tendency (deg/day)
c -----------------------------------
      if( idiabt.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,idiabt+L-1,bi,bj) = qdiag(i,j,idiabt+L-1,bi,bj) +
     .   ( 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
      endif
                                                                  
c Total Diabatic Q-Tendency (g/kg/day)
c ------------------------------------
      if( idiabq.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,idiabq+L-1,bi,bj) =   qdiag(i,j,idiabq+L-1,bi,bj) +
     . ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) * 
     .                                      pinv(i,j)*86400*1000
      enddo
      enddo
      endif
    
c Longwave Heating (deg/day)
c --------------------------
      if (iradlw.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iradlw+l-1,bi,bj) = qdiag(i,j,iradlw+l-1,bi,bj) + 
     . ( 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
      endif
                                                                                
c Longwave Heating Clear-Sky (deg/day)
c ------------------------------------
      if (ilwclr.ne.0) then                                                     
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,ilwclr+l-1,bi,bj) = qdiag(i,j,ilwclr+l-1,bi,bj) +
     . ( 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
      endif
                                                                                
c Solar Radiative Heating (deg/day)
c ---------------------------------
      if (iradsw.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iradsw+l-1,bi,bj) = qdiag(i,j,iradsw+l-1,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
      endif
                                                                                
c Clear Sky Solar Radiative Heating (deg/day)
c -------------------------------------------
      if (iswclr.ne.0) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iswclr+l-1,bi,bj) = qdiag(i,j,iswclr+l-1,bi,bj) +
     .           swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
     .                           pk(i,j,l,bi,bj)*pinv(i,j)*86400
      enddo
      enddo
      endif
                                                                                
c Averaged U-Field (m/sec)
c ------------------------
      if( iuwnd.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) + 
     .                                                 uphy(i,j,L,bi,bj)
      enddo
      enddo
      endif

c Averaged V-Field (m/sec)
c ------------------------
      if( ivwnd.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) + 
     .                                                 vphy(i,j,L,bi,bj)
      enddo
      enddo
      endif

c Averaged T-Field (deg)
c ----------------------
      if( itmpu.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,itmpu+L-1,bi,bj) = qdiag(i,j,itmpu+L-1,bi,bj) + 
     .                               thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
      enddo
      enddo
      endif
                                                                                
c Averaged QQ-Field (m/sec)**2
c ----------------------------
      if( itke.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
      enddo
      enddo
      endif
                                                                                
c Averaged Q-Field (g/kg)
c -----------------------
      if( isphu.ne.0 ) then
      do j=jm1,jm2
      do i=im1,im2
      qdiag(i,j,isphu+L-1,bi,bj) = qdiag(i,j,isphu+L-1,bi,bj) + 
     .                                            sphy(i,j,L,bi,bj)*1000
      enddo
      enddo
      endif
                                                                                
      enddo

      ndiabu = ndiabu + 1
      ndiabv = ndiabv + 1
      ndiabt = ndiabt + 1
      ndiabq = ndiabq + 1
      nradlw = nradlw + 1                                                     
      nlwclr = nlwclr + 1                                                     
      nradsw = nradsw + 1                                                     
      nswclr = nswclr + 1                                                     
      nuwnd  = nuwnd  + 1
      nvwnd  = nvwnd  + 1                           
      ntmpu  = ntmpu  + 1                           
      ntke   = ntke   + 1
      nsphu  = nsphu  + 1

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

c Vertically Averaged Moist-T Increment (K/day)
c ---------------------------------------------
      if( ivdtmoist.ne.0 ) 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
      qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) + 
     .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
      enddo
      enddo
      endif

c Vertically Averaged Turb-T Increment (K/day)
c --------------------------------------------
      if( ivdtturb.ne.0 ) 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
      qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) + 
     .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
      enddo
      enddo
      endif

c Vertically Averaged RADLW Temperature Increment (K/day)
c -------------------------------------------------------
      if( ivdtradlw.ne.0 ) 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
      qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) + 
     .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
      enddo
      enddo
      endif

c Vertically Averaged RADSW Temperature Increment (K/day)
c -------------------------------------------------------
      if( ivdtradsw.ne.0 ) 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
      qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) + 
     . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
      enddo
      enddo
      endif

      nvdtmoist = nvdtmoist + 1
      nvdtturb  = nvdtturb  + 1
      nvdtradlw = nvdtradlw + 1
      nvdtradsw = nvdtradsw + 1
 
      return
      end
1	molod	1.7	C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.6 2004/07/26 18:45:17 molod Exp $
2	molod	1.1	C $Name: $
3	molod	1.4
4	molod	1.6	#include "FIZHI_OPTIONS.h"
5	molod	1.2	subroutine fizhi_step_diag(myThid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
6	molod	1.1	. radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
7			. turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
8	molod	1.5	. lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
9			. im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
10	molod	1.1	C***********************************************************************
11			implicit none
12
13	molod	1.3	#ifdef ALLOW_DIAGNOSTICS
14	molod	1.4	#include "SIZE.h"
15			#include "diagnostics_SIZE.h"
16	molod	1.1	#include "diagnostics.h"
17	molod	1.3	#endif
18	molod	1.1
19	molod	1.5	integer myThid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
20	molod	1.6	_RL p(im2,jm2,Nbi,Nbj)
21	molod	1.7	_RL uphy(im2,jm2,Nrphys,Nbi,Nbj)
22			_RL vphy(im2,jm2,Nrphys,Nbi,Nbj)
23			_RL thphy(im2,jm2,Nrphys,Nbi,Nbj)
24			_RL sphy(im2,jm2,Nrphys,Nbi,Nbj)
25	molod	1.6	_RL qq(im2,jm2,Nrphys),pk(im2,jm2,Nrphys,Nbi,Nbj)
26			_RL dp(im2,jm2,Nrphys,Nbi,Nbj)
27			_RL radswt(im2,jm2,Nbi,Nbj),radswg(im2,jm2,Nbi,Nbj)
28			_RL swgclr(im2,jm2,Nbi,Nbj),osr(im2,jm2,Nbi,Nbj)
29			_RL osrclr(im2,jm2,Nbi,Nbj),st4(im2,jm2,Nbi,Nbj)
30			_RL dst4(im2,jm2,Nbi,Nbj),tgz(im2,jm2,Nbi,Nbj)
31			_RL tg0(im2,jm2,Nbi,Nbj),radlwg(im2,jm2,Nbi,Nbj)
32			_RL lwgclr(im2,jm2,Nbi,Nbj)
33	molod	1.7	_RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
34			_RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
35	molod	1.6	_RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
36			_RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
37	molod	1.7	_RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
38			_RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
39	molod	1.6	_RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
40			_RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
41	molod	1.7	_RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
42			_RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
43	molod	1.6	_RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
44			_RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
45			_RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)
46	molod	1.1
47	molod	1.2	integer i,j,L
48	molod	1.6	_RL pinv(im2,jm2), qbar(im2,jm2)
49	molod	1.1
50			C **********************************************************************
51
52			do j=jm1,jm2
53			do i=im1,im2
54	molod	1.2	pinv(i,j) = 1.0 / p(i,j,bi,bj)
55	molod	1.1	enddo
56			enddo
57	molod	1.2
58	molod	1.1	c Incident Solar Radiation (W/m**2)
59			c ---------------------------------
60			if (iradswt.ne.0) then
61			do j=jm1,jm2
62			do i=im1,im2
63	molod	1.2	qdiag(i,j,iradswt,bi,bj)= qdiag(i,j,iradswt,bi,bj) +
64			. radswt(i,j,bi,bj)
65			enddo
66			enddo
67			endif
68	molod	1.1
69			c Net Solar Radiation at the Ground (W/m**2)
70			c ------------------------------------------
71	molod	1.2	if (iradswg.ne.0) then
72	molod	1.1	do j=jm1,jm2
73			do i=im1,im2
74	molod	1.2	qdiag(i,j,iradswg,bi,bj) = qdiag(i,j,iradswg,bi,bj) +
75			. radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
76			enddo
77			enddo
78			endif
79	molod	1.1
80			c Net Clear Sky Solar Radiation at the Ground (W/m**2)
81			c ----------------------------------------------------
82	molod	1.2	if (iswgclr.ne.0) then
83	molod	1.1	do j=jm1,jm2
84			do i=im1,im2
85	molod	1.2	qdiag(i,j,iswgclr,bi,bj) = qdiag(i,j,iswgclr,bi,bj) +
86			. swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
87			enddo
88			enddo
89			endif
90	molod	1.1
91	molod	1.2	c Outgoing Solar Radiation at top (W/m**2)
92	molod	1.1	c -----------------------------------------
93	molod	1.2	if (iosr.ne.0) then
94	molod	1.1	do j=jm1,jm2
95			do i=im1,im2
96	molod	1.2	qdiag(i,j,iosr,bi,bj) = qdiag(i,j,iosr,bi,bj) +
97			. (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
98			enddo
99			enddo
100			endif
101	molod	1.1
102	molod	1.2	c Outgoing Clear Sky Solar Radiation at top (W/m**2)
103	molod	1.1	c ---------------------------------------------------
104	molod	1.2	if (iosrclr.ne.0) then
105	molod	1.1	do j=jm1,jm2
106			do i=im1,im2
107	molod	1.2	qdiag(i,j,iosrclr,bi,bj) = qdiag(i,j,iosrclr,bi,bj) +
108			. (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)
109			enddo
110			enddo
111			endif
112	molod	1.1
113			c Upward Longwave Flux at the Ground (W/m**2)
114			c -------------------------------------------
115	molod	1.2	if (ilwgup.ne.0) then
116	molod	1.1	do j=jm1,jm2
117			do i=im1,im2
118	molod	1.2	qdiag(i,j,ilwgup,bi,bj) = qdiag(i,j,ilwgup,bi,bj) + st4(i,j,bi,bj)
119			. + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
120			enddo
121			enddo
122			endif
123	molod	1.1
124			c Net Longwave Flux at the Ground (W/m**2)
125			c ----------------------------------------
126	molod	1.2	if (iradlwg.ne.0) then
127	molod	1.1	do j=jm1,jm2
128			do i=im1,im2
129	molod	1.2	qdiag(i,j,iradlwg,bi,bj) = qdiag(i,j,iradlwg,bi,bj) +
130			. radlwg(i,j,bi,bj) +
131			. dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
132			enddo
133			enddo
134			endif
135	molod	1.1
136			c Net Longwave Flux at the Ground Clear Sky (W/m**2)
137			c --------------------------------------------------
138			if (ilwgclr.ne.0) then
139			do j=jm1,jm2
140			do i=im1,im2
141	molod	1.2	qdiag(i,j,ilwgclr,bi,bj) = qdiag(i,j,ilwgclr,bi,bj) +
142			. lwgclr(i,j,bi,bj) +
143			. dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
144	molod	1.1	enddo
145			enddo
146			endif
147	molod	1.2
148			nradswt = nradswt + 1
149			nradswg = nradswg + 1
150			nswgclr = nswgclr + 1
151	molod	1.1	nosr = nosr + 1
152			nosrclr = nosrclr + 1
153			nradlwg = nradlwg + 1
154			nlwgclr = nlwgclr + 1
155			nlwgup = nlwgup + 1
156
157			C **********************************************************************
158			do L=1,Nrphys
159
160			c Total Diabatic U-Tendency (m/sec/day)
161			c -------------------------------------
162			if( idiabu.ne.0 ) then
163			do j=jm1,jm2
164			do i=im1,im2
165	molod	1.2	qdiag(i,j,idiabu+L-1,bi,bj) = qdiag(i,j,idiabu+L-1,bi,bj)
166			. + ( moistu (i,j,L,bi,bj) + turbu(i,j,L,bi,bj) )*86400
167	molod	1.1	enddo
168			enddo
169			endif
170
171			c Total Diabatic V-Tendency (m/sec/day)
172			c -------------------------------------
173			if( idiabv.ne.0 ) then
174			do j=jm1,jm2
175			do i=im1,im2
176	molod	1.2	qdiag(i,j,idiabv+L-1,bi,bj) = qdiag(i,j,idiabv+L-1,bi,bj)
177			. + ( moistv (i,j,L,bi,bj) + turbv(i,j,L,bi,bj) )*86400
178	molod	1.1	enddo
179			enddo
180			endif
181
182			c Total Diabatic T-Tendency (deg/day)
183			c -----------------------------------
184			if( idiabt.ne.0 ) then
185			do j=jm1,jm2
186			do i=im1,im2
187	molod	1.2	qdiag(i,j,idiabt+L-1,bi,bj) = qdiag(i,j,idiabt+L-1,bi,bj) +
188			. ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +
189			. lwdt(i,j,L,bi,bj) +
190			. dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +
191			. swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )
192			. * pk(i,j,L,bi,bj)pinv(i,j)86400
193	molod	1.1	enddo
194			enddo
195			endif
196
197			c Total Diabatic Q-Tendency (g/kg/day)
198			c ------------------------------------
199			if( idiabq.ne.0 ) then
200			do j=jm1,jm2
201			do i=im1,im2
202	molod	1.2	qdiag(i,j,idiabq+L-1,bi,bj) = qdiag(i,j,idiabq+L-1,bi,bj) +
203			. ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) *
204			. pinv(i,j)864001000
205	molod	1.1	enddo
206			enddo
207			endif
208
209	molod	1.2	c Longwave Heating (deg/day)
210			c --------------------------
211			if (iradlw.ne.0) then
212	molod	1.1	do j=jm1,jm2
213			do i=im1,im2
214	molod	1.2	qdiag(i,j,iradlw+l-1,bi,bj) = qdiag(i,j,iradlw+l-1,bi,bj) +
215			. ( lwdt(i,j,l,bi,bj) +
216			. dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
217			. * pk(i,j,l,bi,bj)pinv(i,j)86400
218	molod	1.1	enddo
219			enddo
220			endif
221
222			c Longwave Heating Clear-Sky (deg/day)
223			c ------------------------------------
224			if (ilwclr.ne.0) then
225			do j=jm1,jm2
226			do i=im1,im2
227	molod	1.2	qdiag(i,j,ilwclr+l-1,bi,bj) = qdiag(i,j,ilwclr+l-1,bi,bj) +
228			. ( lwdtclr(i,j,l,bi,bj) +
229			. dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
230			. * pk(i,j,l,bi,bj)pinv(i,j)86400
231			enddo
232			enddo
233			endif
234	molod	1.1
235			c Solar Radiative Heating (deg/day)
236			c ---------------------------------
237	molod	1.2	if (iradsw.ne.0) then
238	molod	1.1	do j=jm1,jm2
239			do i=im1,im2
240	molod	1.2	qdiag(i,j,iradsw+l-1,bi,bj) = qdiag(i,j,iradsw+l-1,bi,bj) +
241			. + swdt(i,j,l,bi,bj)radswt(i,j,bi,bj)
242			. pk(i,j,l,bi,bj)pinv(i,j)86400
243			enddo
244			enddo
245			endif
246	molod	1.1
247			c Clear Sky Solar Radiative Heating (deg/day)
248			c -------------------------------------------
249	molod	1.2	if (iswclr.ne.0) then
250	molod	1.1	do j=jm1,jm2
251			do i=im1,im2
252	molod	1.2	qdiag(i,j,iswclr+l-1,bi,bj) = qdiag(i,j,iswclr+l-1,bi,bj) +
253			. swdtclr(i,j,l,bi,bj)radswt(i,j,bi,bj)
254	molod	1.5	. pk(i,j,l,bi,bj)pinv(i,j)86400
255	molod	1.1	enddo
256			enddo
257			endif
258	molod	1.2
259	molod	1.1	c Averaged U-Field (m/sec)
260			c ------------------------
261			if( iuwnd.ne.0 ) then
262			do j=jm1,jm2
263			do i=im1,im2
264	molod	1.2	qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
265			. uphy(i,j,L,bi,bj)
266	molod	1.1	enddo
267			enddo
268			endif
269
270			c Averaged V-Field (m/sec)
271			c ------------------------
272			if( ivwnd.ne.0 ) then
273			do j=jm1,jm2
274			do i=im1,im2
275	molod	1.2	qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
276			. vphy(i,j,L,bi,bj)
277	molod	1.1	enddo
278			enddo
279			endif
280
281			c Averaged T-Field (deg)
282			c ----------------------
283			if( itmpu.ne.0 ) then
284			do j=jm1,jm2
285			do i=im1,im2
286	molod	1.2	qdiag(i,j,itmpu+L-1,bi,bj) = qdiag(i,j,itmpu+L-1,bi,bj) +
287			. thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
288	molod	1.1	enddo
289			enddo
290			endif
291
292			c Averaged QQ-Field (m/sec)**2
293			c ----------------------------
294			if( itke.ne.0 ) then
295			do j=jm1,jm2
296			do i=im1,im2
297	molod	1.2	qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
298	molod	1.1	enddo
299			enddo
300			endif
301
302			c Averaged Q-Field (g/kg)
303			c -----------------------
304			if( isphu.ne.0 ) then
305			do j=jm1,jm2
306			do i=im1,im2
307	molod	1.2	qdiag(i,j,isphu+L-1,bi,bj) = qdiag(i,j,isphu+L-1,bi,bj) +
308			. sphy(i,j,L,bi,bj)*1000
309	molod	1.1	enddo
310			enddo
311			endif
312
313	molod	1.2	enddo
314	molod	1.1
315			ndiabu = ndiabu + 1
316			ndiabv = ndiabv + 1
317			ndiabt = ndiabt + 1
318	molod	1.2	ndiabq = ndiabq + 1
319	molod	1.1	nradlw = nradlw + 1
320			nlwclr = nlwclr + 1
321			nradsw = nradsw + 1
322			nswclr = nswclr + 1
323			nuwnd = nuwnd + 1
324			nvwnd = nvwnd + 1
325			ntmpu = ntmpu + 1
326			ntke = ntke + 1
327			nsphu = nsphu + 1
328
329			C **********************************************************************
330
331			c Vertically Averaged Moist-T Increment (K/day)
332			c ---------------------------------------------
333			if( ivdtmoist.ne.0 ) then
334			do j=jm1,jm2
335			do i=im1,im2
336			qbar(i,j) = 0.0
337			enddo
338			enddo
339			do L=1,Nrphys
340			do j=jm1,jm2
341			do i=im1,im2
342	molod	1.2	qbar(i,j) = qbar(i,j) +
343			. moistt(i,j,L,bi,bj)pk(i,j,l,bi,bj)dp(i,j,L,bi,bj)
344	molod	1.1	enddo
345			enddo
346			enddo
347			do j=jm1,jm2
348			do i=im1,im2
349	molod	1.2	qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
350	molod	1.5	. qbar(i,j)pinv(i,j)pinv(i,j)*86400
351	molod	1.1	enddo
352			enddo
353			endif
354
355			c Vertically Averaged Turb-T Increment (K/day)
356			c --------------------------------------------
357			if( ivdtturb.ne.0 ) then
358			do j=jm1,jm2
359			do i=im1,im2
360			qbar(i,j) = 0.0
361			enddo
362			enddo
363			do L=1,Nrphys
364			do j=jm1,jm2
365			do i=im1,im2
366	molod	1.2	qbar(i,j) = qbar(i,j) +
367			. turbt(i,j,L,bi,bj)pk(i,j,l,bi,bj)dp(i,j,L,bi,bj)
368	molod	1.1	enddo
369			enddo
370			enddo
371			do j=jm1,jm2
372			do i=im1,im2
373	molod	1.2	qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
374	molod	1.5	. qbar(i,j)pinv(i,j)pinv(i,j)*86400
375	molod	1.1	enddo
376			enddo
377			endif
378
379			c Vertically Averaged RADLW Temperature Increment (K/day)
380			c -------------------------------------------------------
381			if( ivdtradlw.ne.0 ) then
382			do j=jm1,jm2
383			do i=im1,im2
384			qbar(i,j) = 0.0
385			enddo
386			enddo
387			do L=1,Nrphys
388			do j=jm1,jm2
389			do i=im1,im2
390	molod	1.2	qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
391			. dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
392			. pk(i,j,l,bi,bj)dp(i,j,L,bi,bj)
393	molod	1.1	enddo
394			enddo
395			enddo
396			do j=jm1,jm2
397			do i=im1,im2
398	molod	1.2	qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
399	molod	1.5	. qbar(i,j)pinv(i,j)pinv(i,j)*86400
400	molod	1.1	enddo
401			enddo
402			endif
403
404			c Vertically Averaged RADSW Temperature Increment (K/day)
405			c -------------------------------------------------------
406			if( ivdtradsw.ne.0 ) then
407			do j=jm1,jm2
408			do i=im1,im2
409			qbar(i,j) = 0.0
410			enddo
411			enddo
412			do L=1,Nrphys
413			do j=jm1,jm2
414			do i=im1,im2
415	molod	1.2	qbar(i,j) = qbar(i,j) +
416			. swdt(i,j,L,bi,bj)pk(i,j,l,bi,bj)dp(i,j,L,bi,bj)
417	molod	1.1	enddo
418			enddo
419			enddo
420			do j=jm1,jm2
421			do i=im1,im2
422	molod	1.2	qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
423	molod	1.5	. qbar(i,j)radswt(i,j,bi,bj)pinv(i,j)pinv(i,j)86400
424	molod	1.1	enddo
425			enddo
426			endif
427
428			nvdtmoist = nvdtmoist + 1
429			nvdtturb = nvdtturb + 1
430			nvdtradlw = nvdtradlw + 1
431			nvdtradsw = nvdtradsw + 1
432
433			return
434			end