model/src/ini_theta.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_theta.F,v 1.27 2009/06/17 20:00:24 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_THETA
C     !INTERFACE:
      SUBROUTINE INI_THETA( myThid )
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 ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_MNC
#include "MNC_PARAMS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid :: Number of this instance of INI_THETA
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj  :: Tile indices
C     i,j,k  :: Loop counters
      INTEGER bi, bj
      INTEGER i, j, k, localWarnings
      _RL     Tfreezing
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

C--   Initialise temperature field to the vertical reference profile
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           theta(i,j,k,bi,bj) = tRef(k)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      IF ( hydrogThetaFile .NE. ' ' ) THEN
#ifdef ALLOW_MNC
        IF ( useMNC.AND.mnc_read_theta ) THEN
          CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
          CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
          CALL MNC_CW_SET_CITER(hydrogThetaFile, 2, -1, -1, -1, myThid)
          CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
          CALL MNC_CW_RL_R('D',hydrogThetaFile,0,0,'Temp',theta,myThid)
          CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
        ELSE
#endif /*  ALLOW_MNC  */
          CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
#ifdef ALLOW_MNC
        ENDIF
#endif /*  ALLOW_MNC  */
        _EXCH_XYZ_RL(theta,myThid)
      ENDIF

C--   Apply mask and test consistency
      localWarnings=0
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         IF ( maskIniTemp ) THEN
          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
         ENDIF
         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
      IF (localWarnings.NE.0) THEN
       IF ( checkIniTemp ) THEN
        WRITE(msgBuf,'(A,I10,A)')
     &   ' INI_THETA: found', localWarnings,
     &   ' wet grid-pts with theta=0 identically.'
        CALL PRINT_ERROR( msgBuf , myThid)
        WRITE(msgBuf,'(A,A)')
     &  ' If this is intentional, you need to',
     &  ' set checkIniTemp=.false. in "data", namelist PARM05'
        CALL PRINT_ERROR( msgBuf , myThid)
        STOP 'ABNORMAL END: S/R INI_THETA'
       ELSE
        WRITE(msgBuf,'(A,I10,A)')
     &   '** WARNINGS ** INI_THETA: found', localWarnings,
     &   ' wet grid-pts with theta=0 identically.'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
       ENDIF
      ENDIF

C--   Check that there are no values of temperature below freezing point.
      Tfreezing=-1.9 _d 0
      IF ( allowFreezing ) THEN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO k=1,Nr
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            IF (theta(i,j,k,bi,bj) .LT. Tfreezing) THEN
               theta(i,j,k,bi,bj) = Tfreezing
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

      IF (debugMode) THEN
        CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
     &                         Nr, 1, myThid )
      ENDIF

      RETURN
      END
