hs94.cs-32x32x5/code/ini_theta.F

C $Header: /u/gcmpack/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"

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"

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

C     == Functions ==
      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, localWarnings
      _RL     term1,term2,thetaLim,thetaEq
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      _BARRIER

      J = 99+myBxLo(myThid)+nPx*myByLo(myThid)
c     CALL SRAND( J )
c     seed = j

      IF ( hydrogThetaFile .EQ. ' ' ) THEN
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))
     &              *(cos(yC(I,J,bi,bj)*deg2rad)**2)
            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)
c           theta(I,J,K,bi,bj) = tRef(K)
           ENDDO
          ENDDO
         ENDDO
#ifdef ALLOW_ZONAL_FILT
C--   Zonal FFT filter initial conditions
         IF (useZONAL_FILT) THEN
          CALL ZONAL_FILTER(
     U      theta, hFacC,
     I      1, sNy, 1, Nr, bi, bj, 1, myThid)
         ENDIF
#endif /* ALLOW_ZONAL_FILT */
        ENDDO
       ENDDO
      ELSE
       _BEGIN_MASTER( myThid )
       CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
       _END_MASTER(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
         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
      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 )

      CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' ,
     &                       Nr, 1, myThid )

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F,v 1.6 2005/04/15 22:10:13 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	C !DESCRIPTION: \bv
12	C ==========================================================
13	C \| SUBROUTINE INI_THETA
14	C \| o Set model initial temperature field.
15	C ==========================================================
16	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	C ==========================================================
26	C \ev
27
28	C !USES:
29	IMPLICIT NONE
30
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
38	C !INPUT/OUTPUT PARAMETERS:
39	C == Routine arguments ==
40	C myThid - Number of this instance of INI_THETA
41	INTEGER myThid
42
43	C == Functions ==
44	real*8 PORT_RAND
45	real*8 seed
46
47	C !LOCAL VARIABLES:
48	C == Local variables ==
49	C bi,bj - Loop counters
50	C I,J,K
51	INTEGER bi, bj
52	INTEGER I, J, K, localWarnings
53	_RL term1,term2,thetaLim,thetaEq
54	CHARACTER*(MAX_LEN_MBUF) msgBuf
55	CEOP
56
57	_BARRIER
58
59	J = 99+myBxLo(myThid)+nPx*myByLo(myThid)
60	c CALL SRAND( J )
61	c seed = j
62
63	IF ( hydrogThetaFile .EQ. ' ' ) THEN
64	C-- Initialise temperature field to Held & Saurez equilibrium theta
65	DO bj = myByLo(myThid), myByHi(myThid)
66	DO bi = myBxLo(myThid), myBxHi(myThid)
67	DO K=1,Nr
68	thetaLim = 200. _d 0/((rC(K)/atm_po)**atm_kappa)
69	DO J=1,sNy
70	DO I=1,sNx
71	term1=60. _d 0(sin(yC(I,J,bi,bj)deg2rad)**2)
72	term2=10. _d 0*log((rC(K)/atm_po))
73	& (cos(yC(I,J,bi,bj)deg2rad)**2)
74	thetaEq=315. _d 0-term1-term2
75	theta(I,J,K,bi,bj) = MAX( thetaLim, thetaEq )
76	c & + 0.01*(RAND()-0.5)
77	c & + 0.01*(PORT_RAND(seed)-0.5)
78	c theta(I,J,K,bi,bj) = tRef(K)
79	ENDDO
80	ENDDO
81	ENDDO
82	#ifdef ALLOW_ZONAL_FILT
83	C-- Zonal FFT filter initial conditions
84	IF (useZONAL_FILT) THEN
85	CALL ZONAL_FILTER(
86	U theta, hFacC,
87	I 1, sNy, 1, Nr, bi, bj, 1, myThid)
88	ENDIF
89	#endif /* ALLOW_ZONAL_FILT */
90	ENDDO
91	ENDDO
92	ELSE
93	_BEGIN_MASTER( myThid )
94	CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
95	_END_MASTER(myThid)
96	ENDIF
97	C-- Apply mask and test consistency
98	localWarnings=0
99	DO bj = myByLo(myThid), myByHi(myThid)
100	DO bi = myBxLo(myThid), myBxHi(myThid)
101	DO K=1,Nr
102	DO J=1-Oly,sNy+Oly
103	DO I=1-Olx,sNx+Olx
104	IF (maskC(I,J,K,bi,bj).EQ.0.) theta(I,J,K,bi,bj) = 0.
105	ENDDO
106	ENDDO
107	IF ( tRef(k).NE.0. ) THEN
108	DO J=1,sNy
109	DO I=1,sNx
110	IF ( maskC(I,J,K,bi,bj).NE.0.
111	& .AND. theta(I,J,K,bi,bj).EQ.0. ) THEN
112	localWarnings=localWarnings+1
113	ENDIF
114	ENDDO
115	ENDDO
116	ENDIF
117	ENDDO
118	ENDDO
119	ENDDO
120	IF (localWarnings.NE.0) THEN
121	WRITE(msgBuf,'(A,A)')
122	& 'S/R INI_THETA: theta = 0 identically. If this is intentional',
123	& 'you will need to edit ini_theta.F to avoid this safety check'
124	CALL PRINT_ERROR( msgBuf , myThid)
125	STOP 'ABNORMAL END: S/R INI_THETA'
126	ENDIF
127
128	_EXCH_XYZ_R8(theta , myThid )
129
130	CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' ,
131	& Nr, 1, myThid )
132
133	RETURN
134	END