model/src/ini_theta.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_theta.F,v 1.30 2014/08/26 06:15:06 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.
      IF ( checkIniTemp .AND. allowFreezing ) THEN
       Tfreezing=-1.9 _d 0
       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
c     ELSEIF ( allowFreezing ) THEN
c      CALL FREEZE_SURFACE( startTime, nIter0, myThid )
      ENDIF

      IF ( plotLevel.GE.debLevC ) THEN
        CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
     &                         Nr, 1, myThid )
      ENDIF

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