--- MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F	2004/07/19 15:13:07	1.5
+++ MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F	2005/04/15 22:10:13	1.6
@@ -1,24 +1,31 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F,v 1.5 2004/07/19 15:13:07 molod Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F,v 1.6 2005/04/15 22:10:13 jmc Exp $
 C $Name:  $
 
+#include "PACKAGES_CONFIG.h"
 #include "CPP_OPTIONS.h"
 
-CStartOfInterface
+CBOP
+C     !ROUTINE: INI_THETA
+C     !INTERFACE:
       SUBROUTINE INI_THETA( myThid )
-C     /==========================================================\
-C     | SUBROUTINE INI_THETA                                     |
-C     | o Set model initial temperature field.                   |
-C     |==========================================================|
-C     | There are several options for setting the initial        |
-C     | temperature file                                         |
-C     |  1. Inline code                                          |
-C     |  2. Vertical profile ( uniform T in X and Y )            |
-C     |  3. Three-dimensional data from a file. For example from |
-C     |     Levitus or from a checkpoint file from a previous    |
-C     |     integration.                                         |
-C     | In addition to setting the temperature field we also     |
-C     | set the initial temperature tendency term here.          |
-C     \==========================================================/
+C     !DESCRIPTION: \bv
+C     *==========================================================*
+C     | SUBROUTINE INI_THETA                                      
+C     | o Set model initial temperature field.                    
+C     *==========================================================*
+C     | There are several options for setting the initial         
+C     | temperature file                                          
+C     |  1. Inline code                                           
+C     |  2. Vertical profile ( uniform T in X and Y )             
+C     |  3. Three-dimensional data from a file. For example from  
+C     |     Levitus or from a checkpoint file from a previous     
+C     |     integration.                                          
+C     | In addition to setting the temperature field we also      
+C     | set the initial temperature tendency term here.           
+C     *==========================================================*
+C     \ev
+
+C     !USES:
       IMPLICIT NONE
 
 C     === Global variables ===
@@ -28,21 +35,24 @@
 #include "GRID.h"
 #include "DYNVARS.h"
 
+C     !INPUT/OUTPUT PARAMETERS:
 C     == Routine arguments ==
 C     myThid -  Number of this instance of INI_THETA
       INTEGER myThid
-CEndOfInterface
 
 C     == Functions ==
-      Real*8  PORT_RAND
-      real *8 seed
+      real*8  PORT_RAND
+      real*8  seed
 
+C     !LOCAL VARIABLES:
 C     == Local variables ==
 C     bi,bj  - Loop counters
 C     I,J,K
       INTEGER bi, bj
-      INTEGER I, J, K
+      INTEGER I, J, K, localWarnings
       _RL     term1,term2,thetaLim,thetaEq
+      CHARACTER*(MAX_LEN_MBUF) msgBuf
+CEOP
 
       _BARRIER
 
@@ -51,17 +61,17 @@
 c     seed = j
 
       IF ( hydrogThetaFile .EQ. ' ' ) THEN
-C--    Initialise temperature field to Held & Suarez equilibrium theta
+C--    Initialise temperature field to Held & Saurez equilibrium theta
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO K=1,Nr
           thetaLim = 200. _d 0/((rC(K)/atm_po)**atm_kappa)
           DO J=1,sNy
            DO I=1,sNx
-           term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)
-           term2=10. _d 0*log((rC(K)/atm_po))
+            term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)
+            term2=10. _d 0*log((rC(K)/atm_po))
      &              *(cos(yC(I,J,bi,bj)*deg2rad)**2)
-           thetaEq=315. _d 0-term1-term2
+            thetaEq=315. _d 0-term1-term2
             theta(I,J,K,bi,bj) = MAX( thetaLim, thetaEq )
 c    &                          + 0.01*(RAND()-0.5)
 c    &                          + 0.01*(PORT_RAND(seed)-0.5)
@@ -69,49 +79,55 @@
            ENDDO
           ENDDO
          ENDDO
-        ENDDO
-       ENDDO
 #ifdef ALLOW_ZONAL_FILT
 C--   Zonal FFT filter initial conditions
-       DO bj = myByLo(myThid), myByHi(myThid)
-        DO bi = myBxLo(myThid), myBxHi(myThid)
-         DO K=1,Nr
-          DO J=1,sNy
-           CALL ZONAL_FILTER(
+         IF (useZONAL_FILT) THEN
+          CALL ZONAL_FILTER(
      U      theta, hFacC,
-     I      1, sNy, k, k, bi, bj, 1, myThid)
-          ENDDO
-         ENDDO
+     I      1, sNy, 1, Nr, bi, bj, 1, myThid)
+         ENDIF
+#endif /* ALLOW_ZONAL_FILT */
         ENDDO
        ENDDO
-#endif /* INCLUDE_LAT_CIRC_FFT_FILTER_CODE */
       ELSE
        _BEGIN_MASTER( myThid )
        CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
        _END_MASTER(myThid)
       ENDIF
-C     Set initial tendency terms
+C--   Apply mask and test consistency
+      localWarnings=0
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO K=1,Nr
-         DO J=1,sNy
-          DO I=1,sNx
-           gt   (I,J,K,bi,bj) = 0. _d 0
-           gtNM1(I,J,K,bi,bj) = 0. _d 0
-           IF (hFacC(I,J,K,bi,bj).EQ.0.) theta(I,J,K,bi,bj) = 0.
+         DO J=1-Oly,sNy+Oly
+          DO I=1-Olx,sNx+Olx
+           IF (maskC(I,J,K,bi,bj).EQ.0.) theta(I,J,K,bi,bj) = 0.
           ENDDO
          ENDDO
+         IF ( tRef(k).NE.0. ) THEN
+          DO J=1,sNy
+           DO I=1,sNx
+            IF (  maskC(I,J,K,bi,bj).NE.0.
+     &      .AND. theta(I,J,K,bi,bj).EQ.0. ) THEN
+              localWarnings=localWarnings+1
+            ENDIF
+           ENDDO
+          ENDDO
+         ENDIF
         ENDDO
        ENDDO
       ENDDO
-C
-
+      IF (localWarnings.NE.0) THEN
+       WRITE(msgBuf,'(A,A)')
+     &  'S/R INI_THETA: theta = 0 identically. If this is intentional',
+     &  'you will need to edit ini_theta.F to avoid this safety check'
+       CALL PRINT_ERROR( msgBuf , myThid)
+       STOP 'ABNORMAL END: S/R INI_THETA'
+      ENDIF
 
       _EXCH_XYZ_R8(theta , myThid )
-      _EXCH_XYZ_R8(gt , myThid )
-      _EXCH_XYZ_R8(gtNM1 , myThid )
 
-      CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' , 
+      CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' ,
      &                       Nr, 1, myThid )
 
       RETURN