Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:40:24 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
4cd793223d Jean*0002       SUBROUTINE FIZHI_STEP_DIAG(myid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
                0003      &  radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
                0004      &  turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
                0005      &  lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
                0006      &  im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
360ae2c8e1 Andr*0007 C***********************************************************************
4cd793223d Jean*0008       IMPLICIT NONE
360ae2c8e1 Andr*0009 
4cd793223d Jean*0010       INTEGER myid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
a456aa407c Andr*0011       _RL p(im2,jm2,Nbi,Nbj)
07a60b5432 Andr*0012       _RL uphy(im2,jm2,Nrphys)
                0013       _RL vphy(im2,jm2,Nrphys)
                0014       _RL thphy(im2,jm2,Nrphys)
                0015       _RL sphy(im2,jm2,Nrphys)
1144b0a379 Jean*0016       _RL qq(im2,jm2,Nrphys,Nbi,Nbj),pk(im2,jm2,Nrphys,Nbi,Nbj)
a456aa407c Andr*0017       _RL dp(im2,jm2,Nrphys,Nbi,Nbj)
                0018       _RL radswt(im2,jm2,Nbi,Nbj),radswg(im2,jm2,Nbi,Nbj)
                0019       _RL swgclr(im2,jm2,Nbi,Nbj),osr(im2,jm2,Nbi,Nbj)
                0020       _RL osrclr(im2,jm2,Nbi,Nbj),st4(im2,jm2,Nbi,Nbj)
                0021       _RL dst4(im2,jm2,Nbi,Nbj),tgz(im2,jm2,Nbi,Nbj)
                0022       _RL tg0(im2,jm2,Nbi,Nbj),radlwg(im2,jm2,Nbi,Nbj)
                0023       _RL lwgclr(im2,jm2,Nbi,Nbj)
daaf4a8d7f Andr*0024       _RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
                0025       _RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
a456aa407c Andr*0026       _RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
                0027       _RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
daaf4a8d7f Andr*0028       _RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
                0029       _RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
a456aa407c Andr*0030       _RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
                0031       _RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
daaf4a8d7f Andr*0032       _RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
                0033       _RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
a456aa407c Andr*0034       _RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
                0035       _RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
                0036       _RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)
272e3d01f8 Andr*0037 
4cd793223d Jean*0038       INTEGER  i,j,L
f0743b271c Andr*0039       _RL getcon, gravity
8b150b4af9 Andr*0040       _RL pinv(im2,jm2), qbar(im2,jm2),tmpdiag(im2,jm2)
                0041 #ifdef ALLOW_DIAGNOSTICS
4cd793223d Jean*0042       LOGICAL  diagnostics_is_on
                0043       EXTERNAL diagnostics_is_on
8b150b4af9 Andr*0044 #endif
360ae2c8e1 Andr*0045 
4cd793223d Jean*0046 C **********************************************************************
360ae2c8e1 Andr*0047 
09fc20bde8 Andr*0048 #ifdef ALLOW_DIAGNOSTICS
360ae2c8e1 Andr*0049       do j=jm1,jm2
                0050       do i=im1,im2
272e3d01f8 Andr*0051       pinv(i,j) = 1.0 / p(i,j,bi,bj)
360ae2c8e1 Andr*0052       enddo
                0053       enddo
272e3d01f8 Andr*0054 
a5606e7e37 Andr*0055 c Surface Pressure (mb)
                0056 c ---------------------------------
4cd793223d Jean*0057       call diagnostics_fill(p(1,1,bi,bj),'PS      ',0,1,3,bi,bj,myid)
a5606e7e37 Andr*0058 
360ae2c8e1 Andr*0059 c Incident Solar Radiation (W/m**2)
                0060 c ---------------------------------
4cd793223d Jean*0061       call diagnostics_fill(radswt(1,1,bi,bj),'RADSWT  ',
                0062      &                      0,1,3,bi,bj,myid)
8b150b4af9 Andr*0063 
360ae2c8e1 Andr*0064 c Net Solar Radiation at the Ground (W/m**2)
                0065 c ------------------------------------------
