Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

NMM_NEST_UTILS1.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 104.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format