/[MITgcm]/MITgcm/pkg/fizhi/fizhi_step_diag.F
ViewVC logotype

Diff of /MITgcm/pkg/fizhi/fizhi_step_diag.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by molod, Thu Jun 24 15:06:51 2004 UTC revision 1.11 by molod, Sun Aug 29 19:43:43 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3                                                                                          
4  #include "CPP_OPTIONS.h"  #include "FIZHI_OPTIONS.h"
5        subroutine fizhi_step_diag(myThid,p,uphy,vphy,thphy,sphy,qq,        subroutine fizhi_step_diag(myThid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
6       .  radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,       .  radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
7       .  turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,       .  turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
8       .  lwdt,swdt,lwdtclr,swdtclr,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,bi,bj)       .  lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
9         .  im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
10  C***********************************************************************  C***********************************************************************
11        implicit none        implicit none
12    
13    #ifdef ALLOW_DIAGNOSTICS
14    #include "SIZE.h"
15    #include "diagnostics_SIZE.h"
16  #include "diagnostics.h"  #include "diagnostics.h"
17    #endif
18    
19        integer myThid,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,bi,bj        integer myThid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
20        real radswt(im2,jm2)        _RL p(im2,jm2,Nbi,Nbj)
21          _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          _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          _RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
34          _RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
35          _RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
36          _RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
37          _RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
38          _RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
39          _RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
40          _RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
41          _RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
42          _RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
43          _RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
44          _RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
45          _RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)
46    
47        integer  i,j,L,m        integer  i,j,L
48        real getcon        _RL pinv(im2,jm2), qbar(im2,jm2)
       real cp, pstd, tstd, akap, pkstd, thstd, grav, delp  
   
       cp = getcon('CP')  
       pstd = getcon('PSTD')  
       tstd = getcon('TSTD')  
       akap = getcon('KAPPA')  
       pkstd = pstd**akap  
       thstd = tstd/pkstd  
       grav = getcon('GRAVITY')  
49    
50  C **********************************************************************          C **********************************************************************        
 C ****                  Compute 2-D Diagnostics                     ****          
 C **********************************************************************          
                                                                                   
       do j=jm1,jm2  
       do i=im1,im2  
       pinv(i,j) = 1.0 / p(i,j)                                      
       enddo                                                              
       enddo                                                              
51    
52  c Analysis Increment of Surface Pressure (mb/day)  #ifdef ALLOW_DIAGNOSTICS
 c -----------------------------------------------  
       if( ipiau.ne.0 ) then  
53        do j=jm1,jm2        do j=jm1,jm2
54        do i=im1,im2        do i=im1,im2
55        qdiag(i,j,ipiau) = qdiag(i,j,ipiau) + tend%iau%dp(i,j)*86400        pinv(i,j) = 1.0 / p(i,j,bi,bj)
56        enddo        enddo
57        enddo        enddo
58        endif  
                                                                     
