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