1	jmc	1.28	C $Header: /u/gcmpack/MITgcm/model/src/ini_theta.F,v 1.27 2009/06/17 20:00:24 jmc Exp $
2	adcroft	1.15	C $Name: $
3	cnh	1.1
4	jmc	1.27	#include "PACKAGES_CONFIG.h"
5	cnh	1.9	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.16	CBOP
8			C !ROUTINE: INI_THETA
9			C !INTERFACE:
10	cnh	1.1	SUBROUTINE INI_THETA( myThid )
11	cnh	1.16	C !DESCRIPTION: \bv
12			C ==========================================================
13	jmc	1.27	C \| SUBROUTINE INI_THETA
14			C \| o Set model initial temperature field.
15	cnh	1.16	C ==========================================================
16	jmc	1.27	C \| There are several options for setting the initial
17			C \| temperature file
18			C \| 1. Inline code
19			C \| 2. Vertical profile ( uniform T in X and Y )
20			C \| 3. Three-dimensional data from a file. For example from
21			C \| Levitus or from a checkpoint file from a previous
22			C \| integration.
23			C \| In addition to setting the temperature field we also
24			C \| set the initial temperature tendency term here.
25	cnh	1.16	C ==========================================================
26			C \ev
27
28			C !USES:
29	adcroft	1.11	IMPLICIT NONE
30	cnh	1.1	C === Global variables ===
31			#include "SIZE.h"
32			#include "EEPARAMS.h"
33			#include "PARAMS.h"
34			#include "GRID.h"
35			#include "DYNVARS.h"
36	edhill	1.24	#ifdef ALLOW_MNC
37			#include "MNC_PARAMS.h"
38			#endif
39	cnh	1.1
40	cnh	1.16	C !INPUT/OUTPUT PARAMETERS:
41	cnh	1.1	C == Routine arguments ==
42	jmc	1.27	C myThid :: Number of this instance of INI_THETA
43	cnh	1.1	INTEGER myThid
44
45	cnh	1.16	C !LOCAL VARIABLES:
46	cnh	1.1	C == Local variables ==
47	jmc	1.27	C bi,bj :: Tile indices
48	jmc	1.28	C i,j,k :: Loop counters
49	cnh	1.1	INTEGER bi, bj
50	jmc	1.28	INTEGER i, j, k, localWarnings
51	jmc	1.21	_RL Tfreezing
52	adcroft	1.13	CHARACTER*(MAX_LEN_MBUF) msgBuf
53	cnh	1.16	CEOP
54	cnh	1.1
55	adcroft	1.15	C-- Initialise temperature field to the vertical reference profile
56			DO bj = myByLo(myThid), myByHi(myThid)
57			DO bi = myBxLo(myThid), myBxHi(myThid)
58	jmc	1.28	DO k=1,Nr
59			DO j=1-Oly,sNy+Oly
60			DO i=1-Olx,sNx+Olx
61			theta(i,j,k,bi,bj) = tRef(k)
62	cnh	1.1	ENDDO
63			ENDDO
64			ENDDO
65			ENDDO
66	adcroft	1.15	ENDDO
67
68			IF ( hydrogThetaFile .NE. ' ' ) THEN
69	edhill	1.24	#ifdef ALLOW_MNC
70	jmc	1.27	IF ( useMNC.AND.mnc_read_theta ) THEN
71	edhill	1.24	CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
72			CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
73			CALL MNC_CW_SET_CITER(hydrogThetaFile, 2, -1, -1, -1, myThid)
74			CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
75			CALL MNC_CW_RL_R('D',hydrogThetaFile,0,0,'Temp',theta,myThid)
76			CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
77			ELSE
78			#endif /* ALLOW_MNC */
79			CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
80			#ifdef ALLOW_MNC
81			ENDIF
82			#endif /* ALLOW_MNC */
83	jmc	1.26	_EXCH_XYZ_RL(theta,myThid)
84	cnh	1.3	ENDIF
85	adcroft	1.15
86	dimitri	1.19	C-- Apply mask and test consistency
87	adcroft	1.13	localWarnings=0
88	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
89			DO bi = myBxLo(myThid), myBxHi(myThid)
90	jmc	1.28	DO k=1,Nr
91			IF ( maskIniTemp ) THEN
92			DO j=1-Oly,sNy+Oly
93			DO i=1-Olx,sNx+Olx
94			IF (maskC(i,j,k,bi,bj).EQ.0.) theta(i,j,k,bi,bj) = 0.
95			ENDDO
96	jmc	1.21	ENDDO
97	jmc	1.28	ENDIF
98	jmc	1.21	IF ( tRef(k).NE.0. ) THEN
99	jmc	1.28	DO j=1,sNy
100			DO i=1,sNx
101			IF ( maskC(i,j,k,bi,bj).NE.0.
102			& .AND. theta(i,j,k,bi,bj).EQ.0. ) THEN
103	jmc	1.21	localWarnings=localWarnings+1
104	adcroft	1.13	ENDIF
105	jmc	1.21	ENDDO
106	cnh	1.1	ENDDO
107	jmc	1.21	ENDIF
108	cnh	1.1	ENDDO
109			ENDDO
110			ENDDO
111	adcroft	1.13	IF (localWarnings.NE.0) THEN
112	jmc	1.28	IF ( checkIniTemp ) THEN
113			WRITE(msgBuf,'(A,I10,A)')
114			& ' INI_THETA: found', localWarnings,
115			& ' wet grid-pts with theta=0 identically.'
116			CALL PRINT_ERROR( msgBuf , myThid)
117			WRITE(msgBuf,'(A,A)')
118			& ' If this is intentional, you need to',
119			& ' set checkIniTemp=.false. in "data", namelist PARM05'
120			CALL PRINT_ERROR( msgBuf , myThid)
121			STOP 'ABNORMAL END: S/R INI_THETA'
122			ELSE
123			WRITE(msgBuf,'(A,I10,A)')
124			& ' WARNINGS INI_THETA: found', localWarnings,
125			& ' wet grid-pts with theta=0 identically.'
126			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
127			& SQUEEZE_RIGHT, myThid )
128			ENDIF
129	dimitri	1.19	ENDIF
130
131			C-- Check that there are no values of temperature below freezing point.
132	jmc	1.20	Tfreezing=-1.9 _d 0
133	dimitri	1.19	IF ( allowFreezing ) THEN
134			DO bj = myByLo(myThid), myByHi(myThid)
135			DO bi = myBxLo(myThid), myBxHi(myThid)
136	jmc	1.28	DO k=1,Nr
137			DO j=1-Oly,sNy+Oly
138			DO i=1-Olx,sNx+Olx
139			IF (theta(i,j,k,bi,bj) .LT. Tfreezing) THEN
140			theta(i,j,k,bi,bj) = Tfreezing
141	dimitri	1.19	ENDIF
142			ENDDO
143			ENDDO
144			ENDDO
145			ENDDO
146			ENDDO
147	adcroft	1.13	ENDIF
148	cnh	1.5
149	jmc	1.27	IF (debugMode) THEN
150			CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
151			& Nr, 1, myThid )
152			ENDIF
153	cnh	1.5
154	cnh	1.1	RETURN
155			END