Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 650)
+++ /trunk/LMDZ.GENERIC/README	(revision 651)
@@ -667,2 +667,3 @@
 - Added consistency check in inifis
 - Moved 1d water initialization from physiqu to rcm1d
+- All enertests in physiq written in a matricial (F90) way. The rest of physiqu should follow soon
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 650)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 651)
@@ -314,5 +314,6 @@
 
 !     to test energy conservation (RW+JL)
-      real dEtot, dEtots, masse, AtmToSurf_TurbFlux
+      real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
+      real dEtot, dEtots, AtmToSurf_TurbFlux
       real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
       real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
@@ -626,4 +627,6 @@
             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
+	    massarea(ig,l)=mass(ig,l)*area(ig)
          enddo
       enddo
@@ -804,25 +807,10 @@
 ! test energy conservation
          if(enertest)then
-            dEtotSW  = 0.0
-            dEtotLW  = 0.0
-            dEtotsSW = 0.0
-            dEtotsLW = 0.0
-	    dEzRadsw(:,:)=0.
-	    dEzRadlw(:,:)=0.
-            do ig = 1, ngrid
-               do l = 1, nlayer
-                  masse  = (pplev(ig,l) - pplev(ig,l+1))/g
-                  dEtotSW = dEtotSW + cpp*masse*zdtsw(ig,l)*area(ig)
-                  dEtotLW = dEtotLW + cpp*masse*zdtlw(ig,l)*area(ig)
-		  dEzRadsw(ig,l)=cpp*masse*zdtsw(ig,l)
-		  dEzRadlw(ig,l)=cpp*masse*zdtlw(ig,l)
-               enddo
-               dEtotsSW = dEtotsSW + fluxsurf_sw(ig)*(1.-albedo(ig))*area(ig)
-               dEtotsLW = dEtotsLW + emis(ig)*fluxsurf_lw(ig)*area(ig)-zplanck(ig)*area(ig)
-            enddo
-            dEtotSW  = dEtotSW/totarea
-            dEtotLW  = dEtotLW/totarea
-            dEtotsSW = dEtotsSW/totarea
-            dEtotsLW = dEtotsLW/totarea
+            dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
+            dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
+            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
+            dEtotsLW = SUM(fluxsurf_lw(:)*emis(:)*area(:)-zplanck(:))/totarea
+	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
+	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
             print*,'---------------------------------------------------------------'
             print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
@@ -917,27 +905,13 @@
          ! test energy conservation
          if(enertest)then
-            dEtot=0.0
-            dEtots=0.0
-	    dEzdiff(:,:)=0.
-            AtmToSurf_TurbFlux=0.
-            dEdiffs(:)=0.
-            dEdiff(:)=0.
-	    	    
+	    dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
             do ig = 1, ngrid
-               do l = 1, nlayer
-                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
-		  dEzdiff (ig,l)= cpp*masse*(zdtdif(ig,l))
-		  dEdiff(ig)=dEdiff(ig)+dEzdiff (ig,l)
-               enddo
+	       dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
 	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
-               dEdiff(ig)=dEdiff(ig)+ sensibFlux(ig)! subtract flux to the ground
-	       dEtot = dEtot + dEdiff(ig)*area(ig)
-	       dEdiffs(ig)=capcal(ig)*zdtsdif(ig)-zflubid(ig)-sensibFlux(ig)
-               dEtots = dEtots + dEdiffs(ig)*area(ig)
-	       AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
             enddo
-            dEtot  = dEtot/totarea
-            dEtots = dEtots/totarea
-	    AtmToSurf_TurbFlux=AtmToSurf_TurbFlux/totarea
+	    dEtot = SUM(dEdiff(:)*area(:))/totarea
+	    dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
+            dEtots = SUM(dEdiffs(:)*area(:))/totarea
+	    AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
             if (UseTurbDiff) then
 	       print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
@@ -954,43 +928,19 @@
          !-------------------------
 
-
          !-------------------------
          ! test water conservation
          if(watertest.and.water)then
-            dWtot=0.0
-            dWtots=0.0
-            nconsMAX=0.0
+            dWtot = SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap))*ptimestep/totarea
+            dWtots =  SUM(zdqsdif(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
             do ig = 1, ngrid
-
-               vdifcncons(ig)=0.0
-               do l = 1, nlayer
-                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
-
-                  iq    = igcm_h2o_vap
-                  dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
-                  vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
-
-                  iq    = igcm_h2o_ice
-                  dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
-                  vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
-               
-               enddo
-
-               iq     = igcm_h2o_vap
-               dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
-               vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
-
-               iq     = igcm_h2o_ice
-               dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
-               vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
-
-               if(vdifcncons(ig).gt.nconsMAX)then
-                  nconsMAX=vdifcncons(ig)
-               endif
-
-            enddo
-
-            dWtot  = dWtot/totarea
-            dWtots = dWtots/totarea
+	       vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
+	    Enddo
+            dWtot = dWtot + SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice))*ptimestep/totarea
+            dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
+            do ig = 1, ngrid
+	       vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
+	    Enddo	    
+            nconsMAX=MAXVAL(vdifcncons(:))
+
             print*,'---------------------------------------------------------------'
             print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