59  c Incident Solar Radiation (W/m**2)  c Incident Solar Radiation (W/m**2)
60  c ---------------------------------  c ---------------------------------
61        if (iradswt.ne.0) then        if (iradswt.ne.0) then
62        do j=jm1,jm2        do j=jm1,jm2
63        do i=im1,im2        do i=im1,im2
64        qdiag(i,j,iradswt) = qdiag(i,j,iradswt) + radswt(i,j)        qdiag(i,j,iradswt,bi,bj)= qdiag(i,j,iradswt,bi,bj) +
65        enddo                                                                           .                                                 radswt(i,j,bi,bj)
66        enddo                                                                            enddo
67        endif                                                                            enddo
68          endif
69                                                                                                                                                                    
70  c Net Solar Radiation at the Ground (W/m**2)  c Net Solar Radiation at the Ground (W/m**2)
71  c ------------------------------------------  c ------------------------------------------
72        if (iradswg.ne.0) then                                                            if (iradswg.ne.0) then
73        do j=jm1,jm2        do j=jm1,jm2
74        do i=im1,im2        do i=im1,im2
75        qdiag(i,j,iradswg) = qdiag(i,j,iradswg) + coup%sw%radswg(i,j)*radswt(i,j)                qdiag(i,j,iradswg,bi,bj) = qdiag(i,j,iradswg,bi,bj) +
76        enddo                                                                           .                               radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
77        enddo                                                                            enddo
78        endif                                                                            enddo
79          endif
80                                                                                                                                                                    
81  c Net Clear Sky Solar Radiation at the Ground (W/m**2)  c Net Clear Sky Solar Radiation at the Ground (W/m**2)
82  c ----------------------------------------------------  c ----------------------------------------------------
83        if (iswgclr.ne.0) then                                                            if (iswgclr.ne.0) then
84        do j=jm1,jm2        do j=jm1,jm2
85        do i=im1,im2        do i=im1,im2
86        qdiag(i,j,iswgclr) = qdiag(i,j,iswgclr) + coup%sw%swgclr(i,j)*radswt(i,j)                qdiag(i,j,iswgclr,bi,bj) = qdiag(i,j,iswgclr,bi,bj) +
87        enddo                                                                           .                               swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
88        enddo                                                                            enddo
89        endif                                                                            enddo
90          endif
91                                                                                                                                                                    
92  c Outgoing Solar Radiation at Ptop (W/m**2)  c Outgoing Solar Radiation at top (W/m**2)
93  c -----------------------------------------  c -----------------------------------------
94        if (iosr.ne.0) then                                                            if (iosr.ne.0) then
95        do j=jm1,jm2        do j=jm1,jm2
96        do i=im1,im2        do i=im1,im2
97        qdiag(i,j,iosr) = qdiag(i,j,iosr) + (1.0-coup%sw%osr(i,j))*radswt(i,j)                qdiag(i,j,iosr,bi,bj) = qdiag(i,j,iosr,bi,bj) +
98        enddo                                                                           .                            (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
99        enddo                                                                            enddo
100        endif                                                                            enddo
101          endif
102                                                                                                                                                                    
103  c Outgoing Clear Sky Solar Radiation at Ptop (W/m**2)  c Outgoing Clear Sky Solar Radiation at top (W/m**2)
104  c ---------------------------------------------------  c ---------------------------------------------------
105        if (iosrclr.ne.0) then                                                            if (iosrclr.ne.0) then
106        do j=jm1,jm2        do j=jm1,jm2
107        do i=im1,im2        do i=im1,im2
108        qdiag(i,j,iosrclr) = qdiag(i,j,iosrclr) + (1.0-coup%sw%osrclr(i,j))*radswt(i,j)                qdiag(i,j,iosrclr,bi,bj) = qdiag(i,j,iosrclr,bi,bj) +
109        enddo                                                                           .                         (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)        
110        enddo                                                                            enddo
111        endif                                                                            enddo
112          endif
113                                                                                                                                                                    
114  c Upward Longwave Flux at the Ground (W/m**2)  c Upward Longwave Flux at the Ground (W/m**2)
115  c -------------------------------------------  c -------------------------------------------
116        if (ilwgup.ne.0) then                                                            if (ilwgup.ne.0) then
117        do j=jm1,jm2        do j=jm1,jm2
118        do i=im1,im2        do i=im1,im2
119        qdiag(i,j,ilwgup) =  qdiag(i,j,ilwgup) + coup%lw%st4(i,j)        qdiag(i,j,ilwgup,bi,bj) = qdiag(i,j,ilwgup,bi,bj) + st4(i,j,bi,bj)
120       .                   +   coup%lw%dst4(i,j)*(coup%land%tgz(i,j)-coup%lw%tg0(i,j))       .                 + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
121        enddo                                                                            enddo
122        enddo                                                                            enddo
123        endif                                                                            endif
124                                                                                                                                                                    
125  c Net Longwave Flux at the Ground (W/m**2)  c Net Longwave Flux at the Ground (W/m**2)
126  c ----------------------------------------  c ----------------------------------------
127        if (iradlwg.ne.0) then                                                            if (iradlwg.ne.0) then
128        do j=jm1,jm2        do j=jm1,jm2
129        do i=im1,im2        do i=im1,im2
130        qdiag(i,j,iradlwg) =  qdiag(i,j,iradlwg) + coup%lw%radlwg(i,j)        qdiag(i,j,iradlwg,bi,bj) =  qdiag(i,j,iradlwg,bi,bj) +
131       .                   +   coup%lw%dst4(i,j)*(coup%land%tgz(i,j)-coup%lw%tg0(i,j))       .                                              radlwg(i,j,bi,bj) +
132        enddo                                                                           .                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
133        enddo                                                                            enddo
134        endif                                                                            enddo
135          endif
136                                                                                                                                                                    
137  c Net Longwave Flux at the Ground Clear Sky (W/m**2)  c Net Longwave Flux at the Ground Clear Sky (W/m**2)
138  c --------------------------------------------------  c --------------------------------------------------
139        if (ilwgclr.ne.0) then        if (ilwgclr.ne.0) then
140        do j=jm1,jm2        do j=jm1,jm2
141        do i=im1,im2        do i=im1,im2
142        qdiag(i,j,ilwgclr) =  qdiag(i,j,ilwgclr) + coup%lw%lwgclr(i,j)        qdiag(i,j,ilwgclr,bi,bj) = qdiag(i,j,ilwgclr,bi,bj) +
143       .                   +   coup%lw%dst4(i,j)*(coup%land%tgz(i,j)-coup%lw%tg0(i,j))       .                                               lwgclr(i,j,bi,bj) +
144        enddo                                                                           .                   dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
       enddo                                                                      
       endif                                                                      
                                                                                   
 c Total Surface Pressure Tendency (mb/day)  
 c ----------------------------------------  
       if( idpdt.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,idpdt) = qdiag(i,j,idpdt) + (updt%p(i,j)-curr%p(i,j))*fact  
       enddo  
       enddo  
       endif  
                                                                     
 c Averaged P-Field (mb)  
 c ---------------------  
       if( ips.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ips) = qdiag(i,j,ips) + curr%p(i,j)+ptop  
       enddo  
       enddo  
       endif  
   
 c Averaged SLP-Field (mb)  
 c ----------------------  
       if( islp.ne.0 ) then  
       do L=1,Nrphys  
       do j=jm1,jm2  
       do i=im1,im2  
       tmp2(i,j,L) = curr%t(i,j,L) * (1.+0.609*curr%q(i,j,L,1))  
       enddo  
       enddo  
       enddo  
       call slprs ( tmp1(1,1,2),curr%p,ptop,coup%earth%phis_cmp,tmp2,dsig,coup%earth%lw_cmp,im,jm,lm )  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,islp) = qdiag(i,j,islp) + tmp1(i,j,2)  
