| blob | 49e9860f8ba75720ce5dcc5bb5045eb893864574 |
1 SUBROUTINE SETUP SETUP.1
2 C SETUP.2
3 C PURPOSE : INITIALIZES THE CONSTANTS AND VARIABLES WHICH CAN BE SETUP.3
4 C DETERMINED BASED ON THE INPUT NAMELIST VARIABLES. SETUP.4
5 C CALLS : THE INTRINSIC FUNCTIONS SIN, COS, TAN, ALOG, AND ALOG10. SETUP.5
6 C CALLED BY :THE MAIN PROGRAM TERRAIN. SETUP.6
7 C COMMENTS : SETUP.7
8 C THE COMMON BLOCKS /MAPS/, /NESTDMN/, /FUDG/, /TRFUDGE/, /OPTION/, SETUP.8
9 C /LTDATA/, AND /HEADER/ ARE USED IN THIS SUBROUTINE. THERE ARE FIVE SETUP.9
10 C NAMELISTS /MAPBG/, /DOMAIN/, /OPTN/, /FUDGE/, AND /FUDGET/ IN THIS SETUP.10
11 C SUBROUTINE. THIS SUBROUTINE HAS NO DUMMY ARGUMENTS. PLEASE SEE SETUP.11
12 C SECTION 5.1.1 AND 5.2.1 FOR DETAILS. SETUP.12
13 C SETUP.13
14 # include <maps.incl>
15 # include <nestdmn.incl>
16 # include <namelist1.incl>
17 C
18 REAL TRUELAT0
20 DATA CONV/57.29578/, A/6370./ SETUP.56
21 DATA NPROJ/'LAMCON','POLSTR','MERCAT'/ SETUP.57
22 C
23 PRINT *,' ==> CALL SETUP <==' SETUP.60
24 C SETUP.61
26 C ==> READ "MAPBG"
27 READ(15,MAPBG,ERR=55,END=55)
28 c .. Namelist MAPBG available:
29 GO TO 56
30 c .. Namelist MAPBG unavailable:
31 55 rewind (15)
33 C ==> READ "DOMAIN"
34 56 READ(15,DOMAINS,ERR=66,END=66)
35 c .. Namelist DOMAINS available:
36 GO TO 67
37 c .. Namelist DOMAINS unavailable:
38 66 rewind (15)
40 67 CONTINUE
41 C PRINT *, MAPBG
42 C PRINT *, DOMAINS
44 if (iproj.eq.'CYLEQU') then
45 PROJECT='CE'
46 goto 1002
47 endif
48 C SETUP.146
49 C DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ: SETUP.147
50 C SETUP.148
51 XN = -1.0E36 SETUP.150
52 IF(PHIC.LT.0) THEN SETUP.151
53 SIGN=-1. ! SOUTH HEMESPHERE SETUP.152
54 ELSE SETUP.153
55 SIGN=1. ! NORTH HEMESPHERE SETUP.154
56 ENDIF SETUP.155
57 POLE = 90. SETUP.156
59 IF (IPROJ.EQ.NPROJ(1)) THEN SETUP.157
61 print '("Map projection=",a)',NPROJ(1)
62 IF(ABS(TRUELAT1).GT.90.) THEN SETUP.158
63 TRUELAT1=60. SETUP.159
64 TRUELAT2=30. SETUP.160
65 TRUELAT1=SIGN*TRUELAT1 SETUP.161
66 TRUELAT2=SIGN*TRUELAT2 SETUP.162
67 ENDIF SETUP.163
69 ! In case of TRUELAT1 == TRUELAT2 (YRG 09/22/2005):
71 if (abs(TRUELAT1-TRUELAT2)==0.0) then
72 TRUELAT0 = TRUELAT1
73 TRUELAT1 = TRUELAT0 - 1.e-2
74 TRUELAT2 = TRUELAT0 + 1.e-2
75 endif
77 XN = ALOG10(COS(TRUELAT1 / CONV)) - SETUP.164
78 * ALOG10(COS(TRUELAT2 / CONV)) SETUP.165
79 XN = XN/(ALOG10(TAN((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - SETUP.166
80 * ALOG10(TAN((45.0 - SIGN*TRUELAT2/2.0) / CONV))) SETUP.167
82 c TRUELAT1 = TRUELAT0
83 c TRUELAT2 = TRUELAT0
85 PSI1=90.-SIGN*TRUELAT1 SETUP.168
86 PROJECT='LC' SETUP.169
88 ELSE IF (IPROJ.EQ.NPROJ(2)) THEN SETUP.170
90 print '("Map projection=",a)',NPROJ(2)
91 XN = 1.0 SETUP.171
92 IF(ABS(TRUELAT1).GT.90.) THEN SETUP.172
93 TRUELAT1=60. SETUP.173
94 TRUELAT2=0. SETUP.174
95 TRUELAT1=SIGN*TRUELAT1 SETUP.175
96 TRUELAT2=SIGN*TRUELAT2 SETUP.176
97 ENDIF SETUP.177
98 C PSI1 IS THE PSEUDO-LATITUDE SETUP.178
99 PSI1=90.-SIGN*TRUELAT1 SETUP.179
100 PROJECT='ST' SETUP.180
102 ELSE IF (IPROJ.EQ.NPROJ(3)) THEN SETUP.181
104 XN = 0. SETUP.182
105 IF(ABS(TRUELAT1).GT.90.) THEN SETUP.183
106 TRUELAT1=0. SETUP.184
107 TRUELAT2=0. SETUP.185
108 ENDIF SETUP.186
110 C -------------------------------------------------------------------
111 C For Map plotting on Mercator projection, adjusted the
112 C TRUELAT1 and TRUELAT2 to zero by increasing the grid
113 C distance because MAPDRV can only support the TRUELAT=0.0
114 C Y.R. Guo 09/23/2005
116 if (TRUELAT1 .NE. TRUELAT2) THEN
117 TRUELAT1 = max(TRUELAT1,TRUELAT2)
118 TRUELAT2 = max(TRUELAT1,TRUELAT2)
119 endif
120 print '(/"MERCATOR MAP PROJECTION=",a,a,2F8.2)',
121 - IPROJ, " TRUELAT 1,2=",TRUELAT1, TRUELAT2
123 IF(TRUELAT1.NE.0.) THEN SETUP.187
125 do nm = 1, maxnes
126 print '(I2," Original dis=",f10.3)', nm, dis(nm)
127 dis(nm) = dis(nm) / COS(TRUELAT1*3.1415926/180.)
128 print '(" Adjusted DIS=",F10.3)', dis(nm)
129 enddo
130 TRUELAT1=0. SETUP.190
131 TRUELAT2=0. SETUP.190
132 c print '("TRUELAT1 and 2:",2F8.2/)', TRUELAT1, TRUELAT2
134 c PRINT *,'TERRAIN AND MAPDRV ONLY SUPPORT MERCATOR PROJECTION SETUP.188
135 c *AT 0 DEGREE TURE LATITUDE, RESET TRUELAT1=0.' SETUP.189
136 C ---------------------------------------------------------------------
138 ENDIF
139 PSI1=0. SETUP.192
140 PROJECT='ME' SETUP.193
141 END IF SETUP.194
142 C SETUP.195
143 PSI1 = PSI1/CONV SETUP.196
144 IF (PHIC.LT.0.) THEN SETUP.197
145 PSI1 = -PSI1 SETUP.198
146 POLE = -POLE SETUP.199
147 ENDIF SETUP.200
148 C SETUP.201
149 C--------CALCULATE R SETUP.202
150 IF (IPROJ.NE.NPROJ(3)) THEN SETUP.203
151 PSX = (POLE-PHIC)/CONV SETUP.204
152 IF (IPROJ.EQ.NPROJ(1)) THEN SETUP.205
153 CELL = A*SIN(PSI1)/XN SETUP.206
154 CELL2 = (TAN(PSX/2.))/(TAN(PSI1/2.)) SETUP.207
155 ENDIF SETUP.208
156 IF (IPROJ.EQ.NPROJ(2)) THEN SETUP.209
157 CELL = A*SIN(PSX)/XN SETUP.210
158 CELL2 = (1. + COS(PSI1))/(1. + COS(PSX)) SETUP.211
159 ENDIF SETUP.212
160 R = CELL*(CELL2)**XN SETUP.213
161 XCNTR = 0.0 SETUP.214
162 YCNTR = -R SETUP.215
163 ENDIF SETUP.216
164 C-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1 SETUP.217
165 IF (IPROJ.EQ.NPROJ(3)) THEN SETUP.218
166 C2 = A*COS(PSI1) SETUP.219
167 XCNTR = 0.0 SETUP.220
168 PHICTR = PHIC/CONV SETUP.221
169 CELL = COS(PHICTR)/(1.0+SIN(PHICTR)) SETUP.222
170 YCNTR = - C2*ALOG(CELL) SETUP.223
171 ENDIF SETUP.224
172 PRINT 10,XCNTR,YCNTR SETUP.226
173 10 FORMAT(1X,'COARSE GRID CENTER IS AT X =',F8.1,' KM AND Y = ', SETUP.227
174 1 F8.1,' KM ') SETUP.228
175 C SETUP.230
177 1002 continue
179 C CHECK THE COMPATIBILITY OF NEST DOMAINS WITH THE COARSE DOMAINS SETUP.231
180 C AND CALCULATE THE IRATIOS, INORTHS, ISOUTHS, IWESTS AND IEASTS SETUP.232
181 C SETUP.233
182 C A) EXTENDING THE COARSE DOMAIN IF IEXP = 1 SETUP.234
183 C SETUP.235
184 IXEX = NESTIX(1) SETUP.236
185 JXEX = NESTJX(1) SETUP.237
186 IOFFST = 0 SETUP.238
187 JOFFST = 0 SETUP.239
188 IF (IEXP.EQ.1) THEN SETUP.240
189 INCR = INT(AEXP/DIS(1) + 1.001) SETUP.241
190 IXEX = NESTIX(1) + INCR*2 SETUP.242
191 JXEX = NESTJX(1) + INCR*2 SETUP.243
192 IOFFST = INCR SETUP.244
193 JOFFST = INCR SETUP.245
194 C SETUP.252
195 PRINT 20,INCR SETUP.253
196 PRINT 22,IXEX,JXEX SETUP.254
197 20 FORMAT(2X,'$$$$ GRID IS EXPANDED BY ',I3, SETUP.255
198 1 ' GRID POINTS ON EACH EDGE $$$$') SETUP.256
199 22 FORMAT(5X,'EXTENDED IXEX IS ',I3,5X,'EXTENDED JXEX IS ',I3) SETUP.257
200 ENDIF SETUP.258
201 C-----CENTER OF GRID IN THE COARSE DOMAIN SETUP.259
202 CNTRJ0 = FLOAT(JXEX+1)/2. SETUP.260
203 CNTRI0 = FLOAT(IXEX+1)/2. SETUP.261
204 PRINT 21,NESTIX(1),NESTJX(1) SETUP.262
205 21 FORMAT(2X,'COARSE DOMAIN SIZE (IX,JX) =',2I5) SETUP.263
206 C SETUP.264
207 C MIX, MJX ARE USED IN SUB. TFUDGE: SETUP.265
208 MIX = IXEX SETUP.266
209 MJX = JXEX SETUP.267
210 DO 23 M = 1, MAXNES SETUP.268
211 MIX = MAX0(NESTIX(M),MIX) SETUP.269
212 MJX = MAX0(NESTJX(M),MJX) SETUP.270
213 23 CONTINUE SETUP.271
214 PRINT 24, MIX,MJX SETUP.275
215 24 FORMAT(' ++> THE MAXIMUM DIMENSION = (',2I5,') <++') SETUP.276
216 C SETUP.277
217 C CHECK IF POLE IS INSIDE THE DOMAIN OR NOT FOR LAMBERT PROJECTION: SETUP.278
218 HDSIZE = (IXEX-1)*DIS(1)/2. SETUP.279
219 IF (HDSIZE.GT.ABS(YCNTR) .AND. IPROJ.EQ.NPROJ(1))THEN SETUP.280
220 PRINT *,'-------------------------------------------------' SETUP.281
221 PRINT *,'HALF DOMAIN SIZE IN Y-DIRECTION = ',HDSIZE SETUP.282
222 PRINT *,' DISTANCE FROM CNTER TO POLE = ',ABS(YCNTR) SETUP.283
223 PRINT *,'NOT MAKE SENSE WITH THE POLE IS INSIDE THE DOMAIN ' SETUP.284
224 PRINT *,' FOR LAMBERT CONFORMAL PROJECTION!' SETUP.285
225 PRINT *,'=== PLEASE RE-SPECIFY THE CENTER OR DOMAIN SIZE. ===' SETUP.286
226 STOP SETUP.287
227 ENDIF SETUP.288
228 C SETUP.289
229 C B) CALCULATING THE IRATIOS, INORTHS AND IEASTS: SETUP.290
230 C SETUP.291
231 IRATIO(1) = 1 SETUP.292
232 NRATIO(1) = 1 SETUP.293
233 XSOUTH(1) = 1. SETUP.294
234 XWEST(1) = 1. SETUP.295
235 XNORTH(1) = FLOAT(IXEX) SETUP.296
236 XEAST(1) = FLOAT(JXEX) SETUP.297
237 XJC = (XEAST(1) + 1.0)/2. SETUP.298
238 C
239 PRINT 27,XSOUTH(1),XWEST(1),XNORTH(1),XEAST(1),DIS(1), SETUP.300
240 1 IRATIO(1),NRATIO(1) SETUP.301
241 27 FORMAT(1X,'XSOUTH(1)= ',F6.1,2X,'XWEST(1)= ',F6.1,2X, SETUP.302
242 1'XNORTH(1)= ',F6.1,2X,'XEAST(1)= ',F6.1,2X,'DIS(1)= ',F6.1, SETUP.303
243 22X,'IRATIO(1)= ',I3,2X,'NRATIO(1)= ',I3) SETUP.304
244 C SETUP.306
245 MISMATCH = 0 SETUP.307
246 SETUP.308
247 DO 30 NM = 2, MAXNES SETUP.309
248 C SETUP.310
249 C SETUP.311
250 C DOMAINS' CONSISTENCY CHECK: SETUP.312
251 C SETUP.313
252 NMC = NUMC(NM) SETUP.314
253 IF (AMOD((DIS(NMC)+0.09),DIS(NM)) .GT. 0.1) THEN SETUP.315
254 MISMATCH = MISMATCH + 1 SETUP.316
255 PRINT 29, NM,NMC SETUP.317
256 PRINT 31,NM,DIS(NM),NMC,DIS(NMC) SETUP.318
257 29 FORMAT(2X,'DOMAIN ',I2,' HAS INCORRECT GRID SIZE ', SETUP.319
258 1' IT HAS TO BE THE MULTIPLE OF DOMAIN ',I2) SETUP.320
259 31 FORMAT(2X,'DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM', SETUP.321
260 1 ' DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM') SETUP.322
261 GO TO 30 SETUP.323
262 ENDIF SETUP.324
263 IRATIO(NM) = NINT(DIS(NMC)/DIS(NM)) SETUP.325
264 NRATIO(NM) = NINT(DIS(1)/DIS(NM)) SETUP.326
265 C SETUP.327
266 C MAKE SURE THE 4 CORNER POINTS OF NEST DOMAINS ARE ON THE SETUP.328
267 C PREVIOUS DOMAIN GRID-POINTS SETUP.329
268 C SETUP.330
269 IF (MOD((NESTIX(NM)-1),IRATIO(NM)).NE.0) THEN SETUP.331
270 MISMATCH = MISMATCH + 1 SETUP.332
271 IIMXN = (INT(FLOAT(NESTIX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1 SETUP.333
272 PRINT 32,NM,NESTIX(NM),IRATIO(NM),IIMXN SETUP.334
273 32 FORMAT(2X,'NESTIX(',I2,')=',I4,' AND IRATIO=',I2,' DOES NOT MATC SETUP.335
274 1H, YOU MAY SET NESTIX TO ',I4) SETUP.336
275 ENDIF SETUP.337
276 IF (MOD((NESTJX(NM)-1),IRATIO(NM)).NE.0) THEN SETUP.338
277 MISMATCH = MISMATCH + 1 SETUP.339
278 JJMXN = (INT(FLOAT(NESTJX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1 SETUP.340
279 PRINT 33,NM,NESTJX(NM),IRATIO(NM),JJMXN SETUP.341
280 33 FORMAT(2X,'NESTJX(',I2,')=',I4,' AND IRATIO=',I2,' DOES NOT MATC SETUP.342
281 1H, YOU MAY SET NESTJX TO ',I4) SETUP.343
282 ENDIF SETUP.344
283 C SETUP.345
284 C--------REDEFINE LOCATION OF LOWER LEFT CORNER OF FINE MESH (IN TERMS SETUP.346
285 C OF EXTENDED COARSE MESH - DOMAIN 1 INDICES) IF USING EXTENDED G SETUP.347
286 SETUP.348
287 PRINT 34,NM,NESTIX(NM),NESTJX(NM),DIS(NM),NESTI(NM),NESTJ(NM) SETUP.349
288 1, NUMC(NM),IRATIO(NM),NRATIO(NM),IEXP SETUP.350
289 34 FORMAT(/1X,'DOMAIN ',I2,2X,'IX=',I3,2X,'JX=',I3,2X, SETUP.351
290 1 'DS= ',F6.1,2X,'ICNS=',I4,2X,'JCNS=',I4,2X, SETUP.352
291 2 'NUMC= ',I2,2X,'IRATIO= ',I2,2X,'NRATIO= ',I2,2X, SETUP.353
292 3 'IEXP= ',I2) SETUP.354
293 C SETUP.355
294 XDIS = 0.0 SETUP.356
295 YDIS = 0.0 SETUP.357
296 ND1 = NM SETUP.358
297 ND2 = NMC SETUP.359
298 40 CONTINUE SETUP.360
299 XDIS = (NESTI(ND1)-1)*DIS(ND2) + XDIS SETUP.361
300 YDIS = (NESTJ(ND1)-1)*DIS(ND2) + YDIS SETUP.362
301 IF (ND2 .GT. 1) THEN SETUP.363
302 ND1 = ND2 SETUP.364
303 ND2 = NUMC(ND2) SETUP.365
304 GO TO 40 SETUP.366
305 ENDIF SETUP.367
306 C SETUP.368
307 XSOUTH(NM) = XDIS/DIS(1) + FLOAT(IOFFST) + 1 SETUP.369
308 XWEST(NM) = YDIS/DIS(1) + FLOAT(JOFFST) + 1 SETUP.370
309 XNORTH(NM) = XSOUTH(NM) + FLOAT(NESTIX(NM)-1)*DIS(NM)/DIS(1) SETUP.371
310 XEAST(NM) = XWEST(NM) + FLOAT(NESTJX(NM)-1)*DIS(NM)/DIS(1) SETUP.372
311 SETUP.373
312 PRINT 35 SETUP.374
313 PRINT 36,XSOUTH(NM),XWEST(NM),XNORTH(NM),XEAST(NM) SETUP.375
314 35 FORMAT(2X,'COARSE MESH INDICES FOR THE 4 CORNER POINTS ARE') SETUP.376
315 36 FORMAT(2X,'SOUTH=',F6.1,3X,'WEST=',F6.1,3X,'NORTH =',F6.1,3X, SETUP.377
316 1 'EAST = ',F6.1) SETUP.378
317 C SETUP.379
318 30 CONTINUE SETUP.380
319 C SETUP.381
320 IF (MISMATCH.GT.0) THEN SETUP.382
321 PRINT *, 'TERRAIN STOP IN SUBROUTINE SETUP DUE TO INCORRECT NEST SETUP.383
322 1 DOMAIN SET UP' SETUP.384
323 cc STOP 1111 SETUP.385
324 ENDIF SETUP.386
325 C SETUP.387
326 RETURN SETUP.440
327 END SETUP.441