@@ -998,5 +948,4 @@
             print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
             print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
-
 
          endif
@@ -1061,12 +1010,5 @@
          ! test energy conservation
          if(enertest)then
-            dEtot=0.0
-            do ig = 1, ngrid
-               do l = 1, nlayer
-                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                  dEtot = dEtot + cpp*masse*zdtadj(ig,l)*area(ig)
-               enddo
-            enddo
-            dEtot=dEtot/totarea
+            dEtot=cpp*SUM(massarea(:,:)*zdtadj(:,:))/totarea
           print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
          endif
@@ -1076,23 +1018,16 @@
          ! test water conservation
          if(watertest)then
-            dWtot=0.0
+            dWtot = SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap))*ptimestep/totarea
             do ig = 1, ngrid
-               cadjncons(ig)=0.0
-               do l = 1, nlayer
-                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                  
-                  iq    = igcm_h2o_vap
-                  dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
-                  cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
-                  
-                  iq    = igcm_h2o_ice
-                  dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
-                  cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
-                  
-               enddo
-
-            enddo
-            dWtot=dWtot/totarea
+	       cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
+	    Enddo
+            dWtot = dWtot + SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice))*ptimestep/totarea
+            do ig = 1, ngrid
+	       cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
+	    Enddo	    
+            nconsMAX=MAXVAL(cadjncons(:))
+
             print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
+            print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
          endif
          !-------------------------
@@ -1141,15 +1076,6 @@
          ! test energy conservation
          if(enertest)then
-            dEtot=0.0
-            dEtots=0.0
-            do ig = 1, ngrid
-               do l = 1, nlayer
-                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                  dEtot = dEtot + cpp*masse*zdtc(ig,l)*area(ig)
-               enddo
-               dEtots = dEtots + capcal(ig)*zdtsurfc(ig)*area(ig)
-            enddo
-            dEtot=dEtot/totarea
-            dEtots=dEtots/totarea
+            dEtot = cpp*SUM(massarea(:,:)*zdtc(:,:))/totarea
+            dEtots = SUM(capcal(:)*zdtsurfc(:)*area(:))/totarea
             print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
             print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
@@ -1199,32 +1125,13 @@
                ! test energy conservation
                if(enertest)then
-                  dEtot=0.0
-                  madjdE(:)=0.0
+                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
                   do ig = 1, ngrid
-                     do l = 1, nlayer
-                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                        dEtot = dEtot + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
-                        madjdE(ig) = madjdE(ig) + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))
-                     enddo
+                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
                   enddo
-                  dEtot=dEtot/totarea
                   print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
-               endif
-               !-------------------------
-
-               !-------------------------
-               ! test water conservation
-               if(watertest)then
-                  dWtot=0.0
-                  do ig = 1, ngrid
-                     do iq = 1 , nq
-                        do l = 1, nlayer
-                           masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                           dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
-                           dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
-                        enddo
-                     enddo
-                  enddo
-                  dWtot=dWtot/totarea
+		  
+		! test energy conservation
+	          dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
                   print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
                endif
@@ -1264,33 +1171,14 @@
                ! test energy conservation
                if(enertest)then
-                  dEtot=0.0
-                  lscaledE(:)=0.0
+                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
                   do ig = 1, ngrid
-                     do l = 1, nlayer
-                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                        dEtot = dEtot + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
-                        lscaledE(ig) = lscaledE(ig) + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))
-                     enddo
+                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
                   enddo
-                  dEtot=dEtot/totarea
-                  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
-               endif
-               !-------------------------
-
-               !-------------------------
-               ! test water conservation
-               if(watertest)then
-                  dWtot=0.0
-                  do ig = 1, ngrid
-                     do iq = 1 , nq
-                        do l = 1, nlayer
-                           masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                           dWtot = dWtot + masse*dqvaplscale(ig,l)*area(ig)*ptimestep
-                           dWtot = dWtot + masse*dqcldlscale(ig,l)*area(ig)*ptimestep
-                        enddo
-                     enddo
-                  enddo
-                  dWtot=dWtot/totarea
-                  print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
+                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
+		  
+		! test energy conservation
+	          dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
+                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
                endif
                !-------------------------