145        enddo        enddo
146        enddo        enddo
147        endif        endif
148                                                                                    
149        npiau   = npiau   + 1        if( (bi.eq.1) .and. (bj.eq.1) ) then
150        nradswt = nradswt + 1                                                            nradswt = nradswt + 1
151        nradswg = nradswg + 1                                                            nradswg = nradswg + 1
152        nswgclr = nswgclr + 1                                                            nswgclr = nswgclr + 1
153        nosr    = nosr    + 1        nosr    = nosr    + 1
154        nosrclr = nosrclr + 1        nosrclr = nosrclr + 1
155        nradlwg = nradlwg + 1        nradlwg = nradlwg + 1
156        nlwgclr = nlwgclr + 1        nlwgclr = nlwgclr + 1
157        nlwgup  = nlwgup  + 1        nlwgup  = nlwgup  + 1
158        ndpdt   = ndpdt   + 1                                                      endif
       nps     = nps     + 1  
       nslp    = nslp    + 1  
159                                                                                                                                                                    
160  C **********************************************************************          C **********************************************************************        
 C ****                  Compute 3-D Diagnostics                     ****          
 C **********************************************************************          
                                                                                   
161        do L=1,Nrphys        do L=1,Nrphys
162    
       if( itiau.ne.0 .or. ivavetiau.ne.0 .or. idiabt.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       dtiau(i,j,L) = ( tend%iau%dt(i,j,L) + curr%t(i,j,L)*tend%iau%dp(i,j)  
      .               * ( sig(L)*akap*curr%p(i,j)/(sig(L)*curr%p(i,j)+ptop) - 1.0) )  
       enddo  
       enddo  
       endif  
   
       if( iqiau.ne.0 .or. ivaveqiau.ne.0 .or. idiabq.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       dqiau(i,j,L) = tend%iau%dq(i,j,L) - curr%q(i,j,L,1)*tend%iau%dp(i,j)  
       enddo  
       enddo  
       endif  
   
163  c Total Diabatic U-Tendency (m/sec/day)  c Total Diabatic U-Tendency (m/sec/day)
164  c -------------------------------------  c -------------------------------------
165        if( idiabu.ne.0 ) then        if( idiabu.ne.0 ) then
166        do j=jm1,jm2        do j=jm1,jm2
167        do i=im1,im2        do i=im1,im2
168        qdiag(i,j,idiabu+L-1) =   qdiag  (i,j,idiabu+L-1)        qdiag(i,j,idiabu+L-1,bi,bj) = qdiag(i,j,idiabu+L-1,bi,bj)
169       .                      + ( tend%iau%du (i,j,L)       .            + ( moistu (i,j,L,bi,bj) + turbu(i,j,L,bi,bj) )*86400
      .                        + tend%turb%du(i,j,L) )*86400  
170        enddo        enddo
171        enddo        enddo
172        endif        endif
# Line 220  c ------------------------------------- Line 176  c -------------------------------------
176        if( idiabv.ne.0 ) then        if( idiabv.ne.0 ) then
177        do j=jm1,jm2        do j=jm1,jm2
178        do i=im1,im2        do i=im1,im2
179        qdiag(i,j,idiabv+L-1) =   qdiag  (i,j,idiabv+L-1)        qdiag(i,j,idiabv+L-1,bi,bj) = qdiag(i,j,idiabv+L-1,bi,bj)
180       .                      + ( tend%iau%dv (i,j,L)       .            + ( moistv (i,j,L,bi,bj) + turbv(i,j,L,bi,bj) )*86400
      .                        + tend%turb%dv(i,j,L) )*86400  
181        enddo        enddo
182        enddo        enddo
183        endif        endif
# Line 232  c ----------------------------------- Line 187  c -----------------------------------
187        if( idiabt.ne.0 ) then        if( idiabt.ne.0 ) then
188        do j=jm1,jm2        do j=jm1,jm2
189        do i=im1,im2        do i=im1,im2
190        qdiag(i,j,idiabt+L-1) =   qdiag(i,j,idiabt+L-1)        qdiag(i,j,idiabt+L-1,bi,bj) = qdiag(i,j,idiabt+L-1,bi,bj) +
191       .                      + ( dtiau(i,j,L) + tend%turb%dt(i,j,L) + tend%lw%dt(i,j,L)       .   ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +
192       .                        + coup%lw%dlwdtg(i,j,L)*(coup%land%tgc(i,j)-coup%lw%tg0c(i,j))       .      lwdt(i,j,L,bi,bj) +
193       .                        + tend%sw%dt(i,j,L)*radswc(i,j)       .      dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +
194       .                        + tend%moist%dt(i,j,L) )       .      swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )
195       .                        * pkn(i,j,L)*pinv(i,j)*86400       .      * pk(i,j,L,bi,bj)*pinv(i,j)*86400
196        enddo        enddo
197        enddo        enddo
198        endif        endif
# Line 247  c ------------------------------------ Line 202  c ------------------------------------
202        if( idiabq.ne.0 ) then        if( idiabq.ne.0 ) then
203        do j=jm1,jm2        do j=jm1,jm2
204        do i=im1,im2        do i=im1,im2
205        qdiag(i,j,idiabq+L-1) =   qdiag(i,j,idiabq+L-1)        qdiag(i,j,idiabq+L-1,bi,bj) =   qdiag(i,j,idiabq+L-1,bi,bj) +
206       .                      + ( dqiau(i,j,L) + tend%turb%dq(i,j,L,1) + tend%moist%dq(i,j,L,1) )       . ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) *
207       .                      *    pinv(i,j)*86400*1000       .                                      pinv(i,j)*86400*1000
208        enddo        enddo
209        enddo        enddo
210        endif        endif
211            
212  c Analysis U-Wind Increment (m/sec/day)  c Longwave Heating (deg/day)
213  c -------------------------------------  c --------------------------
214        if( iuiau.ne.0 ) then        if (iradlw.ne.0) then
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,iuiau+L-1) = qdiag (i,j,iuiau+L-1) + tend%iau%du(i,j,L)*86400  
       enddo  
       enddo  
       endif  
   
 c Analysis V-Wind Increment (m/sec/day)  
 c -------------------------------------  
       if( iviau.ne.0 ) then  
