Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

NMM_NEST_UTILS1.F @ 3547

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

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

Download in other formats:

Original Format