8b150b4af9 Andr*0066       if(diagnostics_is_on('RADSWG  ',myid) ) then
                0067        do j=jm1,jm2
                0068        do i=im1,im2
                0069         tmpdiag(i,j) = radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
                0070        enddo
                0071        enddo
                0072        call diagnostics_fill(tmpdiag,'RADSWG  ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0073       endif
4cd793223d Jean*0074 
360ae2c8e1 Andr*0075 c Net Clear Sky Solar Radiation at the Ground (W/m**2)
                0076 c ----------------------------------------------------
8b150b4af9 Andr*0077       if(diagnostics_is_on('SWGCLR  ',myid) ) then
                0078        do j=jm1,jm2
                0079        do i=im1,im2
                0080         tmpdiag(i,j) = swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
                0081        enddo
                0082        enddo
                0083        call diagnostics_fill(tmpdiag,'SWGCLR  ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0084       endif
4cd793223d Jean*0085 
272e3d01f8 Andr*0086 c Outgoing Solar Radiation at top (W/m**2)
360ae2c8e1 Andr*0087 c -----------------------------------------
8b150b4af9 Andr*0088       if(diagnostics_is_on('OSR     ',myid) ) then
                0089        do j=jm1,jm2
                0090        do i=im1,im2
                0091         tmpdiag(i,j) = (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
                0092        enddo
                0093        enddo
                0094        call diagnostics_fill(tmpdiag,'OSR     ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0095       endif
4cd793223d Jean*0096 
272e3d01f8 Andr*0097 c Outgoing Clear Sky Solar Radiation at top (W/m**2)
360ae2c8e1 Andr*0098 c ---------------------------------------------------
8b150b4af9 Andr*0099       if(diagnostics_is_on('OSRCLR  ',myid) ) then
                0100        do j=jm1,jm2
                0101        do i=im1,im2
                0102         tmpdiag(i,j) = (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)
                0103        enddo
                0104        enddo
                0105        call diagnostics_fill(tmpdiag,'OSRCLR  ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0106       endif
4cd793223d Jean*0107 
590ca7dd42 Andr*0108 c Planetary Albedo
                0109 c ----------------
                0110       if(diagnostics_is_on('PLALBEDO',myid) ) then
                0111        do j=jm1,jm2
                0112        do i=im1,im2
                0113         if(radswt(i,j,bi,bj).ne.0.) then
                0114          tmpdiag(i,j) = osr(i,j,bi,bj)
                0115         else
                0116          tmpdiag(i,j) = 0.
                0117         endif
                0118        enddo
                0119        enddo
                0120        call diagnostics_fill(tmpdiag,'PLALBEDO',0,1,3,bi,bj,myid)
                0121       endif
4cd793223d Jean*0122 
360ae2c8e1 Andr*0123 c Upward Longwave Flux at the Ground (W/m**2)
                0124 c -------------------------------------------
8b150b4af9 Andr*0125       if(diagnostics_is_on('LWGUP   ',myid) ) then
                0126        do j=jm1,jm2
                0127        do i=im1,im2
4cd793223d Jean*0128         tmpdiag(i,j) = st4(i,j,bi,bj)
                0129      &                 + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
8b150b4af9 Andr*0130        enddo
                0131        enddo
                0132        call diagnostics_fill(tmpdiag,'LWGUP   ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0133       endif
4cd793223d Jean*0134 
360ae2c8e1 Andr*0135 c Net Longwave Flux at the Ground (W/m**2)
                0136 c ----------------------------------------
8b150b4af9 Andr*0137       if(diagnostics_is_on('RADLWG  ',myid) ) then
                0138        do j=jm1,jm2
                0139        do i=im1,im2
4cd793223d Jean*0140         tmpdiag(i,j) = radlwg(i,j,bi,bj) +
                0141      &                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
8b150b4af9 Andr*0142        enddo
                0143        enddo
                0144        call diagnostics_fill(tmpdiag,'RADLWG  ',0,1,3,bi,bj,myid)
272e3d01f8 Andr*0145       endif
4cd793223d Jean*0146 
360ae2c8e1 Andr*0147 c Net Longwave Flux at the Ground Clear Sky (W/m**2)
                0148 c --------------------------------------------------
8b150b4af9 Andr*0149       if(diagnostics_is_on('LWGCLR  ',myid) ) then
                0150        do j=jm1,jm2
                0151        do i=im1,im2
4cd793223d Jean*0152         tmpdiag(i,j) = lwgclr(i,j,bi,bj) +
                0153      &                  dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
8b150b4af9 Andr*0154        enddo
                0155        enddo
                0156        call diagnostics_fill(tmpdiag,'LWGCLR  ',0,1,3,bi,bj,myid)
471f7f004b Andr*0157       endif
4cd793223d Jean*0158 
                0159 C **********************************************************************
360ae2c8e1 Andr*0160       do L=1,Nrphys
                0161 
                0162 c Total Diabatic U-Tendency (m/sec/day)
                0163 c -------------------------------------
8b150b4af9 Andr*0164       if(diagnostics_is_on('DIABU   ',myid) ) then
                0165        do j=jm1,jm2
                0166        do i=im1,im2
                0167         tmpdiag(i,j) = (moistu (i,j,L,bi,bj)+turbu(i,j,L,bi,bj) )*86400
                0168        enddo
                0169        enddo
                0170        call diagnostics_fill(tmpdiag,'DIABU   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0171       endif
8b150b4af9 Andr*0172 
360ae2c8e1 Andr*0173 c Total Diabatic V-Tendency (m/sec/day)
                0174 c -------------------------------------
8b150b4af9 Andr*0175       if(diagnostics_is_on('DIABV   ',myid) ) then
                0176        do j=jm1,jm2
                0177        do i=im1,im2
                0178         tmpdiag(i,j) = (moistv (i,j,L,bi,bj)+turbv(i,j,L,bi,bj) )*86400
                0179        enddo
                0180        enddo
                0181        call diagnostics_fill(tmpdiag,'DIABV   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0182       endif
                0183 
                0184 c Total Diabatic T-Tendency (deg/day)
                0185 c -----------------------------------
8b150b4af9 Andr*0186       if(diagnostics_is_on('DIABT   ',myid) ) then
                0187        do j=jm1,jm2
                0188        do i=im1,im2
4cd793223d Jean*0189         tmpdiag(i,j) =
                0190      &   ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +
                0191      &      lwdt(i,j,L,bi,bj) +
                0192      &      dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +
                0193      &      swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )
                0194      &      * pk(i,j,L,bi,bj)*pinv(i,j)*86400
8b150b4af9 Andr*0195        enddo
                0196        enddo
                0197        call diagnostics_fill(tmpdiag,'DIABT   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0198       endif
8b150b4af9 Andr*0199 
360ae2c8e1 Andr*0200 c Total Diabatic Q-Tendency (g/kg/day)
                0201 c ------------------------------------
8b150b4af9 Andr*0202       if(diagnostics_is_on('DIABQ   ',myid) ) then
                0203        do j=jm1,jm2
                0204        do i=im1,im2
4cd793223d Jean*0205         tmpdiag(i,j) =
                0206      & ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) *
                0207      &                                      pinv(i,j)*86400*1000
8b150b4af9 Andr*0208        enddo
                0209        enddo
                0210        call diagnostics_fill(tmpdiag,'DIABQ   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0211       endif
4cd793223d Jean*0212 
272e3d01f8 Andr*0213 c Longwave Heating (deg/day)
                0214 c --------------------------
8b150b4af9 Andr*0215       if(diagnostics_is_on('RADLW   ',myid) ) then
                0216        do j=jm1,jm2
                0217        do i=im1,im2
4cd793223d Jean*0218         tmpdiag(i,j) =
                0219      & ( lwdt(i,j,l,bi,bj) +
                0220      &            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
                0221      &                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400
8b150b4af9 Andr*0222        enddo
                0223        enddo
                0224        call diagnostics_fill(tmpdiag,'RADLW   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0225       endif
23a6457bc5 Andr*0226 
360ae2c8e1 Andr*0227 c Longwave Heating Clear-Sky (deg/day)
                0228 c ------------------------------------
8b150b4af9 Andr*0229       if(diagnostics_is_on('LWCLR   ',myid) ) then
                0230        do j=jm1,jm2
                0231        do i=im1,im2
4cd793223d Jean*0232         tmpdiag(i,j) =
                0233      & ( lwdtclr(i,j,l,bi,bj) +
                0234      &            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
                0235      &                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400
8b150b4af9 Andr*0236        enddo
                0237        enddo
                0238        call diagnostics_fill(tmpdiag,'LWCLR   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0239       endif
4cd793223d Jean*0240 
272e3d01f8 Andr*0241 c Solar Radiative Heating (deg/day)
                0242 c ---------------------------------
8b150b4af9 Andr*0243       if(diagnostics_is_on('RADSW   ',myid) ) then
                0244        do j=jm1,jm2
                0245        do i=im1,im2
4cd793223d Jean*0246         tmpdiag(i,j) =
                0247      &  + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
                0248      &                   pk(i,j,l,bi,bj)*pinv(i,j)*86400
8b150b4af9 Andr*0249        enddo
                0250        enddo
                0251        call diagnostics_fill(tmpdiag,'RADSW   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0252       endif
4cd793223d Jean*0253 
272e3d01f8 Andr*0254 c Clear Sky Solar Radiative Heating (deg/day)
                0255 c -------------------------------------------
8b150b4af9 Andr*0256       if(diagnostics_is_on('SWCLR   ',myid) ) then
                0257        do j=jm1,jm2
                0258        do i=im1,im2
4cd793223d Jean*0259         tmpdiag(i,j) =
                0260      &  + swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
                0261      &                   pk(i,j,l,bi,bj)*pinv(i,j)*86400
8b150b4af9 Andr*0262        enddo
                0263        enddo
                0264        call diagnostics_fill(tmpdiag,'SWCLR   ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0265       endif
4cd793223d Jean*0266 
360ae2c8e1 Andr*0267 c Averaged U-Field (m/sec)
                0268 c ------------------------
8b150b4af9 Andr*0269       if(diagnostics_is_on('UWND    ',myid) ) then
                0270        do j=jm1,jm2
                0271        do i=im1,im2
07a60b5432 Andr*0272         tmpdiag(i,j) = uphy(i,j,L)
8b150b4af9 Andr*0273        enddo
                0274        enddo
                0275        call diagnostics_fill(tmpdiag,'UWND    ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0276       endif
                0277 
                0278 c Averaged V-Field (m/sec)
                0279 c ------------------------
8b150b4af9 Andr*0280       if(diagnostics_is_on('VWND    ',myid) ) then
                0281        do j=jm1,jm2
                0282        do i=im1,im2
07a60b5432 Andr*0283         tmpdiag(i,j) = vphy(i,j,L)
8b150b4af9 Andr*0284        enddo
                0285        enddo
                0286        call diagnostics_fill(tmpdiag,'VWND    ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0287       endif
                0288 
                0289 c Averaged T-Field (deg)
                0290 c ----------------------
8b150b4af9 Andr*0291       if(diagnostics_is_on('TMPU    ',myid) ) then
                0292        do j=jm1,jm2
                0293        do i=im1,im2
07a60b5432 Andr*0294         tmpdiag(i,j) = thphy(i,j,L)*pk(i,j,L,bi,bj)
8b150b4af9 Andr*0295        enddo
                0296        enddo
                0297        call diagnostics_fill(tmpdiag,'TMPU    ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0298       endif
4cd793223d Jean*0299 
360ae2c8e1 Andr*0300 c Averaged QQ-Field (m/sec)**2
                0301 c ----------------------------
8b150b4af9 Andr*0302       if(diagnostics_is_on('TKE     ',myid) ) then
                0303        do j=jm1,jm2
                0304        do i=im1,im2
                0305         tmpdiag(i,j) = qq(i,j,L,bi,bj)
                0306        enddo
                0307        enddo
                0308        call diagnostics_fill(tmpdiag,'TKE     ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0309       endif
4cd793223d Jean*0310 
360ae2c8e1 Andr*0311 c Averaged Q-Field (g/kg)
                0312 c -----------------------
8b150b4af9 Andr*0313       if(diagnostics_is_on('SPHU    ',myid) ) then
                0314        do j=jm1,jm2
                0315        do i=im1,im2
07a60b5432 Andr*0316         tmpdiag(i,j) = sphy(i,j,L) * 1000.
8b150b4af9 Andr*0317        enddo
                0318        enddo
                0319        call diagnostics_fill(tmpdiag,'SPHU    ',L,1,3,bi,bj,myid)
360ae2c8e1 Andr*0320       endif
4cd793223d Jean*0321 
272e3d01f8 Andr*0322       enddo
360ae2c8e1 Andr*0323 
4cd793223d Jean*0324 C **********************************************************************
360ae2c8e1 Andr*0325 
                0326 c Vertically Averaged Moist-T Increment (K/day)
                0327 c ---------------------------------------------
8b150b4af9 Andr*0328       if(diagnostics_is_on('VDTMOIST',myid) ) then
                0329        do j=jm1,jm2
                0330        do i=im1,im2
                0331        qbar(i,j) = 0.0
                0332        enddo
                0333        enddo
                0334        do L=1,Nrphys
                0335        do j=jm1,jm2
                0336        do i=im1,im2
4cd793223d Jean*0337        qbar(i,j) = qbar(i,j) +
                0338      &             moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
8b150b4af9 Andr*0339        enddo
                0340        enddo
                0341        enddo
                0342        do j=jm1,jm2
                0343        do i=im1,im2
                0344        tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
                0345        enddo
                0346        enddo
                0347        call diagnostics_fill(tmpdiag,'VDTMOIST',0,1,3,bi,bj,myid)
360ae2c8e1 Andr*0348       endif
                0349 
                0350 c Vertically Averaged Turb-T Increment (K/day)
                0351 c --------------------------------------------
8b150b4af9 Andr*0352       if(diagnostics_is_on('VDTTURB ',myid) ) then
                0353        do j=jm1,jm2
                0354        do i=im1,im2
                0355        qbar(i,j) = 0.0
                0356        enddo
                0357        enddo
                0358        do L=1,Nrphys
                0359        do j=jm1,jm2
                0360        do i=im1,im2
4cd793223d Jean*0361        qbar(i,j) = qbar(i,j) +
                0362      &             turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
8b150b4af9 Andr*0363        enddo
                0364        enddo
                0365        enddo
                0366        do j=jm1,jm2
                0367        do i=im1,im2
                0368        tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
                0369        enddo
                0370        enddo
                0371        call diagnostics_fill(tmpdiag,'VDTTURB ',0,1,3,bi,bj,myid)
360ae2c8e1 Andr*0372       endif
                0373 
                0374 c Vertically Averaged RADLW Temperature Increment (K/day)
                0375 c -------------------------------------------------------
8b150b4af9 Andr*0376       if(diagnostics_is_on('VDTRADLW',myid) ) then
                0377        do j=jm1,jm2
                0378        do i=im1,im2
                0379        qbar(i,j) = 0.0
                0380        enddo
                0381        enddo
                0382        do L=1,Nrphys
                0383        do j=jm1,jm2
                0384        do i=im1,im2
                0385         qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
4cd793223d Jean*0386      &  dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
                0387      &             *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
8b150b4af9 Andr*0388        enddo
                0389        enddo
                0390        enddo
                0391        do j=jm1,jm2
                0392        do i=im1,im2
                0393        tmpdiag(i,j) = qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
                0394        enddo
                0395        enddo
                0396        call diagnostics_fill(tmpdiag,'VDTRADLW',0,1,3,bi,bj,myid)
360ae2c8e1 Andr*0397       endif
                0398 
                0399 c Vertically Averaged RADSW Temperature Increment (K/day)
                0400 c -------------------------------------------------------
8b150b4af9 Andr*0401       if(diagnostics_is_on('VDTRADSW',myid) ) then
                0402        do j=jm1,jm2
                0403        do i=im1,im2
                0404        qbar(i,j) = 0.0
                0405        enddo
                0406        enddo
                0407        do L=1,Nrphys
                0408        do j=jm1,jm2
                0409        do i=im1,im2
4cd793223d Jean*0410         qbar(i,j) = qbar(i,j) +
                0411      &             swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
8b150b4af9 Andr*0412        enddo
                0413        enddo
                0414        enddo
                0415        do j=jm1,jm2
                0416        do i=im1,im2
                0417        tmpdiag(i,j) = qbar(i,j) *
4cd793223d Jean*0418      &             radswt(i,j,bi,bj) * pinv(i,j) * pinv(i,j) * 86400
8b150b4af9 Andr*0419        enddo
                0420        enddo
                0421        call diagnostics_fill(tmpdiag,'VDTRADSW',0,1,3,bi,bj,myid)
360ae2c8e1 Andr*0422       endif
                0423 
f0743b271c Andr*0424 c Total Precipitable Water (g/cm^2)
                0425 c ---------------------------------------------
                0426       if(diagnostics_is_on('TPW     ',myid) ) then
                0427        gravity = getcon('GRAVITY')
                0428        do j=jm1,jm2
                0429        do i=im1,im2
                0430        qbar(i,j) = 0.0
                0431        enddo
                0432        enddo
                0433        do L=1,Nrphys
                0434        do j=jm1,jm2
                0435        do i=im1,im2
4cd793223d Jean*0436        qbar(i,j) = qbar(i,j) +
                0437      &             sphy(i,j,L)*dp(i,j,L,bi,bj)
f0743b271c Andr*0438        enddo
                0439        enddo
                0440        enddo
                0441        do j=jm1,jm2
                0442        do i=im1,im2
                0443        tmpdiag(i,j) = qbar(i,j)*10. _d 0 /gravity
                0444        enddo
                0445        enddo
                0446        call diagnostics_fill(tmpdiag,'TPW     ',0,1,3,bi,bj,myid)
                0447       endif
09fc20bde8 Andr*0448 #endif
360ae2c8e1 Andr*0449       return
                0450       end