215        do j=jm1,jm2        do j=jm1,jm2
216        do i=im1,im2        do i=im1,im2
217        qdiag(i,j,iviau+L-1) = qdiag (i,j,iviau+L-1) + tend%iau%dv(i,j,L)*86400        qdiag(i,j,iradlw+l-1,bi,bj) = qdiag(i,j,iradlw+l-1,bi,bj) +
218         . ( lwdt(i,j,l,bi,bj) +
219         .            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
220         .                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400
221        enddo        enddo
222        enddo        enddo
223        endif        endif
224    
 c Analysis Temperature Tendency (deg/day)  
 c ---------------------------------------  
       if( itiau.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,itiau+L-1) = qdiag(i,j,itiau+L-1)  
      .                     + dtiau(i,j,L) * pkn(i,j,L)*pinv(i,j)*86400  
       enddo  
       enddo  
       endif  
                                                                     
 c Analysis Moisture Tendency (g/kg/day)  
 c -------------------------------------  
       if( iqiau.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,iqiau+L-1) = qdiag(i,j,iqiau+L-1)  
      .                     + dqiau(i,j,L) * 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) = qdiag(i,j,iradlw+l-1)  
      .                      + ( tend%lw%dt(i,j,l)  
      .                      +   coup%lw%dlwdtg (i,j,L)*(coup%land%tgc(i,j)-coup%lw%tg0c(i,j)))  
      .                      * pkn(i,j,l)*pinv(i,j)*86400  
       enddo                                                                      
       enddo                                                                      
       endif                                                                      
                                                                                   
