fg/code/ini_theta.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_theta.F,v 1.29 2011/06/08 01:27:59 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-- PRM modif: start
      REAL*8   PORT_RAND
      EXTERNAL PORT_RAND
C-- PRM modif: end

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
C-- PRM modif: start
         DO j=1,sNy
          DO i=1,sNx
            theta(i,j,k,bi,bj) = theta(i,j,k,bi,bj)
     &           + ( PORT_RAND(-1.D0)-0.5 _d 0 )*20. _d -7
          ENDDO
         ENDDO
C-- PRM modif: end
        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-- PRM modif: start
      _EXCH_XYZ_RL(theta,myThid)
C-- PRM modif: end

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 ( debugLevel.GE.debLevC ) THEN
        CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
     &                         Nr, 1, myThid )
      ENDIF

      RETURN
      END
1	jmc	1.1	C $Header: /u/gcmpack/MITgcm/model/src/ini_theta.F,v 1.29 2011/06/08 01:27:59 jmc Exp $
2			C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: INI_THETA
9			C !INTERFACE:
10			SUBROUTINE INI_THETA( myThid )
11
12			C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE INI_THETA
15			C \| o Set model initial temperature field.
16			C ==========================================================
17			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			C ==========================================================
27			C \ev
28
29			C !USES:
30			IMPLICIT NONE
31			C === Global variables ===
32			#include "SIZE.h"
33			#include "EEPARAMS.h"
34			#include "PARAMS.h"
35			#include "GRID.h"
36			#include "DYNVARS.h"
37			#ifdef ALLOW_MNC
38			#include "MNC_PARAMS.h"
39			#endif
40
41			C !INPUT/OUTPUT PARAMETERS:
42			C == Routine arguments ==
43			C myThid :: Number of this instance of INI_THETA
44			INTEGER myThid
45
46			C-- PRM modif: start
47			REAL*8 PORT_RAND
48			EXTERNAL PORT_RAND
49			C-- PRM modif: end
50
51			C !LOCAL VARIABLES:
52			C == Local variables ==
53			C bi,bj :: Tile indices
54			C i,j,k :: Loop counters
55			INTEGER bi, bj
56			INTEGER i, j, k, localWarnings
57			_RL Tfreezing
58			CHARACTER*(MAX_LEN_MBUF) msgBuf
59			CEOP
60
61			C-- Initialise temperature field to the vertical reference profile
62			DO bj = myByLo(myThid), myByHi(myThid)
63			DO bi = myBxLo(myThid), myBxHi(myThid)
64			DO k=1,Nr
65			DO j=1-Oly,sNy+Oly
66			DO i=1-Olx,sNx+Olx
67			theta(i,j,k,bi,bj) = tRef(k)
68			ENDDO
69			ENDDO
70			ENDDO
71			ENDDO
72			ENDDO
73
74			IF ( hydrogThetaFile .NE. ' ' ) THEN
75			#ifdef ALLOW_MNC
76			IF ( useMNC.AND.mnc_read_theta ) THEN
77			CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
78			CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
79			CALL MNC_CW_SET_CITER(hydrogThetaFile, 2, -1, -1, -1, myThid)
80			CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
81			CALL MNC_CW_RL_R('D',hydrogThetaFile,0,0,'Temp',theta,myThid)
82			CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
83			ELSE
84			#endif /* ALLOW_MNC */
85			CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
86			#ifdef ALLOW_MNC
87			ENDIF
88			#endif /* ALLOW_MNC */
89			_EXCH_XYZ_RL(theta,myThid)
90			ENDIF
91
92			C-- Apply mask and test consistency
93			localWarnings=0
94			DO bj = myByLo(myThid), myByHi(myThid)
95			DO bi = myBxLo(myThid), myBxHi(myThid)
96			DO k=1,Nr
97			IF ( maskIniTemp ) THEN
98			DO j=1-Oly,sNy+Oly
99			DO i=1-Olx,sNx+Olx
100			IF (maskC(i,j,k,bi,bj).EQ.0.) theta(i,j,k,bi,bj) = 0.
101			ENDDO
102			ENDDO
103			ENDIF
104			IF ( tRef(k).NE.0. ) THEN
105			DO j=1,sNy
106			DO i=1,sNx
107			IF ( maskC(i,j,k,bi,bj).NE.0.
108			& .AND. theta(i,j,k,bi,bj).EQ.0. ) THEN
109			localWarnings=localWarnings+1
110			ENDIF
111			ENDDO
112			ENDDO
113			ENDIF
114			C-- PRM modif: start
115			DO j=1,sNy
116			DO i=1,sNx
117			theta(i,j,k,bi,bj) = theta(i,j,k,bi,bj)
118			& + ( PORT_RAND(-1.D0)-0.5 _d 0 )*20. _d -7
119			ENDDO
120			ENDDO
121			C-- PRM modif: end
122			ENDDO
123			ENDDO
124			ENDDO
125			IF (localWarnings.NE.0) THEN
126			IF ( checkIniTemp ) THEN
127			WRITE(msgBuf,'(A,I10,A)')
128			& ' INI_THETA: found', localWarnings,
129			& ' wet grid-pts with theta=0 identically.'
130			CALL PRINT_ERROR( msgBuf , myThid)
131			WRITE(msgBuf,'(A,A)')
132			& ' If this is intentional, you need to',
133			& ' set checkIniTemp=.false. in "data", namelist PARM05'
134			CALL PRINT_ERROR( msgBuf , myThid)
135			STOP 'ABNORMAL END: S/R INI_THETA'
136			ELSE
137			WRITE(msgBuf,'(A,I10,A)')
138			& ' WARNINGS INI_THETA: found', localWarnings,
139			& ' wet grid-pts with theta=0 identically.'
140			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
141			& SQUEEZE_RIGHT, myThid )
142			ENDIF
143			ENDIF
144			C-- PRM modif: start
145			_EXCH_XYZ_RL(theta,myThid)
146			C-- PRM modif: end
147
148			C-- Check that there are no values of temperature below freezing point.
149			Tfreezing=-1.9 _d 0
150			IF ( allowFreezing ) THEN
151			DO bj = myByLo(myThid), myByHi(myThid)
152			DO bi = myBxLo(myThid), myBxHi(myThid)
153			DO k=1,Nr
154			DO j=1-Oly,sNy+Oly
155			DO i=1-Olx,sNx+Olx
156			IF (theta(i,j,k,bi,bj) .LT. Tfreezing) THEN
157			theta(i,j,k,bi,bj) = Tfreezing
158			ENDIF
159			ENDDO
160			ENDDO
161			ENDDO
162			ENDDO
163			ENDDO
164			ENDIF
165
166			IF ( debugLevel.GE.debLevC ) THEN
167			CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
168			& Nr, 1, myThid )
169			ENDIF
170
171			RETURN
172			END