@@ -1336,51 +1224,23 @@
 
 
-              !-------------------------
-              ! test energy conservation
-              if(enertest)then
-                 dEtot=0.0
-                 do ig = 1, ngrid
-                    do l = 1, nlayer
-                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                       dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
-                    enddo
-                 enddo
-                 dEtot=dEtot/totarea
+
+               !-------------------------
+               ! test energy conservation
+               if(enertest)then
+                  dEtot=cpp*SUM(massarea(:,:)*zdtrain(:,:))/totarea
                  print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
-
-                 dEtot=0.0
-                 do ig = 1, ngrid
-                    do l = 1, nlayer
-                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                       dItot = dItot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)
-                       dVtot = dVtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)
-                    enddo
-                    dItot = dItot + zdqssnow(ig)*area(ig)
-                    dVtot = dVtot + zdqsrain(ig)*area(ig)
-                 enddo
-                 dEtot=(dItot*RLVTT/cpp + dVtot*RLVTT/cpp)/totarea
-                 print*,'In rain dItot =',dItot*RLVTT/(cpp*totarea),' W m-2'
-                 print*,'In rain dVtot =',dVtot*RLVTT/(cpp*totarea),' W m-2'
+                  dItot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))/totarea*RLVTT/cpp
+		  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
+                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
+		  dVtot = dItot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
+		  dEtot = dItot + dVtot
+                 print*,'In rain dItot =',dItot,' W m-2'
+                 print*,'In rain dVtot =',dVtot,' W m-2'
                  print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
-              endif
-              !-------------------------
-
-              !-------------------------
-              ! test water conservation
-              if(watertest)then
-                 dWtot=0.0
-                 dWtots=0.0
-                 do ig = 1, ngrid
-                    !do iq = 1 , nq
-                       do l = 1, nlayer
-                          masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
-                          dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
-                          dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
-                       enddo
-                    !enddo
-                    dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))*area(ig)*ptimestep
-                 enddo
-                 dWtot=dWtot/totarea
-                 dWtots=dWtots/totarea
+
+               ! test water conservation
+	          dWtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))*ptimestep/totarea
+		  dWtots =  SUM((zdqsrain(:)+zdqssnow(:))*area(:))*ptimestep/totarea
                  print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
                  print*,'In rain surface water change            =',dWtots,' kg m-2'
@@ -1406,16 +1266,7 @@
            ! find qtot
            if(watertest)then
-              dWtot=0.0
-              dWtots=0.0
               iq=3
-              do ig = 1, ngrid
-                 do l = 1, nlayer
-                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
-                    dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
-                    dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
-                 enddo
-              enddo
-              dWtot=dWtot/totarea
-              dWtots=dWtots/totarea
+	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
+	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
               print*,'Before sedim pq  =',dWtot,' kg m-2'
               print*,'Before sedim pdq =',dWtots,' kg m-2'
@@ -1429,16 +1280,7 @@
            ! find qtot
            if(watertest)then
-              dWtot=0.0
-              dWtots=0.0
               iq=3
-              do ig = 1, ngrid
-                 do l = 1, nlayer
-                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
-                    dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
-                    dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
-                 enddo
-              enddo
-              dWtot=dWtot/totarea
-              dWtots=dWtots/totarea
+	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
+	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
               print*,'After sedim pq  =',dWtot,' kg m-2'
               print*,'After sedim pdq =',dWtots,' kg m-2'
@@ -1461,17 +1303,6 @@
            ! test water conservation
            if(watertest)then
-              dWtot=0.0
-              dWtots=0.0
-              do iq=1,nq
-              do ig = 1, ngrid
-                 do l = 1, nlayer
-                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
-                    dWtot = dWtot + masse*zdqsed(ig,l,iq)*area(ig)*ptimestep
-                 enddo
-                 dWtots   = dWtots + zdqssed(ig,iq)*area(ig)*ptimestep
-              enddo
-              enddo
-              dWtot=dWtot/totarea
-              dWtots=dWtots/totarea
+              dWtot = SUM(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))*ptimestep/totarea
+	      dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:))*ptimestep/totarea
               print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
               print*,'In sedim surface ice change             =',dWtots,' kg m-2'
@@ -1507,10 +1338,6 @@
             ! test energy conservation
             if(enertest)then