225  c Longwave Heating Clear-Sky (deg/day)  c Longwave Heating Clear-Sky (deg/day)
226  c ------------------------------------  c ------------------------------------
227        if (ilwclr.ne.0) then                                                            if (ilwclr.ne.0) then                                                    
228        do j=jm1,jm2        do j=jm1,jm2
229        do i=im1,im2        do i=im1,im2
230        qdiag(i,j,ilwclr+l-1) = qdiag(i,j,ilwclr+l-1)        qdiag(i,j,ilwclr+l-1,bi,bj) = qdiag(i,j,ilwclr+l-1,bi,bj) +
231       .                      + ( tend%lw%dtclr(i,j,l)       . ( lwdtclr(i,j,l,bi,bj) +
232       .                      +   coup%lw%dlwdtg(i,j,L)*(coup%land%tgc(i,j)-coup%lw%tg0c(i,j)))       .             dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
233       .                      * pkn(i,j,l)*pinv(i,j)*86400       .                      * 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) = qdiag(i,j,iradsw+l-1)  
      .                      + tend%sw%dt(i,j,l)*radswc(i,j)                          
      .                      * pkn(i,j,l)*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) = qdiag(i,j,iswclr+l-1)  
      .                      + tend%sw%dtclr(i,j,l)*radswc(i,j)                          
      .                      * pkn(i,j,l)*pinv(i,j)*86400  
       enddo                                                                      
       enddo                                                                      
       endif                                                                      
                                                                                   
 c Total U-Tendency (m/sec/day)  
 c ----------------------------  
       if( idudt.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,idudt+L-1) =  qdiag(i,j,idudt+L-1)  
      .                     + ( updt%u(i,j,L)-curr%u(i,j,L) )*fact  
234        enddo        enddo
235        enddo        enddo
236        endif        endif
237                                                                                                                                                      
238  c Total V-Tendency (m/sec/day)  c Solar Radiative Heating (deg/day)
239  c ----------------------------  c ---------------------------------
240        if( idvdt.ne.0 ) then        if (iradsw.ne.0) then
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,idvdt+L-1) =  qdiag(i,j,idvdt+L-1)  
      .                     + ( updt%v(i,j,L)-curr%v(i,j,L) )*fact  
       enddo  
       enddo  
       endif  
   
 c Total T-Tendency (deg/day)  
 c --------------------------  
       if( idtdt.ne.0 ) then  
