/[MITgcm]/MITgcm/model/src/ini_theta.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_theta.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.6.2.2 by cnh, Sat Jun 20 21:04:59 1998 UTC revision 1.31 by jmc, Tue Apr 4 23:22:38 2017 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_THETA
9    C     !INTERFACE:
10        SUBROUTINE INI_THETA( myThid )        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     \==========================================================/  
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 ===  C     === Global variables ===
32  #include "SIZE.h"  #include "SIZE.h"
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "PARAMS.h"  #include "PARAMS.h"
35  #include "GRID.h"  #include "GRID.h"
36  #include "DYNVARS.h"  #include "DYNVARS.h"
37    #ifdef ALLOW_MNC
38    #include "MNC_PARAMS.h"
39    #endif
40    
41    C     !INPUT/OUTPUT PARAMETERS:
42  C     == Routine arguments ==  C     == Routine arguments ==
43  C     myThid -  Number of this instance of INI_THETA  C     myThid :: Number of this instance of INI_THETA
44        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
45    
46    C     !LOCAL VARIABLES:
47  C     == Local variables ==  C     == Local variables ==
48  C     iC, jC - Center of domain  C     bi,bj  :: Tile indices
49  C     iD, jD - Disitance from domain center.  C     i,j,k  :: Loop counters
 C     rad    - Radius of initial patch  
 C     rD     - Radial displacement of point I,J  
 C     iG, jG - Global coordinate index  
 C     bi,bj  - Loop counters  
 C     I,J,K  
       INTEGER iC, jC, iD, jD  
       INTEGER iG, jG  
50        INTEGER bi, bj        INTEGER bi, bj
51        INTEGER  I,  J, K        INTEGER i, j, k, localWarnings
52        REAL rad, rD        _RL     Tfreezing
53          CHARACTER*(MAX_LEN_MBUF) msgBuf
54        _BARRIER  CEOP
55    
56        IF ( hydrogThetaFile .EQ. ' ' ) THEN  C--   Initialise temperature field to the vertical reference profile
57  C--    Example 1        DO bj = myByLo(myThid), myByHi(myThid)
58  C--    Initialise temperature field to a circular patch.         DO bi = myBxLo(myThid), myBxHi(myThid)
59         iC  = Nx/2          DO k=1,Nr
60         jC  = Ny/2           DO j=1-OLy,sNy+OLy
61         rad = MIN(Ny/8,Nx/8)            DO i=1-OLx,sNx+OLx
62         DO bj = myByLo(myThid), myByHi(myThid)             theta(i,j,k,bi,bj) = tRef(k)
         DO bi = myBxLo(myThid), myBxHi(myThid)  
          DO K=1,Nz  
           DO J=1,sNy  
            DO I=1,sNx  
             theta(I,J,K,bi,bj) = 0. _d 0  
             iG = myXGlobalLo-1+(bi-1)*sNx+I  
             jG = myYGlobalLo-1+(bj-1)*sNy+J  
             iD = iG-iC  
             jD = jG-jC  
             rD = SQRT(FLOAT(iD*iD+jD*jD))  
             IF ( rD .LT. rad ) theta(I,J,K,bi,bj) = 1. _d 0  
            ENDDO  
63            ENDDO            ENDDO
64           ENDDO           ENDDO
65          ENDDO          ENDDO
66         ENDDO         ENDDO
67  C--    Example 2        ENDDO
68  C--    Initialise temperature field to the vertical reference profile  
69         DO bj = myByLo(myThid), myByHi(myThid)        IF ( hydrogThetaFile .NE. ' ' ) THEN
70          DO bi = myBxLo(myThid), myBxHi(myThid)  #ifdef ALLOW_MNC
71           DO K=1,Nz          IF ( useMNC.AND.mnc_read_theta ) THEN
72            DO J=1,sNy            CALL MNC_FILE_CLOSE_ALL_MATCHING(hydrogThetaFile, myThid)
73             DO I=1,sNx            CALL MNC_CW_SET_UDIM(hydrogThetaFile, 1, myThid)
74              theta(I,J,K,bi,bj) = tRef(K)            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            _EXCH_XYZ_RL(theta,myThid)
85          ENDIF
86    
87    C--   Apply mask and test consistency
88          localWarnings=0
89          DO bj = myByLo(myThid), myByHi(myThid)
90           DO bi = myBxLo(myThid), myBxHi(myThid)
91            DO k=1,Nr
92             IF ( maskIniTemp ) THEN
93              DO j=1-OLy,sNy+OLy
94               DO i=1-OLx,sNx+OLx
95                IF (maskC(i,j,k,bi,bj).EQ.0.) theta(i,j,k,bi,bj) = 0.
96             ENDDO             ENDDO
97            ENDDO            ENDDO
98           ENDDO           ENDIF
99             IF ( tRef(k).NE.0. ) THEN
100              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                  localWarnings=localWarnings+1
105                ENDIF
106               ENDDO
107              ENDDO
108             ENDIF
109          ENDDO          ENDDO
110         ENDDO         ENDDO
111        ELSE        ENDDO
112         _BEGIN_MASTER( myThid )        IF (localWarnings.NE.0) THEN
113         CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )         IF ( checkIniTemp ) THEN
114         _END_MASTER(myThid)          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         &  ' set checkIniTemp=.FALSE. in "data", namelist PARM05'
121            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        ENDIF        ENDIF
131  C     Set initial tendency terms  
132  C     Note: gTNm1 is set to initial theta to make it consistent with  C--   Check that there are no values of temperature below freezing point.
133  C           the theta* scheme used.        IF ( checkIniTemp .AND. allowFreezing ) THEN
134        DO bj = myByLo(myThid), myByHi(myThid)         Tfreezing=-1.9 _d 0
135         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
136          DO K=1,Nz          DO bi = myBxLo(myThid), myBxHi(myThid)
137           DO J=1,sNy           DO k=1,Nr
138            DO I=1,sNx            DO j=1-OLy,sNy+OLy
139             gt   (I,J,K,bi,bj) = 0. _d 0             DO i=1-OLx,sNx+OLx
140             gtNM1(I,J,K,bi,bj) = 0. _d 0              IF (theta(i,j,k,bi,bj) .LT. Tfreezing) THEN
141                   theta(i,j,k,bi,bj) = Tfreezing
142                ENDIF
143               ENDDO
144            ENDDO            ENDDO
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDDO         ENDDO
148        ENDDO  c     ELSEIF ( allowFreezing ) THEN
149  C  c      CALL FREEZE_SURFACE( startTime, nIter0, myThid )
150        _EXCH_XYZ_R8(theta , myThid )        ENDIF
       _EXCH_XYZ_R8(gt , myThid )  
       _EXCH_XYZ_R8(gtNM1 , myThid )  
151    
152        CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' , Nz, 1, myThid )        IF ( plotLevel.GE.debLevC ) THEN
153            CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature',
154         &                         Nr, 1, myThid )
155          ENDIF
156    
157        RETURN        RETURN
158        END        END
Legend:
Removed from v.1.6.2.2  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22