NMM_NEST_UTILS1.F
References to this file elsewhere.
1 #if (NMM_NEST == 1)
2 !===========================================================================
3 !
4 ! E-GRID NESTING UTILITIES: This is gopal's doing
5 !
6 !===========================================================================
7
8 SUBROUTINE med_nest_egrid_configure ( parent , nest )
9 USE module_domain
10 USE module_configure
11 USE module_timing
12
13 IMPLICIT NONE
14 TYPE(domain) , POINTER :: parent , nest
15 REAL, PARAMETER :: ERAD=6371200.
16 REAL, PARAMETER :: DTR=0.01745329
17 REAL, PARAMETER :: DTAD=1.0
18 REAL, PARAMETER :: CP=1004.6
19 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
20 INTEGER :: IMS,IME,JMS,JME,KMS,KME
21 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
22
23 !----------------------------------------------------------------------------
24 ! PURPOSE:
25 ! - Initialize nested domain configurations including setting up
26 ! wbd,sbd and some other variables and 1D arrays.
27 ! - Note that in order to obtain coincident grid points, which
28 ! is a basic requirement for RSL, WRF infrastructure, we use
29 ! western and southern boundaries of nested domain (nest%wbd0
30 ! and nest%sbd0 derived from the parent domain. In this case
31 ! the nested domain may be considered as a part of the parent
32 ! domain with a higher resolution (telescoping ?).
33 ! - Also note that in this case, the central lat/lons for nested
34 ! domain should coincide with the central lat/lons of the parent,
35 ! although the nested domain NEED NOT be located at the center of
36 ! the domain.
37 !----------------------------------------------------------------------------
38 !
39 ! BASIC TEST FOR PARENT DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
40 ! IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT
41
42 IF(MOD(parent%ed33,2) .NE. 0)THEN
43 CALL wrf_error_fatal("PARENT DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
44 ENDIF
45
46 ! BASIC TEST FOR NESTED DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
47 ! IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT
48
49 IF(MOD(nest%ed33,2) .NE. 0)THEN
50 CALL wrf_error_fatal("NESTED DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
51 ENDIF
52
53 ! Parent grid configuration, including, western and southern boundary
54
55 IDS = parent%sd31
56 IDE = parent%ed31
57 KDS = parent%sd32
58 KDE = parent%ed32
59 JDS = parent%sd33
60 JDE = parent%ed33
61
62 IMS = parent%sm31
63 IME = parent%em31
64 KMS = parent%sm32
65 KME = parent%em32
66 JMS = parent%sm33
67 JME = parent%em33
68
69 ITS = parent%sp31
70 ITE = parent%ep31
71 KTS = parent%sp32
72 KTE = parent%ep32
73 JTS = parent%sp33
74 JTE = parent%ep33
75
76 ! grid configuration
77
78 ! calculate wbd0 and sbd0 only for MOAD i.e. grid with parent_id == 0
79 if (parent%parent_id == 0 ) then ! Dusan's doing
80 parent%wbd0 = -(IDE-2)*parent%dx ! WBD0: in degrees;factor 2 takes care of dummy last column
81 parent%sbd0 = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd
82 end if
83 nest%wbd0 = parent%wbd0 + (nest%i_parent_start-1)*2.*parent%dx + mod(nest%j_parent_start+1,2)*parent%dx
84 nest%sbd0 = parent%sbd0 + (nest%j_parent_start-1)*parent%dy
85 nest%dx = parent%dx/nest%parent_grid_ratio
86 nest%dy = parent%dy/nest%parent_grid_ratio
87
88 write(0,*)" - i_parent_start = ",nest%i_parent_start
89 write(0,*)" - j_parent_start = ",nest%j_parent_start
90 write(0,*)" - parent%wbd0 = ",parent%wbd0
91 write(0,*)" - parent%sbd0 = ",parent%sbd0
92 write(0,*)" - nest%wbd0 = ",nest%wbd0
93 write(0,*)" - nest%sbd0 = ",nest%sbd0
94 write(0,*)" - nest%dx = ",nest%dx
95 write(0,*)" - nest%dy = ",nest%dy
96 !
97 CALL nl_set_dx (nest%id , nest%dx) ! for output purpose
98 CALL nl_set_dy (nest%id , nest%dy) ! for output purpose
99
100 ! set lat-lons; parent set to nested domain
101
102 CALL nl_get_cen_lat (parent%id, parent%cen_lat) ! cen_lat of parent set to nested domain
103 CALL nl_get_cen_lon (parent%id, parent%cen_lon) ! cen_lon of parent set to nested domain
104
105 nest%cen_lat=parent%cen_lat
106 nest%cen_lon=parent%cen_lon
107 !
108 CALL nl_set_cen_lat ( nest%id , nest%cen_lat) ! for output purpose
109 CALL nl_set_cen_lon ( nest%id , nest%cen_lon) ! for output purpose
110
111 write(0,*)" - nest%cen_lat = ",nest%cen_lat
112 write(0,*)" - nest%cen_lon = ",nest%cen_lon
113
114
115 ! soil configuration
116
117 nest%nmm_sldpth = parent%nmm_sldpth
118 nest%nmm_dzsoil = parent%nmm_dzsoil
119 nest%nmm_rtdpth = parent%nmm_rtdpth
120
121 ! numerical set up
122
123 nest%nmm_deta = parent%nmm_deta
124 nest%nmm_aeta = parent%nmm_aeta
125 nest%nmm_etax = parent%nmm_etax
126 nest%nmm_dfl = parent%nmm_dfl
127 nest%nmm_deta1 = parent%nmm_deta1
128 nest%nmm_aeta1 = parent%nmm_aeta1
129 nest%nmm_eta1 = parent%nmm_eta1
130 nest%nmm_deta2 = parent%nmm_deta2
131 nest%nmm_aeta2 = parent%nmm_aeta2
132 nest%nmm_eta2 = parent%nmm_eta2
133 nest%nmm_pdtop = parent%nmm_pdtop
134 nest%nmm_pt = parent%nmm_pt
135 nest%nmm_dfrlg = parent%nmm_dfrlg
136 nest%num_soil_layers = parent%num_soil_layers
137 nest%num_moves = parent%num_moves
138
139 ! Unfortunately, some of the single value constants in used in module_initialize have
140 ! to be defiend here instead of the usual spot in med_initialize_nest_nmm. There
141 ! appears to be a problem in Registry and related code in this area.
142 !
143 ! state logical upstrm - dyn_nmm - - -
144
145
146 nest%nmm_dlmd = nest%dx
147 nest%nmm_dphd = nest%dy
148 nest%nmm_dy_nmm = erad*(nest%nmm_dphd*dtr)
149 nest%nmm_cpgfv = -nest%dt/(48.*nest%nmm_dy_nmm)
150 nest%nmm_en = nest%dt/( 4.*nest%nmm_dy_nmm)*dtad
151 nest%nmm_ent = nest%dt/(16.*nest%nmm_dy_nmm)*dtad
152 nest%nmm_f4d = -.5*nest%dt*dtad
153 nest%nmm_f4q = -nest%dt*dtad
154 nest%nmm_ef4t = .5*nest%dt/cp
155
156 ! Other output configurations that will make grads happy
157
158 CALL nl_get_truelat1 (parent%id, parent%truelat1 )
159 CALL nl_get_truelat2 (parent%id, parent%truelat2 )
160 CALL nl_get_map_proj (parent%id, parent%map_proj )
161 CALL nl_get_iswater (parent%id, parent%iswater )
162
163 nest%truelat1=parent%truelat1
164 nest%truelat2=parent%truelat2
165 nest%map_proj=parent%map_proj
166 nest%iswater=parent%iswater
167
168 CALL nl_set_truelat1(nest%id, nest%truelat1)
169 CALL nl_set_truelat2(nest%id, nest%truelat2)
170 CALL nl_set_map_proj(nest%id, nest%map_proj)
171 CALL nl_set_iswater(nest%id, nest%iswater)
172
173 ! physics and other configurations
174 ! CALL nl_get_iswater (parent%id, nest%iswater) ! iswater is just based on parents
175 ! CALL nl_get_bl_surface_physics (nest%id, nest%bl_surface_physics )
176 ! CALL nl_get_num_soil_layers( parent%num_soil_layers )
177 ! CALL nl_get_real_data_init_type (parent%real_data_init_type)
178
179 END SUBROUTINE med_nest_egrid_configure
180
181 SUBROUTINE med_construct_egrid_weights ( parent , nest )
182 USE module_domain
183 USE module_configure
184 USE module_timing
185
186 IMPLICIT NONE
187 TYPE(domain) , POINTER :: parent , nest
188 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
189 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
190 INTEGER :: IMS,IME,JMS,JME,KMS,KME
191 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
192 INTEGER :: I,J,II,JJ,NII,NJJ
193 REAL :: parent_CLAT,parent_CLON,parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
194 REAL :: nest_WBD,nest_SBD,nest_DLMD,nest_DPHD
195 REAL :: SW_LATD,SW_LOND
196 REAL :: ADDSUM1,ADDSUM2
197 REAL :: xr,zr,xc
198 !-----------------------------------------------------------------------------------------------------------
199 ! PURPOSE:
200 ! - Initialize lat-lons and determine weights
201 !
202 !----------------------------------------------------------------------------------------------------------
203
204 ! First obtain central latitude and longitude for the parent domain
205
206 CALL nl_get_cen_lat (parent%ID, parent_CLAT)
207 CALL nl_get_cen_lon (parent%ID, parent_CLON)
208
209 ! Parent grid configuration, including, western and southern boundary
210
211 IDS = parent%sd31
212 IDE = parent%ed31
213 KDS = parent%sd32
214 KDE = parent%ed32
215 JDS = parent%sd33
216 JDE = parent%ed33
217
218 IMS = parent%sm31
219 IME = parent%em31
220 KMS = parent%sm32
221 KME = parent%em32
222 JMS = parent%sm33
223 JME = parent%em33
224
225 ITS = parent%sp31
226 ITE = parent%ep31
227 KTS = parent%sp32
228 KTE = parent%ep32
229 JTS = parent%sp33
230 JTE = parent%ep33
231 !
232 parent_DLMD = parent%dx ! DLMD: dlamda in degrees
233 parent_DPHD = parent%dy ! DPHD: dphi in degrees
234 parent_WBD = parent%wbd0
235 parent_SBD = parent%sbd0
236
237 ! Now compute Geodetic lat/lon (Positive East) of parent grid in degrees
238
239 CALL EARTH_LATLON ( parent%nmm_HLAT,parent%nmm_HLON,parent%nmm_VLAT,parent%nmm_VLON, & !output
240 parent_DLMD,parent_DPHD,parent_WBD,parent_SBD, & !inputs
241 parent_CLAT,parent_CLON, &
242 IDS,IDE,JDS,JDE,KDS,KDE, &
243 IMS,IME,JMS,JME,KMS,KME, &
244 ITS,ITE,JTS,JTE,KTS,KTE )
245
246 ! Nested grid configuration, including, western and southern boundary
247
248 IDS = nest%sd31
249 IDE = nest%ed31
250 KDS = nest%sd32
251 KDE = nest%ed32
252 JDS = nest%sd33
253 JDE = nest%ed33
254
255 IMS = nest%sm31
256 IME = nest%em31
257 KMS = nest%sm32
258 KME = nest%em32
259 JMS = nest%sm33
260 JME = nest%em33
261
262 ITS = nest%sp31
263 ITE = nest%ep31
264 KTS = nest%sp32
265 KTE = nest%ep32
266 JTS = nest%sp33
267 JTE = nest%ep33
268 !
269 nest_DLMD = nest%dx
270 nest_DPHD = nest%dy
271 nest_WBD = nest%wbd0
272 nest_SBD = nest%sbd0
273
274 !
275 ! Now compute Geodetic lat/lon (Positive East) of nest in degrees, with the same central lat-lon
276 ! as the parent grid
277 !
278
279 CALL EARTH_LATLON ( nest%nmm_HLAT,nest%nmm_HLON,nest%nmm_VLAT,nest%nmm_VLON, & ! output
280 nest_DLMD,nest_DPHD,nest_WBD,nest_SBD, & ! nest inputs
281 parent_CLAT,parent_CLON, & ! parent central lat/lon
282 IDS,IDE,JDS,JDE,KDS,KDE, & ! nested domain dimension
283 IMS,IME,JMS,JME,KMS,KME, &
284 ITS,ITE,JTS,JTE,KTS,KTE )
285
286 ! Determine the weights of nested grid h points nearest to H points of parent domain
287
288 CALL G2T2H( nest%nmm_IIH,nest%nmm_JJH, & ! output grid index on nested grid
289 nest%nmm_HBWGT1,nest%nmm_HBWGT2, & ! output weights on the nested grid
290 nest%nmm_HBWGT3,nest%nmm_HBWGT4, &
291 nest%nmm_HLAT,nest%nmm_HLON, & ! target (nest) input lat lon in degrees
292 parent_DLMD,parent_DPHD,parent_WBD,parent_SBD, & ! parent res, western and south boundaries
293 parent_CLAT,parent_CLON, & ! parent central lat,lon, all in degrees
294 parent%ed31,parent%ed33, & ! parent imax and jmax
295 IDS,IDE,JDS,JDE,KDS,KDE, & !
296 IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
297 ITS,ITE,JTS,JTE,KTS,KTE ) !
298
299
300 ! Determine the weights of nested grid v points nearest to V points of parent domain
301
302 CALL G2T2V( nest%nmm_IIV,nest%nmm_JJV, & ! output grid index on nested grid
303 nest%nmm_VBWGT1,nest%nmm_VBWGT2, & ! output weights on the nested grid
304 nest%nmm_VBWGT3,nest%nmm_VBWGT4, &
305 nest%nmm_VLAT,nest%nmm_VLON, & ! target (nest) input lat lon in degrees
306 parent_DLMD,parent_DPHD,parent_WBD,parent_SBD, & ! parent res, western and south boundaries
307 parent_CLAT,parent_CLON, & ! parent central lat,lon, all in degrees
308 parent%ed31,parent%ed33, & ! parent imax and jmax
309 IDS,IDE,JDS,JDE,KDS,KDE, & !
310 IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
311 ITS,ITE,JTS,JTE,KTS,KTE ) !
312
313 !*** CHECK WEIGHTS AT MASS AND VELOCITY POINTS
314
315 CALL WEIGTS_CHECK(nest%nmm_HBWGT1,nest%nmm_HBWGT2, & ! output weights on the nested grid
316 nest%nmm_HBWGT3,nest%nmm_HBWGT4, &
317 nest%nmm_VBWGT1,nest%nmm_VBWGT2, & ! output weights on the nested grid
318 nest%nmm_VBWGT3,nest%nmm_VBWGT4, &
319 IDS,IDE,JDS,JDE,KDS,KDE, & !
320 IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
321 ITS,ITE,JTS,JTE,KTS,KTE )
322
323 !*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
324 !
325 CALL BOUNDS_CHECK( nest%nmm_IIH,nest%nmm_JJH,nest%nmm_IIV,nest%nmm_JJV, &
326 nest%i_parent_start,nest%j_parent_start,nest%shw, &
327 IDS,IDE,JDS,JDE,KDS,KDE, & !
328 IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
329 ITS,ITE,JTS,JTE,KTS,KTE )
330
331 !------------------------------------------------------------------------------------------
332
333 END SUBROUTINE med_construct_egrid_weights
334
335 !======================================================================================
336 !
337 ! compute earth lat-lons for parent and the nest before interpolations
338 !------------------------------------------------------------------------------
339
340 SUBROUTINE EARTH_LATLON ( HLAT,HLON,VLAT,VLON, & !Earth lat,lon at H and V points
341 DLMD1,DPHD1,WBD1,SBD1, & !input res,west & south boundaries,
342 CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees
343 IDS,IDE,JDS,JDE,KDS,KDE, &
344 IMS,IME,JMS,JME,KMS,KME, &
345 ITS,ITE,JTS,JTE,KTS,KTE )
346 !============================================================================
347 !
348 IMPLICIT NONE
349 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
350 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
351 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
352 REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
353 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
354 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HLAT,HLON,VLAT,VLON
355
356 ! local
357
358 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
359 INTEGER :: I,J
360 REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
361 REAL(KIND=KNUM) :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR
362 REAL(KIND=KNUM) :: STPH,CTPH,STPV,CTPV,PI_2
363 REAL(KIND=KNUM) :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV
364 REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV
365 !-------------------------------------------------------------------------
366
367 !
368 PI_2 = ACOS(0.)
369 DTR = PI_2/90.
370 WB = WBD1 * DTR ! WB: western boundary in radians
371 SB = SBD1 * DTR ! SB: southern boundary in radians
372 DLM = DLMD1 * DTR ! DLM: dlamda in radians
373 DPH = DPHD1 * DTR ! DPH: dphi in radians
374 TDLM = DLM + DLM ! TDLM: 2.0*dlamda
375 TDPH = DPH + DPH ! TDPH: 2.0*DPH
376
377 ! For earth lat lon only
378
379 TPH0 = CENTRAL_LAT*DTR ! TPH0: central lat in radians
380 STPH0 = SIN(TPH0)
381 CTPH0 = COS(TPH0)
382
383 ! WRITE(0,*) 'WB,SB,DLM,DPH,DTR: ',WBD1,SBD1,DLM,DPH,DTR
384 ! WRITE(0,*) 'IMS,IME,JMS,JME,KMS,KME',IMS,IME,JMS,JME,KMS,KME
385 ! WRITE(0,*) 'IDS,IDE,JDS,JDE,KDS,KDE',IDS,IDE,JDS,JDE,KDS,KDE
386 ! WRITE(0,*) 'ITS,ITE,JTS,JTE,KTS,KTE',ITS,ITE,JTS,JTE,KTS,KTE
387
388 ! .H
389 DO J = JTS,MIN(JTE,JDE-1) ! H./ This loop takes care of zig-zag
390 ! ! \.H starting points along j
391 TLMH0 = WB - TDLM + MOD(J+1,2) * DLM ! ./ TLMH (rotated lats at H points)
392 TLMV0 = WB - TDLM + MOD(J,2) * DLM ! H (//ly for V points)
393 TPHH = SB + (J-1)*DPH ! TPHH (rotated lons at H points) are simple trans.
394 TPHV = TPHH ! TPHV (rotated lons at V points) are simple trans.
395 STPH = SIN(TPHH)
396 CTPH = COS(TPHH)
397 STPV = SIN(TPHV)
398 CTPV = COS(TPHV)
399
400 ! .H
401 DO I = ITS,MIN(ITE,IDE-1) ! /
402 TLMH = TLMH0 + I*TDLM ! \.H .U .H
403 ! !H./ ----><----
404 SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH) ! DLM + DLM
405 GLATH(I,J)=ASIN(SPHH) ! GLATH: Earth Lat in radians
406 CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) &
407 - TAN(GLATH(I,J))*TAN(TPH0)
408 IF(CLMH .GT. 1.) CLMH = 1.0
409 IF(CLMH .LT. -1.) CLMH = -1.0
410 FACTH = 1.
411 IF(TLMH .GT. 0.) FACTH = -1.
412 GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH)
413
414 ENDDO
415
416 DO I = ITS,MIN(ITE,IDE-1)
417 TLMV = TLMV0 + I*TDLM
418 SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV)
419 GLATV(I,J) = ASIN(SPHV)
420 CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) &
421 - TAN(GLATV(I,J))*TAN(TPH0)
422 IF(CLMV .GT. 1.) CLMV = 1.
423 IF(CLMV .LT. -1.) CLMV = -1.
424 FACTV = 1.
425 IF(TLMV .GT. 0.) FACTV = -1.
426 GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV)
427
428 ENDDO
429
430 ENDDO
431
432 ! Conversion to degrees (may not be required, eventually)
433
434 DO J = JTS, MIN(JTE,JDE-1)
435 DO I = ITS, MIN(ITE,IDE-1)
436 HLAT(I,J) = GLATH(I,J) / DTR
437 HLON(I,J)= -GLONH(I,J)/DTR
438 IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J) - 360.
439 IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360.
440 !
441 VLAT(I,J) = GLATV(I,J) / DTR
442 VLON(I,J) = -GLONV(I,J) / DTR
443 IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J) - 360.
444 IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360.
445
446 ENDDO
447 ENDDO
448
449 END SUBROUTINE EARTH_LATLON
450
451 !-----------------------------------------------------------------------------
452
453 SUBROUTINE G2T2H( IIH,JJH, & ! output grid index and weights
454 HBWGT1,HBWGT2, & ! output weights in terms of parent grid
455 HBWGT3,HBWGT4, &
456 HLAT,HLON, & ! target (nest) input lat lon in degrees
457 DLMD1,DPHD1,WBD1,SBD1, & ! parent res, west and south boundaries
458 CENTRAL_LAT,CENTRAL_LON, & ! parent central lat,lon, all in degrees
459 P_IDE,P_JDE, & ! parent imax and jmax
460 IDS,IDE,JDS,JDE,KDS,KDE, & ! target (nest) dIMEnsions
461 IMS,IME,JMS,JME,KMS,KME, &
462 ITS,ITE,JTS,JTE,KTS,KTE )
463
464 !
465 !*** Tom Black - Initial Version
466 !*** Gopal - Revised Version for WRF (includes coincident grid points)
467 !***
468 !*** GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
469 !*** AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE
470 !*** INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE
471 !*** h POINTS OF THE NESTED DOMAIN
472 !
473 !============================================================================
474 !
475 IMPLICIT NONE
476 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
477 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
478 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
479 INTEGER, INTENT(IN ) :: P_IDE,P_JDE
480 REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
481 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
482 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HLAT,HLON
483 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
484 INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
485
486 ! local
487
488 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
489 INTEGER :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
490 INTEGER :: NROW,NCOL,KROWS
491 REAL(KIND=KNUM) :: X,Y,Z,TLAT,TLON
492 REAL(KIND=KNUM) :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
493 REAL(KIND=KNUM) :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE
494 REAL(KIND=KNUM) :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
495 REAL(KIND=KNUM) :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3 ! Q
496 REAL(KIND=KNUM) :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
497 REAL(KIND=KNUM) :: DTEMP
498 REAL , DIMENSION(IMS:IME,JMS:JME) :: TLATHX,TLONHX
499 INTEGER, DIMENSION(IMS:IME,JMS:JME) :: KOUTB
500 !-------------------------------------------------------------------------------
501
502 IMT=2*P_IDE-2 ! parent i dIMEnsions
503 JMT=P_JDE/2 ! parent j dIMEnsions
504 PI_2=ACOS(0.)
505 D2R=PI_2/90.
506 R2D=1./D2R
507 DPH=DPHD1*D2R
508 DLM=DLMD1*D2R
509 TPH0= CENTRAL_LAT*D2R
510 TLM0=-CENTRAL_LON*D2R ! NOTE THE MINUS HERE
511 WB=WBD1*D2R ! CONVERT NESTED GRID H POINTS FROM GEODETIC
512 SB=SBD1*D2R
513 SLP0=DPHD1/DLMD1
514 DSLP0=NINT(R2D*ATAN(SLP0))
515 DS1=SQRT(DPH*DPH+DLM*DLM) ! Q
516 AN1=ASIN(DLM/DS1)
517 AN2=ASIN(DPH/DS1)
518
519 DO J = JTS,MIN(JTE,JDE-1)
520 DO I = ITS,MIN(ITE,IDE-1)
521
522 !***
523 !*** LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND
524 !*** DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
525 !*** CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED
526 !*** COORDINATE ON THE PARENT GRID
527 !
528
529 GLAT=HLAT(I,J)*D2R
530 GLON= (360. - HLON(I,J))*D2R
531 X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
532 Y=-COS(GLAT)*SIN(GLON-TLM0)
533 Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
534 TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
535 TLON=R2D*ATAN(Y/X)
536
537 ! ROW=TLAT/DPHD1+JMT ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
538 ! COL=TLON/DLMD1+P_IDE-1 ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
539
540 ROW=(TLAT-SBD1)/DPHD1+1 ! Dusan's doing
541 COL=(TLON-WBD1)/DLMD1+1 ! Dusan's doing
542
543 NROW=INT(ROW + 0.001) ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
544 NCOL=INT(COL + 0.001)
545 TLAT=TLAT*D2R
546 TLON=TLON*D2R
547
548 !***
549 !***
550 !*** FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT
551 !***
552 !*** V H
553 !***
554 !***
555 !*** h
556 !*** H V
557 !***
558 !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
559 !***
560 IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR. &
561 MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN
562 TLAT1=(NROW-JMT)*DPH
563 TLAT2=TLAT1+DPH
564 TLON1=(NCOL-(P_IDE-1))*DLM
565 TLON2=TLON1+DLM
566 DLM1=TLON-TLON1
567 DLM2=TLON-TLON2
568 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
569 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
570 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
571 D1=ACOS(DTEMP)
572 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
573 D2=ACOS(DTEMP)
574 IF(D1.GT.D2)THEN
575 NROW=NROW+1 ! FIND THE NEAREST H ROW
576 NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
577 ENDIF
578 ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
579 ELSE
580 !***
581 !*** NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT
582 !***
583 !*** H V
584 !***
585 !***
586 !*** h
587 !*** V H
588 !***
589 !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
590 !***
591 !***
592 TLAT1=(NROW+1-JMT)*DPH
593 TLAT2=TLAT1-DPH
594 TLON1=(NCOL-(P_IDE-1))*DLM
595 TLON2=TLON1+DLM
596 DLM1=TLON-TLON1
597 DLM2=TLON-TLON2
598 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
599 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
600 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
601 D1=ACOS(DTEMP)
602 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
603 D2=ACOS(DTEMP)
604 IF(D1.LT.D2)THEN
605 NROW=NROW+1 ! FIND THE NEAREST H ROW
606 ELSE
607 NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
608 ENDIF
609 ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
610 ENDIF
611
612 KROWS=((NROW-1)/2)*IMT
613 IF(MOD(NROW,2).EQ.1)THEN
614 K=KROWS+(NCOL+1)/2
615 ELSE
616 K=KROWS+P_IDE-1+NCOL/2
617 ENDIF
618
619 !***
620 !*** WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
621 !*** NEAREST TO THE CENTER K AS SEEN BELOW. WE MUST FIND
622 !*** WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS
623 !*** A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION.
624 !***
625 !**
626 !*** H
627 !***
628 !***
629 !***
630 !*** H V H
631 !***
632 !***
633 !*** h
634 !*** H V H V H
635 !***
636 !***
637 !***
638 !*** H V H
639 !***
640 !***
641 !***
642 !*** H
643 !***
644 !***
645 !*** FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H.
646 !***
647 N2R=K/IMT
648 MK=MOD(K,IMT)
649 !
650 IF(MK.EQ.0)THEN
651 TLATHC=SB+(2*N2R-1)*DPH
652 ELSE
653 TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH
654 ENDIF
655 !
656 IF(MK.LE.(P_IDE-1))THEN
657 TLONHC=WB+2*(MK-1)*DLM
658 ELSE
659 TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM
660 ENDIF
661
662 !
663 !*** EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
664 !*** DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE
665 !*** CAREFUL HERE
666 !
667
668 IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON
669 IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT
670 DENOM=(TLON-TLONHC)
671 !
672 !***
673 !***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
674 !***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
675 !***
676 !*** COINCIDENT CONDITIONS
677
678 IF(DENOM.EQ.0.0)THEN
679
680 IF(TLATHC.EQ.TLAT)THEN
681 KOUTB(I,J)=K
682 IIH(I,J) = NCOL
683 JJH(I,J) = NROW
684 TLATHX(I,J)=TLATHC
685 TLONHX(I,J)=TLONHC
686 HBWGT1(I,J)=1.0
687 HBWGT2(I,J)=0.0
688 HBWGT3(I,J)=0.0
689 HBWGT4(I,J)=0.0
690 ! WRITE(60,*)'TRIVIAL SOLUTION'
691 ELSE ! SAME LONGITUDE BUT DIFFERENT LATS
692 !
693 IF(TLATHC .GT. TLAT)THEN ! NESTED POINT SOUTH OF PARENT
694 KOUTB(I,J)=K-(P_IDE-1)
695 IIH(I,J) = NCOL-1
696 JJH(I,J) = NROW-1
697 TLATHX(I,J)=TLATHC-DPH
698 TLONHX(I,J)=TLONHC-DLM
699 ! WRITE(60,*)'VANISHING SLOPE, -ve: TLATHC-DPH, TLONHC-DLM'
700 ELSE ! NESTED POINT NORTH OF PARENT
701 KOUTB(I,J)=K+(P_IDE-1)-1
702 IIH(I,J) = NCOL-1
703 JJH(I,J) = NROW+1
704 TLATHX(I,J)=TLATHC+DPH
705 TLONHX(I,J)=TLONHC-DLM
706 ! WRITE(60,*)'VANISHING SLOPE, +ve: TLATHC+DPH, TLONHC-DLM'
707 ENDIF
708 !***
709 !***
710 !*** 4
711 !***
712 !*** h
713 !*** 1 2
714 !***
715 !*** 3
716 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
717
718 TLATO=TLATHX(I,J)
719 TLONO=TLONHX(I,J)
720 DLM1=TLON-TLONO
721 DLA1=TLAT-TLATO ! Q
722 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
723 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
724 !
725 TLATO=TLATHX(I,J)
726 TLONO=TLONHX(I,J)+2.*DLM
727 DLM2=TLON-TLONO
728 DLA2=TLAT-TLATO ! Q
729 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
730 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
731 !
732 TLATO=TLATHX(I,J)-DPH
733 TLONO=TLONHX(I,J)+DLM
734 DLM3=TLON-TLONO
735 DLA3=TLAT-TLATO ! Q
736 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
737 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
738
739 TLATO=TLATHX(I,J)+DPH
740 TLONO=TLONHX(I,J)+DLM
741 DLM4=TLON-TLONO
742 DLA4=TLAT-TLATO ! Q
743 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
744 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
745
746
747 ! THE BILINEAR WEIGHTS
748 !***
749 !***
750 AN3=ATAN2(DLA1,DLM1) ! Q
751 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
752 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
753 R1=R1/DS1
754 S1=S1/DS1
755 DL1I=(1.-R1)*(1.-S1)
756 DL2I=R1*S1
757 DL3I=R1*(1.-S1)
758 DL4I=(1.-R1)*S1
759 !
760 HBWGT1(I,J)=DL1I
761 HBWGT2(I,J)=DL2I
762 HBWGT3(I,J)=DL3I
763 HBWGT4(I,J)=DL4I
764 !
765 ENDIF
766
767 ELSE
768 !
769 !*** NON-COINCIDENT POINTS
770 !
771 SLOPE=(TLAT-TLATHC)/DENOM
772 DSLOPE=NINT(R2D*ATAN(SLOPE))
773
774 IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
775 IF(TLON.GT.TLONHC)THEN
776 IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
777 KOUTB(I,J)=K
778 IIH(I,J) = NCOL
779 JJH(I,J) = NROW
780 TLATHX(I,J)=TLATHC
781 TLONHX(I,J)=TLONHC
782 ! WRITE(60,*)'HERE WE GO1: TLATHC, TLONHC'
783 ELSE
784 IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
785 KOUTB(I,J)=K-1
786 IIH(I,J) = NCOL-2
787 JJH(I,J) = NROW
788 TLATHX(I,J)=TLATHC
789 TLONHX(I,J)=TLONHC -2.*DLM
790 ! WRITE(60,*)'HERE WE GO2: TLATHC, TLONHC -2.*DLM'
791 ENDIF
792
793 !
794 ELSEIF(DSLOPE.GT.DSLP0)THEN
795 IF(TLON.GT.TLONHC)THEN
796 IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
797 KOUTB(I,J)=K+(P_IDE-1)-1
798 IIH(I,J) = NCOL-1
799 JJH(I,J) = NROW+1
800 TLATHX(I,J)=TLATHC+DPH
801 TLONHX(I,J)=TLONHC-DLM
802 ! WRITE(60,*)'HERE WE GO3: TLATHC+DPH, TLONHC-DLM'
803 ELSE
804 IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
805 KOUTB(I,J)=K-(P_IDE-1)
806 IIH(I,J) = NCOL-1
807 JJH(I,J) = NROW-1
808 TLATHX(I,J)=TLATHC-DPH
809 TLONHX(I,J)=TLONHC-DLM
810 ! WRITE(60,*)'HERE WE GO4: TLATHC-DPH, TLONHC-DLM'
811 ENDIF
812
813 !
814 ELSEIF(DSLOPE.LT.-DSLP0)THEN
815 IF(TLON.GT.TLONHC)THEN
816 IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
817 KOUTB(I,J)=K-(P_IDE-1)
818 IIH(I,J) = NCOL-1
819 JJH(I,J) = NROW-1
820 TLATHX(I,J)=TLATHC-DPH
821 TLONHX(I,J)=TLONHC-DLM
822 ! WRITE(60,*)'HERE WE GO5: TLATHC-DPH, TLONHC-DLM'
823 ELSE
824 IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
825 KOUTB(I,J)=K+(P_IDE-1)-1
826 IIH(I,J) = NCOL-1
827 JJH(I,J) = NROW+1
828 TLATHX(I,J)=TLATHC+DPH
829 TLONHX(I,J)=TLONHC-DLM
830 ! WRITE(60,*)'HERE WE GO6: TLATHC+DPH, TLONHC-DLM'
831 ENDIF
832 ENDIF
833
834 !
835 !*** NOW WE WILL MOVE AS FOLLOWS:
836 !***
837 !***
838 !*** 4
839 !***
840 !***
841 !***
842 !*** h
843 !*** 1 2
844 !***
845 !***
846 !***
847 !***
848 !*** 3
849 !***
850 !***
851 !***
852 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
853
854 TLATO=TLATHX(I,J)
855 TLONO=TLONHX(I,J)
856 DLM1=TLON-TLONO
857 DLA1=TLAT-TLATO ! Q
858 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
859 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
860 !
861 TLATO=TLATHX(I,J) ! redundant computations
862 TLONO=TLONHX(I,J)+2.*DLM
863 DLM2=TLON-TLONO
864 DLA2=TLAT-TLATO ! Q
865 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
866 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
867 !
868 TLATO=TLATHX(I,J)-DPH
869 TLONO=TLONHX(I,J)+DLM
870 DLM3=TLON-TLONO
871 DLA3=TLAT-TLATO ! Q
872 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
873 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
874 !
875 TLATO=TLATHX(I,J)+DPH
876 TLONO=TLONHX(I,J)+DLM
877 DLM4=TLON-TLONO
878 DLA4=TLAT-TLATO ! Q
879 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
880 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
881
882 ! THE BILINEAR WEIGHTS
883 !***
884 AN3=ATAN2(DLA1,DLM1) ! Q
885 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
886 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
887 R1=R1/DS1
888 S1=S1/DS1
889 DL1I=(1.-R1)*(1.-S1)
890 DL2I=R1*S1
891 DL3I=R1*(1.-S1)
892 DL4I=(1.-R1)*S1
893 !
894 HBWGT1(I,J)=DL1I
895 HBWGT2(I,J)=DL2I
896 HBWGT3(I,J)=DL3I
897 HBWGT4(I,J)=DL4I
898 !
899 ENDIF
900
901 !
902 !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX
903 !
904 IIH(I,J)=NINT(0.5*IIH(I,J))
905
906 HBWGT1(I,J)=MAX(HBWGT1(I,J),0.0) ! all weights must be GE zero (postive def)
907 HBWGT2(I,J)=MAX(HBWGT2(I,J),0.0) ! all weights must be GE zero (postive def)
908 HBWGT3(I,J)=MAX(HBWGT3(I,J),0.0) ! all weights must be GE zero (postive def)
909 HBWGT4(I,J)=MAX(HBWGT4(I,J),0.0) ! all weights must be GE zero (postive def)
910
911 ! write(0,105)"H WEIGHTS:",I,J,HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J), &
912 ! HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J),IIH(i,j),JJH(i,j)
913 ! 105 format(a,2i4,5f7.3,2i4)
914
915 ENDDO
916 ENDDO
917
918
919 RETURN
920 END SUBROUTINE G2T2H
921 !========================================================================================
922
923
924 SUBROUTINE G2T2V( IIV,JJV, & ! output grid index and weights
925 VBWGT1,VBWGT2, & ! output weights in terms of parent grid
926 VBWGT3,VBWGT4, &
927 VLAT,VLON, & ! target (nest) input lat lon in degrees
928 DLMD1,DPHD1,WBD1,SBD1, & ! parent res, west and south boundaries
929 CENTRAL_LAT,CENTRAL_LON, & ! parent central lat,lon, all in degrees
930 P_IDE,P_JDE, & ! parent imax and jmax
931 IDS,IDE,JDS,JDE,KDS,KDE, & ! target (nest) dIMEnsions
932 IMS,IME,JMS,JME,KMS,KME, &
933 ITS,ITE,JTS,JTE,KTS,KTE )
934
935 !
936 !*** Tom Black - Initial Version
937 !*** Gopal - Revised Version for WRF (includes coincIDEnt grid points)
938 !***
939 !*** GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
940 !*** AND THE NESTED GRID LAT-LONS AT v POINTS, THIS ROUTINE FIRST LOCATES THE
941 !*** INDICES,IIV,JJV, OF THE PARENT DOMAIN'S v POINTS THAT LIES CLOSEST TO THE
942 !*** v POINTS OF THE NESTED DOMAIN
943 !
944 !============================================================================
945
946 IMPLICIT NONE
947 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
948 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
949 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
950 INTEGER, INTENT(IN ) :: P_IDE,P_JDE
951 REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
952 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
953 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: VLAT,VLON
954 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
955 INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
956
957 ! local
958
959 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13) ! (6) single precision
960 INTEGER :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
961 INTEGER :: NROW,NCOL,KROWS
962 REAL(KIND=KNUM) :: X,Y,Z,TLAT,TLON
963 REAL(KIND=KNUM) :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
964 REAL(KIND=KNUM) :: ROW,COL,SLP0,TLATVC,TLONVC,DENOM,SLOPE
965 REAL(KIND=KNUM) :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
966 REAL(KIND=KNUM) :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3 ! Q
967 REAL(KIND=KNUM) :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
968 REAL(KIND=KNUM) :: DTEMP
969 REAL , DIMENSION(IMS:IME,JMS:JME) :: TLATVX,TLONVX
970 INTEGER, DIMENSION(IMS:IME,JMS:JME) :: KOUTB
971 !-------------------------------------------------------------------------------------
972
973 IMT=2*P_IDE-2 ! parent i dIMEnsions
974 JMT=P_JDE/2 ! parent j dIMEnsions
975 PI_2=ACOS(0.)
976 D2R=PI_2/90.
977 R2D=1./D2R
978 DPH=DPHD1*D2R
979 DLM=DLMD1*D2R
980 TPH0= CENTRAL_LAT*D2R
981 TLM0=-CENTRAL_LON*D2R ! NOTE THE MINUS HERE
982 WB=WBD1*D2R ! DEGREES TO RADIANS
983 SB=SBD1*D2R
984 SLP0=DPHD1/DLMD1
985 DSLP0=NINT(R2D*ATAN(SLP0))
986 DS1=SQRT(DPH*DPH+DLM*DLM) ! Q
987 AN1=ASIN(DLM/DS1)
988 AN2=ASIN(DPH/DS1)
989
990 DO J = JTS,MIN(JTE,JDE-1)
991 DO I = ITS,MIN(ITE,IDE-1)
992 !***
993 !*** LOCATE TARGET v POINTS (VLAT AND VLON) ON THE PARENT DOMAIN AND
994 !*** DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
995 !*** CONVERT NESTED GRID v POINTS FROM GEODETIC TO TRANSFORMED
996 !*** COORDINATE ON THE PARENT GRID
997 !
998
999 GLAT=VLAT(I,J)*D2R
1000 GLON=(360. - VLON(I,J))*D2R
1001 X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
1002 Y=-COS(GLAT)*SIN(GLON-TLM0)
1003 Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
1004 TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
1005 TLON=R2D*ATAN(Y/X)
1006
1007 ! ROW=TLAT/DPHD1+JMT ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
1008 ! COL=TLON/DLMD1+P_IDE-1 ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
1009
1010 ROW=(TLAT-SBD1)/DPHD1+1 ! Dusan's doing
1011 COL=(TLON-WBD1)/DLMD1+1 ! Dusan's doing
1012
1013 NROW=INT(ROW + 0.001) ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
1014 NCOL=INT(COL + 0.001)
1015 TLAT=TLAT*D2R
1016 TLON=TLON*D2R
1017
1018 !***
1019 !***
1020 !*** FIRST CONSIDER THE SITUATION WHERE THE POINT v IS AT
1021 !***
1022 !*** H V
1023 !***
1024 !***
1025 !*** v
1026 !*** V H
1027 !***
1028 !*** THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
1029 !***
1030
1031 IF(MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.1.OR. &
1032 MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.0)THEN
1033 TLAT1=(NROW-JMT)*DPH
1034 TLAT2=TLAT1+DPH
1035 TLON1=(NCOL-(P_IDE-1))*DLM
1036 TLON2=TLON1+DLM
1037 DLM1=TLON-TLON1
1038 DLM2=TLON-TLON2
1039 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1040 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1041 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1042 D1=ACOS(DTEMP)
1043 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1044 D2=ACOS(DTEMP)
1045 IF(D1.GT.D2)THEN
1046 NROW=NROW+1 ! FIND THE NEAREST V ROW
1047 NCOL=NCOL+1 ! FIND THE NEAREST V COLUMN
1048 ENDIF
1049 ! WRITE(61,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
1050 ELSE
1051
1052 !***
1053 !*** NOW CONSIDER THE SITUATION WHERE THE POINT v IS AT
1054 !***
1055 !*** V H
1056 !***
1057 !***
1058 !*** v
1059 !*** H V
1060 !***
1061 !*** THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
1062 !***
1063 TLAT1=(NROW+1-JMT)*DPH
1064 TLAT2=TLAT1-DPH
1065 TLON1=(NCOL-(P_IDE-1))*DLM
1066 TLON2=TLON1+DLM
1067 DLM1=TLON-TLON1
1068 DLM2=TLON-TLON2
1069 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1070 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1071 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1072 D1=ACOS(DTEMP)
1073 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1074 D2=ACOS(DTEMP)
1075 IF(D1.LT.D2)THEN
1076 NROW=NROW+1 ! FIND THE NEAREST H ROW
1077 ELSE
1078 NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
1079 ENDIF
1080 ! WRITE(61,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
1081
1082 ENDIF
1083
1084 KROWS=((NROW-1)/2)*IMT
1085 IF(MOD(NROW,2).EQ.1)THEN
1086 K=KROWS+NCOL/2
1087 ELSE
1088 K=KROWS+P_IDE-2+(NCOL+1)/2 ! check this one should this not be P_IDE-2 ????
1089 ENDIF
1090
1091 !***
1092 !*** WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
1093 !*** NEAREST TO THE CENTER K AS SEEN BELOW. WE MUST FIND
1094 !*** WHICH OF THE FOUR V-BOXES (OF WHICH THIS V POINT IS
1095 !*** A VERTEX) SURROUNDS THE INNER GRID v POINT IN QUESTION.
1096 !***
1097 !***
1098 !*** V
1099 !***
1100 !***
1101 !***
1102 !*** V H V
1103 !***
1104 !***
1105 !*** v
1106 !*** V H V H V
1107 !***
1108 !***
1109 !***
1110 !*** V H V
1111 !***
1112 !***
1113 !***
1114 !*** V
1115 !***
1116 !***
1117 !*** FIND THE SLOPE OF THE LINE CONNECTING v AND THE CENTER V.
1118 !***
1119 N2R=K/IMT
1120 MK=MOD(K,IMT)
1121 !
1122 IF(MK.EQ.0)THEN
1123 TLATVC=SB+(2*N2R-1)*DPH
1124 ELSE
1125 TLATVC=SB+(2*N2R+MK/(P_IDE-1))*DPH
1126 ENDIF
1127 !
1128 IF(MK.LE.(P_IDE-1)-1)THEN
1129 TLONVC=WB+(2*MK-1)*DLM
1130 ELSE
1131 TLONVC=WB+2*(MK-(P_IDE-1))*DLM
1132 ENDIF
1133
1134 !
1135 !*** EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
1136 !*** DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE V BOXES, WE NEED TO BE
1137 !*** CAREFUL HERE
1138 !
1139 IF(ABS(TLON-TLONVC) .LE. 1.E-4)TLONVC=TLON
1140 IF(ABS(TLAT-TLATVC) .LE. 1.E-4)TLATVC=TLAT
1141 DENOM=(TLON-TLONVC)
1142 !
1143 !***
1144 !***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
1145 !***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
1146 !***
1147 !*** COINCIDENT CONDITIONS
1148
1149 IF(DENOM.EQ.0.0)THEN
1150
1151 IF(TLATVC.EQ.TLAT)THEN
1152 KOUTB(I,J)=K
1153 IIV(I,J) = NCOL
1154 JJV(I,J) = NROW
1155 TLATVX(I,J)=TLATVC
1156 TLONVX(I,J)=TLONVC
1157 VBWGT1(I,J)=1.0
1158 VBWGT2(I,J)=0.0
1159 VBWGT3(I,J)=0.0
1160 VBWGT4(I,J)=0.0
1161 ! WRITE(61,*)'TRIVIAL SOLUTION'
1162 ELSE ! SAME LONGITUDE BUT DIFFERENT LATS
1163
1164 IF(TLATVC .GT. TLAT)THEN ! NESTED POINT SOUTH OF PARENT
1165 KOUTB(I,J)=K-(P_IDE-1)
1166 IIV(I,J) = NCOL-1
1167 JJV(I,J) = NROW-1
1168 TLATVX(I,J)=TLATVC-DPH
1169 TLONVX(I,J)=TLONVC-DLM
1170 ! WRITE(61,*)'VANISHING SLOPE, -ve: TLATVC-DPH, TLONVC-DLM'
1171 ELSE ! NESTED POINT NORTH OF PARENT
1172 KOUTB(I,J)=K+(P_IDE-1)-1
1173 IIV(I,J) = NCOL-1
1174 JJV(I,J) = NROW+1
1175 TLATVX(I,J)=TLATVC+DPH
1176 TLONVX(I,J)=TLONVC-DLM
1177 ! WRITE(61,*)'VANISHING SLOPE, +ve: TLATVC+DPH, TLONVC-DLM'
1178 ENDIF
1179
1180 !***
1181 !***
1182 !*** 4
1183 !***
1184 !*** v
1185 !*** 1 2
1186 !***
1187 !*** 3
1188 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
1189
1190 TLATO=TLATVX(I,J)
1191 TLONO=TLONVX(I,J)
1192 DLM1=TLON-TLONO
1193 DLA1=TLAT-TLATO ! Q
1194 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
1195 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
1196 !
1197 TLATO=TLATVX(I,J)
1198 TLONO=TLONVX(I,J)+2.*DLM
1199 DLM2=TLON-TLONO
1200 DLA2=TLAT-TLATO ! Q
1201 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
1202 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
1203
1204 TLATO=TLATVX(I,J)-DPH
1205 TLONO=TLONVX(I,J)+DLM
1206 DLM3=TLON-TLONO
1207 DLA3=TLAT-TLATO ! Q
1208 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
1209 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
1210
1211 TLATO=TLATVX(I,J)+DPH
1212 TLONO=TLONVX(I,J)+DLM
1213 DLM4=TLON-TLONO
1214 DLA4=TLAT-TLATO ! Q
1215 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
1216 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
1217
1218 ! THE BILINEAR WEIGHTS
1219 !***
1220 AN3=ATAN2(DLA1,DLM1) ! Q
1221 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
1222 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
1223 R1=R1/DS1
1224 S1=S1/DS1
1225 DL1I=(1.-R1)*(1.-S1)
1226 DL2I=R1*S1
1227 DL3I=R1*(1.-S1)
1228 DL4I=(1.-R1)*S1
1229 !
1230 VBWGT1(I,J)=DL1I
1231 VBWGT2(I,J)=DL2I
1232 VBWGT3(I,J)=DL3I
1233 VBWGT4(I,J)=DL4I
1234
1235 ENDIF
1236
1237 ELSE
1238
1239 !
1240 !*** NON-COINCIDENT POINTS
1241 !
1242 SLOPE=(TLAT-TLATVC)/DENOM
1243 DSLOPE=NINT(R2D*ATAN(SLOPE))
1244
1245 IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
1246 IF(TLON.GT.TLONVC)THEN
1247 IF(TLONVC.GE.-WB-DLM)CALL wrf_error_fatal("1V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1248 KOUTB(I,J)=K
1249 IIV(I,J)=NCOL
1250 JJV(I,J)=NROW
1251 TLATVX(I,J)=TLATVC
1252 TLONVX(I,J)=TLONVC
1253 ! WRITE(61,*)'HERE WE GO1: TLATHC, TLONHC'
1254 ELSE
1255 IF(TLONVC.LE.WB+DLM)CALL wrf_error_fatal("2V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1256 KOUTB(I,J)=K-1
1257 IIV(I,J) = NCOL-2
1258 JJV(I,J) = NROW
1259 TLATVX(I,J)=TLATVC
1260 TLONVX(I,J)=TLONVC-2.*DLM
1261 ! WRITE(61,*)'HERE WE GO2: TLATHC, TLONHC -2.*DLM'
1262 ENDIF
1263
1264 ELSEIF(DSLOPE.GT.DSLP0)THEN
1265 IF(TLON.GT.TLONVC)THEN
1266 IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("3V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1267 KOUTB(I,J)=K+(P_IDE-1)-1
1268 IIV(I,J) = NCOL-1
1269 JJV(I,J) = NROW+1
1270 TLATVX(I,J)=TLATVC+DPH
1271 TLONVX(I,J)=TLONVC-DLM
1272 ! WRITE(61,*)'HERE WE GO3: TLATHC+DPH, TLONHC-DLM'
1273 ELSE
1274 IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("4V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1275 KOUTB(I,J)=K-(P_IDE-1)
1276 IIV(I,J) = NCOL-1
1277 JJV(I,J) = NROW-1
1278 TLATVX(I,J)=TLATVC-DPH
1279 TLONVX(I,J)=TLONVC-DLM
1280 ! WRITE(61,*)'HERE WE GO4: TLATHC-DPH, TLONHC-DLM'
1281 ENDIF
1282
1283 ELSEIF(DSLOPE.LT.-DSLP0)THEN
1284 IF(TLON.GT.TLONVC)THEN
1285 IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("5V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1286 KOUTB(I,J)=K-(P_IDE-1)
1287 IIV(I,J) = NCOL-1
1288 JJV(I,J) = NROW-1
1289 TLATVX(I,J)=TLATVC-DPH
1290 TLONVX(I,J)=TLONVC-DLM
1291 ! WRITE(61,*)'HERE WE GO5: TLATHC-DPH, TLONHC-DLM'
1292 ELSE
1293 IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("6V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1294 KOUTB(I,J)=K+(P_IDE-1)-1
1295 IIV(I,J) = NCOL-1
1296 JJV(I,J) = NROW+1
1297 TLATVX(I,J)=TLATVC+DPH
1298 TLONVX(I,J)=TLONVC-DLM
1299 ! WRITE(61,*)'HERE WE GO6: TLATHC+DPH, TLONHC-DLM'
1300 ENDIF
1301 ENDIF
1302 !
1303 !*** NOW WE WILL MOVE AS FOLLOWS:
1304 !***
1305 !***
1306 !*** 4
1307 !***
1308 !***
1309 !***
1310 !*** v
1311 !*** 1 2
1312 !***
1313 !***
1314 !***
1315 !***
1316 !*** 3
1317 !***
1318 !***
1319 !***
1320 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM v TO EACH VERTEX
1321
1322 TLATO=TLATVX(I,J)
1323 TLONO=TLONVX(I,J)
1324 DLM1=TLON-TLONO
1325 DLA1=TLAT-TLATO ! Q
1326 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
1327 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
1328 !
1329 TLATO=TLATVX(I,J)
1330 TLONO=TLONVX(I,J)+2.*DLM
1331 DLM2=TLON-TLONO
1332 DLA2=TLAT-TLATO ! Q
1333 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
1334 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
1335 !
1336 TLATO=TLATVX(I,J)-DPH
1337 TLONO=TLONVX(I,J)+DLM
1338 DLM3=TLON-TLONO
1339 DLA3=TLAT-TLATO ! Q
1340 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
1341 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
1342 !
1343 TLATO=TLATVX(I,J)+DPH
1344 TLONO=TLONVX(I,J)+DLM
1345 DLM4=TLON-TLONO
1346 DLA4=TLAT-TLATO ! Q
1347 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
1348 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
1349
1350 ! THE BILINEAR WEIGHTS
1351 !***
1352 AN3=ATAN2(DLA1,DLM1) ! Q
1353 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
1354 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
1355 R1=R1/DS1
1356 S1=S1/DS1
1357 DL1I=(1.-R1)*(1.-S1)
1358 DL2I=R1*S1
1359 DL3I=R1*(1.-S1)
1360 DL4I=(1.-R1)*S1
1361 !
1362 VBWGT1(I,J)=DL1I
1363 VBWGT2(I,J)=DL2I
1364 VBWGT3(I,J)=DL3I
1365 VBWGT4(I,J)=DL4I
1366
1367 ENDIF
1368
1369 !
1370 !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX
1371 !
1372 IIV(I,J)=NINT(0.5*IIV(I,J))
1373
1374 VBWGT1(I,J)=MAX(VBWGT1(I,J),0.0) ! all weights must be GE zero (postive def)
1375 VBWGT2(I,J)=MAX(VBWGT2(I,J),0.0) ! all weights must be GE zero (postive def)
1376 VBWGT3(I,J)=MAX(VBWGT3(I,J),0.0) ! all weights must be GE zero (postive def)
1377 VBWGT4(I,J)=MAX(VBWGT4(I,J),0.0) ! all weights must be GE zero (postive def)
1378
1379 ! WRITE(61,*)I,J,VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J), &
1380 ! VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J),IIV(i,j),JJV(i,j)
1381
1382 ENDDO
1383 ENDDO
1384
1385 RETURN
1386 END SUBROUTINE G2T2V
1387
1388 !------------------------------------------------------------------------------
1389 !
1390 SUBROUTINE WEIGTS_CHECK(HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1391 VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1392 IDS,IDE,JDS,JDE,KDS,KDE, &
1393 IMS,IME,JMS,JME,KMS,KME, &
1394 ITS,ITE,JTS,JTE,KTS,KTE )
1395
1396 IMPLICIT NONE
1397 INTEGER, INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1398 IMS,IME,JMS,JME,KMS,KME, &
1399 ITS,ITE,JTS,JTE,KTS,KTE
1400 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1401 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
1402
1403 ! local
1404
1405 REAL , PARAMETER :: EPSI=1.0E-3
1406 INTEGER :: I,J
1407 REAL :: ADDSUM
1408
1409 !-------------------------------------------------------------------------------------
1410
1411 ! DUE TO THE NEED FOR HALO EXCHANGES IN PARALLEL RUNS ONE HAS TO ENSURE CONSISTENT
1412 ! USAGE OF NUMBER OF PROCESSORS BEFORE ANY FURTHER COMPUTATIONS. WE INTRODUCE THIS
1413 ! CHECK FIRST
1414
1415 IF((ITE-ITS) .LE. 5 .OR. (JTE-JTS) .LE. 5)THEN
1416 WRITE(0,*)'ITE-ITS=',ITE-ITS,'JTE-JTS=',JTE-JTS
1417 CALL wrf_error_fatal ('NESTED DOMAIN:PLEASE OPTIMIZE THE NUMBER OF PROCESSES; TRY SQUARES OF NUMBERS')
1418 ENDIF
1419
1420 !
1421 ! NOW CHECK WEIGHTS
1422 !
1423
1424 ADDSUM=0.
1425 DO J = JTS, MIN(JTE,JDE-1)
1426 DO I = ITS, MIN(ITE,IDE-1)
1427 ADDSUM=HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J)
1428 IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
1429 WRITE(0,*)'I=',I,'J=',J,'WEIGHTS=',HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J),1-ADDSUM
1430 CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT MASS POINTS')
1431 ENDIF
1432 ENDDO
1433 ENDDO
1434
1435 ADDSUM=0.
1436 DO J = JTS, MIN(JTE,JDE-1)
1437 DO I = ITS, MIN(ITE,IDE-1)
1438 ADDSUM=VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J)
1439 IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
1440 WRITE(0,*)'I=',I,'J=',J,'WEIGHTS=',VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J),1-ADDSUM
1441 CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT VELOCITY POINTS')
1442 ENDIF
1443 ENDDO
1444 ENDDO
1445
1446 END SUBROUTINE WEIGTS_CHECK
1447
1448 !-----------------------------------------------------------------------------------
1449
1450 SUBROUTINE BOUNDS_CHECK( IIH,JJH,IIV,JJV, &
1451 IPOS,JPOS,SHW, &
1452 IDS,IDE,JDS,JDE,KDS,KDE, & !
1453 IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
1454 ITS,ITE,JTS,JTE,KTS,KTE )
1455
1456 IMPLICIT NONE
1457 INTEGER, INTENT(IN) :: IPOS,JPOS,SHW, &
1458 IDS,IDE,JDS,JDE,KDS,KDE, &
1459 IMS,IME,JMS,JME,KMS,KME, &
1460 ITS,ITE,JTS,JTE,KTS,KTE
1461
1462 INTEGER, DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IIH,JJH,IIV,JJV
1463
1464 ! local variables
1465
1466 INTEGER :: I,J
1467
1468 !*** Gopal - Initial version
1469 !***
1470 !*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
1471 !
1472 !============================================================================
1473
1474 IF(IPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs X-BOUNDARY')
1475 IF(JPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs Y-BOUNDARY')
1476
1477 DO J = JTS, MIN(JTE,JDE-1)
1478 DO I = ITS, MIN(ITE,IDE-1)
1479 IF(IIH(I,J) .EQ. 0)CALL wrf_error_fatal ('IIH=0: SOMETHING IS WRONG')
1480 IF(JJH(I,J) .EQ. 0)CALL wrf_error_fatal ('JJH=0: SOMETHING IS WRONG')
1481 ENDDO
1482 ENDDO
1483
1484 DO J = JTS, MIN(JTE,JDE-1)
1485 DO I = ITS, MIN(ITE,IDE-1)
1486 IF(IIH(I,J) .LT. (IPOS-SHW) .OR. JJH(I,J) .LT. (JPOS-SHW) .OR. &
1487 IIV(I,J) .LT. (IPOS-SHW) .OR. JJV(I,J) .LT. (JPOS-SHW))THEN
1488 WRITE(0,*)I,J,IIH(I,J),IPOS,JJH(I,J),JPOS,SHW
1489 WRITE(0,*)I,J,IIV(I,J),IPOS,JJV(I,J),JPOS,SHW
1490 CALL wrf_error_fatal ('CHECK NESTED DOMAIN BOUNDS: TRY INCREASING STENCIL WIDTH')
1491 ENDIF
1492 ENDDO
1493 ENDDO
1494
1495 END SUBROUTINE BOUNDS_CHECK
1496
1497 !==========================================================================================
1498
1499
1500 SUBROUTINE BASE_STATE_PARENT ( Z3d,Q3d,T3d,PSTD, &
1501 PINT,T,Q,CWM, &
1502 FIS,QS,PD,PDTOP,PTOP, &
1503 ETA1,ETA2, &
1504 DETA1,DETA2, &
1505 IDS,IDE,JDS,JDE,KDS,KDE, &
1506 IMS,IME,JMS,JME,KMS,KME, &
1507 ITS,ITE,JTS,JTE,KTS,KTE )
1508 !
1509
1510 USE MODULE_MODEL_CONSTANTS
1511 IMPLICIT NONE
1512 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
1513 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
1514 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
1515 REAL, INTENT(IN ) :: PDTOP,PTOP
1516 REAL, DIMENSION(KMS:KME), INTENT(IN) :: ETA1,ETA2,DETA1,DETA2
1517 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: FIS,PD,QS
1518 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(IN) :: PINT,T,Q,CWM
1519 REAL, DIMENSION(KMS:KME), INTENT(OUT):: PSTD
1520 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(OUT):: Z3d,Q3d,T3d
1521
1522 ! local
1523
1524 INTEGER,PARAMETER :: JTB=134
1525 INTEGER :: I,J,K,ILOC,JLOC
1526 REAL, PARAMETER :: LAPSR=6.5E-3, GI=1./G,D608=0.608
1527 REAL, PARAMETER :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
1528 REAL, PARAMETER :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
1529 REAL, PARAMETER :: P_REF=103000.
1530 REAL :: A,B,APELP,RTOPP,DZ,ZMID
1531 REAL, DIMENSION(IMS:IME,JMS:JME) :: SLP,TSFC,ZMSLP
1532 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: Z3d_IN
1533 REAL,DIMENSION(JTB) :: PIN,ZIN,Y2,PIO,ZOUT,DUM1,DUM2
1534 REAL,DIMENSION(JTB) :: QIN,QOUT,TIN,TOUT
1535 !--------------------------------------------------------------------------------------
1536
1537 ! CLEAN Z3D ARRAY FIRST
1538
1539 DO J = JTS, MIN(JTE,JDE-1)
1540 DO K=KDS,KDE
1541 DO I = ITS, MIN(ITE,IDE-1)
1542 Z3d(I,K,J)=0.0
1543 T3d(I,K,J)=0.0
1544 Q3d(I,K,J)=0.0
1545 ENDDO
1546 ENDDO
1547 ENDDO
1548
1549
1550 ! DETERMINE THE HEIGHTS ON THE PARENT DOMAIN
1551
1552 DO J = JTS, MIN(JTE,JDE-1)
1553 DO I = ITS, MIN(ITE,IDE-1)
1554 Z3d_IN(I,1,J)=FIS(I,J)*GI
1555 ENDDO
1556 ENDDO
1557
1558 DO J = JTS, MIN(JTE,JDE-1)
1559 DO K = KDS,KDE-1
1560 DO I = ITS, MIN(ITE,IDE-1)
1561 APELP = (PINT(I,K+1,J)+PINT(I,K,J))
1562 ! RTOPP = TRG*T(I,K,J)*(1.0+Q(I,K,J)*P608-CWM(I,K,J))/APELP
1563 RTOPP = TRG*T(I,K,J)*(1.0+Q(I,K,J)*P608)/APELP
1564 DZ = RTOPP*(DETA1(K)*PDTOP+DETA2(K)*PD(I,J)) ! (RTv/P_TOT)*D(P_HYDRO)
1565 Z3d_IN(I,K+1,J) = Z3d_IN(I,K,J) + DZ
1566 ! IF(I==2 .AND. J==2)WRITE(0,*)'INSIDE BASE_STATE',K,T(I,K,J)
1567 ENDDO
1568 ENDDO
1569 ENDDO
1570
1571
1572 ! CONSTRUCT STANDARD ISOBARIC SURFACES
1573
1574 DO K=KDS,KDE ! target points in model interface levels (pint)
1575 PSTD(K) = ETA1(K)*PDTOP + ETA2(K)*(P_REF -PDTOP - PTOP) + PTOP
1576 ENDDO
1577
1578 ! DETERMINE THE MSLP USE THAT TO CREATE HEIGHTS AT 1000. mb LEVEL. THESE HEIGHTS
1579 ! MAY ONLY BE USED IN VERTICAL INTERPOLATION TO ISOBARIC SURFACES WHICH ARE LOCATED
1580 ! BELOW GROUND LEVEL.
1581
1582 DO J = JTS, MIN(JTE,JDE-1)
1583 DO I = ITS, MIN(ITE,IDE-1)
1584 TSFC(I,J) = T(I,1,J)*(1.+D608*Q(I,1,J)) + LAPSR*(Z3d_IN(I,1,J)+Z3d_IN(I,2,J))*0.5
1585 A = LAPSR*Z3d_IN(I,1,J)/TSFC(I,J)
1586 SLP(I,J) = PINT(I,1,J)*(1-A)**COEF2 ! sea level pressure
1587 B = (PSTD(1)/SLP(I,J))**COEF3
1588 ZMSLP(I,J)= TSFC(I,J)*LAPSI*(1.0 - B) ! Height at 1000. mb level
1589 ENDDO
1590 ENDDO
1591
1592 ! INTERPOLATE Z3d_IN TO STANDARD PRESSURE INTERFACES. FOR LEVELS BELOW
1593 ! GROUND USE ZMSLP(I,J)
1594
1595 DO J = JTS, MIN(JTE,JDE-1)
1596 DO I = ITS, MIN(ITE,IDE-1)
1597 !
1598 ! clean local array before use of spline
1599
1600 PIN=0.;ZIN=0.;Y2=0;PIO=0.;ZOUT=0.;DUM1=0.;DUM2=0.
1601
1602 DO K=KDS,KDE ! inputs at model interfaces
1603 PIN(K) = PINT(I,KDE-K+1,J)
1604 ZIN(K) = Z3d_IN(I,KDE-K+1,J)
1605 ENDDO
1606
1607 IF(PINT(I,1,J) .LE. PSTD(1))THEN
1608 PIN(KDE) = PSTD(1)
1609 ZIN(KDE) = ZMSLP(I,J)
1610 ENDIF
1611 !
1612 Y2(1 )=0.
1613 Y2(KDE)=0.
1614 !
1615 DO K=KDS,KDE
1616 PIO(K)=PSTD(K)
1617 ENDDO
1618 !
1619 CALL SPLINE1(I,J,JTB,KDE,PIN,ZIN,Y2,KDE,PIO,ZOUT,DUM1,DUM2) ! interpolate
1620 !
1621
1622 DO K=KDS,KDE ! inputs at model interfaces
1623 Z3d(I,K,J)=ZOUT(K)
1624 ENDDO
1625
1626 ENDDO
1627 ENDDO
1628 !
1629 ! INTERPOLATE TEMPERATURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW
1630 ! GROUND USE A LAPSE RATE ATMOSPHERE
1631 !
1632 DO J = JTS, MIN(JTE,JDE-1)
1633 DO I = ITS, MIN(ITE,IDE-1)
1634 !
1635 ! clean local array before use of spline or linear interpolation
1636 !
1637 PIN=0.;TIN=0.;Y2=0;PIO=0.;TOUT=0.;DUM1=0.;DUM2=0.
1638
1639 DO K=KDS+1,KDE ! inputs at model levels
1640 PIN(K-1) = EXP((ALOG(PINT(I,KDE-K+1,J))+ALOG(PINT(I,KDE-K+2,J)))*0.5)
1641 TIN(K-1) = T(I,KDE-K+1,J)
1642 ENDDO
1643
1644 IF(PINT(I,1,J) .LE. PSTD(1))THEN
1645 PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
1646 ZMID = 0.5*(Z3d_IN(I,1,J)+Z3d_IN(I,2,J))
1647 TIN(KDE-1) = T(I,1,J) + LAPSR*(ZMID-ZMSLP(I,J))
1648 ENDIF
1649 !
1650 Y2(1 )=0.
1651 Y2(KDE-1)=0.
1652 !
1653 DO K=KDS,KDE-1
1654 PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
1655 ENDDO
1656
1657 CALL SPLINE1(I,J,JTB,KDE-1,PIN,TIN,Y2,KDE-1,PIO,TOUT,DUM1,DUM2) ! interpolate
1658
1659
1660 DO K=KDS,KDE-1 ! inputs at model levels
1661 T3d(I,K,J)=TOUT(K)
1662 ENDDO
1663
1664 ENDDO
1665 ENDDO
1666
1667 !
1668 ! INTERPOLATE MOISTURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW
1669 ! GROUND USE THE SURFACE MOISTURE
1670 !
1671 DO J = JTS, MIN(JTE,JDE-1)
1672 DO I = ITS, MIN(ITE,IDE-1)
1673 !
1674 ! clean local array before use of spline or linear interpolation
1675
1676
1677 PIN=0.;QIN=0.;Y2=0;PIO=0.;QOUT=0.;DUM1=0.;DUM2=0.
1678
1679 DO K=KDS+1,KDE ! inputs at model levels
1680 PIN(K-1) = EXP((ALOG(PINT(I,KDE-K+1,J))+ALOG(PINT(I,KDE-K+2,J)))*0.5)
1681 QIN(K-1) = Q(I,KDE-K+1,J)
1682 ENDDO
1683
1684 IF(PINT(I,1,J) .LE. PSTD(1))THEN
1685 PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
1686 ! QIN(KDE-1) = QS(I,J)
1687 ENDIF
1688
1689 Y2(1 )=0.
1690 Y2(KDE-1)=0.
1691 !
1692 DO K=KDS,KDE-1
1693 PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
1694 ENDDO
1695
1696 CALL SPLINE1(I,J,JTB,KDE-1,PIN,QIN,Y2,KDE-1,PIO,QOUT,DUM1,DUM2) ! interpolate
1697
1698 DO K=KDS,KDE-1 ! inputs at model levels
1699 Q3d(I,K,J)=QOUT(K)
1700 ENDDO
1701
1702 ENDDO
1703 ENDDO
1704
1705 END SUBROUTINE BASE_STATE_PARENT
1706 !=============================================================================
1707 SUBROUTINE SPLINE1(I,J,JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
1708 !
1709 ! ******************************************************************
1710 ! * *
1711 ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
1712 ! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
1713 ! * *
1714 ! * PROGRAMER Z. JANJIC *
1715 ! * *
1716 ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
1717 ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
1718 ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
1719 ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
1720 ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
1721 ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
1722 ! * SPECIFIED. *
1723 ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
1724 ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
1725 ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
1726 ! * AND LE XOLD(NOLD). *
1727 ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
1728 ! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
1729 ! * *
1730 ! ******************************************************************
1731 !---------------------------------------------------------------------
1732 IMPLICIT NONE
1733 !---------------------------------------------------------------------
1734 INTEGER,INTENT(IN) :: I,J,JTBX,NNEW,NOLD
1735 REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
1736 REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
1737 REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
1738 !
1739 INTEGER :: II,JJ,K,K1,K2,KOLD,NOLDM1
1740 REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
1741 ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
1742 !---------------------------------------------------------------------
1743
1744 ! debug
1745
1746 II=9999 !67 !35 !50 !4
1747 JJ=9999 !31 !73 !115 !192
1748 IF(I.eq.II.and.J.eq.JJ)THEN
1749 WRITE(0,*)'DEBUG in SPLINE1:HSO= ',xnew(1:nold)
1750 DO K=1,NOLD
1751 WRITE(0,*)'DEBUG in SPLINE1:L,ZETAI,PINTI= ' &
1752 ,K,YOLD(K),XOLD(K)
1753 ENDDO
1754 ENDIF
1755
1756 !
1757 NOLDM1=NOLD-1
1758 !
1759 DXL=XOLD(2)-XOLD(1)
1760 DXR=XOLD(3)-XOLD(2)
1761 DYDXL=(YOLD(2)-YOLD(1))/DXL
1762 DYDXR=(YOLD(3)-YOLD(2))/DXR
1763 RTDXC=0.5/(DXL+DXR)
1764 !
1765 P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
1766 Q(1)=-RTDXC*DXR
1767 !
1768 IF(NOLD.EQ.3)GO TO 150
1769 !---------------------------------------------------------------------
1770 K=3
1771 !
1772 100 DXL=DXR
1773 DYDXL=DYDXR
1774 DXR=XOLD(K+1)-XOLD(K)
1775 DYDXR=(YOLD(K+1)-YOLD(K))/DXR
1776 DXC=DXL+DXR
1777 DEN=1./(DXL*Q(K-2)+DXC+DXC)
1778 !
1779 P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
1780 Q(K-1)=-DEN*DXR
1781 !
1782 K=K+1
1783 IF(K.LT.NOLD)GO TO 100
1784 !-----------------------------------------------------------------------
1785 150 K=NOLDM1
1786 !
1787 200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
1788 !
1789 K=K-1
1790 IF(K.GT.1)GO TO 200
1791 !-----------------------------------------------------------------------
1792 K1=1
1793 !
1794 300 XK=XNEW(K1)
1795 !
1796 DO 400 K2=2,NOLD
1797 !
1798 IF(XOLD(K2).GT.XK)THEN
1799 KOLD=K2-1
1800 GO TO 450
1801 ENDIF
1802 !
1803 400 CONTINUE
1804 !
1805 YNEW(K1)=YOLD(NOLD)
1806 GO TO 600
1807 !
1808 450 IF(K1.EQ.1)GO TO 500
1809 IF(K.EQ.KOLD)GO TO 550
1810 !
1811 500 K=KOLD
1812 !
1813 Y2K=Y2(K)
1814 Y2KP1=Y2(K+1)
1815 DX=XOLD(K+1)-XOLD(K)
1816 RDX=1./DX
1817 !
1818 AK=.1666667*RDX*(Y2KP1-Y2K)
1819 BK=0.5*Y2K
1820 CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
1821 !
1822 550 X=XK-XOLD(K)
1823 XSQ=X*X
1824 !
1825 YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)
1826
1827 ! debug
1828
1829 if(i.eq.ii.and.j.eq.jj)then
1830 write(0,*) 'DEBUG:: k1,xnew(k1),ynew(k1): ', k1,xnew(k1),ynew(k1)
1831 endif
1832
1833 !
1834 600 K1=K1+1
1835 IF(K1.LE.NNEW)GO TO 300
1836
1837 RETURN
1838 END SUBROUTINE SPLINE1
1839 !---------------------------------------------------------------------
1840
1841 SUBROUTINE NEST_TERRAIN ( nest, config_flags )
1842
1843 USE module_domain
1844 USE module_configure
1845 USE module_timing
1846
1847 USE wrfsi_static
1848
1849 IMPLICIT NONE
1850
1851 TYPE(domain) , POINTER :: nest
1852 TYPE(grid_config_rec_type) , INTENT(IN) :: config_flags
1853
1854 !
1855 ! Local variables
1856 !
1857
1858 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
1859 INTEGER :: ids,ide,jds,jde,kds,kde
1860 INTEGER :: ims,ime,jms,jme,kms,kme
1861 INTEGER :: its,ite,jts,jte,kts,kte
1862 INTEGER :: i_parent_start, j_parent_start
1863 INTEGER :: parent_grid_ratio
1864 INTEGER :: n,i,j,ii,jj,nnxp,nnyp
1865 INTEGER :: i_start,j_start,level
1866 REAL, ALLOCATABLE, DIMENSION(:,:) :: data1 ! for highres topo
1867 REAL, ALLOCATABLE, DIMENSION(:,:) :: avc_big, lnd_big, lah_big, loh_big
1868 REAL, ALLOCATABLE, DIMENSION(:,:) :: avc_nest, lnd_nest, lah_nest, loh_nest
1869 INTEGER :: im_big, jm_big, i_add
1870 INTEGER :: im, jm
1871 CHARACTER(LEN=6) :: nestpath
1872
1873 integer :: input_type
1874 character(len=128) :: input_fname
1875 character (len=32) :: cname
1876 integer :: ndim
1877 character (len=3) :: memorder
1878 character (len=32) :: stagger
1879 integer, dimension(3) :: domain_start, domain_end
1880 integer :: wrftype
1881 character (len=128), dimension(3) :: dimnames
1882
1883 integer :: istatus
1884 integer :: handle
1885 integer :: comm_1, comm_2
1886
1887 real, allocatable, dimension(:,:,:) :: real_domain
1888
1889 character (len=10), dimension(4) :: name = (/ "XLAT_M ", &
1890 "XLONG_M ", &
1891 "LANDMASK ", &
1892 "HGT_M " /)
1893
1894 integer, parameter :: IO_BIN=1, IO_NET=2
1895
1896 integer :: io_form_input
1897
1898 write(0,*)"in NEST_TERRAIN config_flags%io_form_input = ", config_flags%io_form_input
1899 write(0,*)"in NEST_TERRAIN config_flags%auxinput1_inname = ", config_flags%auxinput1_inname
1900 io_form_input = config_flags%io_form_input
1901 if (config_flags%auxinput1_inname(1:7) == "met_nmm") then
1902 input_type = 2
1903 else
1904 input_type = 1
1905 end if
1906
1907 !----------------------------------------------------------------------------------
1908
1909 IDS = nest%sd31
1910 IDE = nest%ed31
1911 KDS = nest%sd32
1912 KDE = nest%ed32
1913 JDS = nest%sd33
1914 JDE = nest%ed33
1915
1916 IMS = nest%sm31
1917 IME = nest%em31
1918 KMS = nest%sm32
1919 KME = nest%em32
1920 JMS = nest%sm33
1921 JME = nest%em33
1922
1923 ITS = nest%sp31
1924 ITE = nest%ep31
1925 KTS = nest%sp32
1926 KTE = nest%ep32
1927 JTS = nest%sp33
1928 JTE = nest%ep33
1929
1930 i_parent_start = nest%i_parent_start
1931 j_parent_start = nest%j_parent_start
1932 parent_grid_ratio = nest%parent_grid_ratio
1933
1934 NNXP=IDE-1
1935 NNYP=JDE-1
1936
1937 ALLOCATE(DATA1(1:NNXP,1:NNYP))
1938 !
1939 !
1940 !--- Read in high resolution topography
1941 !
1942 IF ( wrf_dm_on_monitor() ) THEN ! first assign a status
1943 !
1944 ! This part of the code is Dusan's doing. Extended by gopal for multiple nest (Feb 19,2005)
1945 !
1946 call find_ijstart_level (nest,i_start,j_start,level)
1947 write(0,*)" nest%id =", nest%id , " i_start,j_start,level =", i_start,j_start,level
1948
1949 write(nestpath,"(a4,i1,a1)") 'nest',level,'/'
1950
1951 if ( level > 0 ) then
1952
1953 if (input_type == 1) then
1954 !
1955 ! SI version of the static file
1956 !
1957 CALL get_wrfsi_static_dims(nestpath, im_big, jm_big)
1958 ALLOCATE (avc_big(im_big,jm_big))
1959 ALLOCATE (lnd_big(im_big,jm_big))
1960 ALLOCATE (lah_big(im_big,jm_big))
1961 ALLOCATE (loh_big(im_big,jm_big))
1962 CALL get_wrfsi_static_2d(nestpath, 'avc', avc_big)
1963 CALL get_wrfsi_static_2d(nestpath, 'lnd', lnd_big)
1964 CALL get_wrfsi_static_2d(nestpath, 'lah', lah_big)
1965 CALL get_wrfsi_static_2d(nestpath, 'loh', loh_big)
1966
1967 else if (input_type == 2) then
1968 !
1969 ! WPS version of the static file
1970 !
1971
1972 #ifdef INTIO
1973 if (io_form_input == IO_BIN) write(input_fname,"(A,I2.2,A)") "geo_nmm_nest.l",level,".int"
1974 #endif
1975 #ifdef NETCDF
1976 if (io_form_input == IO_NET) write(input_fname,"(A,I2.2,A)") "geo_nmm_nest.l",level,".nc"
1977 #endif
1978
1979 comm_1 = 1
1980 comm_2 = 1
1981
1982 #ifdef INTIO
1983 if (io_form_input == IO_BIN) &
1984 call ext_int_open_for_read(trim(input_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
1985 #endif
1986 #ifdef NETCDF
1987 if (io_form_input == IO_NET) &
1988 call ext_ncd_open_for_read(trim(input_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
1989 #endif
1990 if (istatus /= 0) CALL wrf_error_fatal('NEST_TERRAIN error after ext_XXX_open_for_read '//trim(input_fname))
1991
1992
1993 do n=1,4
1994
1995 cname = name(n)
1996
1997 domain_start = 1
1998 domain_end = 1
1999 #ifdef INTIO
2000 if (io_form_input == IO_BIN) &
2001 call ext_int_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus)
2002 #endif
2003 #ifdef NETCDF
2004 if (io_form_input == IO_NET) &
2005 call ext_ncd_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus)
2006 #endif
2007 print *, "istatus=", istatus
2008 print *, "ndim=", ndim
2009 print *, "memorder=", memorder
2010 print *, "stagger=", stagger
2011 print *, "domain_start=", domain_start
2012 print *, "domain_end=", domain_end
2013 print *, "wrftype=", wrftype
2014
2015
2016 if (allocated(real_domain)) deallocate(real_domain)
2017 allocate(real_domain(domain_start(1):domain_end(1), domain_start(2):domain_end(2), domain_start(3):domain_end(3)))
2018
2019 #ifdef INTIO
2020 if (io_form_input == IO_BIN) then
2021 call ext_int_read_field(handle, '0000-00-00_00:00:00', cname, real_domain, wrftype, &
2022 1, 1, 0, memorder, stagger, &
2023 dimnames, domain_start, domain_end, domain_start, domain_end, &
2024 domain_start, domain_end, istatus)
2025 end if
2026 #endif
2027 #ifdef NETCDF
2028 if (io_form_input == IO_NET) then
2029 call ext_ncd_read_field(handle, '0000-00-00_00:00:00', cname, real_domain, wrftype, &
2030 1, 1, 0, memorder, stagger, &
2031 dimnames, domain_start, domain_end, domain_start, domain_end, &
2032 domain_start, domain_end, istatus)
2033 end if
2034 #endif
2035 print *, "istatus=", istatus
2036
2037 im_big = domain_end(1)
2038 jm_big = domain_end(2)
2039 if (cname(1:10) == "XLAT_M ") then
2040 ALLOCATE (lah_big(im_big,jm_big))
2041 do j=1,jm_big
2042 do i=1,im_big
2043 lah_big(i,j) = real_domain(i,j,1)
2044 end do
2045 end do
2046 else if (cname(1:10) == "XLONG_M ") then
2047 ALLOCATE (loh_big(im_big,jm_big))
2048 do j=1,jm_big
2049 do i=1,im_big
2050 loh_big(i,j) = real_domain(i,j,1)
2051 end do
2052 end do
2053 else if (cname(1:10) == "LANDMASK ") then
2054 ALLOCATE (lnd_big(im_big,jm_big))
2055 do j=1,jm_big
2056 do i=1,im_big
2057 lnd_big(i,j) = real_domain(i,j,1)
2058 end do
2059 end do
2060 else if (cname(1:10) == "HGT_M ") then
2061 ALLOCATE (avc_big(im_big,jm_big))
2062 do j=1,jm_big
2063 do i=1,im_big
2064 avc_big(i,j) = real_domain(i,j,1)
2065 end do
2066 end do
2067 end if
2068
2069 end do
2070
2071 #ifdef INTIO
2072 if (io_form_input == IO_BIN) then
2073 call ext_int_ioclose(handle, istatus)
2074 end if
2075 #endif
2076 #ifdef NETCDF
2077 if (io_form_input == IO_NET) then
2078 call ext_ncd_ioclose(handle, istatus)
2079 end if
2080 #endif
2081
2082 else
2083 CALL wrf_error_fatal('NEST_TERRAIN wrong input_type')
2084 end if
2085
2086 else
2087 CALL wrf_error_fatal('this routine NEST_TERRAIN should nou be called for top-level domain')
2088 end if
2089
2090 ! select subdomain from big fine grid
2091
2092 im = NNXP
2093 jm = NNYP
2094
2095 ALLOCATE (avc_nest(im,jm))
2096 ALLOCATE (lnd_nest(im,jm))
2097 ALLOCATE (lah_nest(im,jm))
2098 ALLOCATE (loh_nest(im,jm))
2099
2100 i_add = mod(j_start+1,2)
2101 DO j=1,jm
2102 DO i=1,im
2103 avc_nest(i,j) = avc_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2104 lnd_nest(i,j) = lnd_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2105 lah_nest(i,j) = lah_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2106 loh_nest(i,j) = loh_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2107 END DO
2108 END DO
2109
2110 WRITE(0,*)'SOME MATCHING TEST i_parent_start, j_parent_start',i_parent_start,j_parent_start
2111 WRITE(0,*)'WRFSI LAT COMPUTED LAT'
2112 WRITE(0,*)lah_nest(1,1),nest%nmm_hlat(1,1)
2113 WRITE(0,*)'WRFSI LON COMPUTED LON'
2114 WRITE(0,*)loh_nest(1,1),nest%nmm_hlon(1,1)
2115
2116 IF(ABS(lah_nest(1,1)-nest%nmm_hlat(1,1)) .GE. 0.5 .OR. &
2117 ABS(loh_nest(1,1)-nest%nmm_hlon(1,1)) .GE. 0.5)THEN
2118 WRITE(0,*)'CHECK WRFSI CONFIGURATION AND INPUT HIGH RESOLUTION TOPOGRAPHY AND/OR GRID RATIO'
2119 CALL wrf_error_fatal('LATLON MISMATCH: ERROR READING static FILE FOR THE NEST')
2120 ENDIF
2121
2122 call smdhld(im,jm,avc_nest,1.0-lnd_nest,12,12)
2123
2124 !-------------4-point averaging of mountains along inner boundary-------
2125
2126 do i=1,im-1
2127 avc_nest(i,2)=0.25*(avc_nest(i,1)+avc_nest(i+1,1)+ &
2128 & avc_nest(i,3)+avc_nest(i+1,3))
2129 enddo
2130
2131 do i=1,im-1
2132 avc_nest(i,jm-1)=0.25*(avc_nest(i,jm-2)+avc_nest(i+1,jm-2)+ &
2133 & avc_nest(i,jm)+avc_nest(i+1,jm))
2134 enddo
2135
2136 do j=4,jm-3,2
2137 avc_nest(1,j)=0.25*(avc_nest(1,j-1)+avc_nest(2,j-1)+ &
2138 & avc_nest(1,j+1)+avc_nest(2,j+1))
2139 enddo
2140
2141 do j=4,jm-3,2
2142 avc_nest(im,j)=0.25*(avc_nest(im-1,j-1)+avc_nest(im,j-1)+ &
2143 & avc_nest(im-1,j+1)+avc_nest(im,j+1))
2144 enddo
2145
2146 DO J = 1,NNYP
2147 DO I = 1,NNXP
2148 DATA1(I,J) = 9.81*avc_nest(I,J)
2149 ENDDO
2150 ENDDO
2151
2152 DEALLOCATE (avc_big,lnd_big)
2153 DEALLOCATE (avc_nest,lnd_nest)
2154 !
2155 ENDIF
2156
2157 CALL wrf_dm_bcast_bytes (DATA1,NNXP*NNYP*RWORDSIZE)
2158
2159 DO J=JDS,JDE
2160 DO I =IDS,IDE
2161 IF(I.GE.ITS .AND. I .LE. MIN(ide-1,ite) .AND. J.GE.JTS .AND. J .LE. MIN(jde-1,jte))THEN
2162 nest%nmm_hres_fis(I,J)=DATA1(I,J)
2163 ENDIF
2164 ENDDO
2165 ENDDO
2166
2167 DEALLOCATE(DATA1)
2168 WRITE(0,*)'end of NEST_TERRAIN'
2169
2170 END SUBROUTINE NEST_TERRAIN
2171 !===========================================================================================
2172
2173
2174 SUBROUTINE med_init_domain_constants_nmm ( parent, nest) !, config_flags)
2175 ! Driver layer
2176 USE module_domain
2177 USE module_configure
2178 USE module_timing
2179 IMPLICIT NONE
2180 TYPE(domain) , POINTER :: parent, nest, grid
2181 !
2182 #ifdef DEREF_KLUDGE
2183 INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
2184 INTEGER :: sm31x, em31x, sm32x, em32x, sm33x, em33x
2185 INTEGER :: sm31y, em31y, sm32y, em32y, sm33y, em33y
2186 #endif
2187 !
2188 INTERFACE
2189 SUBROUTINE med_initialize_nest_nmm ( grid &
2190 !
2191 # include <nmm_dummy_args.inc>
2192 !
2193 )
2194 USE module_domain
2195 USE module_configure
2196 USE module_timing
2197 IMPLICIT NONE
2198 TYPE(domain) , POINTER :: grid
2199 #include <nmm_dummy_decl.inc>
2200 END SUBROUTINE med_initialize_nest_nmm
2201 END INTERFACE
2202
2203 !------------------------------------------------------------------------------
2204 ! PURPOSE:
2205 ! - initialize some data, mainly 2D & 3D nmm arrays very similar to
2206 ! those done in ./dyn_nmm/module_initialize_real.F
2207 !-----------------------------------------------------------------------------
2208 !
2209
2210 grid => nest
2211
2212 #ifdef DEREF_KLUDGE
2213 sm31 = grid%sm31
2214 em31 = grid%em31
2215 sm32 = grid%sm32
2216 em32 = grid%em32
2217 sm33 = grid%sm33
2218 em33 = grid%em33
2219 sm31x = grid%sm31x
2220 em31x = grid%em31x
2221 sm32x = grid%sm32x
2222 em32x = grid%em32x
2223 sm33x = grid%sm33x
2224 em33x = grid%em33x
2225 sm31y = grid%sm31y
2226 em31y = grid%em31y
2227 sm32y = grid%sm32y
2228 em32y = grid%em32y
2229 sm33y = grid%sm33y
2230 em33y = grid%em33y
2231 #endif
2232
2233 CALL med_initialize_nest_nmm( grid &
2234 !
2235 # include <nmm_actual_args.inc>
2236 !
2237 )
2238
2239 END SUBROUTINE med_init_domain_constants_nmm
2240
2241 SUBROUTINE med_initialize_nest_nmm( grid &
2242 !
2243 # include <nmm_dummy_args.inc>
2244 !
2245 )
2246
2247 USE module_domain
2248 USE module_configure
2249 USE module_timing
2250 IMPLICIT NONE
2251
2252 ! Local domain indices and counters.
2253
2254 INTEGER :: ids, ide, jds, jde, kds, kde, &
2255 ims, ime, jms, jme, kms, kme, &
2256 its, ite, jts, jte, kts, kte, &
2257 i, j, k, nnxp, nnyp
2258
2259 TYPE(domain) , POINTER :: grid
2260
2261 ! Local data
2262
2263 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
2264 INTEGER :: KHH,KVH,JAM,JA,IHL, IHH, L
2265 INTEGER :: II,JJ,ISRCH,ISUM
2266 INTEGER, ALLOCATABLE, DIMENSION(:) :: KHL2,KVL2,KHH2,KVH2,KHLA,KHHA,KVLA,KVHA
2267 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
2268 !
2269 REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
2270 REAL(KIND=KNUM) :: STPH,CTPH,TDLM,TDPH,FP,TPH,TLM,TLM0
2271 REAL :: TPH0D,TLM0D,ANBI,TSPH,DTAD,DTCF,DT
2272 REAL :: ACDT,CDDAMP,DXP
2273 REAL :: WBD,SBD,WBI,SBI,EBI
2274 REAL :: DY_NMM0
2275 REAL :: RSNOW,SNOFAC
2276 REAL, ALLOCATABLE, DIMENSION(:) :: DXJ,WPDARJ,CPGFUJ,CURVJ, &
2277 FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
2278 HDACJ,DDMPUJ,DDMPVJ
2279 !
2280 REAL, PARAMETER:: SALP=2.60
2281 REAL, PARAMETER:: SNUP=0.040
2282 REAL, PARAMETER:: W_NMM=0.08
2283 REAL, PARAMETER:: COAC=0.75
2284 REAL, PARAMETER:: CODAMP=6.4
2285 REAL, PARAMETER:: TWOM=.00014584
2286 REAL, PARAMETER:: CP=1004.6
2287 REAL, PARAMETER:: DFC=1.0
2288 REAL, PARAMETER:: DDFC=1.0
2289 REAL, PARAMETER:: ROI=916.6
2290 REAL, PARAMETER:: R=287.04
2291 REAL, PARAMETER:: CI=2060.0
2292 REAL, PARAMETER:: ROS=1500.
2293 REAL, PARAMETER:: CS=1339.2
2294 REAL, PARAMETER:: DS=0.050
2295 REAL, PARAMETER:: AKS=.0000005
2296 REAL, PARAMETER:: DZG=2.85
2297 REAL, PARAMETER:: DI=.1000
2298 REAL, PARAMETER:: AKI=0.000001075
2299 REAL, PARAMETER:: DZI=2.0
2300 REAL, PARAMETER:: THL=210.
2301 REAL, PARAMETER:: PLQ=70000.
2302 REAL, PARAMETER:: ERAD=6371200.
2303 REAL, PARAMETER:: DTR=0.01745329
2304
2305 ! Definitions of dummy arguments to solve
2306 #include <nmm_dummy_decl.inc>
2307
2308 #ifdef DEREF_KLUDGE
2309 INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
2310 INTEGER :: sm31x, em31x, sm32x, em32x, sm33x, em33x
2311 INTEGER :: sm31y, em31y, sm32y, em32y, sm33y, em33y
2312 #endif
2313
2314 #ifdef DEREF_KLUDGE
2315 sm31 = grid%sm31
2316 em31 = grid%em31
2317 sm32 = grid%sm32
2318 em32 = grid%em32
2319 sm33 = grid%sm33
2320 em33 = grid%em33
2321 sm31x = grid%sm31x
2322 em31x = grid%em31x
2323 sm32x = grid%sm32x
2324 em32x = grid%em32x
2325 sm33x = grid%sm33x
2326 em33x = grid%em33x
2327 sm31y = grid%sm31y
2328 em31y = grid%em31y
2329 sm32y = grid%sm32y
2330 em32y = grid%em32y
2331 sm33y = grid%sm33y
2332 em33y = grid%em33y
2333 #endif
2334
2335 #define COPY_IN
2336 #include <nmm_scalar_derefs.inc>
2337 #ifdef DM_PARALLEL
2338 # include <nmm_data_calls.inc>
2339 #endif
2340
2341 CALL get_ijk_from_grid ( grid , &
2342 ids, ide, jds, jde, kds, kde, &
2343 ims, ime, jms, jme, kms, kme, &
2344 its, ite, jts, jte, kts, kte )
2345
2346
2347 !=================================================================================
2348 !
2349 !
2350
2351 DT=grid%dt !float(TIME_STEP)/parent_time_step_ratio
2352 NNXP=min(ITE,IDE-1)
2353 NNYP=min(JTE,JDE-1)
2354 JAM=6+2*((JDE-1)-10) ! this should be the fix instead of JAM=6+2*(NNYP-10)
2355
2356 WRITE(0,*)'TIME STEP ON DOMAIN',grid%id,'==',dt
2357
2358 !
2359 ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
2360 ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
2361 ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),FADJ(JTS:NNYP))
2362 ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
2363 ALLOCATE(KHLA(JAM),KHHA(JAM))
2364 ALLOCATE(KVLA(JAM),KVHA(JAM))
2365
2366 ! INITIALIZE SOME LAND/WATER SURFACE DATA ON THE BASIS OF INPUTS: SM, XICE, WEASD,
2367 ! INTERPOLATED FROM MOTHER (WRFSI) DOMAIN. THIS PART OF THE CODE HAS TO BE REVISITED
2368 ! LATER ON
2369
2370 DO J = JTS, MIN(JTE,JDE-1)
2371 DO I = ITS, MIN(ITE,IDE-1)
2372 IF(SM(I,J).GT.0.9) THEN ! OVER WATER SURFACE
2373 !
2374 IF (XICE(I,J) .gt. 0)THEN ! XICE: SI INPUT ON PARENT, INTERPOLATED ONTO NEST
2375 SI(I,J)=1.0 ! INITIALIZE SI BASED ON XICE FROM INTERPOLATED INPUT
2376 ENDIF
2377 !
2378 EPSR(I,J)= 0.97 ! VALID OVER SEA SURFACE
2379 GFFC(I,J)= 0.
2380 ALBEDO(I,J)=.06
2381 ALBASE(I,J)=.06
2382 !
2383 IF(SI (I,J) .GT. 0.)THEN ! VALID OVER SEA-ICE
2384 SM(I,J)=0.
2385 SI(I,J)=0. !
2386 SICE(I,J)=1.
2387 GFFC(I,J)=0. ! just leave zero as irrelevant
2388 ALBEDO(I,J)=.60 ! DEFINE ALBEDO
2389 ALBASE(I,J)=.60
2390 ENDIF
2391 !
2392 ELSE ! OVER LAND SURFACE
2393 !
2394 SI(I,J)=5.0*WEASD(I,J)/1000. ! SNOW WATER EQ (mm) OBTAINED FROM PARENT (SI) IS INTERPOLATED
2395 EPSR(I,J)=1.0 ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
2396 GFFC(I,J)=0.0 ! just leave zero as irrelevant
2397 SICE(I,J)=0. ! SEA ICE
2398 SNO(I,J)=SI(I,J)*.20 ! LAND-SNOW COVER
2399 !
2400 ENDIF
2401 !
2402 ENDDO
2403 ENDDO
2404
2405 ! This may just be a fix and may need some Registry related changes, later on
2406
2407 DO J = JTS, MIN(JTE,JDE-1)
2408 DO I = ITS, MIN(ITE,IDE-1)
2409 VEGFRA(I,J)=VEGFRC(I,J)
2410 ENDDO
2411 ENDDO
2412
2413 ! DETERMINE ALBEDO OVER LAND ON THE BASIS OF INPUTS: SM, ALBASE, MXSNAL & VEGFRA
2414 ! INTERPOLATED FROM MOTHER (WRFSI) DOMAIN
2415
2416
2417 DO J = JTS, MIN(JTE,JDE-1)
2418 DO I = ITS, MIN(ITE,IDE-1)
2419
2420 IF(SM(I,J).LT.0.9.AND.SICE(I,J).LT.0.9) THEN
2421 !
2422 IF ( (SNO(I,J) .EQ. 0.0) .OR. & ! SNOWFREE ALBEDO
2423 (ALBASE(I,J) .GE. MXSNAL(I,J) ) ) THEN
2424 ALBEDO(I,J) = ALBASE(I,J)
2425 ELSE
2426 IF (SNO(I,J) .LT. SNUP) THEN ! MODIFY ALBEDO IF SNOWCOVER:
2427 RSNOW = SNO(I,J)/SNUP ! BELOW SNOWDEPTH THRESHOLD
2428 SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
2429 ELSE
2430 SNOFAC = 1.0 ! ABOVE SNOWDEPTH THRESHOLD
2431 ENDIF
2432 ALBEDO(I,J) = ALBASE(I,J) &
2433 + (1.0-VEGFRA(I,J))*SNOFAC*(MXSNAL(I,J)-ALBASE(I,J))
2434 ENDIF
2435 !
2436 END IF
2437
2438 SI(I,J)=5.0*WEASD(I,J)
2439 SNO(I,J)=WEASD(I,J)
2440 ! this block probably superfluous. Meant to guarantee land/sea agreement
2441
2442 IF (SM(I,J) .gt. 0.5)THEN
2443 landmask(I,J)=0.0
2444 ELSE
2445 landmask(I,J)=1.0
2446 ENDIF
2447
2448 IF (SICE(I,J) .eq. 1.0) then !!!! change vegtyp and sltyp to fit seaice (desireable??)
2449 ISLTYP(I,J)=16
2450 IVGTYP(I,J)=24
2451 ENDIF
2452
2453 ENDDO
2454 ENDDO
2455
2456 ! Check land water interface
2457
2458 DO J = JTS, MIN(JTE,JDE-1)
2459 DO I = ITS,MIN(ITE,IDE-1)
2460 IF(SM(I,J).GT.0.9 .AND. VEGFRA(I,J) .NE. 0) THEN
2461 WRITE(20,*)'PROBLEM AT THE LAND-WATER INTERFACE:',I,J,SM(I-1,J),VEGFRA(I-1,j),SM(I,J),VEGFRA(I,J)
2462 ENDIF
2463 !
2464 IF(SM(I,J).GT.0.9 .AND. NMM_TSK(I,J) .NE. 0) THEN
2465 WRITE(20,*)'PROBLEM AT THE LAND-WATER INTERFACE:',I,J,SM(I-1,J),NMM_TSK(I-1,J),SM(I,J),NMM_TSK(I,J)
2466 ENDIF
2467 ENDDO
2468 ENDDO
2469
2470
2471 ! hardwire root depth for time being
2472
2473 RTDPTH=0.
2474 RTDPTH(1)=0.1
2475 RTDPTH(2)=0.3
2476 RTDPTH(3)=0.6
2477
2478 ! hardwire soil depth for time being
2479
2480 SLDPTH=0.
2481 SLDPTH(1)=0.1
2482 SLDPTH(2)=0.3
2483 SLDPTH(3)=0.6
2484 SLDPTH(4)=1.0
2485
2486 !----------- END OF LAND SURFACE INITIALIZATION -------------------------------------
2487 !
2488 ! INITIALIZE 3D HEIGHT MASK AND VELOCITY FIELDS (HTM AND VTM),
2489 ! AND LOWEST ABV GROUND LEVEL (LMH AND LMV) AND RECIPROCAL
2490 ! ETAS (RES) OVER THE NESTED DOMAIN
2491
2492
2493 DO J = JTS, MIN(JTE,JDE-1)
2494 DO K = KTS,KTE
2495 DO I = ITS, MIN(ITE,IDE-1)
2496 HTM(I,K,J)=1.0
2497 VTM(I,K,J)=1.0
2498 ENDDO
2499 ENDDO
2500 ENDDO
2501
2502 DO J = JTS, MIN(JTE,JDE-1)
2503 DO I = ITS, MIN(ITE,IDE-1)
2504 LMH(I,J)= KME-1 ! note the flipping for start_domain_nmm.F
2505 LMV(I,J)= KME-1 ! this is consistent with Tom's version
2506 RES(I,J)=1.
2507 ENDDO
2508 ENDDO
2509
2510 ! INITIALIZE 2D BOUNDARY MASKS
2511
2512 !! HBM2:
2513
2514 HBM2=0.
2515 DO J = JTS, MIN(JTE,JDE-1)
2516 DO I = ITS, MIN(ITE,IDE-1)
2517 IF((J .GE. 3 .and. J .LE. (JDE-1)-2) .AND. &
2518 (I .GE. 2 .and. I .LE. (IDE-1)-2+MOD(J,2))) THEN
2519 HBM2(I,J)=1.
2520 ENDIF
2521 ENDDO
2522 ENDDO
2523
2524 !! HBM3:
2525
2526 HBM3=0.
2527 DO J=JTS,MIN(JTE,JDE-1)
2528 IHWG(J)=mod(J+1,2)-1
2529 IF (J .ge. 4 .and. J .le. (JDE-1)-3) THEN
2530 IHL=(IDS+1)-IHWG(J)
2531 IHH=(IDE-1)-2
2532 DO I=ITS,MIN(ITE,IDE-1)
2533 IF (I .ge. IHL .and. I .le. IHH) HBM3(I,J)=1.
2534 ENDDO
2535 ENDIF
2536 ENDDO
2537
2538 !! VBM2
2539
2540 VBM2=0.
2541 DO J=JTS,MIN(JTE,JDE-1)
2542 DO I=ITS,MIN(ITE,IDE-1)
2543 IF((J .ge. 3 .and. J .le. (JDE-1)-2) .AND. &
2544 (I .ge. 2 .and. I .le. (IDE-1)-1-MOD(J,2))) THEN
2545 VBM2(I,J)=1.
2546 ENDIF
2547 ENDDO
2548 ENDDO
2549
2550 !! VBM3
2551
2552 VBM3=0.
2553 DO J=JTS,MIN(JTE,JDE-1)
2554 DO I=ITS,MIN(ITE,IDE-1)
2555 IF((J .ge. 4 .and. J .le. (JDE-1)-3) .AND. &
2556 (I .ge. 3-MOD(J,2) .and. I .le. (IDE-1)-2)) THEN
2557 VBM3(I,J)=1.
2558 ENDIF
2559 ENDDO
2560 ENDDO
2561
2562 TPH0D = grid%CEN_LAT
2563 TLM0D = grid%CEN_LON
2564 TPH0 = TPH0D*DTR
2565 WBD = grid%WBD0 ! gopal's doing: may use Registry WBD0 now
2566 WB = WBD*DTR
2567 SBD = grid%SBD0 ! gopal's doing: may use Registry SBD0 now
2568 SB = SBD*DTR
2569 DLM = DLMD*DTR ! input now from med_nest_egrid_configure
2570 DPH = DPHD*DTR ! input now from med_nest_egrid_configure
2571 TDLM = DLM+DLM
2572 TDPH = DPH+DPH
2573 WBI = WB+TDLM
2574 SBI = SB+TDPH
2575 EBI = WB+((ide-1)-2)*TDLM ! gopal's doing: check this for nested domain
2576 ANBI = SB+((jde-1)-3)*DPH ! gopal's doing: check this for nested domain
2577 STPH0 = SIN(TPH0)
2578 CTPH0 = COS(TPH0)
2579 TSPH = 3600./grid%DT
2580 DTAD = 1.0
2581 DTCF = 4.0
2582 DY_NMM0= DY_NMM ! ERAD*DPH; input now from med_nest_egrid_configure
2583
2584 ! CORIOLIS PARAMETER (There appears to be some roundoff in computing TLM & STPH and other terms,
2585 ! in the nested domain. The problem needs to be revisited
2586
2587 DO J=JTS,MIN(JTE,JDE-1)
2588 TLM0=WB-TDLM+MOD(J,2)*DLM ! remember this is a wind point
2589 TPH =SB+float(J-1)*DPH
2590 STPH=SIN(TPH)
2591 CTPH=COS(TPH)
2592 DO I=ITS,MIN(ITE,IDE-1)
2593 TLM=TLM0 + I*TDLM
2594 FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
2595 F(I,J)=0.5*grid%DT*FP
2596 ENDDO
2597 ENDDO
2598
2599
2600 DO J=JTS,MIN(JTE,JDE-1)
2601 KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
2602 KVL2(J)=(IDE-1)*(J-1)-J/2+2
2603 KHH2(J)=(IDE-1)*J-J/2-1
2604 KVH2(J)=(IDE-1)*J-(J+1)/2-1
2605 ENDDO
2606
2607
2608 TPH=SB-DPH
2609 DO J=JTS,MIN(JTE,JDE-1)
2610 TPH=SB+float(J-1)*DPH
2611 DXP=ERAD*DLM*COS(TPH)
2612 DXJ(J)=DXP
2613 WPDARJ(J)=-W_NMM*((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)/ &
2614 (grid%DT*32.*DXP*DY_NMM0)
2615 CPGFUJ(J)=-grid%DT/(48.*DXP)
2616 CURVJ(J)=.5*grid%DT*TAN(TPH)/ERAD
2617 FCPJ(J)=grid%DT/(CP*192.*DXP*DY_NMM0)
2618 FDIVJ(J)=1./(12.*DXP*DY_NMM0)
2619 FADJ(J)=-grid%DT/(48.*DXP*DY_NMM0)*DTAD
2620 ACDT=grid%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)
2621 CDDAMP=CODAMP*ACDT
2622 HDACJ(J)=COAC*ACDT/(4.*DXP*DY_NMM0)
2623 DDMPUJ(J)=CDDAMP/DXP
2624 DDMPVJ(J)=CDDAMP/DY_NMM0
2625 ENDDO
2626
2627 ! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
2628
2629 WRITE(0,*)'NEW CHANGE',F4D,EF4T,F4Q
2630
2631 DO L=KDS,KDE-1
2632 RDETA(L)=1./DETA(L)
2633 F4Q2(L)=-.25*grid%DT*DTAD/DETA(L)
2634 ENDDO
2635
2636 DO J=JTS,MIN(JTE,JDE-1)
2637 DO I=ITS,MIN(ITE,IDE-1)
2638 DX_NMM(I,J)=DXJ(J)
2639 WPDAR(I,J)=WPDARJ(J)*HBM2(I,J)
2640 CPGFU(I,J)=CPGFUJ(J)*VBM2(I,J)
2641 CURV(I,J)=CURVJ(J)*VBM2(I,J)
2642 FCP(I,J)=FCPJ(J)*HBM2(I,J)
2643 FDIV(I,J)=FDIVJ(J)*HBM2(I,J)
2644 FAD(I,J)=FADJ(J)
2645 HDACV(I,J)=HDACJ(J)*VBM2(I,J)
2646 HDAC(I,J)=HDACJ(J)*1.25*HBM2(I,J)
2647 ENDDO
2648 ENDDO
2649
2650 DO J=JTS, MIN(JTE,JDE-1)
2651 IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
2652 KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
2653 DO I=ITS,MIN(ITE,IDE-1)
2654 IF (I .ge. 2 .and. I .le. KHH) THEN
2655 HDAC(I,J)=HDAC(I,J)* DFC
2656 ENDIF
2657 ENDDO
2658 ELSE
2659 KHH=2+MOD(J,2)
2660 DO I=ITS,MIN(ITE,IDE-1)
2661 IF (I .ge. 2 .and. I .le. KHH) THEN
2662 HDAC(I,J)=HDAC(I,J)* DFC
2663 ENDIF
2664 ENDDO
2665 KHH=(IDE-1)-2+MOD(J,2)
2666
2667 DO I=ITS,MIN(ITE,IDE-1)
2668 IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
2669 HDAC(I,J)=HDAC(I,J)* DFC
2670 ENDIF
2671 ENDDO
2672 ENDIF
2673 ENDDO
2674
2675 DO J=JTS,min(JTE,JDE-1)
2676 DO I=ITS,min(ITE,IDE-1)
2677 DDMPU(I,J)=DDMPUJ(J)*VBM2(I,J)
2678 DDMPV(I,J)=DDMPVJ(J)*VBM2(I,J)
2679 HDACV(I,J)=HDACV(I,J)*VBM2(I,J)
2680 ENDDO
2681 ENDDO
2682
2683 ! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
2684
2685 DO J=JTS,MIN(JTE,JDE-1)
2686 IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
2687 KVH=(IDE-1)-1-MOD(J,2)
2688 DO I=ITS,MIN(ITE,IDE-1)
2689 IF (I .ge. 2 .and. I .le. KVH) THEN
2690 DDMPU(I,J)=DDMPU(I,J)*DDFC
2691 DDMPV(I,J)=DDMPV(I,J)*DDFC
2692 HDACV(I,J)=HDACV(I,J)*DFC
2693 ENDIF
2694 ENDDO
2695 ELSE
2696 KVH=3-MOD(J,2)
2697 DO I=ITS,MIN(ITE,IDE-1)
2698 IF (I .ge. 2 .and. I .le. KVH) THEN
2699 DDMPU(I,J)=DDMPU(I,J)*DDFC
2700 DDMPV(I,J)=DDMPV(I,J)*DDFC
2701 HDACV(I,J)=HDACV(I,J)*DFC
2702 ENDIF
2703 ENDDO
2704 KVH=(IDE-1)-1-MOD(J,2)
2705 DO I=ITS,MIN(ITE,IDE-1)
2706 IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
2707 DDMPU(I,J)=DDMPU(I,J)*DDFC
2708 DDMPV(I,J)=DDMPV(I,J)*DDFC
2709 HDACV(I,J)=HDACV(I,J)*DFC
2710 ENDIF
2711 ENDDO
2712 ENDIF
2713 ENDDO
2714
2715 ! This one was left over for nested domain
2716
2717 DO J = JTS, MIN(JTE,JDE-1)
2718 DO I = ITS, MIN(ITE,IDE-1)
2719 GLAT(I,J)=HLAT(I,J)*DTR
2720 GLON(I,J)=HLON(I,J)*DTR
2721 ENDDO
2722 ENDDO
2723
2724 !! compute EMT, EM on global domain, and only on task 0.
2725
2726 ! IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
2727
2728 ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
2729 write(0,*) 'FIGURING OUT EMJ, EMTJ ', JDS, JDE-1
2730 DO J=JDS,JDE-1
2731 TPH=SB+float(J-1)*DPH
2732 DXP=ERAD*DLM*COS(TPH)
2733 EMJ(J)= grid%DT/( 4.*DXP)*DTAD
2734 EMTJ(J)=grid%DT/(16.*DXP)*DTAD
2735 ! write(0,*) 'J, EMTJ(J): ', J, EMTJ(J)
2736 ENDDO
2737
2738 JA=0
2739 DO 161 J=3,5
2740 JA=JA+1
2741 KHLA(JA)=2
2742 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
2743 161 EMT(JA)=EMTJ(J)
2744 DO 162 J=(JDE-1)-4,(JDE-2)-2
2745 JA=JA+1
2746 KHLA(JA)=2
2747 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
2748 162 EMT(JA)=EMTJ(J)
2749 DO 163 J=6,(JDE-1)-5
2750 JA=JA+1
2751 KHLA(JA)=2
2752 KHHA(JA)=2+MOD(J,2)
2753 163 EMT(JA)=EMTJ(J)
2754 DO 164 J=6,(JDE-1)-5
2755 JA=JA+1
2756 KHLA(JA)=(IDE-1)-2
2757 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
2758 164 EMT(JA)=EMTJ(J)
2759
2760 ! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
2761
2762 JA=0
2763 DO 171 J=3,5
2764 JA=JA+1
2765 KVLA(JA)=2
2766 KVHA(JA)=(IDE-1)-1-MOD(J,2)
2767 171 EM(JA)=EMJ(J)
2768 DO 172 J=(JDE-1)-4,(JDE-2)-2
2769 JA=JA+1
2770 KVLA(JA)=2
2771 KVHA(JA)=(IDE-1)-1-MOD(J,2)
2772 172 EM(JA)=EMJ(J)
2773 DO 173 J=6,(JDE-1)-5
2774 JA=JA+1
2775 KVLA(JA)=2
2776 KVHA(JA)=2+MOD(J+1,2)
2777 173 EM(JA)=EMJ(J)
2778 DO 174 J=6,(JDE-1)-5
2779 JA=JA+1
2780 KVLA(JA)=(IDE-1)-2
2781 KVHA(JA)=(IDE-1)-1-MOD(J,2)
2782 174 EM(JA)=EMJ(J)
2783
2784 ! ENDIF ! wrf_dm_on_monitor
2785
2786 !! must be a better place to put this, but will eliminate "phantom"
2787 !! wind points here (no wind point on eastern boundary of odd numbered rows)
2788 !!
2789 ! phantom
2790 IF (ABS(IDE-1-ITE) .eq. 1 ) THEN ! |
2791 WRITE(0,*)'zero phantom winds' ! H [x] H V
2792 DO K=KDS,KDE-1 !
2793 DO J=JDS,JDE-1,2 ! V [H] V H
2794 IF (J .ge. JTS .and. J .le. JTE) THEN !
2795 U(IDE-1,K,J)=0. ! H [x] H V
2796 V(IDE-1,K,J)=0. ! ------ ------
2797 ENDIF ! ide-1 ide
2798 ENDDO ! NMM/SI WRF
2799 ENDDO ! domain domain
2800 ENDIF ! (dummy)
2801
2802
2803 ! just a test for gravity waves
2804
2805 ! PD=62000.
2806 ! U=0.0
2807 ! V=0.0
2808 ! T=300.
2809 ! Q=0.0
2810 ! Q2=0.0
2811 ! CWM=0.0
2812 ! FIS=0.0
2813
2814 ! testx
2815 ! DO J = JTS, MIN(JTE,JDE-1)
2816 ! DO K = KTS,KTE
2817 ! DO I = ITS, MIN(ITE,IDE-1)
2818 ! SM(I,J)=I
2819 ! U(I,K,J)=J
2820 ! ENDDO
2821 ! ENDDO
2822 ! ENDDO
2823 !
2824
2825 ! deallocs
2826
2827 DEALLOCATE(KHL2,KVL2,KHH2,KVH2)
2828 DEALLOCATE(DXJ,WPDARJ,CPGFUJ,CURVJ)
2829 DEALLOCATE(FCPJ,FDIVJ,FADJ)
2830 DEALLOCATE(HDACJ,DDMPUJ,DDMPVJ)
2831 DEALLOCATE(KHLA,KHHA)
2832 DEALLOCATE(KVLA,KVHA)
2833
2834
2835 END SUBROUTINE med_initialize_nest_nmm
2836 !======================================================================
2837
2838 subroutine smdhld(ime,jme,h,s,lines,nsmud)
2839 dimension ihw(jme),ihe(jme)
2840 dimension h(ime,jme),s(ime,jme) &
2841 & ,hbms(ime,jme),hne(ime,jme),hse(ime,jme)
2842 !-----------------------------------------------------------------------
2843 do j=1,jme
2844 ihw(j)=-mod(j,2)
2845 ihe(j)=ihw(j)+1
2846 enddo
2847 !-----------------------------------------------------------------------
2848
2849 do j=1,jme
2850 do i=1,ime
2851 hbms(i,j)=1.-s(i,j)
2852 enddo
2853 enddo
2854 !
2855 jmelin=jme-lines+1
2856 ibas=lines/2
2857 m2l=mod(lines,2)
2858 !
2859 do j=lines,jmelin
2860 ihl=ibas+mod(j,2)+m2l*mod(j+1,2)
2861 ihh=ime-ibas-m2l*mod(j+1,2)
2862
2863 ! write(6,*) 'no smooth limits for J: ', J, 'are ', ihl,ihh
2864 !
2865 do i=ihl,ihh
2866 hbms(i,j)=0.
2867 enddo
2868 enddo
2869 !-----------------------------------------------------------------------
2870 do ks=1,nsmud
2871
2872 write(6,*) 'H(1,1): ', h(1,1)
2873 write(6,*) 'H(3,1): ', h(1,1)
2874 !-----------------------------------------------------------------------
2875 do j=1,jme-1
2876 do i=1,ime-1
2877 hne(i,j)=h(i+ihe(j),j+1)-h(i,j)
2878 enddo
2879 enddo
2880 do j=2,jme
2881 do i=1,ime-1
2882 hse(i,j)=h(i+ihe(j),j-1)-h(i,j)
2883 enddo
2884 enddo
2885 !
2886 do j=2,jme-1
2887 do i=1+mod(j,2),ime-1
2888 h(i,j)=(hne(i,j)-hne(i+ihw(j),j-1) &
2889 & +hse(i,j)-hse(i+ihw(j),j+1))*hbms(i,j)*0.125+h(i,j)
2890 enddo
2891 enddo
2892 !-----------------------------------------------------------------------
2893
2894 !!! smooth around boundary somehow?
2895
2896 ! special treatment for four corners
2897
2898 if (hbms(1,1) .eq. 1) then
2899 h(1,1)=0.75*h(1,1)+0.125*h(1+ihe(1),2)+ &
2900 & 0.0625*(h(2,1)+h(1,3))
2901 endif
2902
2903 if (hbms(ime,1) .eq. 1) then
2904 h(ime,1)=0.75*h(ime,1)+0.125*h(ime+ihw(1),2)+ &
2905 & 0.0625*(h(ime-1,1)+h(ime,3))
2906 endif
2907
2908 if (hbms(1,jme) .eq. 1) then
2909 h(1,jme)=0.75*h(1,jme)+0.125*h(1+ihe(jme),jme-1)+ &
2910 & 0.0625*(h(2,jme)+h(1,jme-2))
2911 endif
2912
2913 if (hbms(ime,jme) .eq. 1) then
2914 h(ime,jme)=0.75*h(ime,jme)+0.125*h(ime+ihw(jme),jme-1)+ &
2915 & 0.0625*(h(ime-1,jme)+h(ime,jme-2))
2916 endif
2917
2918
2919 ! S bound
2920
2921 J=1
2922 do I=2,ime-1
2923 if (hbms(I,J) .eq. 1) then
2924 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
2925 endif
2926 enddo
2927
2928 ! N bound
2929
2930 J=JME
2931 do I=2,ime-1
2932 if (hbms(I,J) .eq. 1) then
2933 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
2934 endif
2935 enddo
2936
2937 ! W bound
2938
2939 I=1
2940 do J=3,jme-2
2941 if (hbms(I,J) .eq. 1) then
2942 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
2943 endif
2944 enddo
2945
2946 ! E bound
2947
2948 I=IME
2949 do J=3,jme-2
2950 if (hbms(I,J) .eq. 1) then
2951 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
2952 endif
2953 enddo
2954
2955
2956 enddo ! end ks loop
2957
2958 !! (light touch) with 5-point filter over untouched interior?
2959
2960 ! do ks=1,5
2961 ! do J=lines-1,jme-(lines-1)
2962 ! do I=lines-1,ime-(lines-1)
2963 ! if (s(I,J) .eq. 0 .and.
2964 ! & h(I,J) .gt. h(i+ihw(J),J+1) .and.
2965 ! & h(I,J) .gt. h(I+ihe(J),J+1) .and.
2966 ! & h(I,J) .gt. h(i+ihw(J),J-1) .and.
2967 ! & h(I,J) .gt. h(I+ihe(J),J-1)) then
2968 ! write(6,*) 'smoothing topo at I,J...', I,J,H(I,J)
2969 ! h(I,J)=h(I,J)+0.125*( h(i+ihw(J),J+1) + h(I+ihe(J),J+1) +
2970 ! & h(i+ihw(J),J-1) + h(I+ihe(J),J-1) -
2971 ! & 4*h(I,J) )
2972 ! write(6,*) 'post smoothing val', ks,H(I,J)
2973 ! endif
2974 ! enddo
2975 ! enddo
2976 ! enddo
2977
2978 !-----------------------------------------------------------------------
2979 return
2980 end subroutine smdhld
2981
2982 !--------------------------------------------------------------------------------------
2983 #if 0
2984 SUBROUTINE initial_nest_pivot ( parent , nest, iloc, jloc )
2985
2986 !==========================================================================================
2987 !
2988 ! This program produces i_start and j_start for the nested domain depending on the
2989 ! central lat-lon of the storm.
2990 !
2991 !==========================================================================================
2992
2993 USE module_domain
2994 USE module_configure
2995 USE module_timing
2996 USE module_dm
2997
2998 IMPLICIT NONE
2999 TYPE(domain) , POINTER :: parent , nest
3000 INTEGER, INTENT(OUT) :: ILOC,JLOC
3001 INTEGER :: IMS,IME,JMS,JME,KMS,KME
3002 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3003 INTEGER :: IMS,IME,JMS,JME,KMS,KME
3004 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
3005 INTEGER :: NIDE,NJDE ! nest dimension
3006 INTEGER :: I,J,ITER,IDUM,JDUM
3007 REAL :: ALAT,ALON,DIFF1,DIFF2,ERR
3008 REAL :: parent_CLAT,parent_CLON,parent_SLAT,parent_SLON
3009 REAL :: parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
3010 !========================================================================================
3011
3012 ! First obtain central latitude and longitude for the parent domain
3013
3014 CALL nl_get_cen_lat (parent%ID, parent_CLAT)
3015 CALL nl_get_cen_lon (parent%ID, parent_CLON)
3016 ! CALL nl_get_storm_lat (parent%ID, parent_SLAT)
3017 ! CALL nl_get_storm_lon (parent%ID, parent_SLON)
3018
3019 ! Parent grid configuration, including, western and southern boundary
3020
3021 IDS = parent%sd31
3022 IDE = parent%ed31
3023 KDS = parent%sd32
3024 KDE = parent%ed32
3025 JDS = parent%sd33
3026 JDE = parent%ed33
3027
3028 IMS = parent%sm31
3029 IME = parent%em31
3030 KMS = parent%sm32
3031 KME = parent%em32
3032 JMS = parent%sm33
3033 JME = parent%em33
3034
3035 ITS = parent%sp31
3036 ITE = parent%ep31
3037 KTS = parent%sp32
3038 KTE = parent%ep32
3039 JTS = parent%sp33
3040 JTE = parent%ep33
3041
3042 NIDE = nest%ed31
3043 NJDE = nest%ed33
3044
3045 parent_DLMD = parent%dx ! DLMD: dlamda in degrees
3046 parent_DPHD = parent%dy ! DPHD: dphi in degrees
3047 parent_WBD = -(IDE-2)*parent%dx ! WBD0: in deg;factor 2 takes care of dummy last column
3048 parent_SBD = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd
3049 ALAT = parent_SLAT - 0.5*(NJDE-2)*parent_DPHD/nest%parent_grid_ratio
3050 ALON = parent_SLON - 1.0*(NIDE-2)*parent_DLMD/nest%parent_grid_ratio
3051
3052 ! WRITE(0,*)'ALAT AND ALON=',ALAT,ALON
3053
3054 CALL EARTH_LATLON ( parent%nmm_HLAT,parent%nmm_HLON,parent%nmm_VLAT,parent%nmm_VLON, & !output
3055 parent_DLMD,parent_DPHD,parent_WBD,parent_SBD, & !inputs
3056 parent_CLAT,parent_CLON, &
3057 IDS,IDE,JDS,JDE,KDS,KDE, &
3058 IMS,IME,JMS,JME,KMS,KME, &
3059 ITS,ITE,JTS,JTE,KTS,KTE )
3060
3061 ! start iteration
3062
3063 ILOC=-99
3064 JLOC=-99
3065 ERR=0.1
3066 ITER=1
3067 100 CONTINUE
3068
3069 DO J = JTS,min(JTE,JDE-1)
3070 DO I = ITS,min(ITE,IDE-1)
3071 DIFF1 = ABS(ALAT - parent%nmm_HLAT(I,J))
3072 DIFF2 = ABS(ALON - parent%nmm_HLON(I,J))
3073 IF(DIFF1 .LE. ERR .AND. DIFF2 .LE. ERR)THEN
3074 ILOC=I
3075 JLOC=J
3076 ! WRITE(0,*)'ITERATED',ERR,ITER,I,J,parent%nmm_HLAT(I,J),ALAT,parent%nmm_HLON(I,J),ALON
3077 ENDIF
3078 ENDDO
3079 ENDDO
3080
3081 CALL wrf_dm_maxval_integer ( ILOC, idum, jdum )
3082 CALL wrf_dm_maxval_integer ( JLOC, idum, jdum )
3083
3084 IF(ILOC .EQ. -99 .AND. JLOC .EQ. -99)THEN
3085 ERR=ERR+0.1
3086 ITER=ITER+1
3087 IF(ITER .LE. 100)GO TO 100
3088 ENDIF
3089
3090 IF(ILOC .NE. -99 .AND. JLOC .NE. -99)THEN
3091 WRITE(0,*)'NOTE: I_PARENT_START AND J_PARENT_START FOUND FOR THE NESTED DOMAIN CONFIGURATION AT ITER=',ITER
3092 WRITE(0,*)'istart=',ILOC
3093 WRITE(0,*)'jstart=',JLOC
3094 ELSE
3095 ILOC=IDE/3
3096 JLOC=JDE/3
3097 !
3098 WRITE(0,*)'WARNING: COULD NOT LOCATE I_PARENT_START AND J_PARENT_START FROM INPUT STORM INFO'
3099 WRITE(0,*)'ISTART=',IDE/3
3100 WRITE(0,*)'JSTART=',JDE/3
3101 ENDIF
3102
3103 RETURN
3104 END SUBROUTINE initial_nest_pivot
3105
3106 !============================================================================================
3107 #endif
3108
3109 LOGICAL FUNCTION nmm_cd_feedback_mask( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
3110 INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
3111 LOGICAL, INTENT(IN) :: xstag, ystag
3112
3113 INTEGER ioff, joff
3114
3115 ioff = 0 ; joff = 0
3116 IF ( xstag ) ioff = 1
3117 IF ( ystag ) joff = 1
3118
3119 nmm_cd_feedback_mask = ( pig .ge. ips_save+1 .and. &
3120 pjg .ge. jps_save+1 .and. &
3121 pig .le. ipe_save-1 +ioff .and. &
3122 pjg .le. jpe_save-1 +joff )
3123
3124 END FUNCTION nmm_cd_feedback_mask
3125
3126 !----------------------------------------------------------------------------
3127 #else
3128 SUBROUTINE stub_nmm_nest_stub
3129 END SUBROUTINE stub_nmm_nest_stub
3130 #endif
3131
3132 RECURSIVE SUBROUTINE find_ijstart_level ( grid, i_start, j_start, level )
3133
3134 ! Dusan Jovic
3135
3136 USE module_domain
3137
3138 IMPLICIT NONE
3139
3140 ! Input data.
3141
3142 TYPE(domain) :: grid
3143 INTEGER, INTENT (OUT) :: i_start, j_start, level
3144 INTEGER :: iadd
3145
3146 if (grid%parent_id == 0 ) then
3147 i_start = 1
3148 j_start = 1
3149 level = 0
3150 else
3151 call find_ijstart_level ( grid%parents(1)%ptr, i_start, j_start, level )
3152 if (level > 0) then
3153 iadd = (i_start-1)*3
3154 if ( mod(j_start,2).ne.0 .and. mod(grid%j_parent_start,2).ne.0 ) iadd = iadd - 1
3155 if ( mod(j_start,2).eq.0 .and. mod(grid%j_parent_start,2).eq.0 ) iadd = iadd + 2
3156 else
3157 iadd = -mod(grid%j_parent_start,2)
3158 end if
3159 i_start = iadd + grid%i_parent_start*3 - 1
3160 j_start = ( (j_start-1) + (grid%j_parent_start-1) ) * 3 + 1
3161 level = level + 1
3162 end if
3163
3164 END SUBROUTINE find_ijstart_level