241        do j=jm1,jm2        do j=jm1,jm2
242        do i=im1,im2        do i=im1,im2
243        qdiag(i,j,idtdt+L-1) =  qdiag(i,j,idtdt+L-1)        qdiag(i,j,iradsw+l-1,bi,bj) = qdiag(i,j,iradsw+l-1,bi,bj) +
244       .                     + ( updt%t(i,j,L)*pknp1(i,j,L) - curr%t(i,j,L)*pkn(i,j,L) )*fact       .  + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
245         .                   pk(i,j,l,bi,bj)*pinv(i,j)*86400
246        enddo        enddo
247        enddo        enddo
248        endif        endif
249                                                                                    
250  c Total Q-Tendency (g/kg/day)  c Clear Sky Solar Radiative Heating (deg/day)
251  c ---------------------------  c -------------------------------------------
252        if( idqdt.ne.0 ) then        if (iswclr.ne.0) then
253        do j=jm1,jm2        do j=jm1,jm2
254        do i=im1,im2        do i=im1,im2
255        qdiag(i,j,idqdt+L-1) =  qdiag(i,j,idqdt+L-1)        qdiag(i,j,iswclr+l-1,bi,bj) = qdiag(i,j,iswclr+l-1,bi,bj) +
256       .                     + ( updt%q(i,j,L,1)-curr%q(i,j,L,1) )*fact*1000       .           swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
257         .                           pk(i,j,l,bi,bj)*pinv(i,j)*86400
258        enddo        enddo
259        enddo        enddo
260        endif        endif
261                                                                                    
262  c Averaged U-Field (m/sec)  c Averaged U-Field (m/sec)
263  c ------------------------  c ------------------------
264        if( iuwnd.ne.0 ) then        if( iuwnd.ne.0 ) then
265        do j=jm1,jm2        do j=jm1,jm2
266        do i=im1,im2        do i=im1,im2
267        qdiag(i,j,iuwnd+L-1) = qdiag(i,j,iuwnd+L-1) + curr%u(i,j,L)        qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
268         .                                                 uphy(i,j,L,bi,bj)
269        enddo        enddo
270        enddo        enddo
271        endif        endif
# Line 405  c ------------------------ Line 275  c ------------------------
275        if( ivwnd.ne.0 ) then        if( ivwnd.ne.0 ) then
276        do j=jm1,jm2        do j=jm1,jm2
277        do i=im1,im2        do i=im1,im2
278        qdiag(i,j,ivwnd+L-1) = qdiag(i,j,ivwnd+L-1) + curr%v(i,j,L)        qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
279         .                                                 vphy(i,j,L,bi,bj)
280        enddo        enddo
281        enddo        enddo
282        endif        endif
# Line 415  c ---------------------- Line 286  c ----------------------
286        if( itmpu.ne.0 ) then        if( itmpu.ne.0 ) then
287        do j=jm1,jm2        do j=jm1,jm2
288        do i=im1,im2        do i=im1,im2
289        qdiag(i,j,itmpu+L-1) = qdiag(i,j,itmpu+L-1) + curr%t(i,j,L)*pkn(i,j,L)        qdiag(i,j,itmpu+L-1,bi,bj) = qdiag(i,j,itmpu+L-1,bi,bj) +
290         .                               thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
291        enddo        enddo
292        enddo        enddo
293        endif        endif
# Line 425  c ---------------------------- Line 297  c ----------------------------
297        if( itke.ne.0 ) then        if( itke.ne.0 ) then
298        do j=jm1,jm2        do j=jm1,jm2
299        do i=im1,im2        do i=im1,im2
300        qdiag(i,j,itke+L-1) = qdiag(i,j,itke+L-1) + coup%turb%qq(i,j,L)        qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
301        enddo        enddo
302        enddo        enddo
303        endif        endif
# Line 435  c ----------------------- Line 307  c -----------------------
307        if( isphu.ne.0 ) then        if( isphu.ne.0 ) then
308        do j=jm1,jm2        do j=jm1,jm2
309        do i=im1,im2        do i=im1,im2
310        qdiag(i,j,isphu+L-1) = qdiag(i,j,isphu+L-1) + curr%q(i,j,L,1)*1000        qdiag(i,j,isphu+L-1,bi,bj) = qdiag(i,j,isphu+L-1,bi,bj) +
311         .                                            sphy(i,j,L,bi,bj)*1000
312        enddo        enddo
313        enddo        enddo
314        endif        endif
315                                                                                                                                                                    
316        enddo   ! End Level Loop        enddo
317    
318          if( (bi.eq.1) .and. (bj.eq.1) ) then
319    
320        ndiabu = ndiabu + 1        ndiabu = ndiabu + 1
321        ndiabv = ndiabv + 1        ndiabv = ndiabv + 1
322        ndiabt = ndiabt + 1        ndiabt = ndiabt + 1
323        nuiau  = nuiau  + 1        ndiabq = ndiabq + 1
       nviau  = nviau  + 1  
       ntiau  = ntiau  + 1  
324        nradlw = nradlw + 1                                                            nradlw = nradlw + 1                                                    
325        nlwclr = nlwclr + 1                                                            nlwclr = nlwclr + 1                                                    
326        nradsw = nradsw + 1                                                            nradsw = nradsw + 1                                                    
327        nswclr = nswclr + 1                                                            nswclr = nswclr + 1                                                    
       ndudt  = ndudt  + 1                                                
       ndvdt  = ndvdt  + 1                                                
       ndtdt  = ndtdt  + 1                                                
       ndiabq = ndiabq + 1  
       nqiau  = nqiau  + 1  
       ndqdt  = ndqdt  + 1                                                