-               dEtot=0.0
-               do ig = 1, ngrid
-                  dEtots = dEtots + capcal(ig)*zdtsurf_hyd(ig)*area(ig)
-               enddo
-               dEtot=dEtot/totarea
-               print*,'In hydrol atmospheric energy change     =',dEtot,' W m-2'
+               dEtots = SUM(area(:)*capcal(:)*zdtsurf_hyd(:))/totarea
+               print*,'In hydrol surface energy change     =',dEtots,' W m-2'
             endif
             !-------------------------
@@ -1519,15 +1346,7 @@
             ! test water conservation
             if(watertest)then
-               dWtots=0.0
-               do ig = 1, ngrid
-                  dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)*area(ig)*ptimestep
-               enddo
-               dWtots=dWtots/totarea
+	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
                print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
-               dWtots=0.0
-               do ig = 1, ngrid
-                  dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)*area(ig)*ptimestep
-               enddo
-               dWtots=dWtots/totarea
+	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
                print*,'In hydrol surface water change          =',dWtots,' kg m-2'
                print*,'---------------------------------------------------------------'
@@ -1648,17 +1467,9 @@
 !     ---------------------------------------------------------
 
-      Ts1 = 0.0
-      Ts2 = 99999.9
-      Ts3 = 0.0
-      TsS = 0.0 ! mean temperature at bottom soil layer
-      do ig=1,ngrid
-         Ts1 = Ts1 + area(ig)*tsurf(ig)
-         Ts2 = min(Ts2,tsurf(ig))  
-         Ts3 = max(Ts3,tsurf(ig))
-         TsS = TsS + area(ig)*tsoil(ig,nsoilmx)
-      end do
-      Ts1=Ts1/totarea
-      TsS=TsS/totarea
+      Ts1 = SUM(area(:)*tsurf(:))/totarea
+      Ts2 = MINVAL(tsurf(:))
+      Ts3 = MAXVAL(tsurf(:))
       if(callsoil)then
+         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
          print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
          print*,Ts1,Ts2,Ts3,TsS
@@ -1674,16 +1485,10 @@
       if(corrk)then
 
-         ISR = 0.0
-         ASR = 0.0
-         OLR = 0.0
-         GND = 0.0
-         DYN = 0.0
-         do ig=1,ngrid 
-            ISR = ISR + area(ig)*fluxtop_dn(ig) 
-            ASR = ASR + area(ig)*fluxabs_sw(ig)
-            OLR = OLR + area(ig)*fluxtop_lw(ig) 
-            GND = GND + area(ig)*fluxgrd(ig) 
-            DYN = DYN + area(ig)*fluxdyn(ig) 
-            
+         ISR = SUM(area(:)*fluxtop_dn(:))/totarea
+         ASR = SUM(area(:)*fluxabs_sw(:))/totarea
+         OLR = SUM(area(:)*fluxtop_lw(:))/totarea
+         GND = SUM(area(:)*fluxgrd(:))/totarea
+         DYN = SUM(area(:)*fluxdyn(:))/totarea
+         do ig=1,ngrid             
             if(fluxtop_dn(ig).lt.0.0)then
                print*,'fluxtop_dn has gone crazy'
@@ -1702,10 +1507,10 @@
 
          print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
-         print*, ISR/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
+         print*, ISR,ASR,OLR,GND,DYN
          
          if(enertest)then
-            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
-            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
-            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
+            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
+            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
+            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
          endif
 
@@ -1714,5 +1519,5 @@
                ! to record global radiative balance
                open(92,file="rad_bal.out",form='formatted',position='append')
-               write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
+               write(92,*) zday,ISR,ASR,OLR
                close(92)
                open(93,file="tem_bal.out",form='formatted',position='append')
@@ -1723,4 +1528,5 @@
 
       endif
+
 
 !     ------------------------------------------------------------------
@@ -1777,25 +1583,14 @@
       if(water)then
 
-         icesrf = 0.0
-         liqsrf = 0.0
-         icecol = 0.0
-         vapcol = 0.0
-
-         h2otot = 0.0
-         do ig=1,ngrid
-
-            icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
-            liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
-            icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
-            vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)
-
-            h2otot = h2otot + area(ig)*                       &
-            (qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap)      &
-            +qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
-         end do
-
-         print*,' Total water amount [kg m^-2]: ',h2otot/totarea
+         icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
+         liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
+         icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
+         vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
+
+         h2otot = icesrf + liqsrf + icecol + vapcol
+
+         print*,' Total water amount [kg m^-2]: ',h2otot
          print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
-         print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
+         print*, icesrf,liqsrf,icecol,vapcol
 
          if(meanOLR)then
@@ -1803,5 +1598,5 @@
                ! to record global water balance
                open(98,file="h2o_bal.out",form='formatted',position='append')
-               write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
+               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
                close(98)
             endif
