-
Notifications
You must be signed in to change notification settings - Fork 237
/
ini_theta.F
130 lines (119 loc) · 3.89 KB
/
ini_theta.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#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 ==
c real*8 PORT_RAND
c 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
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(1-OLx,1-OLy,1,bi,bj),
I hFacC(1-OLx,1-OLy,1,bi,bj),
I 1, sNy, Nr, bi, bj, 1, myThid )
ENDIF
#endif /* ALLOW_ZONAL_FILT */
ENDDO
ENDDO
ELSE
CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, 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_RL(theta , myThid )
IF (debugMode) THEN
CALL PLOT_FIELD_XYZRL( theta, 'Initial Temperature' ,
& Nr, 1, myThid )
ENDIF
RETURN
END