328        nuwnd  = nuwnd  + 1        nuwnd  = nuwnd  + 1
329        nvwnd  = nvwnd  + 1                                  nvwnd  = nvwnd  + 1                          
330        ntmpu  = ntmpu  + 1                                  ntmpu  = ntmpu  + 1                          
331        ntke   = ntke   + 1        ntke   = ntke   + 1
332        nsphu  = nsphu  + 1        nsphu  = nsphu  + 1
333    
 C **********************************************************************          
 C ****         Compute Vertically Integrated Diagnostics            ****          
 C **********************************************************************          
   
 c Compute Moisture and Temperature Convergence Diagnostic  
 c -------------------------------------------------------  
       if( ivaveut.ne.0 .or. ivavevt.ne.0  .or.  
      .    ivaveuq.ne.0 .or. ivavevq.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       vintut(i,j) = 0.0  
       vintvt(i,j) = 0.0  
       vintuq(i,j) = 0.0  
       vintvq(i,j) = 0.0  
       enddo  
       enddo  
   
       call ctoa   ( curr%t,curr%t,dlam,dphi,im,jm,lm,0,grid%lattice )  
       call ctoa   ( curr%q,curr%q,dlam,dphi,im,jm,lm,0,grid%lattice )  
       call ctoa_winds ( curr%u,curr%v,tmp1,tmp2,dlam,dphi,im,jm,lm,grid%lattice )  
       do L=1,Nrphys  
       do j=jm1,jm2  
       do i=im1,im2  
       vintut(i,j) = vintut(i,j) + tmp1(i,j,L)*curr%t(i,j,L)  *dsig(L)* pkn(i,j,L)  
       vintvt(i,j) = vintvt(i,j) + tmp2(i,j,L)*curr%t(i,j,L)  *dsig(L)* pkn(i,j,L)  
       vintuq(i,j) = vintuq(i,j) + tmp1(i,j,L)*curr%q(i,j,L,1)*dsig(L)  
       vintvq(i,j) = vintvq(i,j) + tmp2(i,j,L)*curr%q(i,j,L,1)*dsig(L)  
       enddo  
       enddo  
       enddo  
   
       if( ivaveut.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivaveut) = qdiag(i,j,ivaveut) + vintut(i,j)  
       enddo  
       enddo  
       endif  
       if( ivavevt.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivavevt) = qdiag(i,j,ivavevt) + vintvt(i,j)  
       enddo  
       enddo  
       endif  
       if( ivaveuq.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivaveuq) = qdiag(i,j,ivaveuq) + vintuq(i,j)*1000  
       enddo  
       enddo  
       endif  
       if( ivavevq.ne.0 ) then  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivavevq) = qdiag(i,j,ivavevq) + vintvq(i,j)*1000  
       enddo  
       enddo  
334        endif        endif
335    
336        endif    ! End Convergence Diagnostic  C **********************************************************************        
   
 c Total Precipitable Water (gm/cm**2)  
 c -----------------------------------  
       if( itpw.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) + curr%q(i,j,L,1)*dsig(L)  
       enddo  
       enddo  
       enddo  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,itpw) = qdiag(i,j,itpw) + qbar(i,j)*10*curr%p(i,j)/grav  
       enddo  
       enddo  
       endif  
                                                                                   
 c Total Precipitable Analysis Increment (mm/day)  
 c ----------------------------------------------  
       if( ivaveqiau.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) + dqiau(i,j,L)*dsig(L)  
       enddo  
       enddo  
       enddo  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivaveqiau) = qdiag(i,j,ivaveqiau) + qbar(i,j)*(100*86400/grav)  
       enddo  
       enddo  
       endif  
   
 c Vertically Averaged Analysis Temperature Increment (K/day)  
 c ----------------------------------------------------------  
       if( ivavetiau.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) + dtiau(i,j,L)*pkn(i,j,l)*dsig(L)  
       enddo  
       enddo  
       enddo  
       do j=jm1,jm2  
       do i=im1,im2  
       qdiag(i,j,ivavetiau) = qdiag(i,j,ivavetiau) + qbar(i,j)*pinv(i,j)*86400  
       enddo  
       enddo  
       endif  
