-
Notifications
You must be signed in to change notification settings - Fork 237
/
ini_sigma_hfac.F
240 lines (227 loc) · 8.23 KB
/
ini_sigma_hfac.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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
c#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: INI_SIGMA_HFAC
C !INTERFACE:
SUBROUTINE INI_SIGMA_HFAC( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INI_SIGMA_HFAC
C | o Initialise grid factors when using Sigma coordiante
C *==========================================================*
C | These arrays are used throughout the code and describe
C | fractional height factors.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: Number of this instance of INI_SIGMA_HFAC
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C bi, bj :: tile indices
C i, j, k :: Loop counters
C rEmpty :: empty column r-position
C rFullDepth :: maximum depth of a full column
C tmpFld :: Temporary array used to compute & write Total Depth
C min_hFac :: actual minimum of cell-centered hFac
C msgBuf :: Informational/error message buffer
INTEGER bi, bj
INTEGER i, j, k
_RS rEmpty
_RL rFullDepth
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL min_hFac
_RL hFactmp
CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
C r(ij,k,t) = rLow(ij) + aHybSigm(k)*[rF(1)-rF(Nr+1)]
C + bHybSigm(k)*[eta(ij,t)+Ro_surf(ij) - rLow(ij)]
IF ( usingPCoords ) rEmpty = rF(Nr+1)
IF ( usingZCoords ) rEmpty = rF(1)
rFullDepth = rF(1)-rF(Nr+1)
C--- Calculate partial-cell factor hFacC :
min_hFac = 1.
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
C- Remove column (mask=0) thinner than hFacMin*rFullDepth
C ensures hFac > hFacMin (assuming we use pure Sigma)
C Note: because of unfortunate hFacMin default value (=1) (would produce
C unexpected empty column), for now, use hFacInf instead of hFacMin
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
tmpFld(i,j) = Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj)
c IF ( tmpFld(i,j).LT.hFacMin*rFullDepth )
IF ( tmpFld(i,j).LT.hFacInf*rFullDepth )
& tmpFld(i,j) = 0. _d 0
ENDDO
ENDDO
c#ifdef ALLOW_SHELFICE
C-- Would need a specific call here similar to SHELFICE_UPDATE_MASKS
c IF ( useShelfIce ) THEN
c ENDIF
c#endif /* ALLOW_SHELFICE */
C- Set (or reset) other 2-D cell-centered fields
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF ( tmpFld(i,j).GT.0. _d 0 ) THEN
kSurfC (i,j,bi,bj) = 1
kLowC (i,j,bi,bj) = Nr
maskInC(i,j,bi,bj) = 1.
recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpFld(i,j)
ELSE
kSurfC (i,j,bi,bj) = Nr+1
kLowC (i,j,bi,bj) = 0
maskInC(i,j,bi,bj) = 0.
recip_Rcol(i,j,bi,bj) = 0. _d 0
Ro_surf(i,j,bi,bj) = rEmpty
R_low(i,j,bi,bj) = rEmpty
ENDIF
ENDDO
ENDDO
C- Set 3-D hFacC
DO k=1, Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF ( maskInC(i,j,bi,bj).NE.0. _d 0 ) THEN
hFactmp = ( dAHybSigF(k)*rFullDepth
& + dBHybSigF(k)*tmpFld(i,j)
& )*recip_drF(k)
hFacC(i,j,k,bi,bj) = hFactmp
min_hFac = MIN( min_hFac, hFactmp )
ELSE
hFacC(i,j,k,bi,bj) = 0.
ENDIF
ENDDO
ENDDO
ENDDO
C- end bi,bj loops.
ENDDO
ENDDO
WRITE(msgBuf,'(A,1PE14.6)')
& 'S/R INI_SIGMA_HFAC: minimum hFacC=', min_hFac
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
c CALL PLOT_FIELD_XYRS(R_low,
c & 'Model R_low (ini_masks_etc)', 1, myThid)
c CALL PLOT_FIELD_XYRS(Ro_surf,
c & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
C-- Set Western & Southern fields (at U and V points)
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
C- set 2-D mask and rLow & reference rSurf at Western & Southern edges
i = 1-OlX
DO j=1-OLy,sNy+OLy
rSurfW(i,j,bi,bj) = rEmpty
rLowW (i,j,bi,bj) = rEmpty
maskInW(i,j,bi,bj)= 0.
ENDDO
j = 1-OLy
DO i=1-OLx,sNx+OLx
rSurfS(i,j,bi,bj) = rEmpty
rLowS (i,j,bi,bj) = rEmpty
maskInS(i,j,bi,bj)= 0.
ENDDO
DO j=1-OLy,sNy+OLy
DO i=2-OLx,sNx+OLx
maskInW(i,j,bi,bj)= maskInC(i-1,j,bi,bj)*maskInC(i,j,bi,bj)
rSurfW(i,j,bi,bj) =
& ( Ro_surf(i-1,j,bi,bj)
& + Ro_surf( i, j,bi,bj) )*0.5 _d 0
rLowW(i,j,bi,bj) =
& ( R_low(i-1,j,bi,bj)
& + R_low( i, j,bi,bj) )*0.5 _d 0
c rSurfW(i,j,bi,bj) =
c & ( Ro_surf(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
c & + Ro_surf( i, j,bi,bj)*rA( i, j,bi,bj)
c & )*recip_rAw(i,j,bi,bj)*0.5 _d 0
c rLowW(i,j,bi,bj) =
c & ( R_low(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
c & + R_low( i, j,bi,bj)*rA( i, j,bi,bj)
c & )*recip_rAw(i,j,bi,bj)*0.5 _d 0
IF ( maskInW(i,j,bi,bj).EQ.0. ) THEN
rSurfW(i,j,bi,bj) = rEmpty
rLowW (i,j,bi,bj) = rEmpty
ENDIF
ENDDO
ENDDO
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
maskInS(i,j,bi,bj)= maskInC(i,j-1,bi,bj)*maskInC(i,j,bi,bj)
rSurfS(i,j,bi,bj) =
& ( Ro_surf(i,j-1,bi,bj)
& + Ro_surf(i, j, bi,bj) )*0.5 _d 0
rLowS(i,j,bi,bj) =
& ( R_low(i,j-1,bi,bj)
& + R_low(i, j, bi,bj) )*0.5 _d 0
c rSurfS(i,j,bi,bj) =
c & ( Ro_surf(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
c & + Ro_surf(i, j, bi,bj)*rA(i, j, bi,bj)
c & )*recip_rAs(i,j,bi,bj)*0.5 _d 0
c rLowS(i,j,bi,bj) =
c & ( R_low(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
c & + R_low(i, j, bi,bj)*rA(i, j, bi,bj)
c & )*recip_rAs(i,j,bi,bj)*0.5 _d 0
IF ( maskInS(i,j,bi,bj).EQ.0. ) THEN
rSurfS(i,j,bi,bj) = rEmpty
rLowS (i,j,bi,bj) = rEmpty
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid )
CALL EXCH_UV_XY_RS( maskInW, maskInS, .FALSE., myThid )
C- Set hFacW and hFacS (at U and V points)
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
hFactmp =
& ( dAHybSigF(k)*rFullDepth
& + dBHybSigF(k)*( rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj) )
& )*recip_drF(k)
hFacW(i,j,k,bi,bj) = hFactmp*maskInW(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
DO k=1, Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hFactmp =
& ( dAHybSigF(k)*rFullDepth
& + dBHybSigF(k)*( rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj) )
& )*recip_drF(k)
hFacS(i,j,k,bi,bj) = hFactmp*maskInS(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
C- Set surface k index for interface W & S (U & V points)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
kSurfW(i,j,bi,bj) = Nr+1
kSurfS(i,j,bi,bj) = Nr+1
IF ( maskInW(i,j,bi,bj).NE.0. ) kSurfW(i,j,bi,bj) = 1
IF ( maskInS(i,j,bi,bj).NE.0. ) kSurfS(i,j,bi,bj) = 1
ENDDO
ENDDO
C- end bi,bj loops.
ENDDO
ENDDO
C-- Additional closing of Western and Southern grid-cell edges: for example,
C a) might add some "thin walls" in specific location
C b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
C new: location now reccorded as kSurfW/S = Nr+2
CALL ADD_WALLS2MASKS( rEmpty, myThid )
RETURN
END