--- MITgcm/pkg/diagnostics/diagnostics_fill_state.F	2004/04/16 17:50:43	1.5
+++ MITgcm/pkg/diagnostics/diagnostics_fill_state.F	2004/07/06 03:55:53	1.10
@@ -1,18 +1,23 @@
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagnostics_fill_state.F,v 1.10 2004/07/06 03:55:53 edhill Exp $
+C $Name:  $
+
+#include "DIAG_OPTIONS.h"
+
        subroutine diagnostics_fill_state(myThid)
+
        implicit none
-#include "PACKAGES_CONFIG.h"
 #include "SIZE.h"
 #include "EEPARAMS.h"
-#include "CPP_OPTIONS.h"
 #include "GRID.h"
 #include "DYNVARS.h"
-# ifdef ALLOW_PTRACERS
-#  include "PTRACERS.h"
-# endif
 
       integer myThid 
-      _RL dummy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nr,Nsx,Nsy)
+
+#ifdef ALLOW_DIAGNOSTICS
+      _RL dummy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
+      _RL dummys(1-OLx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
       integer i,j,K,bi,bj
+      integer km1
 
         call fill_diagnostics(myThid,'ETAN    ',0,1,0,1,1,etaN)
 
@@ -20,12 +25,12 @@
         do bi = myBxLo(myThid), myBxHi(myThid)
           do j = 1,sNy
           do i = 1,sNx
-           dummy(i,j,1,bi,bj) = etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
+           dummys(i,j,bi,bj) = etaN(i,j,bi,bj)*etaN(i,j,bi,bj)
           enddo
           enddo
         enddo
         enddo
-        call fill_diagnostics(myThid,'ETANSQ  ',0,1,0,1,1,dummy)
+        call fill_diagnostics(myThid,'ETANSQ  ',0,1,0,1,1,dummys)
 
         call fill_diagnostics(myThid,'UVEL    ',0,Nr,0,1,1,uVel)
         call fill_diagnostics(myThid,'VVEL    ',0,Nr,0,1,1,vVel)
@@ -116,7 +121,8 @@
          do K=1,Nr
           do j = 1,sNy
           do i = 1,sNx
-           dummy(i,j,K,bi,bj) = uVel(i,j,K,bi,bj)*theta(i,j,K,bi,bj)
+           dummy(i,j,K,bi,bj) = uVel(i,j,K,bi,bj)*
+     &         0.5*(theta(i,j,K,bi,bj)+theta(i-1,j,K,bi,bj))
           enddo
           enddo
          enddo
@@ -129,7 +135,8 @@
          do K=1,Nr
           do j = 1,sNy
           do i = 1,sNx
-           dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*theta(i,j,K,bi,bj)
+           dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*
+     &         0.5*(theta(i,j,K,bi,bj)+theta(i,j-1,K,bi,bj))
           enddo
           enddo
          enddo
@@ -138,28 +145,18 @@
         call fill_diagnostics(myThid,'VVELTH  ',0,Nr,0,1,1,dummy)
 
         do bj = myByLo(myThid), myByHi(myThid)
-        do bi = myBxLo(myThid), myBxHi(myThid)
-         do K=2,Nr
-          do j = 1,sNy
-          do i = 1,sNx
-           dummy(i,j,K,bi,bj) = 0.
-          enddo
+         do bi = myBxLo(myThid), myBxHi(myThid)
+          do K=1,Nr
+           km1 = MAX(k-1,1)
+           do j = 1,sNy
+            do i = 1,sNx
+             dummy(i,j,K,bi,bj) = wVel(i,j,K,bi,bj)*0.5*
+     &                 (theta(i,j,K,bi,bj)+theta(i,j,km1,bi,bj))
+            enddo
+           enddo
           enddo
          enddo
         enddo
-        enddo
-        do bj = myByLo(myThid), myByHi(myThid)
-        do bi = myBxLo(myThid), myBxHi(myThid)
-         do K=2,Nr
-          do j = 1,sNy
-          do i = 1,sNx
-           dummy(i,j,K,bi,bj) = wVel(i,j,K,bi,bj)*0.5*
-     .                 (theta(i,j,K,bi,bj)+theta(i,j,K-1,bi,bj))
-          enddo
-          enddo
-         enddo
-        enddo
-        enddo
         call fill_diagnostics(myThid,'WVELTH  ',0,Nr,0,1,1,dummy)
 
         do bj = myByLo(myThid), myByHi(myThid)
@@ -167,7 +164,8 @@
          do K=1,Nr
           do j = 1,sNy
           do i = 1,sNx
-           dummy(i,j,K,bi,bj) = uVel(i,j,K,bi,bj)*salt(i,j,K,bi,bj)
+           dummy(i,j,K,bi,bj) = uVel(i,j,K,bi,bj)*
+     &         0.5*(salt(i,j,K,bi,bj)+salt(i-1,j,K,bi,bj))
           enddo
           enddo
          enddo
@@ -180,7 +178,8 @@
          do K=1,Nr
           do j = 1,sNy
           do i = 1,sNx
-           dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*salt(i,j,K,bi,bj)
+           dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*
+     &         0.5*(salt(i,j,K,bi,bj)+salt(i,j-1,K,bi,bj))
           enddo
           enddo
          enddo
@@ -189,17 +188,18 @@
         call fill_diagnostics(myThid,'VVELSLT ',0,Nr,0,1,1,dummy)
 
         do bj = myByLo(myThid), myByHi(myThid)
-        do bi = myBxLo(myThid), myBxHi(myThid)
-         do K=1,Nr
-          do j = 1,sNy
-          do i = 1,sNx
-           dummy(i,j,K,bi,bj) = wVel(i,j,K,bi,bj)*0.5*
-     .                 (salt(i,j,K,bi,bj)+salt(i,j,K-1,bi,bj))
-          enddo
+         do bi = myBxLo(myThid), myBxHi(myThid)
+          do K=1,Nr
+           km1 = MAX(k-1,1)
+           do j = 1,sNy
+            do i = 1,sNx
+             dummy(i,j,K,bi,bj) = wVel(i,j,K,bi,bj)*0.5*
+     &                 (salt(i,j,K,bi,bj)+salt(i,j,km1,bi,bj))
+            enddo
+           enddo
           enddo
          enddo
         enddo
-        enddo
         call fill_diagnostics(myThid,'WVELSLT ',0,Nr,0,1,1,dummy)
 
         do bj = myByLo(myThid), myByHi(myThid)
@@ -213,7 +213,7 @@
          enddo
         enddo
         enddo
-        call fill_diagnostics(myThid,'UVELMASS',0,Nr,0,1,1,salt)
+        call fill_diagnostics(myThid,'UVELMASS',0,Nr,0,1,1,dummy)
 
         do bj = myByLo(myThid), myByHi(myThid)
         do bi = myBxLo(myThid), myBxHi(myThid)
@@ -250,7 +250,7 @@
           do i = 1,sNx
            dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*
      .         0.5*(theta(i,j,K,bi,bj)+theta(i,j-1,K,bi,bj))
-     .                                * hFacW(i,j,K,bi,bj)
+     .                                * hFacS(i,j,K,bi,bj)
           enddo
           enddo
          enddo
@@ -280,7 +280,7 @@
           do i = 1,sNx
            dummy(i,j,K,bi,bj) = vVel(i,j,K,bi,bj)*
      .         0.5*(salt(i,j,K,bi,bj)+salt(i,j-1,K,bi,bj))
-     .                                * hFacW(i,j,K,bi,bj)
+     .                                * hFacS(i,j,K,bi,bj)
           enddo
           enddo
          enddo
@@ -289,5 +289,7 @@
         call fill_diagnostics(myThid,'VSLTMASS',0,Nr,0,1,1,dummy)
 
 
+#endif /* ALLOW_DIAGNOSTICS */
+
       return
       end