337    
338  c Vertically Averaged Moist-T Increment (K/day)  c Vertically Averaged Moist-T Increment (K/day)
339  c ---------------------------------------------  c ---------------------------------------------
# Line 603  c -------------------------------------- Line 346  c --------------------------------------
346        do L=1,Nrphys        do L=1,Nrphys
347        do j=jm1,jm2        do j=jm1,jm2
348        do i=im1,im2        do i=im1,im2
349        qbar(i,j) = qbar(i,j) + tend%moist%dt(i,j,L)*pkn(i,j,l)*dsig(L)        qbar(i,j) = qbar(i,j) +
350         .             moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
351        enddo        enddo
352        enddo        enddo
353        enddo        enddo
354        do j=jm1,jm2        do j=jm1,jm2
355        do i=im1,im2        do i=im1,im2
356        qdiag(i,j,ivdtmoist) = qdiag(i,j,ivdtmoist) + qbar(i,j)*pinv(i,j)*86400        qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
357         .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
358        enddo        enddo
359        enddo        enddo
360        endif        endif
# Line 625  c -------------------------------------- Line 370  c --------------------------------------
370        do L=1,Nrphys        do L=1,Nrphys
371        do j=jm1,jm2        do j=jm1,jm2
372        do i=im1,im2        do i=im1,im2
373        qbar(i,j) = qbar(i,j) + tend%turb%dt(i,j,L)*pkn(i,j,l)*dsig(L)        qbar(i,j) = qbar(i,j) +
374         .             turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
375        enddo        enddo
376        enddo        enddo
377        enddo        enddo
378        do j=jm1,jm2        do j=jm1,jm2
379        do i=im1,im2        do i=im1,im2
380        qdiag(i,j,ivdtturb) = qdiag(i,j,ivdtturb) + qbar(i,j)*pinv(i,j)*86400        qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
381         .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
382        enddo        enddo
383        enddo        enddo
384        endif        endif
# Line 647  c -------------------------------------- Line 394  c --------------------------------------
394        do L=1,Nrphys        do L=1,Nrphys
395        do j=jm1,jm2        do j=jm1,jm2
396        do i=im1,im2        do i=im1,im2
397        qbar(i,j) = qbar(i,j) + ( tend%lw%dt(i,j,L)        qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
398       .          + coup%lw%dlwdtg(i,j,L)*(coup%land%tgc(i,j)-coup%lw%tg0c(i,j)) )*pkn(i,j,l)*dsig(L)       .  dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
399         .             *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
400        enddo        enddo
401        enddo        enddo
402        enddo        enddo
403        do j=jm1,jm2        do j=jm1,jm2
404        do i=im1,im2        do i=im1,im2
405        qdiag(i,j,ivdtradlw) = qdiag(i,j,ivdtradlw) + qbar(i,j)*pinv(i,j)*86400        qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
406         .      qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
407        enddo        enddo
408        enddo        enddo
409        endif        endif
# Line 670  c -------------------------------------- Line 419  c --------------------------------------
419        do L=1,Nrphys        do L=1,Nrphys
420        do j=jm1,jm2        do j=jm1,jm2
421        do i=im1,im2        do i=im1,im2
422        qbar(i,j) = qbar(i,j) + tend%sw%dt(i,j,L)*pkn(i,j,l)*dsig(L)        qbar(i,j) = qbar(i,j) +
423         .             swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
424        enddo        enddo
425        enddo        enddo
426        enddo        enddo
427        do j=jm1,jm2        do j=jm1,jm2
428        do i=im1,im2        do i=im1,im2
429        qdiag(i,j,ivdtradsw) = qdiag(i,j,ivdtradsw) + qbar(i,j)*radswc(i,j)*pinv(i,j)*86400        qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
430         . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
431        enddo        enddo
432        enddo        enddo
433        endif        endif
434    
435        nvaveut   = nvaveut   + 1        if( (bi.eq.1) .and. (bj.eq.1) ) then
       nvavevt   = nvavevt   + 1  
       nvaveuq   = nvaveuq   + 1  
       nvavevq   = nvavevq   + 1  
       ntpw      = ntpw      + 1  
       nvaveqiau = nvaveqiau + 1  
       nvavetiau = nvavetiau + 1  
436        nvdtmoist = nvdtmoist + 1        nvdtmoist = nvdtmoist + 1
437        nvdtturb  = nvdtturb  + 1        nvdtturb  = nvdtturb  + 1
438        nvdtradlw = nvdtradlw + 1        nvdtradlw = nvdtradlw + 1
439        nvdtradsw = nvdtradsw + 1        nvdtradsw = nvdtradsw + 1
440          endif
441    
442  C *****************************************************************      #endif
 C ****                  Release Workspace                      ****      
 C *****************************************************************      
   
       deallocate (  dlam )  
       deallocate (  dphi )  
       deallocate (  dsig )  
       deallocate (   sig )  
       deallocate (  sige )  
   
       deallocate (  pinv  )  
       deallocate (  tmp1  )  
       deallocate (  tmp2  )  
       deallocate (   qbar )  
       deallocate ( vintuq )  
       deallocate ( vintvq )  
       deallocate ( vintut )  
       deallocate ( vintvt )  
       deallocate (    pkn )  
       deallocate (  pknp1 )  
       deallocate (  dtiau )  
       deallocate (  dqiau )  
   
       call timeend (' step_diag')  
443        return        return
444        end        end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22