source: lmdz_wrf/WRFV3/dyn_nmm/NMM_NEST_UTILS1.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 143.1 KB
Line 
1#if (NMM_NEST == 1)
2!===========================================================================
3!
4! E-GRID NESTING UTILITIES: This is gopal's doing
5!
6!===========================================================================
7
8SUBROUTINE 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 CHARACTER(LEN=255)                 :: message
23
24!----------------------------------------------------------------------------
25!  PURPOSE:
26!         - Initialize nested domain configurations including setting up
27!           wbd,sbd and some other variables and 1D arrays.
28!         - Note that in order to obtain coincident grid points, which 
29!           is a basic requirement for RSL, WRF infrastructure, we use
30!           western and southern boundaries of nested domain (nest%wbd0
31!           and nest%sbd0 derived from the parent domain. In this case
32!           the nested domain may be considered as a part of the parent
33!           domain with a higher resolution (telescoping ?).
34!         - Also note that in this case, the central lat/lons for nested
35!           domain should coincide with the central lat/lons of the parent,
36!           although the nested domain NEED NOT be located at the center of
37!           the domain.
38!----------------------------------------------------------------------------
39!
40!   BASIC TEST FOR PARENT DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
41!   IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT
42
43    IF(MOD(parent%ed32,2) .NE. 0)THEN
44     CALL wrf_error_fatal("PARENT DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
45    ENDIF
46
47!   BASIC TEST FOR NESTED DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
48!   IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT
49
50    IF(MOD(nest%ed32,2) .NE. 0)THEN
51     CALL wrf_error_fatal("NESTED DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
52    ENDIF
53
54!   Parent grid configuration, including, western and southern boundary
55
56    IDS = parent%sd31
57    IDE = parent%ed31
58    JDS = parent%sd32
59    JDE = parent%ed32
60    KDS = parent%sd33
61    KDE = parent%ed33
62
63    IMS = parent%sm31
64    IME = parent%em31
65    JMS = parent%sm32
66    JME = parent%em32
67    KMS = parent%sm33
68    KME = parent%em33
69
70    ITS = parent%sp31
71    ITE = parent%ep31
72    JTS = parent%sp32
73    JTE = parent%ep32
74    KTS = parent%sp33
75    KTE = parent%ep33
76
77!   grid configuration
78
79    ! calculate wbd0 and sbd0 only for MOAD i.e. grid with parent_id == 0
80    if (parent%parent_id == 0 ) then        ! Dusan's doing
81       parent%wbd0 = -(IDE-2)*parent%dx     ! WBD0: in degrees;factor 2 takes care of dummy last column
82       parent%sbd0 = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd 
83    end if
84    nest%wbd0   = parent%wbd0 + (nest%i_parent_start-1)*2.*parent%dx + mod(nest%j_parent_start+1,2)*parent%dx
85    nest%sbd0   = parent%sbd0 + (nest%j_parent_start-1)*parent%dy
86    nest%dx     = parent%dx/nest%parent_grid_ratio
87    nest%dy     = parent%dy/nest%parent_grid_ratio
88
89    write(message,*)" - i_parent_start = ",nest%i_parent_start
90    CALL wrf_message(trim(message))
91    write(message,*)" - j_parent_start = ",nest%j_parent_start
92    CALL wrf_message(trim(message))
93    write(message,*)" - parent%wbd0    = ",parent%wbd0
94    CALL wrf_message(trim(message))
95    write(message,*)" - parent%sbd0    = ",parent%sbd0
96    CALL wrf_message(trim(message))
97    write(message,*)" - nest%wbd0      = ",nest%wbd0
98    CALL wrf_message(trim(message))
99    write(message,*)" - nest%sbd0      = ",nest%sbd0
100    CALL wrf_message(trim(message))
101    write(message,*)" - nest%dx        = ",nest%dx
102    CALL wrf_message(trim(message))
103    write(message,*)" - nest%dy        = ",nest%dy
104    CALL wrf_message(trim(message))
105!
106    CALL nl_set_dx (nest%id , nest%dx)   ! for output purpose
107    CALL nl_set_dy (nest%id , nest%dy)   ! for output purpose
108
109!   set lat-lons; parent set to nested domain
110
111    CALL nl_get_cen_lat (parent%id, parent%cen_lat) ! cen_lat of parent set to nested domain
112    CALL nl_get_cen_lon (parent%id, parent%cen_lon) ! cen_lon of parent set to nested domain
113
114    nest%cen_lat=parent%cen_lat
115    nest%cen_lon=parent%cen_lon
116!
117    CALL nl_set_cen_lat ( nest%id , nest%cen_lat)  ! for output purpose
118    CALL nl_set_cen_lon ( nest%id , nest%cen_lon)  ! for output purpose
119
120    write(message,*)" - nest%cen_lat   = ",nest%cen_lat
121    CALL wrf_message(trim(message))
122    write(message,*)" - nest%cen_lon   = ",nest%cen_lon
123    CALL wrf_message(trim(message))
124
125
126!   soil configuration
127
128#ifdef HWRF
129!zhang
130    if ( .not.nest%analysis ) then
131#endif
132    nest%sldpth  = parent%sldpth
133    nest%dzsoil  = parent%dzsoil
134    nest%rtdpth  = parent%rtdpth
135#ifdef HWRF
136    endif
137#endif
138
139!   numerical set up
140
141    nest%deta        = parent%deta
142    nest%aeta        = parent%aeta
143    nest%etax        = parent%etax
144    nest%dfl         = parent%dfl
145    nest%deta1       = parent%deta1
146    nest%aeta1       = parent%aeta1
147    nest%eta1        = parent%eta1
148    nest%deta2       = parent%deta2
149    nest%aeta2       = parent%aeta2
150    nest%eta2        = parent%eta2
151    nest%pdtop       = parent%pdtop
152    nest%pt          = parent%pt
153    nest%dfrlg       = parent%dfrlg
154    nest%num_soil_layers = parent%num_soil_layers
155    nest%num_moves       = parent%num_moves
156
157! Unfortunately, some of the single value constants in used in module_initialize have
158! to be defiend here instead of the usual spot in med_initialize_nest_nmm. There
159! appears to be a problem in Registry and related code in this area.
160!
161! state  logical upstrm   -      dyn_nmm     -      -      -
162
163
164    nest%dlmd   = nest%dx
165    nest%dphd   = nest%dy
166    nest%dy_nmm = erad*(nest%dphd*dtr)
167    nest%cpgfv  = -nest%dt/(48.*nest%dy_nmm)
168    nest%en     = nest%dt/( 4.*nest%dy_nmm)*dtad
169    nest%ent    = nest%dt/(16.*nest%dy_nmm)*dtad
170    nest%f4d    = -.5*nest%dt*dtad
171    nest%f4q    = -nest%dt*dtad
172    nest%ef4t   = .5*nest%dt/cp
173
174!  Other output configurations that will make grads happy
175
176   CALL nl_get_truelat1 (parent%id, parent%truelat1 )
177   CALL nl_get_truelat2 (parent%id, parent%truelat2 )
178#ifdef HWRF
179! bao : to make the restart output identical at the restart initial time for stand_lon
180   CALL nl_get_stand_lon (parent%id, parent%stand_lon )
181#endif
182   CALL nl_get_map_proj (parent%id, parent%map_proj )
183   CALL nl_get_iswater (parent%id, parent%iswater )
184
185   nest%truelat1=parent%truelat1
186   nest%truelat2=parent%truelat2
187!bao
188   nest%stand_lon=parent%stand_lon
189!bao
190   nest%map_proj=parent%map_proj
191   nest%iswater=parent%iswater
192
193   CALL nl_set_truelat1(nest%id, nest%truelat1)
194   CALL nl_set_truelat2(nest%id, nest%truelat2)
195!bao
196   CALL nl_set_stand_lon(nest%id, nest%stand_lon)
197!bao
198   CALL nl_set_map_proj(nest%id, nest%map_proj)
199   CALL nl_set_iswater(nest%id, nest%iswater)
200
201!   physics and other configurations
202!   CALL nl_get_iswater (parent%id, nest%iswater) ! iswater is just based on parents
203!   CALL nl_get_bl_surface_physics (nest%id,  nest%bl_surface_physics )
204!   CALL nl_get_num_soil_layers( parent%num_soil_layers )
205!   CALL nl_get_real_data_init_type (parent%real_data_init_type)
206
207END SUBROUTINE med_nest_egrid_configure
208
209SUBROUTINE med_construct_egrid_weights ( parent , nest )
210 USE module_domain
211 USE module_configure
212 USE module_timing
213
214 IMPLICIT NONE
215 TYPE(domain) , POINTER             :: parent , nest
216 LOGICAL, EXTERNAL                  :: wrf_dm_on_monitor
217 INTEGER                            :: IDS,IDE,JDS,JDE,KDS,KDE
218 INTEGER                            :: IMS,IME,JMS,JME,KMS,KME
219 INTEGER                            :: ITS,ITE,JTS,JTE,KTS,KTE
220 INTEGER                            :: I,J,II,JJ,NII,NJJ
221 REAL                               :: parent_CLAT,parent_CLON,parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
222 REAL                               :: nest_WBD,nest_SBD,nest_DLMD,nest_DPHD
223 REAL                               :: SW_LATD,SW_LOND
224 REAL                               :: ADDSUM1,ADDSUM2
225 REAL                               :: xr,zr,xc
226!-----------------------------------------------------------------------------------------------------------
227!   PURPOSE:
228!           - Initialize lat-lons and determine weights
229!
230!----------------------------------------------------------------------------------------------------------
231
232!   First obtain central latitude and longitude for the parent domain
233
234    CALL nl_get_cen_lat (parent%ID, parent_CLAT)
235    CALL nl_get_cen_lon (parent%ID, parent_CLON)
236
237!   Parent grid configuration, including, western and southern boundary
238
239    IDS = parent%sd31
240    IDE = parent%ed31
241    JDS = parent%sd32
242    JDE = parent%ed32
243    KDS = parent%sd33
244    KDE = parent%ed33
245
246    IMS = parent%sm31
247    IME = parent%em31
248    JMS = parent%sm32
249    JME = parent%em32
250    KMS = parent%sm33
251    KME = parent%em33
252
253    ITS  = parent%sp31
254    ITE  = parent%ep31
255    JTS  = parent%sp32
256    JTE  = parent%ep32
257    KTS  = parent%sp33
258    KTE  = parent%ep33
259!
260    parent_DLMD = parent%dx                ! DLMD: dlamda in degrees
261    parent_DPHD = parent%dy                ! DPHD: dphi in degrees
262    parent_WBD  = parent%wbd0
263    parent_SBD  = parent%sbd0
264
265!   Now compute Geodetic lat/lon (Positive East) of parent grid in degrees
266
267    CALL EARTH_LATLON ( parent%HLAT,parent%HLON,parent%VLAT,parent%VLON,  & !output
268                        parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,                    & !inputs
269                        parent_CLAT,parent_CLON,                                          &
270                        IDS,IDE,JDS,JDE,KDS,KDE,                                          &
271                        IMS,IME,JMS,JME,KMS,KME,                                          &
272                        ITS,ITE,JTS,JTE,KTS,KTE                                           )
273
274!   Nested grid configuration, including, western and southern boundary
275
276    IDS = nest%sd31
277    IDE = nest%ed31
278    JDS = nest%sd32
279    JDE = nest%ed32
280    KDS = nest%sd33
281    KDE = nest%ed33
282
283    IMS = nest%sm31
284    IME = nest%em31
285    JMS = nest%sm32
286    JME = nest%em32
287    KMS = nest%sm33
288    KME = nest%em33
289
290    ITS  = nest%sp31
291    ITE  = nest%ep31
292    JTS  = nest%sp32
293    JTE  = nest%ep32
294    KTS  = nest%sp33
295    KTE  = nest%ep33
296!
297    nest_DLMD = nest%dx
298    nest_DPHD = nest%dy
299    nest_WBD  = nest%wbd0
300    nest_SBD  = nest%sbd0
301
302!
303!   Now compute Geodetic lat/lon (Positive East) of nest in degrees, with the same central lat-lon
304!   as the parent grid
305!
306
307    CALL EARTH_LATLON ( nest%HLAT,nest%HLON,nest%VLAT,nest%VLON, & ! output
308                        nest_DLMD,nest_DPHD,nest_WBD,nest_SBD,                   & ! nest inputs
309                        parent_CLAT,parent_CLON,                                 & ! parent central lat/lon
310                        IDS,IDE,JDS,JDE,KDS,KDE,                                 & ! nested domain dimension
311                        IMS,IME,JMS,JME,KMS,KME,                                 &
312                        ITS,ITE,JTS,JTE,KTS,KTE                                  )
313
314!   Determine the weights of nested grid h points nearest to H points of parent domain
315
316#ifdef HWRFX
317    CALL G2T2H_new(    nest%IIH,nest%JJH,                            & ! output grid index in parent grid
318                       nest%HBWGT1,nest%HBWGT2,                      & ! output weights in terms of parent grid
319                       nest%HBWGT3,nest%HBWGT4,                      &
320                       nest%I_PARENT_START,nest%J_PARENT_START,      & ! nest start I, J in parent domain
321                       3,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
322                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
323                       IMS,IME,JMS,JME,KMS,KME,            &
324                       ITS,ITE,JTS,JTE,KTS,KTE      )
325#else
326    CALL G2T2H( nest%IIH,nest%JJH,                       & ! output grid index on nested grid
327                nest%HBWGT1,nest%HBWGT2,                 & ! output weights on the nested grid
328                nest%HBWGT3,nest%HBWGT4,                 &
329                nest%HLAT,nest%HLON,                     & ! target (nest) input lat lon in degrees
330                parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,   & ! parent res, western and south boundaries
331                parent_CLAT,parent_CLON,                         & ! parent central lat,lon, all in degrees
332                parent%ed31,parent%ed32,                         & ! parent imax and jmax
333                IDS,IDE,JDS,JDE,KDS,KDE,                         & !
334                IMS,IME,JMS,JME,KMS,KME,                         & ! nested grid configuration
335                ITS,ITE,JTS,JTE,KTS,KTE                          ) !
336#endif
337
338
339!   Determine the weights of nested grid v points nearest to V points of parent domain
340
341#ifdef HWRFX
342    CALL G2T2V_new(    nest%IIV,nest%JJV,                            & ! output grid index in parent grid
343                       nest%VBWGT1,nest%VBWGT2,                      & ! output weights in terms of parent grid
344                       nest%VBWGT3,nest%VBWGT4,                      &
345                       nest%I_PARENT_START,nest%J_PARENT_START,      & ! nest start I, J in parent domain
346                       3,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
347                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
348                       IMS,IME,JMS,JME,KMS,KME,            &
349                       ITS,ITE,JTS,JTE,KTS,KTE      )
350#else
351    CALL G2T2V( nest%IIV,nest%JJV,                       & ! output grid index on nested grid
352                nest%VBWGT1,nest%VBWGT2,                 & ! output weights on the nested grid
353                nest%VBWGT3,nest%VBWGT4,                 &
354                nest%VLAT,nest%VLON,                     & ! target (nest) input lat lon in degrees
355                parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,   & ! parent res, western and south boundaries
356                parent_CLAT,parent_CLON,                         & ! parent central lat,lon, all in degrees
357                parent%ed31,parent%ed32,                         & ! parent imax and jmax
358                IDS,IDE,JDS,JDE,KDS,KDE,                         & !
359                IMS,IME,JMS,JME,KMS,KME,                         & ! nested grid configuration
360                ITS,ITE,JTS,JTE,KTS,KTE                          ) !
361#endif
362
363!*** CHECK WEIGHTS AT MASS AND VELOCITY POINTS
364
365     CALL WEIGTS_CHECK(nest%HBWGT1,nest%HBWGT2,          & ! output weights on the nested grid
366                       nest%HBWGT3,nest%HBWGT4,          &
367                       nest%VBWGT1,nest%VBWGT2,          & ! output weights on the nested grid
368                       nest%VBWGT3,nest%VBWGT4,          &
369                       IDS,IDE,JDS,JDE,KDS,KDE,                  & !
370                       IMS,IME,JMS,JME,KMS,KME,                  & ! nested grid configuration
371                       ITS,ITE,JTS,JTE,KTS,KTE                   )
372
373!*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
374!
375    CALL BOUNDS_CHECK( nest%IIH,nest%JJH,nest%IIV,nest%JJV,   &
376                       nest%i_parent_start,nest%j_parent_start,nest%shw,      &
377                       IDS,IDE,JDS,JDE,KDS,KDE,                               & !
378                       IMS,IME,JMS,JME,KMS,KME,                               & ! nested grid configuration
379                       ITS,ITE,JTS,JTE,KTS,KTE                                )
380
381!------------------------------------------------------------------------------------------
382
383END SUBROUTINE med_construct_egrid_weights
384
385!======================================================================================
386!
387! compute earth lat-lons for parent and the nest before interpolations
388!------------------------------------------------------------------------------
389
390SUBROUTINE EARTH_LATLON ( HLAT,HLON,VLAT,VLON,     & !Earth lat,lon at H and V points
391                          DLMD1,DPHD1,WBD1,SBD1,   & !input res,west & south boundaries,
392                          CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees   
393                          IDS,IDE,JDS,JDE,KDS,KDE, & 
394                          IMS,IME,JMS,JME,KMS,KME, &
395                          ITS,ITE,JTS,JTE,KTS,KTE  )
396!============================================================================
397!
398 IMPLICIT NONE
399 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
400 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
401 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE 
402 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
403 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
404 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT)        :: HLAT,HLON,VLAT,VLON
405
406! local
407
408 INTEGER,PARAMETER                           :: KNUM=SELECTED_REAL_KIND(13)
409 INTEGER                                     :: I,J
410 REAL(KIND=KNUM)                             :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
411 REAL(KIND=KNUM)                             :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR
412 REAL(KIND=KNUM)                             :: STPH,CTPH,STPV,CTPV,PI_2
413 REAL(KIND=KNUM)                             :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV
414 REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV
415!-------------------------------------------------------------------------
416
417!
418      PI_2 = ACOS(0.)
419      DTR  = PI_2/90.
420      WB   = WBD1 * DTR                 ! WB:   western boundary in radians
421      SB   = SBD1 * DTR                 ! SB:   southern boundary in radians
422      DLM  = DLMD1 * DTR                ! DLM:  dlamda in radians
423      DPH  = DPHD1 * DTR                ! DPH:  dphi   in radians
424      TDLM = DLM + DLM                  ! TDLM: 2.0*dlamda
425      TDPH = DPH + DPH                  ! TDPH: 2.0*DPH
426
427!     For earth lat lon only
428
429      TPH0  = CENTRAL_LAT*DTR                ! TPH0: central lat in radians
430      STPH0 = SIN(TPH0)
431      CTPH0 = COS(TPH0)
432
433                                                !    .H
434      DO J = JTS,MIN(JTE,JDE-1)                 ! H./    This loop takes care of zig-zag
435!                                               !   \.H  starting points along j 
436         TLMH0 = WB - TDLM + MOD(J+1,2) * DLM   !  ./    TLMH (rotated lats at H points)
437         TLMV0 = WB - TDLM + MOD(J,2) * DLM     !  H     (//ly for V points)
438         TPHH = SB + (J-1)*DPH                  !   TPHH (rotated lons at H points) are simple trans.
439         TPHV = TPHH                            !   TPHV (rotated lons at V points) are simple trans.
440         STPH = SIN(TPHH)
441         CTPH = COS(TPHH)
442         STPV = SIN(TPHV)
443         CTPV = COS(TPHV)
444
445                                                              !   .H
446         DO I = ITS,MIN(ITE,IDE-1)                            !  /
447           TLMH = TLMH0 + I*TDLM                              !  \.H   .U   .H
448!                                                             !H./ ----><----
449           SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH)     !     DLM + DLM
450           GLATH(I,J)=ASIN(SPHH)                              ! GLATH: Earth Lat in radians
451           CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) &
452                - TAN(GLATH(I,J))*TAN(TPH0)
453           IF(CLMH .GT. 1.) CLMH = 1.0
454           IF(CLMH .LT. -1.) CLMH = -1.0
455           FACTH = 1.
456           IF(TLMH .GT. 0.) FACTH = -1.
457           GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH)
458
459         ENDDO                                   
460
461         DO I = ITS,MIN(ITE,IDE-1)
462           TLMV = TLMV0 + I*TDLM
463           SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV)
464           GLATV(I,J) = ASIN(SPHV)
465           CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) &
466                - TAN(GLATV(I,J))*TAN(TPH0)
467           IF(CLMV .GT. 1.) CLMV = 1.
468           IF(CLMV .LT. -1.) CLMV = -1.
469           FACTV = 1.
470           IF(TLMV .GT. 0.) FACTV = -1.
471           GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV)
472
473         ENDDO
474
475      ENDDO
476
477!     Conversion to degrees (may not be required, eventually)
478
479      DO J = JTS, MIN(JTE,JDE-1)
480       DO I = ITS, MIN(ITE,IDE-1)
481          HLAT(I,J) = GLATH(I,J) / DTR
482          HLON(I,J)= -GLONH(I,J)/DTR
483          IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J)  - 360.
484          IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360.
485!
486          VLAT(I,J) = GLATV(I,J) / DTR
487          VLON(I,J) = -GLONV(I,J) / DTR
488          IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J)  - 360.
489          IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360.
490
491       ENDDO
492      ENDDO
493
494END SUBROUTINE EARTH_LATLON
495
496!----------------------------------------------------------------------------- 
497
498  SUBROUTINE G2T2H( IIH,JJH,                     & ! output grid index and weights
499                    HBWGT1,HBWGT2,               & ! output weights in terms of parent grid
500                    HBWGT3,HBWGT4,               &
501                    HLAT,HLON,                   & ! target (nest) input lat lon in degrees
502                    DLMD1,DPHD1,WBD1,SBD1,       & ! parent res, west and south boundaries
503                    CENTRAL_LAT,CENTRAL_LON,     & ! parent central lat,lon, all in degrees
504                    P_IDE,P_JDE,                 & ! parent imax and jmax
505                    IDS,IDE,JDS,JDE,KDS,KDE,     & ! target (nest) dIMEnsions
506                    IMS,IME,JMS,JME,KMS,KME,     &
507                    ITS,ITE,JTS,JTE,KTS,KTE      )
508
509!
510!***  Tom Black   - Initial Version
511!***  Gopal       - Revised Version for WRF (includes coincident grid points)
512!***
513!***  GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
514!***  AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE
515!***  INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE
516!***  h POINTS OF THE NESTED DOMAIN   
517!
518!============================================================================
519!
520 IMPLICIT NONE
521 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
522 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
523 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
524 INTEGER,    INTENT(IN   )                            :: P_IDE,P_JDE
525 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
526 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
527 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(IN)   :: HLAT,HLON
528 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
529 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
530
531! local
532
533 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
534 INTEGER           :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
535 INTEGER           :: NROW,NCOL,KROWS
536 REAL(KIND=KNUM)   :: X,Y,Z,TLAT,TLON
537 REAL(KIND=KNUM)   :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB               
538 REAL(KIND=KNUM)   :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE
539 REAL(KIND=KNUM)   :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
540 REAL(KIND=KNUM)   :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3                    ! Q
541 REAL(KIND=KNUM)   :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
542 REAL(KIND=KNUM)   :: DTEMP
543 REAL   , DIMENSION(IMS:IME,JMS:JME)    :: TLATHX,TLONHX
544 INTEGER, DIMENSION(IMS:IME,JMS:JME)    :: KOUTB
545!-------------------------------------------------------------------------------
546
547  IMT=2*P_IDE-2             ! parent i dIMEnsions
548  JMT=P_JDE/2               ! parent j dIMEnsions
549  PI_2=ACOS(0.)
550  D2R=PI_2/90.
551  R2D=1./D2R
552  DPH=DPHD1*D2R
553  DLM=DLMD1*D2R
554  TPH0= CENTRAL_LAT*D2R
555  TLM0=-CENTRAL_LON*D2R        ! NOTE THE MINUS HERE
556  WB=WBD1*D2R                   ! CONVERT NESTED GRID H POINTS FROM GEODETIC
557  SB=SBD1*D2R
558  SLP0=DPHD1/DLMD1
559  DSLP0=NINT(R2D*ATAN(SLP0))
560  DS1=SQRT(DPH*DPH+DLM*DLM)    ! Q
561  AN1=ASIN(DLM/DS1)
562  AN2=ASIN(DPH/DS1)
563
564  DO J =  JTS,MIN(JTE,JDE-1)
565    DO I = ITS,MIN(ITE,IDE-1)
566
567!***
568!***  LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND
569!***  DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
570!***  CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED
571!***  COORDINATE ON THE PARENT GRID
572!
573
574      GLAT=HLAT(I,J)*D2R
575      GLON= (360. - HLON(I,J))*D2R
576      X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
577      Y=-COS(GLAT)*SIN(GLON-TLM0)
578      Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
579      TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
580      TLON=R2D*ATAN(Y/X)
581
582!      ROW=TLAT/DPHD1+JMT         ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
583!      COL=TLON/DLMD1+P_IDE-1     ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
584
585      ROW=(TLAT-SBD1)/DPHD1+1     ! Dusan's doing
586      COL=(TLON-WBD1)/DLMD1+1     ! Dusan's doing
587
588      NROW=INT(ROW + 0.001)     ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
589      NCOL=INT(COL + 0.001)     
590      TLAT=TLAT*D2R
591      TLON=TLON*D2R
592
593!*** 
594!***
595!***  FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT
596!***
597!***              V      H
598!***
599!***
600!***                 H           
601!***              H      V
602!***
603!***  THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
604!***
605      IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR.     &     
606         MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN       
607           TLAT1=(NROW-JMT)*DPH                           
608           TLAT2=TLAT1+DPH                             
609           TLON1=(NCOL-(P_IDE-1))*DLM                                   
610           TLON2=TLON1+DLM                                 
611           DLM1=TLON-TLON1                                 
612           DLM2=TLON-TLON2
613!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
614!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
615           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
616           D1=ACOS(DTEMP)
617           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
618           D2=ACOS(DTEMP)
619            IF(D1.GT.D2)THEN
620             NROW=NROW+1                    ! FIND THE NEAREST H ROW
621             NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
622            ENDIF
623      ELSE
624!***
625!***  NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT
626!***
627!***              H      V
628!***
629!***
630!***                 H
631!***              V      H
632!***
633!***  THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
634!***
635!***
636           TLAT1=(NROW+1-JMT)*DPH
637           TLAT2=TLAT1-DPH
638           TLON1=(NCOL-(P_IDE-1))*DLM
639           TLON2=TLON1+DLM
640           DLM1=TLON-TLON1
641           DLM2=TLON-TLON2
642!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
643!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
644           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
645           D1=ACOS(DTEMP)
646           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
647           D2=ACOS(DTEMP)
648             IF(D1.LT.D2)THEN
649              NROW=NROW+1                    ! FIND THE NEAREST H ROW
650             ELSE
651              NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
652             ENDIF
653      ENDIF
654
655      KROWS=((NROW-1)/2)*IMT
656      IF(MOD(NROW,2).EQ.1)THEN
657        K=KROWS+(NCOL+1)/2
658      ELSE
659        K=KROWS+P_IDE-1+NCOL/2
660      ENDIF
661
662!***
663!***  WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
664!***  NEAREST TO THE CENTER K AS SEEN BELOW.  WE MUST FIND
665!***  WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS
666!***  A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION.
667!***
668!**
669!***                  H
670!***
671!***
672!***
673!***            H     V     H
674!***
675!***
676!***               H
677!***      H     V     H     V     H
678!***
679!***
680!***
681!***            H     V     H
682!***
683!***
684!***
685!***                  H
686!***
687!***
688!***  FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H.
689!***
690    N2R=K/IMT
691    MK=MOD(K,IMT)
692!
693    IF(MK.EQ.0)THEN
694      TLATHC=SB+(2*N2R-1)*DPH
695    ELSE
696      TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH
697    ENDIF
698!
699    IF(MK.LE.(P_IDE-1))THEN
700      TLONHC=WB+2*(MK-1)*DLM
701    ELSE
702      TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM
703    ENDIF
704 
705!
706!***  EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
707!***  DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE
708!***  CAREFUL HERE     
709!
710
711    IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON
712    IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT
713    DENOM=(TLON-TLONHC)
714!
715!***
716!***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
717!***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
718!***
719!*** COINCIDENT CONDITIONS
720
721    IF(DENOM.EQ.0.0)THEN
722
723      IF(TLATHC.EQ.TLAT)THEN
724        KOUTB(I,J)=K
725        IIH(I,J) = NCOL
726        JJH(I,J) = NROW
727        TLATHX(I,J)=TLATHC
728        TLONHX(I,J)=TLONHC
729        HBWGT1(I,J)=1.0
730        HBWGT2(I,J)=0.0
731        HBWGT3(I,J)=0.0
732        HBWGT4(I,J)=0.0
733      ELSE                                      ! SAME LONGITUDE BUT DIFFERENT LATS
734!
735         IF(TLATHC .GT. TLAT)THEN      ! NESTED POINT SOUTH OF PARENT
736          KOUTB(I,J)=K-(P_IDE-1)
737          IIH(I,J) = NCOL-1
738          JJH(I,J) = NROW-1
739          TLATHX(I,J)=TLATHC-DPH
740          TLONHX(I,J)=TLONHC-DLM
741         ELSE                                   ! NESTED POINT NORTH OF PARENT
742          KOUTB(I,J)=K+(P_IDE-1)-1
743          IIH(I,J) = NCOL-1
744          JJH(I,J) = NROW+1
745          TLATHX(I,J)=TLATHC+DPH
746          TLONHX(I,J)=TLONHC-DLM
747         ENDIF
748!***
749!***
750!***                  4
751!***
752!***                  H
753!***             1         2
754!***
755!***                  3
756!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
757
758       TLATO=TLATHX(I,J)
759       TLONO=TLONHX(I,J)
760       DLM1=TLON-TLONO
761       DLA1=TLAT-TLATO                                               ! Q
762!      DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
763       DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
764!
765       TLATO=TLATHX(I,J)
766       TLONO=TLONHX(I,J)+2.*DLM
767       DLM2=TLON-TLONO
768       DLA2=TLAT-TLATO                                               ! Q
769!      DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
770       DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
771!
772       TLATO=TLATHX(I,J)-DPH
773       TLONO=TLONHX(I,J)+DLM
774       DLM3=TLON-TLONO
775       DLA3=TLAT-TLATO                                               ! Q
776!      DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
777       DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
778 
779       TLATO=TLATHX(I,J)+DPH
780       TLONO=TLONHX(I,J)+DLM
781       DLM4=TLON-TLONO
782       DLA4=TLAT-TLATO                                               ! Q
783!      DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
784       DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q
785
786
787!      THE BILINEAR WEIGHTS
788!***
789!***
790       AN3=ATAN2(DLA1,DLM1)                                          ! Q
791       R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
792       S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
793       R1=R1/DS1
794       S1=S1/DS1
795       DL1I=(1.-R1)*(1.-S1)
796       DL2I=R1*S1
797       DL3I=R1*(1.-S1)
798       DL4I=(1.-R1)*S1
799!
800       HBWGT1(I,J)=DL1I
801       HBWGT2(I,J)=DL2I
802       HBWGT3(I,J)=DL3I
803       HBWGT4(I,J)=DL4I
804!
805      ENDIF
806
807    ELSE
808!
809!*** NON-COINCIDENT POINTS   
810!
811      SLOPE=(TLAT-TLATHC)/DENOM
812      DSLOPE=NINT(R2D*ATAN(SLOPE))
813
814      IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
815        IF(TLON.GT.TLONHC)THEN
816          IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
817          KOUTB(I,J)=K
818          IIH(I,J) = NCOL
819          JJH(I,J) = NROW
820          TLATHX(I,J)=TLATHC
821          TLONHX(I,J)=TLONHC
822        ELSE
823          IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
824          KOUTB(I,J)=K-1
825          IIH(I,J) = NCOL-2
826          JJH(I,J) = NROW
827          TLATHX(I,J)=TLATHC
828          TLONHX(I,J)=TLONHC -2.*DLM
829        ENDIF
830
831!
832      ELSEIF(DSLOPE.GT.DSLP0)THEN
833        IF(TLON.GT.TLONHC)THEN
834          IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
835          KOUTB(I,J)=K+(P_IDE-1)-1
836          IIH(I,J) = NCOL-1
837          JJH(I,J) = NROW+1
838          TLATHX(I,J)=TLATHC+DPH
839          TLONHX(I,J)=TLONHC-DLM
840        ELSE
841          IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
842          KOUTB(I,J)=K-(P_IDE-1)
843          IIH(I,J) = NCOL-1
844          JJH(I,J) = NROW-1
845          TLATHX(I,J)=TLATHC-DPH
846          TLONHX(I,J)=TLONHC-DLM
847        ENDIF
848
849!
850      ELSEIF(DSLOPE.LT.-DSLP0)THEN
851        IF(TLON.GT.TLONHC)THEN
852          IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
853          KOUTB(I,J)=K-(P_IDE-1)
854          IIH(I,J) = NCOL-1
855          JJH(I,J) = NROW-1
856          TLATHX(I,J)=TLATHC-DPH
857          TLONHX(I,J)=TLONHC-DLM
858        ELSE
859          IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
860          KOUTB(I,J)=K+(P_IDE-1)-1
861          IIH(I,J) = NCOL-1
862          JJH(I,J) = NROW+1
863          TLATHX(I,J)=TLATHC+DPH
864          TLONHX(I,J)=TLONHC-DLM
865        ENDIF
866      ENDIF
867
868!
869!***  NOW WE WILL MOVE AS FOLLOWS:
870!***
871!***
872!***                      4
873!***
874!***
875!*** 
876!***                   H
877!***             1                 2
878!***
879!***
880!***
881!***
882!***                      3
883!***
884!***
885!***
886!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
887     
888      TLATO=TLATHX(I,J)
889      TLONO=TLONHX(I,J)
890      DLM1=TLON-TLONO
891      DLA1=TLAT-TLATO                                               ! Q
892!     DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
893      DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
894!
895      TLATO=TLATHX(I,J)                                             ! redundant computations
896      TLONO=TLONHX(I,J)+2.*DLM
897      DLM2=TLON-TLONO
898      DLA2=TLAT-TLATO                                               ! Q
899!     DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
900      DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
901!
902      TLATO=TLATHX(I,J)-DPH
903      TLONO=TLONHX(I,J)+DLM
904      DLM3=TLON-TLONO
905      DLA3=TLAT-TLATO                                               ! Q
906!     DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
907      DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
908!
909      TLATO=TLATHX(I,J)+DPH
910      TLONO=TLONHX(I,J)+DLM
911      DLM4=TLON-TLONO
912      DLA4=TLAT-TLATO                                               ! Q
913!     DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
914      DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q
915
916!     THE BILINEAR WEIGHTS
917!***
918      AN3=ATAN2(DLA1,DLM1)                                          ! Q
919      R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
920      S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
921      R1=R1/DS1
922      S1=S1/DS1
923      DL1I=(1.-R1)*(1.-S1)
924      DL2I=R1*S1
925      DL3I=R1*(1.-S1)
926      DL4I=(1.-R1)*S1
927!
928      HBWGT1(I,J)=DL1I
929      HBWGT2(I,J)=DL2I
930      HBWGT3(I,J)=DL3I
931      HBWGT4(I,J)=DL4I
932!
933    ENDIF
934
935!
936!***  FINALLY STORE IIH IN TERMS OF E-GRID INDEX
937!
938      IIH(I,J)=NINT(0.5*IIH(I,J))
939
940      HBWGT1(I,J)=MAX(HBWGT1(I,J),0.0)   ! all weights must be GE zero (postive def)
941      HBWGT2(I,J)=MAX(HBWGT2(I,J),0.0)   ! all weights must be GE zero (postive def)
942      HBWGT3(I,J)=MAX(HBWGT3(I,J),0.0)   ! all weights must be GE zero (postive def)
943      HBWGT4(I,J)=MAX(HBWGT4(I,J),0.0)   ! all weights must be GE zero (postive def)
944
945!      write(message,105)"H WEIGHTS:",I,J,HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J), &
946!                               HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J),IIH(i,j),JJH(i,j)
947!      CALL wrf_message(trim(message))
948! 105  format(a,2i4,5f7.3,2i4)
949
950   ENDDO
951  ENDDO
952
953
954  RETURN
955  END SUBROUTINE G2T2H
956!========================================================================================
957
958
959  SUBROUTINE G2T2V( IIV,JJV,                     & ! output grid index and weights
960                    VBWGT1,VBWGT2,               & ! output weights in terms of parent grid
961                    VBWGT3,VBWGT4,               &
962                    VLAT,VLON,                   & ! target (nest) input lat lon in degrees
963                    DLMD1,DPHD1,WBD1,SBD1,       & ! parent res, west and south boundaries
964                    CENTRAL_LAT,CENTRAL_LON,     & ! parent central lat,lon, all in degrees
965                    P_IDE,P_JDE,                 & ! parent imax and jmax
966                    IDS,IDE,JDS,JDE,KDS,KDE,     & ! target (nest) dIMEnsions
967                    IMS,IME,JMS,JME,KMS,KME,     &
968                    ITS,ITE,JTS,JTE,KTS,KTE      )
969
970!
971!***  Tom Black   - Initial Version
972!***  Gopal       - Revised Version for WRF (includes coincIDEnt grid points)
973!***
974!***  GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
975!***  AND THE NESTED GRID LAT-LONS AT V POINTS, THIS ROUTINE FIRST LOCATES THE
976!***  INDICES,IIV,JJV, OF THE PARENT DOMAIN'S V POINTS THAT LIES CLOSEST TO THE
977!***  V POINTS OF THE NESTED DOMAIN
978!
979!============================================================================
980
981 IMPLICIT NONE
982 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
983 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
984 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
985 INTEGER,    INTENT(IN   )                            :: P_IDE,P_JDE
986 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
987 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
988 REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(IN)    :: VLAT,VLON
989 REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(OUT)   :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
990 INTEGER, DIMENSION(IMS:IME,JMS:JME),   INTENT(OUT)   :: IIV,JJV
991
992! local
993
994 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)     ! (6) single precision
995 INTEGER           :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
996 INTEGER           :: NROW,NCOL,KROWS
997 REAL(KIND=KNUM)   :: X,Y,Z,TLAT,TLON
998 REAL(KIND=KNUM)   :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
999 REAL(KIND=KNUM)   :: ROW,COL,SLP0,TLATVC,TLONVC,DENOM,SLOPE
1000 REAL(KIND=KNUM)   :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
1001 REAL(KIND=KNUM)   :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3                    ! Q
1002 REAL(KIND=KNUM)   :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
1003 REAL(KIND=KNUM)   :: DTEMP
1004 REAL  , DIMENSION(IMS:IME,JMS:JME)      :: TLATVX,TLONVX
1005 INTEGER, DIMENSION(IMS:IME,JMS:JME)     :: KOUTB
1006!-------------------------------------------------------------------------------------
1007
1008  IMT=2*P_IDE-2             ! parent i dIMEnsions
1009  JMT=P_JDE/2               ! parent j dIMEnsions
1010  PI_2=ACOS(0.)
1011  D2R=PI_2/90.
1012  R2D=1./D2R
1013  DPH=DPHD1*D2R
1014  DLM=DLMD1*D2R
1015  TPH0= CENTRAL_LAT*D2R
1016  TLM0=-CENTRAL_LON*D2R        ! NOTE THE MINUS HERE
1017  WB=WBD1*D2R                   ! DEGREES TO RADIANS
1018  SB=SBD1*D2R
1019  SLP0=DPHD1/DLMD1
1020  DSLP0=NINT(R2D*ATAN(SLP0))
1021  DS1=SQRT(DPH*DPH+DLM*DLM)    ! Q
1022  AN1=ASIN(DLM/DS1)
1023  AN2=ASIN(DPH/DS1)
1024
1025  DO J =  JTS,MIN(JTE,JDE-1)
1026    DO I = ITS,MIN(ITE,IDE-1)
1027!***
1028!***  LOCATE TARGET V POINTS (VLAT AND VLON) ON THE PARENT DOMAIN AND
1029!***  DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
1030!***  CONVERT NESTED GRID V POINTS FROM GEODETIC TO TRANSFORMED
1031!***  COORDINATE ON THE PARENT GRID
1032!
1033
1034      GLAT=VLAT(I,J)*D2R
1035      GLON=(360. - VLON(I,J))*D2R
1036      X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
1037      Y=-COS(GLAT)*SIN(GLON-TLM0)
1038      Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
1039      TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
1040      TLON=R2D*ATAN(Y/X)
1041
1042!      ROW=TLAT/DPHD1+JMT         ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
1043!      COL=TLON/DLMD1+P_IDE-1     ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
1044
1045      ROW=(TLAT-SBD1)/DPHD1+1     ! Dusan's doing
1046      COL=(TLON-WBD1)/DLMD1+1     ! Dusan's doing
1047
1048      NROW=INT(ROW + 0.001)     ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
1049      NCOL=INT(COL + 0.001)
1050      TLAT=TLAT*D2R
1051      TLON=TLON*D2R
1052
1053!***
1054!***
1055!***  FIRST CONSIDER THE SITUATION WHERE THE POINT V IS AT
1056!***
1057!***              H      V
1058!***
1059!***
1060!***                 V
1061!***              V      H
1062!***
1063!***  THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
1064!***
1065
1066      IF(MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.1.OR.     &
1067         MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.0)THEN
1068           TLAT1=(NROW-JMT)*DPH
1069           TLAT2=TLAT1+DPH
1070           TLON1=(NCOL-(P_IDE-1))*DLM
1071           TLON2=TLON1+DLM
1072           DLM1=TLON-TLON1
1073           DLM2=TLON-TLON2
1074!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1075!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1076           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1077           D1=ACOS(DTEMP)
1078           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1079           D2=ACOS(DTEMP)
1080            IF(D1.GT.D2)THEN
1081             NROW=NROW+1                    ! FIND THE NEAREST V ROW
1082             NCOL=NCOL+1                    ! FIND THE NEAREST V COLUMN
1083            ENDIF
1084      ELSE
1085
1086!***
1087!***  NOW CONSIDER THE SITUATION WHERE THE POINT V IS AT
1088!***
1089!***              V      H
1090!***
1091!***
1092!***                 V
1093!***              H      V
1094!***
1095!*** THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
1096!***
1097           TLAT1=(NROW+1-JMT)*DPH
1098           TLAT2=TLAT1-DPH
1099           TLON1=(NCOL-(P_IDE-1))*DLM
1100           TLON2=TLON1+DLM
1101           DLM1=TLON-TLON1
1102           DLM2=TLON-TLON2
1103!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1104!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1105           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
1106           D1=ACOS(DTEMP)
1107           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
1108           D2=ACOS(DTEMP)
1109             IF(D1.LT.D2)THEN
1110              NROW=NROW+1                    ! FIND THE NEAREST H ROW
1111             ELSE
1112              NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
1113             ENDIF
1114
1115      ENDIF
1116
1117      KROWS=((NROW-1)/2)*IMT
1118      IF(MOD(NROW,2).EQ.1)THEN
1119        K=KROWS+NCOL/2
1120      ELSE
1121        K=KROWS+P_IDE-2+(NCOL+1)/2     ! check this one should this not be P_IDE-2 ????
1122      ENDIF
1123
1124!***
1125!***  WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
1126!***  NEAREST TO THE CENTER K AS SEEN BELOW.  WE MUST FIND
1127!***  WHICH OF THE FOUR v-BOXES (OF WHICH THIS V POINT IS
1128!***  A VERTEX) SURROUNDS THE INNER GRID V POINT IN QUESTION.
1129!***
1130!***
1131!***                  V
1132!***
1133!***
1134!***
1135!***            V     H     V
1136!***
1137!***
1138!***               V
1139!***      V     H     V     H     V
1140!***
1141!***
1142!***
1143!***            V     H     V
1144!***
1145!***
1146!***
1147!***                  V
1148!***
1149!***
1150!***  FIND THE SLOPE OF THE LINE CONNECTING V AND THE CENTER v.
1151!***
1152      N2R=K/IMT
1153      MK=MOD(K,IMT)
1154!
1155      IF(MK.EQ.0)THEN
1156        TLATVC=SB+(2*N2R-1)*DPH
1157      ELSE
1158        TLATVC=SB+(2*N2R+MK/(P_IDE-1))*DPH
1159      ENDIF
1160!
1161      IF(MK.LE.(P_IDE-1)-1)THEN
1162        TLONVC=WB+(2*MK-1)*DLM
1163      ELSE
1164        TLONVC=WB+2*(MK-(P_IDE-1))*DLM
1165      ENDIF
1166
1167!
1168!***  EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
1169!***  DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE V BOXES, WE NEED TO BE
1170!***  CAREFUL HERE
1171!
1172       IF(ABS(TLON-TLONVC) .LE. 1.E-4)TLONVC=TLON
1173       IF(ABS(TLAT-TLATVC) .LE. 1.E-4)TLATVC=TLAT
1174       DENOM=(TLON-TLONVC)
1175!
1176!***
1177!***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
1178!***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
1179!***
1180!*** COINCIDENT CONDITIONS
1181
1182     IF(DENOM.EQ.0.0)THEN
1183
1184       IF(TLATVC.EQ.TLAT)THEN
1185         KOUTB(I,J)=K
1186         IIV(I,J) = NCOL
1187         JJV(I,J) = NROW
1188         TLATVX(I,J)=TLATVC
1189         TLONVX(I,J)=TLONVC
1190         VBWGT1(I,J)=1.0
1191         VBWGT2(I,J)=0.0
1192         VBWGT3(I,J)=0.0
1193         VBWGT4(I,J)=0.0
1194       ELSE                              ! SAME LONGITUDE BUT DIFFERENT LATS
1195         
1196         IF(TLATVC .GT. TLAT)THEN      ! NESTED POINT SOUTH OF PARENT
1197          KOUTB(I,J)=K-(P_IDE-1)
1198          IIV(I,J) = NCOL-1
1199          JJV(I,J) = NROW-1
1200          TLATVX(I,J)=TLATVC-DPH
1201          TLONVX(I,J)=TLONVC-DLM
1202         ELSE                                   ! NESTED POINT NORTH OF PARENT
1203          KOUTB(I,J)=K+(P_IDE-1)-1
1204          IIV(I,J) = NCOL-1
1205          JJV(I,J) = NROW+1
1206          TLATVX(I,J)=TLATVC+DPH
1207          TLONVX(I,J)=TLONVC-DLM
1208         ENDIF
1209
1210!***
1211!***
1212!***                  4
1213!***
1214!***                  V
1215!***             1         2
1216!***
1217!***                  3
1218!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
1219
1220       TLATO=TLATVX(I,J)
1221       TLONO=TLONVX(I,J)
1222       DLM1=TLON-TLONO
1223       DLA1=TLAT-TLATO                                               ! Q
1224!      DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
1225       DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
1226!
1227       TLATO=TLATVX(I,J)
1228       TLONO=TLONVX(I,J)+2.*DLM
1229       DLM2=TLON-TLONO
1230       DLA2=TLAT-TLATO                                               ! Q
1231!      DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
1232       DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
1233
1234       TLATO=TLATVX(I,J)-DPH
1235       TLONO=TLONVX(I,J)+DLM
1236       DLM3=TLON-TLONO
1237       DLA3=TLAT-TLATO                                               ! Q
1238!      DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
1239       DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
1240
1241       TLATO=TLATVX(I,J)+DPH
1242       TLONO=TLONVX(I,J)+DLM
1243       DLM4=TLON-TLONO
1244       DLA4=TLAT-TLATO                                               ! Q
1245!      DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
1246       DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q
1247                 
1248!      THE BILINEAR WEIGHTS
1249!***
1250       AN3=ATAN2(DLA1,DLM1)                                          ! Q
1251       R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
1252       S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
1253       R1=R1/DS1
1254       S1=S1/DS1
1255       DL1I=(1.-R1)*(1.-S1)
1256       DL2I=R1*S1
1257       DL3I=R1*(1.-S1)
1258       DL4I=(1.-R1)*S1
1259!
1260       VBWGT1(I,J)=DL1I
1261       VBWGT2(I,J)=DL2I
1262       VBWGT3(I,J)=DL3I
1263       VBWGT4(I,J)=DL4I     
1264
1265      ENDIF
1266
1267    ELSE
1268
1269!
1270!*** NON-COINCIDENT POINTS
1271!
1272      SLOPE=(TLAT-TLATVC)/DENOM
1273      DSLOPE=NINT(R2D*ATAN(SLOPE))
1274
1275      IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
1276        IF(TLON.GT.TLONVC)THEN
1277          IF(TLONVC.GE.-WB-DLM)CALL wrf_error_fatal("1V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1278          KOUTB(I,J)=K
1279          IIV(I,J)=NCOL
1280          JJV(I,J)=NROW
1281          TLATVX(I,J)=TLATVC
1282          TLONVX(I,J)=TLONVC
1283        ELSE
1284          IF(TLONVC.LE.WB+DLM)CALL wrf_error_fatal("2V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1285          KOUTB(I,J)=K-1
1286          IIV(I,J) = NCOL-2
1287          JJV(I,J) = NROW
1288          TLATVX(I,J)=TLATVC
1289          TLONVX(I,J)=TLONVC-2.*DLM
1290        ENDIF
1291 
1292      ELSEIF(DSLOPE.GT.DSLP0)THEN
1293        IF(TLON.GT.TLONVC)THEN
1294          IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("3V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1295          KOUTB(I,J)=K+(P_IDE-1)-1
1296          IIV(I,J) = NCOL-1
1297          JJV(I,J) = NROW+1
1298          TLATVX(I,J)=TLATVC+DPH
1299          TLONVX(I,J)=TLONVC-DLM
1300        ELSE
1301          IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("4V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1302          KOUTB(I,J)=K-(P_IDE-1)
1303          IIV(I,J) = NCOL-1
1304          JJV(I,J) = NROW-1
1305          TLATVX(I,J)=TLATVC-DPH
1306          TLONVX(I,J)=TLONVC-DLM
1307        ENDIF
1308 
1309      ELSEIF(DSLOPE.LT.-DSLP0)THEN
1310        IF(TLON.GT.TLONVC)THEN
1311          IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("5V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1312          KOUTB(I,J)=K-(P_IDE-1)
1313          IIV(I,J) = NCOL-1
1314          JJV(I,J) = NROW-1
1315          TLATVX(I,J)=TLATVC-DPH
1316          TLONVX(I,J)=TLONVC-DLM
1317        ELSE
1318          IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("6V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
1319          KOUTB(I,J)=K+(P_IDE-1)-1
1320          IIV(I,J) = NCOL-1
1321          JJV(I,J) = NROW+1
1322          TLATVX(I,J)=TLATVC+DPH
1323          TLONVX(I,J)=TLONVC-DLM
1324        ENDIF
1325      ENDIF
1326!
1327!***  NOW WE WILL MOVE AS FOLLOWS:
1328!***
1329!***
1330!***                      4
1331!***
1332!***
1333!***
1334!***                   V
1335!***             1                 2
1336!***
1337!***
1338!***
1339!***
1340!***                      3
1341!***
1342!***
1343!***
1344!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM V TO EACH VERTEX
1345
1346      TLATO=TLATVX(I,J)
1347      TLONO=TLONVX(I,J)
1348      DLM1=TLON-TLONO
1349      DLA1=TLAT-TLATO                                               ! Q
1350!     DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
1351      DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
1352!
1353      TLATO=TLATVX(I,J)
1354      TLONO=TLONVX(I,J)+2.*DLM
1355      DLM2=TLON-TLONO
1356      DLA2=TLAT-TLATO                                               ! Q
1357!     DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
1358      DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
1359!
1360      TLATO=TLATVX(I,J)-DPH
1361      TLONO=TLONVX(I,J)+DLM
1362      DLM3=TLON-TLONO
1363      DLA3=TLAT-TLATO                                               ! Q
1364!     DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
1365      DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
1366!
1367      TLATO=TLATVX(I,J)+DPH
1368      TLONO=TLONVX(I,J)+DLM
1369      DLM4=TLON-TLONO
1370      DLA4=TLAT-TLATO                                               ! Q
1371!     DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
1372      DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q
1373 
1374!     THE BILINEAR WEIGHTS
1375!***
1376      AN3=ATAN2(DLA1,DLM1)                                          ! Q
1377      R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
1378      S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
1379      R1=R1/DS1
1380      S1=S1/DS1
1381      DL1I=(1.-R1)*(1.-S1)
1382      DL2I=R1*S1
1383      DL3I=R1*(1.-S1)
1384      DL4I=(1.-R1)*S1
1385!
1386      VBWGT1(I,J)=DL1I
1387      VBWGT2(I,J)=DL2I
1388      VBWGT3(I,J)=DL3I
1389      VBWGT4(I,J)=DL4I
1390
1391     ENDIF
1392
1393!
1394!***  FINALLY STORE IIH IN TERMS OF E-GRID INDEX
1395!
1396      IIV(I,J)=NINT(0.5*IIV(I,J))
1397
1398      VBWGT1(I,J)=MAX(VBWGT1(I,J),0.0)   ! all weights must be GE zero (postive def)
1399      VBWGT2(I,J)=MAX(VBWGT2(I,J),0.0)   ! all weights must be GE zero (postive def)
1400      VBWGT3(I,J)=MAX(VBWGT3(I,J),0.0)   ! all weights must be GE zero (postive def)
1401      VBWGT4(I,J)=MAX(VBWGT4(I,J),0.0)   ! all weights must be GE zero (postive def)
1402
1403    ENDDO
1404  ENDDO
1405
1406 RETURN
1407 END SUBROUTINE G2T2V
1408
1409!----------------------------------------------------------------------------- 
1410
1411 SUBROUTINE G2T2H_new( IIH,JJH,                            & ! output grid index and weights
1412                       HBWGT1,HBWGT2,                      & ! output weights in terms of parent grid
1413                       HBWGT3,HBWGT4,                      &
1414                       I_PARENT_START,J_PARENT_START,      & ! nest start I and J in parent domain 
1415                       RATIO,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
1416                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
1417                       IMS,IME,JMS,JME,KMS,KME,            &
1418                       ITS,ITE,JTS,JTE,KTS,KTE      )
1419!
1420!*** XUEJIN ZHANG --- Initial version (09/08/2008)
1421!*** XUEJIN ZHANG --- Modified for parallel purpose (09/10/2009)
1422!*** XUEJIN ZHANG --- Modified for parallel purpose (09/29/2009)
1423!Function: Bilnear interpolation weight and indexing for E-grid
1424!
1425 IMPLICIT NONE
1426 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1427 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1428 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1429 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1430 INTEGER,    INTENT(IN   )                            :: RATIO
1431 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1432 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1433!LOCAL VARIABLES
1434 INTEGER                                              :: I,J
1435 INTEGER                                              :: JP
1436
1437
1438!*** ARRANGEMENT OF 4 VERTEXES FROM PARENT DOMAIN
1439!***
1440!***                  4
1441!***
1442!***                  h
1443!***             1         2
1444!***
1445!***
1446!***                  3
1447!
1448!*************************************************************
1449!***  POINT (i,j) SPANS 9 NESTED POINTS
1450!***  A VERTEX IN THE NEST DOMAIN AND INDEXING FOR PARENT DOMAIN
1451!***
1452!***
1453!***                  H
1454!***               
1455!***               h     h
1456!***              / \
1457!***            h     h     h
1458!***           / \   / \
1459!***         H     h     h     H 
1460!***          \   / \   /   
1461!***            h     h     h
1462!***             \   /
1463!***               h     h
1464!***     
1465!***                  H
1466!***                 
1467!***
1468!***
1469!*************************************************************
1470!*** MOVING NEST ALWAYS STARTING FROM MASS GRID H
1471!*** PLEASE REFER TO E-GRID ARRANGEMENT FIGURE FOR MASS POINTS
1472
1473
1474 DO J=JTS,MIN(JTE,JDE-1)
1475 DO I=ITS,MIN(ITE,IDE-1)
1476    JP = MOD(J,RATIO*2)
1477    SELECT CASE(JP)
1478    CASE ( 1 )
1479      CALL SUB1H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1480                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1481                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1482                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1483                IMS,IME,JMS,JME,KMS,KME,                 &
1484                ITS,ITE,JTS,JTE,KTS,KTE      )
1485    CASE ( 2 )
1486      CALL SUB2H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &   
1487                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1488                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1489                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1490                IMS,IME,JMS,JME,KMS,KME,                 &
1491                ITS,ITE,JTS,JTE,KTS,KTE      )
1492    CASE ( 3 )
1493      CALL SUB3H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1494                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1495                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1496                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1497                IMS,IME,JMS,JME,KMS,KME,                 &
1498                ITS,ITE,JTS,JTE,KTS,KTE      )
1499    CASE ( 4 )
1500      CALL SUB4H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &   
1501                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1502                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1503                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1504                IMS,IME,JMS,JME,KMS,KME,                 &
1505                ITS,ITE,JTS,JTE,KTS,KTE      )
1506    CASE ( 5 )
1507      CALL SUB5H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1508                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1509                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1510                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1511                IMS,IME,JMS,JME,KMS,KME,                 &
1512                ITS,ITE,JTS,JTE,KTS,KTE      )
1513    CASE ( 0 )
1514      CALL SUB6H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &   
1515                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1516                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1517                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1518                IMS,IME,JMS,JME,KMS,KME,                 &
1519                ITS,ITE,JTS,JTE,KTS,KTE      )
1520    END SELECT
1521       write(0,105)"NEW H WEIGHTS:",I,J,IIH(i,j),JJH(i,j),HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J), &
1522                                HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J)
1523  105  format(a,4i4,5f7.3)
1524 END DO
1525 END DO
1526
1527 RETURN
1528 END SUBROUTINE G2T2H_new
1529 SUBROUTINE SUB1H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1530                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1531                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1532                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1533                 IMS,IME,JMS,JME,KMS,KME,                 &
1534                 ITS,ITE,JTS,JTE,KTS,KTE      )
1535 IMPLICIT NONE
1536 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1537 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1538 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1539 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1540 INTEGER,    INTENT(IN   )                            :: RATIO
1541 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1542 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1543!LOCAL VARIABLES
1544 INTEGER                                              :: I,J
1545 INTEGER                                              :: IP
1546
1547  IP = MOD(I,RATIO)
1548  SELECT CASE(IP)
1549   CASE ( 1 )
1550    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1551    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1552    HBWGT1(I,J) = 1.0
1553    HBWGT2(I,J) = 0.0
1554    HBWGT3(I,J) = 0.0
1555    HBWGT4(I,J) = 0.0
1556   CASE ( 2 )
1557    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1558    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1559    HBWGT1(I,J) = 4.0/9.0
1560    HBWGT2(I,J) = 1.0/9.0
1561    HBWGT3(I,J) = 2.0/9.0
1562    HBWGT4(I,J) = 2.0/9.0
1563   CASE ( 0 )
1564    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1565    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1566    HBWGT1(I,J) = 1.0/9.0
1567    HBWGT2(I,J) = 4.0/9.0
1568    HBWGT3(I,J) = 2.0/9.0
1569    HBWGT4(I,J) = 2.0/9.0
1570   END SELECT
1571 RETURN
1572 END SUBROUTINE SUB1H
1573 SUBROUTINE SUB2H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1574                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1575                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1576                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1577                 IMS,IME,JMS,JME,KMS,KME,                 &
1578                 ITS,ITE,JTS,JTE,KTS,KTE      )
1579 IMPLICIT NONE
1580 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1581 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1582 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1583 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1584 INTEGER,    INTENT(IN   )                            :: RATIO
1585 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1586 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1587!LOCAL VARIABLES
1588 INTEGER                                              :: I,J
1589 INTEGER                                              :: IP
1590
1591  IP = MOD(I,RATIO)
1592  SELECT CASE(IP)
1593   CASE ( 1 )
1594    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1595    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1596    HBWGT1(I,J) = 2.0/3.0
1597    HBWGT2(I,J) = 0.0
1598    HBWGT3(I,J) = 0.0
1599    HBWGT4(I,J) = 1.0/3.0
1600   CASE ( 2 )
1601    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1602    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1603    HBWGT1(I,J) = 2.0/9.0
1604    HBWGT2(I,J) = 2.0/9.0
1605    HBWGT3(I,J) = 1.0/9.0
1606    HBWGT4(I,J) = 4.0/9.0
1607   CASE ( 0 )
1608    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1609    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)+1
1610    HBWGT1(I,J) = 1.0/3.0
1611    HBWGT2(I,J) = 0.0
1612    HBWGT3(I,J) = 2.0/3.0
1613    HBWGT4(I,J) = 0.0
1614   END SELECT
1615 RETURN
1616 END SUBROUTINE SUB2H
1617 SUBROUTINE SUB3H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1618                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1619                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1620                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1621                 IMS,IME,JMS,JME,KMS,KME,                 &
1622                 ITS,ITE,JTS,JTE,KTS,KTE      )
1623 IMPLICIT NONE
1624 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1625 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1626 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1627 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1628 INTEGER,    INTENT(IN   )                            :: RATIO
1629 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1630 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1631!LOCAL VARIABLES
1632 INTEGER                                              :: I,J
1633 INTEGER                                              :: IP
1634
1635  IP = MOD(I,RATIO)
1636  SELECT CASE(IP)
1637   CASE ( 1 )
1638    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)-1
1639    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)+1
1640    HBWGT1(I,J) = 2.0/9.0
1641    HBWGT2(I,J) = 2.0/9.0
1642    HBWGT3(I,J) = 4.0/9.0
1643    HBWGT4(I,J) = 1.0/9.0
1644   CASE ( 2 )
1645    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1646    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1647    HBWGT1(I,J) = 1.0/3.0
1648    HBWGT2(I,J) = 0.0
1649    HBWGT3(I,J) = 0.0
1650    HBWGT4(I,J) = 2.0/3.0
1651   CASE ( 0 )
1652    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1653    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)+1
1654    HBWGT1(I,J) = 2.0/3.0
1655    HBWGT2(I,J) = 0.0
1656    HBWGT3(I,J) = 1.0/3.0
1657    HBWGT4(I,J) = 0.0
1658   END SELECT
1659 RETURN
1660 END SUBROUTINE SUB3H
1661 SUBROUTINE SUB4H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1662                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1663                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1664                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1665                 IMS,IME,JMS,JME,KMS,KME,                 &
1666                 ITS,ITE,JTS,JTE,KTS,KTE      )
1667 IMPLICIT NONE
1668 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1669 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1670 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1671 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1672 INTEGER,    INTENT(IN   )                            :: RATIO
1673 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1674 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1675!LOCAL VARIABLES
1676 INTEGER                                              :: I,J
1677 INTEGER                                              :: IP
1678
1679  IP = MOD(I,RATIO)
1680  SELECT CASE(IP)
1681   CASE ( 1 )
1682    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)-1
1683    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1684    HBWGT1(I,J) = 1.0/9.0
1685    HBWGT2(I,J) = 4.0/9.0
1686    HBWGT3(I,J) = 2.0/9.0
1687    HBWGT4(I,J) = 2.0/9.0
1688   CASE ( 2 )
1689    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1690    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1691    HBWGT1(I,J) = 1.0
1692    HBWGT2(I,J) = 0.0
1693    HBWGT3(I,J) = 0.0
1694    HBWGT4(I,J) = 0.0
1695   CASE ( 0 )
1696    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1697    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1698    HBWGT1(I,J) = 4.0/9.0
1699    HBWGT2(I,J) = 1.0/9.0
1700    HBWGT3(I,J) = 2.0/9.0
1701    HBWGT4(I,J) = 2.0/9.0
1702   END SELECT
1703 RETURN
1704 END SUBROUTINE SUB4H
1705 SUBROUTINE SUB5H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1706                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1707                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1708                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1709                 IMS,IME,JMS,JME,KMS,KME,                 &
1710                 ITS,ITE,JTS,JTE,KTS,KTE      )
1711 IMPLICIT NONE
1712 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1713 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1714 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1715 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1716 INTEGER,    INTENT(IN   )                            :: RATIO
1717 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1718 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1719!LOCAL VARIABLES
1720 INTEGER                                              :: I,J
1721 INTEGER                                              :: IP
1722
1723  IP = MOD(I,RATIO)
1724  SELECT CASE(IP)
1725   CASE ( 1 )
1726    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)-1
1727    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1728    HBWGT1(I,J) = 2.0/9.0
1729    HBWGT2(I,J) = 2.0/9.0
1730    HBWGT3(I,J) = 1.0/9.0
1731    HBWGT4(I,J) = 4.0/9.0
1732   CASE ( 2 )
1733    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1734    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)+1
1735    HBWGT1(I,J) = 1.0/3.0
1736    HBWGT2(I,J) = 0.0
1737    HBWGT3(I,J) = 2.0/3.0
1738    HBWGT4(I,J) = 0.0
1739   CASE ( 0 )
1740    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1741    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1742    HBWGT1(I,J) = 2.0/3.0
1743    HBWGT2(I,J) = 0.0
1744    HBWGT3(I,J) = 0.0
1745    HBWGT4(I,J) = 1.0/3.0
1746   END SELECT
1747 RETURN
1748 END SUBROUTINE SUB5H
1749 SUBROUTINE SUB6H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
1750                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1751                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1752                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1753                 IMS,IME,JMS,JME,KMS,KME,                 &
1754                 ITS,ITE,JTS,JTE,KTS,KTE      )
1755 IMPLICIT NONE
1756 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1757 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1758 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1759 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1760 INTEGER,    INTENT(IN   )                            :: RATIO
1761 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
1762 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH
1763!LOCAL VARIABLES
1764 INTEGER                                              :: I,J
1765 INTEGER                                              :: IP
1766
1767  IP = MOD(I,RATIO)
1768  SELECT CASE(IP)
1769   CASE ( 1 )
1770    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1771    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
1772    HBWGT1(I,J) = 2.0/3.0
1773    HBWGT2(I,J) = 0.0
1774    HBWGT3(I,J) = 1.0/3.0
1775    HBWGT4(I,J) = 0.0
1776   CASE ( 2 )
1777    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1778    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
1779    HBWGT1(I,J) = 2.0/9.0
1780    HBWGT2(I,J) = 2.0/9.0
1781    HBWGT3(I,J) = 4.0/9.0
1782    HBWGT4(I,J) = 1.0/9.0
1783   CASE ( 0 )
1784    IIH(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1785    JJH(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1786    HBWGT1(I,J) = 1.0/3.0
1787    HBWGT2(I,J) = 0.0
1788    HBWGT3(I,J) = 0.0
1789    HBWGT4(I,J) = 2.0/3.0
1790   END SELECT
1791 RETURN
1792 END SUBROUTINE SUB6H
1793
1794!----------------------------------------------------------------------------- 
1795
1796 SUBROUTINE G2T2V_new( IIV,JJV,                            & ! output grid index and weights
1797                       VBWGT1,VBWGT2,                      & ! output weights in terms of parent grid
1798                       VBWGT3,VBWGT4,                      &
1799                       I_PARENT_START,J_PARENT_START,      & ! nest start I and J in parent domain 
1800                       RATIO,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
1801                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
1802                       IMS,IME,JMS,JME,KMS,KME,            &
1803                       ITS,ITE,JTS,JTE,KTS,KTE      )
1804!
1805!*** XUEJIN ZHANG --- Initial version (09/08/2008)
1806!*** XUEJIN ZHANG --- Modified for parallel purpose (09/10/2009)
1807!*** XUEJIN ZHANG --- Modified for parallel purpose (09/29/2009)
1808!Function: Bilnear interpolation weight and indexing for E-grid
1809!
1810 IMPLICIT NONE
1811 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1812 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1813 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1814 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1815 INTEGER,    INTENT(IN   )                            :: RATIO
1816 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
1817 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
1818!LOCAL VARIABLES
1819 INTEGER                                              :: I,J
1820 INTEGER                                              :: JP
1821
1822!*** ARRANGEMENT OF 4 VERTEXES FROM PARENT DOMAIN
1823!***
1824!***                  4
1825!***
1826!***                  v
1827!***             1         2
1828!***
1829!***
1830!***                  3
1831!
1832!*************************************************************
1833!***  POINT (i,j) SPANS 9 NESTED POINTS
1834!***  A VERTEX IN THE NEST DOMAIN AND INDEXING FOR PARENT DOMAIN
1835!***
1836!***
1837!***                  V
1838!***               
1839!***               v  h  v
1840!***              / \
1841!***            v  h  v  h  v
1842!***           / \   / \
1843!***         V  H  v  h  v  h  V  H
1844!***          \   / \   /   
1845!***            v  h  v  h  v
1846!***             \   /
1847!***               v  h  v
1848!***     
1849!***                  V
1850!***                 
1851!***
1852!***
1853!*************************************************************
1854!*** MOVING NEST ALWAYS STARTING FROM MASS GRID H
1855!*** PLEASE REFER TO E-GRID ARRANGEMENT FIGURE FOR WIND POINTS
1856
1857 DO J=JTS,MIN(JTE,JDE-1)
1858 DO I=ITS,MIN(ITE,IDE-1)
1859    JP = MOD(J,RATIO*2)
1860    SELECT CASE(JP)
1861    CASE ( 1 )
1862      CALL SUB1V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1863                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1864                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1865                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1866                IMS,IME,JMS,JME,KMS,KME,                 &
1867                ITS,ITE,JTS,JTE,KTS,KTE      )
1868    CASE ( 2 )
1869      CALL SUB2V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &   
1870                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1871                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1872                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1873                IMS,IME,JMS,JME,KMS,KME,                 &
1874                ITS,ITE,JTS,JTE,KTS,KTE      )
1875    CASE ( 3 )
1876      CALL SUB3V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1877                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1878                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1879                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1880                IMS,IME,JMS,JME,KMS,KME,                 &
1881                ITS,ITE,JTS,JTE,KTS,KTE      )
1882    CASE ( 4 )
1883      CALL SUB4V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &   
1884                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1885                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1886                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1887                IMS,IME,JMS,JME,KMS,KME,                 &
1888                ITS,ITE,JTS,JTE,KTS,KTE      )
1889    CASE ( 5 )
1890      CALL SUB5V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1891                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1892                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1893                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1894                IMS,IME,JMS,JME,KMS,KME,                 &
1895                ITS,ITE,JTS,JTE,KTS,KTE      )
1896    CASE ( 0 )
1897      CALL SUB6V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &   
1898                I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1899                RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1900                IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1901                IMS,IME,JMS,JME,KMS,KME,                 &
1902                ITS,ITE,JTS,JTE,KTS,KTE      )
1903    END SELECT
1904!      WRITE(0,105)'NEW V WEIGHTS:',I,J,VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J), &
1905!                    VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J),IIV(i,j),JJV(i,j)
1906! 105  format(a,2i4,5f7.3,2i4)
1907 END DO
1908 END DO
1909
1910 RETURN
1911 END SUBROUTINE G2T2V_new
1912 SUBROUTINE SUB1V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1913                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1914                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1915                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1916                 IMS,IME,JMS,JME,KMS,KME,                 &
1917                 ITS,ITE,JTS,JTE,KTS,KTE      )
1918 IMPLICIT NONE
1919 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1920 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1921 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1922 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1923 INTEGER,    INTENT(IN   )                            :: RATIO
1924 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
1925 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
1926!LOCAL VARIABLES
1927 INTEGER                                              :: I,J
1928 INTEGER                                              :: IP
1929  IP = MOD(I,RATIO)
1930  SELECT CASE(IP)
1931   CASE ( 1 )
1932    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO) - 1
1933    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1934    VBWGT1(I,J) = 1.0/9.0
1935    VBWGT2(I,J) = 4.0/9.0
1936    VBWGT3(I,J) = 2.0/9.0
1937    VBWGT4(I,J) = 2.0/9.0
1938   CASE ( 2 )
1939    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1940    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1941    VBWGT1(I,J) = 1.0
1942    VBWGT2(I,J) = 0.0
1943    VBWGT3(I,J) = 0.0
1944    VBWGT4(I,J) = 0.0
1945   CASE ( 0 )
1946    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1947    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1948    VBWGT1(I,J) = 4.0/9.0
1949    VBWGT2(I,J) = 1.0/9.0
1950    VBWGT3(I,J) = 2.0/9.0
1951    VBWGT4(I,J) = 2.0/9.0
1952   END SELECT
1953 RETURN
1954 END SUBROUTINE SUB1V
1955 SUBROUTINE SUB2V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1956                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
1957                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
1958                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
1959                 IMS,IME,JMS,JME,KMS,KME,                 &
1960                 ITS,ITE,JTS,JTE,KTS,KTE      )
1961 IMPLICIT NONE
1962 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
1963 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
1964 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
1965 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
1966 INTEGER,    INTENT(IN   )                            :: RATIO
1967 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
1968 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
1969!LOCAL VARIABLES
1970 INTEGER                                              :: I,J
1971 INTEGER                                              :: IP
1972  IP = MOD(I,RATIO)
1973  SELECT CASE(IP)
1974   CASE ( 1 )
1975    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO) - 1
1976    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1977    VBWGT1(I,J) = 2.0/9.0
1978    VBWGT2(I,J) = 2.0/9.0
1979    VBWGT3(I,J) = 1.0/9.0
1980    VBWGT4(I,J) = 4.0/9.0
1981   CASE ( 2 )
1982    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1983    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
1984    VBWGT1(I,J) = 1.0/3.0
1985    VBWGT2(I,J) = 0.0
1986    VBWGT3(I,J) = 2.0/3.0
1987    VBWGT4(I,J) = 0.0
1988   CASE ( 0 )
1989    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
1990    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
1991    VBWGT1(I,J) = 2.0/3.0
1992    VBWGT2(I,J) = 0.0
1993    VBWGT3(I,J) = 0.0
1994    VBWGT4(I,J) = 1.0/3.0
1995   END SELECT
1996 RETURN
1997 END SUBROUTINE SUB2V
1998 SUBROUTINE SUB3V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
1999                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
2000                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
2001                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
2002                 IMS,IME,JMS,JME,KMS,KME,                 &
2003                 ITS,ITE,JTS,JTE,KTS,KTE      )
2004 IMPLICIT NONE
2005 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
2006 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
2007 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
2008 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
2009 INTEGER,    INTENT(IN   )                            :: RATIO
2010 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
2011 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
2012!LOCAL VARIABLES
2013 INTEGER                                              :: I,J
2014 INTEGER                                              :: IP
2015  IP = MOD(I,RATIO)
2016  SELECT CASE(IP)
2017   CASE ( 1 )
2018    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2019    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
2020    VBWGT1(I,J) = 2.0/3.0
2021    VBWGT2(I,J) = 0.0
2022    VBWGT3(I,J) = 1.0/3.0
2023    VBWGT4(I,J) = 0.0
2024   CASE ( 2 )
2025    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2026    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
2027    VBWGT1(I,J) = 2.0/9.0
2028    VBWGT2(I,J) = 2.0/9.0
2029    VBWGT3(I,J) = 4.0/9.0
2030    VBWGT4(I,J) = 1.0/9.0
2031   CASE ( 0 )
2032    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2033    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2034    VBWGT1(I,J) = 1.0/3.0
2035    VBWGT2(I,J) = 0.0
2036    VBWGT3(I,J) = 0.0
2037    VBWGT4(I,J) = 2.0/3.0
2038   END SELECT
2039 RETURN
2040 END SUBROUTINE SUB3V
2041 SUBROUTINE SUB4V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
2042                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
2043                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
2044                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
2045                 IMS,IME,JMS,JME,KMS,KME,                 &
2046                 ITS,ITE,JTS,JTE,KTS,KTE      )
2047 IMPLICIT NONE
2048 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
2049 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
2050 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
2051 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
2052 INTEGER,    INTENT(IN   )                            :: RATIO
2053 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
2054 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
2055!LOCAL VARIABLES
2056 INTEGER                                              :: I,J
2057 INTEGER                                              :: IP
2058  IP = MOD(I,RATIO)
2059  SELECT CASE(IP)
2060   CASE ( 1 )
2061    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2062    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2063    VBWGT1(I,J) = 1.0
2064    VBWGT2(I,J) = 0.0
2065    VBWGT3(I,J) = 0.0
2066    VBWGT4(I,J) = 0.0
2067   CASE ( 2 )
2068    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2069    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2070    VBWGT1(I,J) = 4.0/9.0
2071    VBWGT2(I,J) = 1.0/9.0
2072    VBWGT3(I,J) = 2.0/9.0
2073    VBWGT4(I,J) = 2.0/9.0
2074   CASE ( 0 )
2075    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2076    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2077    VBWGT1(I,J) = 1.0/9.0
2078    VBWGT2(I,J) = 4.0/9.0
2079    VBWGT3(I,J) = 2.0/9.0
2080    VBWGT4(I,J) = 2.0/9.0
2081   END SELECT
2082 RETURN
2083 END SUBROUTINE SUB4V
2084 SUBROUTINE SUB5V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
2085                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
2086                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
2087                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
2088                 IMS,IME,JMS,JME,KMS,KME,                 &
2089                 ITS,ITE,JTS,JTE,KTS,KTE      )
2090 IMPLICIT NONE
2091 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
2092 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
2093 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
2094 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
2095 INTEGER,    INTENT(IN   )                            :: RATIO
2096 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
2097 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
2098!LOCAL VARIABLES
2099 INTEGER                                              :: I,J
2100 INTEGER                                              :: IP
2101  IP = MOD(I,RATIO)
2102  SELECT CASE(IP)
2103   CASE ( 1 )
2104    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2105    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2106    VBWGT1(I,J) = 2.0/3.0
2107    VBWGT2(I,J) = 0.0
2108    VBWGT3(I,J) = 0.0
2109    VBWGT4(I,J) = 1.0/3.0
2110   CASE ( 2 )
2111    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2112    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2113    VBWGT1(I,J) = 2.0/9.0
2114    VBWGT2(I,J) = 2.0/9.0
2115    VBWGT3(I,J) = 1.0/9.0
2116    VBWGT4(I,J) = 4.0/9.0
2117   CASE ( 0 )
2118    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2119    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
2120    VBWGT1(I,J) = 1.0/3.0
2121    VBWGT2(I,J) = 0.0
2122    VBWGT3(I,J) = 2.0/3.0
2123    VBWGT4(I,J) = 0.0
2124   END SELECT
2125 RETURN
2126 END SUBROUTINE SUB5V
2127 SUBROUTINE SUB6V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
2128                 I_PARENT_START,J_PARENT_START,           & ! nest start I and J in parent domain 
2129                 RATIO,                                   & ! Ratio of parent and child grid ( always = 3 for NMM)
2130                 IDS,IDE,JDS,JDE,KDS,KDE,                 & ! target (nest) dimensions
2131                 IMS,IME,JMS,JME,KMS,KME,                 &
2132                 ITS,ITE,JTS,JTE,KTS,KTE      )
2133 IMPLICIT NONE
2134 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
2135 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
2136 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
2137 INTEGER,    INTENT(IN   )                            :: I_PARENT_START,J_PARENT_START
2138 INTEGER,    INTENT(IN   )                            :: RATIO
2139 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
2140 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIV,JJV
2141!LOCAL VARIABLES
2142 INTEGER                                              :: I,J
2143 INTEGER                                              :: IP
2144  IP = MOD(I,RATIO)
2145  SELECT CASE(IP)
2146   CASE ( 1 )
2147    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO) - 1
2148    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
2149    VBWGT1(I,J) = 2.0/9.0
2150    VBWGT2(I,J) = 2.0/9.0
2151    VBWGT3(I,J) = 4.0/9.0
2152    VBWGT4(I,J) = 1.0/9.0
2153   CASE ( 2 )
2154    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2155    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO)
2156    VBWGT1(I,J) = 1.0/3.0
2157    VBWGT2(I,J) = 0.0
2158    VBWGT3(I,J) = 0.0
2159    VBWGT4(I,J) = 2.0/3.0
2160   CASE ( 0 )
2161    IIV(I,J)    = I_PARENT_START + INT((I-1)/RATIO)
2162    JJV(I,J)    = J_PARENT_START + INT((J-1)/RATIO) + 1
2163    VBWGT1(I,J) = 2.0/3.0
2164    VBWGT2(I,J) = 0.0
2165    VBWGT3(I,J) = 1.0/3.0
2166    VBWGT4(I,J) = 0.0
2167   END SELECT
2168 RETURN
2169 END SUBROUTINE SUB6V
2170
2171
2172!------------------------------------------------------------------------------
2173!
2174SUBROUTINE WEIGTS_CHECK(HBWGT1,HBWGT2,HBWGT3,HBWGT4,       &
2175                        VBWGT1,VBWGT2,VBWGT3,VBWGT4,       &
2176                        IDS,IDE,JDS,JDE,KDS,KDE,           &
2177                        IMS,IME,JMS,JME,KMS,KME,           &
2178                        ITS,ITE,JTS,JTE,KTS,KTE            )
2179
2180  IMPLICIT NONE
2181  INTEGER, INTENT(IN)                                 :: IDS,IDE,JDS,JDE,KDS,KDE,  &
2182                                                         IMS,IME,JMS,JME,KMS,KME,  &
2183                                                         ITS,ITE,JTS,JTE,KTS,KTE
2184  REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(IN)   :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
2185  REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(IN)   :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
2186
2187! local
2188
2189  REAL , PARAMETER :: EPSI=1.0E-3
2190  INTEGER          :: I,J
2191  REAL             :: ADDSUM
2192 CHARACTER(LEN=255):: message
2193
2194!-------------------------------------------------------------------------------------
2195
2196! DUE TO THE NEED FOR HALO EXCHANGES IN PARALLEL RUNS ONE HAS TO ENSURE CONSISTENT
2197! USAGE OF NUMBER OF PROCESSORS BEFORE ANY FURTHER COMPUTATIONS. WE INTRODUCE THIS
2198! CHECK FIRST
2199
2200  IF((ITE-ITS) .LE. 5 .OR. (JTE-JTS) .LE. 5)THEN
2201   WRITE(message,*)'ITE-ITS=',ITE-ITS,'JTE-JTS=',JTE-JTS
2202   CALL wrf_message(trim(message))
2203   CALL wrf_error_fatal ('NESTED DOMAIN:PLEASE OPTIMIZE THE NUMBER OF PROCESSES; TRY SQUARES OF NUMBERS')
2204  ENDIF
2205
2206
2207! NOW CHECK WEIGHTS
2208!
2209
2210  ADDSUM=0.
2211  DO J = JTS, MIN(JTE,JDE-1)
2212   DO I = ITS, MIN(ITE,IDE-1)
2213      ADDSUM=HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J)
2214      IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
2215       WRITE(message,*)'I=',I,'J=',J,'WEIGHTS=',HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J),1-ADDSUM
2216       CALL wrf_message(trim(message))
2217       CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT MASS POINTS')
2218      ENDIF
2219   ENDDO
2220  ENDDO
2221
2222  ADDSUM=0.
2223  DO J = JTS, MIN(JTE,JDE-1)
2224   DO I = ITS, MIN(ITE,IDE-1)
2225      ADDSUM=VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J)
2226      IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
2227       WRITE(message,*)'I=',I,'J=',J,'WEIGHTS=',VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J),1-ADDSUM
2228       CALL wrf_message(trim(message))
2229       CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT VELOCITY POINTS')
2230      ENDIF
2231   ENDDO
2232  ENDDO
2233
2234END SUBROUTINE WEIGTS_CHECK
2235
2236!-----------------------------------------------------------------------------------
2237
2238SUBROUTINE BOUNDS_CHECK( IIH,JJH,IIV,JJV,          &
2239                         IPOS,JPOS,SHW,            &
2240                         IDS,IDE,JDS,JDE,KDS,KDE,  & !
2241                         IMS,IME,JMS,JME,KMS,KME,  & ! nested grid configuration
2242                         ITS,ITE,JTS,JTE,KTS,KTE   )
2243
2244 IMPLICIT NONE
2245 INTEGER, INTENT(IN) :: IPOS,JPOS,SHW,            &
2246                        IDS,IDE,JDS,JDE,KDS,KDE,  &
2247                        IMS,IME,JMS,JME,KMS,KME,  &
2248                        ITS,ITE,JTS,JTE,KTS,KTE   
2249
2250 INTEGER, DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IIH,JJH,IIV,JJV
2251
2252! local variables
2253
2254 INTEGER :: I,J
2255 CHARACTER(LEN=255)                 :: message
2256
2257!***  Gopal       - Initial version
2258!***
2259!*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
2260!
2261!============================================================================
2262
2263  IF(IPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs X-BOUNDARY')
2264  IF(JPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs Y-BOUNDARY')
2265
2266  DO J = JTS, MIN(JTE,JDE-1)
2267   DO I = ITS, MIN(ITE,IDE-1)
2268      IF(IIH(I,J) .EQ. 0)CALL wrf_error_fatal ('IIH=0: SOMETHING IS WRONG')
2269      IF(JJH(I,J) .EQ. 0)CALL wrf_error_fatal ('JJH=0: SOMETHING IS WRONG')
2270   ENDDO
2271  ENDDO
2272
2273  DO J = JTS, MIN(JTE,JDE-1)
2274   DO I = ITS, MIN(ITE,IDE-1)
2275      IF(IIH(I,J) .LT. (IPOS-SHW) .OR. JJH(I,J) .LT. (JPOS-SHW) .OR.   &
2276         IIV(I,J) .LT. (IPOS-SHW) .OR. JJV(I,J) .LT. (JPOS-SHW))THEN
2277         WRITE(message,*)I,J,IIH(I,J),IPOS,JJH(I,J),JPOS,SHW
2278         CALL wrf_message(trim(message))
2279         WRITE(message,*)I,J,IIV(I,J),IPOS,JJV(I,J),JPOS,SHW
2280         CALL wrf_message(trim(message))
2281         CALL wrf_error_fatal ('CHECK NESTED DOMAIN BOUNDS: TRY INCREASING STENCIL WIDTH')
2282      ENDIF
2283   ENDDO
2284  ENDDO
2285
2286END SUBROUTINE BOUNDS_CHECK
2287
2288!==========================================================================================
2289
2290
2291SUBROUTINE BASE_STATE_PARENT ( Z3d,Q3d,T3d,PSTD,        &
2292                               PINT,T,Q,CWM,            &
2293                               FIS,QS,PD,PDTOP,PTOP,    &
2294                               ETA1,ETA2,               &
2295                               DETA1,DETA2,             &
2296                               IDS,IDE,JDS,JDE,KDS,KDE, &
2297                               IMS,IME,JMS,JME,KMS,KME, &
2298                               ITS,ITE,JTS,JTE,KTS,KTE  )
2299!
2300
2301 USE MODULE_MODEL_CONSTANTS
2302 IMPLICIT NONE
2303 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
2304 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
2305 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
2306 REAL,       INTENT(IN   )                            :: PDTOP,PTOP
2307 REAL, DIMENSION(KMS:KME),                 INTENT(IN) :: ETA1,ETA2,DETA1,DETA2
2308 REAL, DIMENSION(IMS:IME,JMS:JME),         INTENT(IN) :: FIS,PD,QS
2309 REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME), INTENT(IN) :: PINT,T,Q,CWM
2310 REAL, DIMENSION(KMS:KME),                 INTENT(OUT):: PSTD
2311 REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME), INTENT(OUT):: Z3d,Q3d,T3d
2312
2313! local
2314
2315 INTEGER,PARAMETER                                    :: JTB=134
2316 INTEGER                                              :: I,J,K,ILOC,JLOC
2317 REAL, PARAMETER                                      :: LAPSR=6.5E-3, GI=1./G,D608=0.608
2318 REAL, PARAMETER                                      :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
2319 REAL, PARAMETER                                      :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
2320 REAL, PARAMETER                                      :: P_REF=103000.
2321 REAL                                                 :: A,B,APELP,RTOPP,DZ,ZMID
2322 REAL, DIMENSION(IMS:IME,JMS:JME)                     :: SLP,TSFC,ZMSLP
2323 REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME)             :: Z3d_IN
2324 REAL,DIMENSION(JTB)                                  :: PIN,ZIN,Y2,PIO,ZOUT,DUM1,DUM2
2325 REAL,DIMENSION(JTB)                                  :: QIN,QOUT,TIN,TOUT
2326!--------------------------------------------------------------------------------------
2327
2328!  CLEAN Z3D ARRAY FIRST
2329
2330    DO K=KDS,KDE
2331     DO J = JTS, MIN(JTE,JDE-1)
2332      DO I = ITS, MIN(ITE,IDE-1)
2333       Z3d(I,J,K)=0.0
2334       T3d(I,J,K)=0.0
2335       Q3d(I,J,K)=0.0
2336      ENDDO
2337     ENDDO
2338    ENDDO
2339
2340
2341!  DETERMINE THE HEIGHTS ON THE PARENT DOMAIN
2342
2343    DO J = JTS, MIN(JTE,JDE-1)
2344      DO I = ITS, MIN(ITE,IDE-1)
2345       Z3d_IN(I,J,1)=FIS(I,J)*GI
2346      ENDDO
2347    ENDDO
2348
2349    DO K = KDS,KDE-1
2350     DO J = JTS, MIN(JTE,JDE-1)
2351      DO I = ITS, MIN(ITE,IDE-1)
2352        APELP    = (PINT(I,J,K+1)+PINT(I,J,K))
2353!       RTOPP    = TRG*T(I,J,K)*(1.0+Q(I,J,K)*P608-CWM(I,J,K))/APELP
2354        RTOPP    = TRG*T(I,J,K)*(1.0+Q(I,J,K)*P608)/APELP
2355        DZ       = RTOPP*(DETA1(K)*PDTOP+DETA2(K)*PD(I,J))   ! (RTv/P_TOT)*D(P_HYDRO)
2356        Z3d_IN(I,J,K+1) = Z3d_IN(I,J,K) + DZ
2357      ENDDO
2358     ENDDO
2359    ENDDO
2360
2361
2362!  CONSTRUCT STANDARD ISOBARIC SURFACES
2363
2364    DO K=KDS,KDE                         ! target points in model interface levels (pint)
2365       PSTD(K) = ETA1(K)*PDTOP + ETA2(K)*(P_REF -PDTOP - PTOP) + PTOP
2366    ENDDO
2367
2368!   DETERMINE THE MSLP USE THAT TO CREATE HEIGHTS AT 1000. mb LEVEL. THESE HEIGHTS 
2369!   MAY ONLY BE USED IN VERTICAL INTERPOLATION TO ISOBARIC SURFACES WHICH ARE LOCATED
2370!   BELOW GROUND LEVEL.
2371
2372    DO J = JTS, MIN(JTE,JDE-1)
2373      DO I = ITS, MIN(ITE,IDE-1)
2374        TSFC(I,J) = T(I,J,1)*(1.+D608*Q(I,J,1)) + LAPSR*(Z3d_IN(I,J,1)+Z3d_IN(I,J,2))*0.5
2375        A         = LAPSR*Z3d_IN(I,J,1)/TSFC(I,J)
2376        SLP(I,J)  = PINT(I,J,1)*(1-A)**COEF2    ! sea level pressure
2377        B         = (PSTD(1)/SLP(I,J))**COEF3
2378        ZMSLP(I,J)= TSFC(I,J)*LAPSI*(1.0 - B)   ! Height at 1000. mb level
2379      ENDDO
2380    ENDDO
2381
2382!   INTERPOLATE Z3d_IN TO STANDARD PRESSURE INTERFACES. FOR LEVELS BELOW
2383!   GROUND USE ZMSLP(I,J)
2384
2385    DO J = JTS, MIN(JTE,JDE-1)
2386      DO I = ITS, MIN(ITE,IDE-1)
2387!   
2388!     clean local array before use of spline
2389
2390      PIN=0.;ZIN=0.;Y2=0;PIO=0.;ZOUT=0.;DUM1=0.;DUM2=0.
2391
2392       DO K=KDS,KDE                           ! inputs at model interfaces
2393         PIN(K) = PINT(I,J,KDE-K+1)
2394         ZIN(K) = Z3d_IN(I,J,KDE-K+1)
2395       ENDDO
2396
2397       IF(PINT(I,J,1) .LE. PSTD(1))THEN
2398          PIN(KDE) = PSTD(1)
2399          ZIN(KDE) = ZMSLP(I,J)
2400       ENDIF
2401!
2402       Y2(1  )=0.
2403       Y2(KDE)=0.
2404!
2405       DO K=KDS,KDE
2406          PIO(K)=PSTD(K)
2407       ENDDO
2408!
2409       CALL SPLINE1(I,J,JTB,KDE,PIN,ZIN,Y2,KDE,PIO,ZOUT,DUM1,DUM2)  ! interpolate
2410!
2411
2412       DO K=KDS,KDE                           ! inputs at model interfaces
2413         Z3d(I,J,K)=ZOUT(K)
2414       ENDDO
2415
2416      ENDDO
2417    ENDDO
2418!
2419!   INTERPOLATE TEMPERATURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW 
2420!   GROUND USE A LAPSE RATE ATMOSPHERE
2421!
2422    DO J = JTS, MIN(JTE,JDE-1)
2423      DO I = ITS, MIN(ITE,IDE-1)
2424
2425!     clean local array before use of spline or linear interpolation
2426!
2427      PIN=0.;TIN=0.;Y2=0;PIO=0.;TOUT=0.;DUM1=0.;DUM2=0.
2428
2429       DO K=KDS+1,KDE                           ! inputs at model levels
2430         PIN(K-1) = EXP((ALOG(PINT(I,J,KDE-K+1))+ALOG(PINT(I,J,KDE-K+2)))*0.5)
2431         TIN(K-1) = T(I,J,KDE-K+1)
2432       ENDDO
2433
2434       IF(PINT(I,J,1) .LE. PSTD(1))THEN
2435         PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
2436         ZMID     = 0.5*(Z3d_IN(I,J,1)+Z3d_IN(I,J,2))
2437         TIN(KDE-1) = T(I,J,1) + LAPSR*(ZMID-ZMSLP(I,J))
2438       ENDIF
2439!
2440       Y2(1    )=0.
2441       Y2(KDE-1)=0.
2442!
2443       DO K=KDS,KDE-1
2444          PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
2445       ENDDO
2446
2447       CALL SPLINE1(I,J,JTB,KDE-1,PIN,TIN,Y2,KDE-1,PIO,TOUT,DUM1,DUM2)  ! interpolate
2448
2449
2450       DO K=KDS,KDE-1                           ! inputs at model levels
2451         T3d(I,J,K)=TOUT(K)
2452       ENDDO
2453
2454      ENDDO
2455    ENDDO
2456
2457!
2458!   INTERPOLATE MOISTURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW
2459!   GROUND USE THE SURFACE MOISTURE
2460!
2461    DO J = JTS, MIN(JTE,JDE-1)
2462      DO I = ITS, MIN(ITE,IDE-1)
2463!
2464!     clean local array before use of spline or linear interpolation
2465
2466
2467      PIN=0.;QIN=0.;Y2=0;PIO=0.;QOUT=0.;DUM1=0.;DUM2=0.
2468
2469       DO K=KDS+1,KDE                           ! inputs at model levels
2470         PIN(K-1) = EXP((ALOG(PINT(I,J,KDE-K+1))+ALOG(PINT(I,J,KDE-K+2)))*0.5)
2471         QIN(K-1) = Q(I,J,KDE-K+1)
2472       ENDDO
2473
2474       IF(PINT(I,J,1) .LE. PSTD(1))THEN
2475          PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
2476!         QIN(KDE-1) =  QS(I,J)
2477       ENDIF
2478
2479       Y2(1    )=0.
2480       Y2(KDE-1)=0.
2481!
2482       DO K=KDS,KDE-1
2483          PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
2484       ENDDO
2485
2486       CALL SPLINE1(I,J,JTB,KDE-1,PIN,QIN,Y2,KDE-1,PIO,QOUT,DUM1,DUM2)  ! interpolate
2487
2488       DO K=KDS,KDE-1                          ! inputs at model levels
2489         Q3d(I,J,K)=QOUT(K)
2490       ENDDO
2491
2492      ENDDO
2493    ENDDO
2494
2495END SUBROUTINE BASE_STATE_PARENT
2496!=============================================================================
2497      SUBROUTINE SPLINE1(I,J,JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
2498!
2499!   ******************************************************************
2500!   *                                                                *
2501!   *  THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE        *
2502!   *  PROGRAMED FOR A SMALL SCALAR MACHINE.                         *
2503!   *                                                                *
2504!   *  PROGRAMER Z. JANJIC                                           *
2505!   *                                                                *
2506!   *  NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION.  MUST BE GE 3. *
2507!   *  XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
2508!   *         FUNCTION ARE GIVEN.  MUST BE IN ASCENDING ORDER.       *
2509!   *  YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD.   *
2510!   *  Y2   - THE SECOND DERIVATIVES AT THE POINTS XOLD.  IF NATURAL *
2511!   *         SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE      *
2512!   *         SPECIFIED.                                             *
2513!   *  NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED.     *
2514!   *  XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
2515!   *         FUNCTION ARE CALCULATED.  XNEW(K) MUST BE GE XOLD(1)   *
2516!   *         AND LE XOLD(NOLD).                                     *
2517!   *  YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED.           *
2518!   *  P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2.                *
2519!   *                                                                *
2520!   ******************************************************************
2521!---------------------------------------------------------------------
2522      IMPLICIT NONE
2523!---------------------------------------------------------------------
2524      INTEGER,INTENT(IN) :: I,J,JTBX,NNEW,NOLD
2525      REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
2526      REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
2527      REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
2528!
2529      INTEGER :: II,JJ,K,K1,K2,KOLD,NOLDM1
2530      REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR                 &
2531             ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
2532      CHARACTER(LEN=255) :: message
2533!---------------------------------------------------------------------
2534
2535!     debug
2536
2537      II=9999 !67 !35 !50   !4
2538      JJ=9999 !31 !73 !115  !192
2539      IF(I.eq.II.and.J.eq.JJ)THEN
2540        WRITE(message,*)'DEBUG in SPLINE1:HSO= ',xnew(1:nold)
2541        CALL wrf_debug(1,trim(message))
2542        DO K=1,NOLD
2543         WRITE(message,*)'DEBUG in SPLINE1:L,ZETAI,PINTI= ' &
2544                        ,K,YOLD(K),XOLD(K)
2545         CALL wrf_debug(1,trim(message))
2546        ENDDO
2547      ENDIF
2548
2549!
2550      NOLDM1=NOLD-1
2551!
2552      DXL=XOLD(2)-XOLD(1)
2553      DXR=XOLD(3)-XOLD(2)
2554      DYDXL=(YOLD(2)-YOLD(1))/DXL
2555      DYDXR=(YOLD(3)-YOLD(2))/DXR
2556      RTDXC=0.5/(DXL+DXR)
2557!
2558      P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
2559      Q(1)=-RTDXC*DXR
2560!
2561      IF(NOLD.EQ.3)GO TO 150
2562!---------------------------------------------------------------------
2563      K=3
2564!
2565  100 DXL=DXR
2566      DYDXL=DYDXR
2567      DXR=XOLD(K+1)-XOLD(K)
2568      DYDXR=(YOLD(K+1)-YOLD(K))/DXR
2569      DXC=DXL+DXR
2570      DEN=1./(DXL*Q(K-2)+DXC+DXC)
2571!
2572      P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
2573      Q(K-1)=-DEN*DXR
2574!
2575      K=K+1
2576      IF(K.LT.NOLD)GO TO 100
2577!-----------------------------------------------------------------------
2578  150 K=NOLDM1
2579!
2580  200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
2581!
2582      K=K-1
2583      IF(K.GT.1)GO TO 200
2584!-----------------------------------------------------------------------
2585      K1=1
2586!
2587  300 XK=XNEW(K1)
2588!
2589      DO 400 K2=2,NOLD
2590!
2591      IF(XOLD(K2).GT.XK)THEN
2592        KOLD=K2-1
2593        GO TO 450
2594      ENDIF
2595!
2596  400 CONTINUE
2597!
2598      YNEW(K1)=YOLD(NOLD)
2599      GO TO 600
2600!
2601  450 IF(K1.EQ.1)GO TO 500
2602      IF(K.EQ.KOLD)GO TO 550
2603!
2604  500 K=KOLD
2605!
2606      Y2K=Y2(K)
2607      Y2KP1=Y2(K+1)
2608      DX=XOLD(K+1)-XOLD(K)
2609      RDX=1./DX
2610!
2611      AK=.1666667*RDX*(Y2KP1-Y2K)
2612      BK=0.5*Y2K
2613      CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
2614!
2615  550 X=XK-XOLD(K)
2616      XSQ=X*X
2617!
2618      YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)
2619
2620!  debug
2621
2622      if(i.eq.ii.and.j.eq.jj)then
2623        write(message,*) 'DEBUG:: k1,xnew(k1),ynew(k1): ', k1,xnew(k1),ynew(k1)
2624        CALL wrf_debug(1,trim(message))
2625      endif
2626
2627!
2628  600 K1=K1+1
2629      IF(K1.LE.NNEW)GO TO 300
2630
2631      RETURN
2632      END SUBROUTINE SPLINE1
2633!---------------------------------------------------------------------
2634
2635SUBROUTINE NEST_TERRAIN ( nest, config_flags )
2636
2637 USE module_domain
2638 USE module_configure
2639 USE module_timing
2640
2641 USE wrfsi_static
2642
2643 IMPLICIT NONE
2644
2645 TYPE(domain) , POINTER                        :: nest
2646 TYPE(grid_config_rec_type) , INTENT(IN)       :: config_flags
2647
2648!
2649! Local variables
2650!
2651
2652 LOGICAL, EXTERNAL                 :: wrf_dm_on_monitor
2653 INTEGER                           :: ids,ide,jds,jde,kds,kde
2654 INTEGER                           :: ims,ime,jms,jme,kms,kme
2655 INTEGER                           :: its,ite,jts,jte,kts,kte
2656 INTEGER                           :: i_parent_start, j_parent_start
2657 INTEGER                           :: parent_grid_ratio
2658 INTEGER                           :: n,i,j,ii,jj,nnxp,nnyp
2659 INTEGER                           :: i_start,j_start,level
2660 REAL, ALLOCATABLE, DIMENSION(:,:) :: data1     ! for highres topo
2661 REAL, ALLOCATABLE, DIMENSION(:,:) :: avc_big, lnd_big, lah_big, loh_big
2662 REAL, ALLOCATABLE, DIMENSION(:,:) :: avc_nest, lnd_nest, lah_nest, loh_nest
2663 INTEGER                           :: im_big, jm_big, i_add
2664 INTEGER                           :: im, jm
2665 CHARACTER(LEN=6)                  :: nestpath
2666
2667 integer                           :: input_type
2668 character(len=128)                :: input_fname
2669 character (len=32)                :: cname
2670 integer                           :: ndim
2671 character (len=3)                 :: memorder
2672 character (len=32)                :: stagger
2673 integer, dimension(3)             :: domain_start, domain_end
2674 integer                           :: wrftype
2675 character (len=128), dimension(3) :: dimnames
2676
2677 integer :: istatus
2678 integer :: handle
2679 integer :: comm_1, comm_2
2680
2681 real, allocatable, dimension(:,:,:) :: real_domain
2682
2683 character (len=10), dimension(24)  :: name = (/ "XLAT_M    ", &
2684                                                "XLONG_M   ", &
2685                                                "XLAT_V    ", &
2686                                                "XLONG_V   ", &
2687                                                "E         ", &
2688                                                "F         ", &
2689                                                "LANDMASK  ", &
2690                                                "LANDUSEF  ", &
2691                                                "LU_INDEX  ", &
2692                                                "HCNVX     ", &
2693                                                "HSTDV     ", &
2694                                                "HASYW     ", &
2695                                                "HASYS     ", &
2696                                                "HASYSW    ", &
2697                                                "HASYNW    ", &
2698                                                "HLENW     ", &
2699                                                "HLENS     ", &
2700                                                "HLENSW    ", &
2701                                                "HLENNW    ", &
2702                                                "HANIS     ", &
2703                                                "HSLOP     ", &
2704                                                "HANGL     ", &
2705                                                "HZMAX     ", &
2706                                                "HGT_M     " /)
2707
2708 integer, parameter :: IO_BIN=1, IO_NET=2
2709
2710 integer :: io_form_input
2711
2712 CHARACTER(LEN=512) :: message
2713
2714 write(message,*)"in NEST_TERRAIN config_flags%io_form_input = ", config_flags%io_form_input
2715 CALL wrf_message(trim(message))
2716 write(message,*)"in NEST_TERRAIN config_flags%auxinput1_inname = ", config_flags%auxinput1_inname
2717 CALL wrf_message(trim(message))
2718
2719 io_form_input = config_flags%io_form_input
2720 if (config_flags%auxinput1_inname(1:7) == "met_nmm") then
2721    input_type = 2
2722 else
2723    input_type = 1
2724 end if
2725
2726!----------------------------------------------------------------------------------
2727
2728      IDS = nest%sd31
2729      IDE = nest%ed31
2730      JDS = nest%sd32
2731      JDE = nest%ed32
2732      KDS = nest%sd33
2733      KDE = nest%ed33
2734
2735      IMS = nest%sm31
2736      IME = nest%em31
2737      JMS = nest%sm32
2738      JME = nest%em32
2739      KMS = nest%sm33
2740      KME = nest%em33
2741
2742      ITS = nest%sp31
2743      ITE = nest%ep31
2744      JTS = nest%sp32
2745      JTE = nest%ep32
2746      KTS = nest%sp33
2747      KTE = nest%ep33
2748
2749      i_parent_start = nest%i_parent_start
2750      j_parent_start = nest%j_parent_start
2751      parent_grid_ratio = nest%parent_grid_ratio
2752
2753      NNXP=IDE-1
2754      NNYP=JDE-1
2755
2756      ALLOCATE(DATA1(1:NNXP,1:NNYP))
2757!
2758!
2759!--- Read in high resolution topography
2760!
2761      IF ( wrf_dm_on_monitor() ) THEN                    ! first assign a status
2762!
2763!      This part of the code is Dusan's doing. Extended by gopal for multiple nest (Feb 19,2005)
2764!
2765       call find_ijstart_level (nest,i_start,j_start,level)
2766       write(message,*)" nest%id =", nest%id , " i_start,j_start,level =", i_start,j_start,level
2767       CALL wrf_message(trim(message))
2768
2769       write(nestpath,"(a4,i1,a1)") 'nest',level,'/'
2770
2771       if ( level > 0 ) then
2772
2773          if (input_type == 1) then
2774!
2775!        si version of the static file
2776!
2777             CALL get_wrfsi_static_dims(nestpath, im_big, jm_big)
2778             ALLOCATE (avc_big(im_big,jm_big))
2779             ALLOCATE (lnd_big(im_big,jm_big))
2780             ALLOCATE (lah_big(im_big,jm_big))
2781             ALLOCATE (loh_big(im_big,jm_big))
2782             CALL get_wrfsi_static_2d(nestpath, 'avc', avc_big)
2783             CALL get_wrfsi_static_2d(nestpath, 'lnd', lnd_big)
2784             CALL get_wrfsi_static_2d(nestpath, 'lah', lah_big)
2785             CALL get_wrfsi_static_2d(nestpath, 'loh', loh_big)
2786
2787          else if (input_type == 2) then
2788!
2789!        WPS version of the static file
2790!
2791
2792#ifdef INTIO
2793             if (io_form_input == IO_BIN) write(input_fname,"(A,I2.2,A)") "geo_nmm_nest.l",level,".int"
2794#endif
2795#ifdef NETCDF
2796             if (io_form_input == IO_NET) write(input_fname,"(A,I2.2,A)") "geo_nmm_nest.l",level,".nc"
2797#endif
2798
2799             comm_1 = 1
2800             comm_2 = 1
2801
2802#ifdef INTIO
2803             if (io_form_input == IO_BIN) &
2804                call ext_int_open_for_read(trim(input_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
2805#endif
2806#ifdef NETCDF
2807             if (io_form_input == IO_NET) &
2808                call ext_ncd_open_for_read(trim(input_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
2809#endif
2810             if (istatus /= 0) CALL wrf_error_fatal('NEST_TERRAIN error after ext_XXX_open_for_read '//trim(input_fname))
2811
2812
2813             do n=1,24
2814
2815             cname = name(n)
2816
2817             domain_start = 1
2818             domain_end = 1
2819#ifdef INTIO
2820             if (io_form_input == IO_BIN) &
2821                call ext_int_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus)
2822#endif
2823#ifdef NETCDF
2824             if (io_form_input == IO_NET) &
2825                call ext_ncd_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus)
2826#endif
2827             print *, "cname=", cname
2828             print *, "istatus=", istatus
2829             print *, "ndim=", ndim
2830             print *, "memorder=", memorder
2831             print *, "stagger=", stagger
2832             print *, "domain_start=", domain_start
2833             print *, "domain_end=", domain_end
2834             print *, "wrftype=", wrftype
2835
2836
2837             if (allocated(real_domain)) deallocate(real_domain)
2838             allocate(real_domain(domain_start(1):domain_end(1), domain_start(2):domain_end(2), domain_start(3):domain_end(3)))
2839
2840#ifdef INTIO
2841             if (io_form_input == IO_BIN) then
2842                call ext_int_read_field(handle, '0000-00-00_00:00:00', cname, real_domain, wrftype, &
2843                                        1, 1, 0, memorder, stagger, &
2844                                        dimnames, domain_start, domain_end, domain_start, domain_end, &
2845                                        domain_start, domain_end, istatus)
2846             end if
2847#endif
2848#ifdef NETCDF
2849             if (io_form_input == IO_NET) then
2850                call ext_ncd_read_field(handle, '0000-00-00_00:00:00', cname, real_domain, wrftype, &
2851                                        1, 1, 0, memorder, stagger, &
2852                                        dimnames, domain_start, domain_end, domain_start, domain_end, &
2853                                        domain_start, domain_end, istatus)
2854             end if
2855#endif
2856             print *, "istatus=", istatus
2857
2858             im_big = domain_end(1)
2859             jm_big = domain_end(2)
2860             if (cname(1:10) == "XLAT_M    ") then
2861                ALLOCATE (lah_big(im_big,jm_big))
2862                do j=1,jm_big
2863                do i=1,im_big
2864                   lah_big(i,j) = real_domain(i,j,1)
2865                end do
2866                end do
2867             else if (cname(1:10) == "XLONG_M   ") then
2868                ALLOCATE (loh_big(im_big,jm_big))
2869                do j=1,jm_big
2870                do i=1,im_big
2871                   loh_big(i,j) = real_domain(i,j,1)
2872                end do
2873                end do
2874             else if (cname(1:10) == "LANDMASK  ") then
2875                ALLOCATE (lnd_big(im_big,jm_big))
2876                do j=1,jm_big
2877                do i=1,im_big
2878                   lnd_big(i,j) = real_domain(i,j,1)
2879                end do
2880                end do
2881             else if (cname(1:10) == "HGT_M     ") then
2882                ALLOCATE (avc_big(im_big,jm_big))
2883                do j=1,jm_big
2884                do i=1,im_big
2885                   avc_big(i,j) = real_domain(i,j,1)
2886                end do
2887                end do
2888             end if
2889
2890             end do
2891
2892#ifdef INTIO
2893             if (io_form_input == IO_BIN) then
2894                call ext_int_ioclose(handle, istatus)
2895             end if
2896#endif
2897#ifdef NETCDF
2898             if (io_form_input == IO_NET) then
2899                call ext_ncd_ioclose(handle, istatus)
2900             end if
2901#endif
2902
2903          else
2904             CALL wrf_error_fatal('NEST_TERRAIN wrong input_type')
2905          end if
2906
2907       else
2908          CALL wrf_error_fatal('this routine NEST_TERRAIN should not be called for top-level domain')
2909       end if
2910
2911! select subdomain from big fine grid
2912
2913        im = NNXP
2914        jm = NNYP
2915
2916        ALLOCATE (avc_nest(im,jm))
2917        ALLOCATE (lnd_nest(im,jm))
2918        ALLOCATE (lah_nest(im,jm))
2919        ALLOCATE (loh_nest(im,jm))
2920
2921        i_add = mod(j_start+1,2)
2922        DO j=1,jm
2923        DO i=1,im
2924           avc_nest(i,j) = avc_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2925           lnd_nest(i,j) = lnd_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2926           lah_nest(i,j) = lah_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2927           loh_nest(i,j) = loh_big(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
2928        END DO
2929        END DO
2930
2931        WRITE(message,*)'SOME MATCHING TEST i_parent_start, j_parent_start',i_parent_start,j_parent_start
2932       CALL wrf_message(trim(message))
2933       CALL wrf_message('WRFSI LAT      COMPUTED LAT')
2934        WRITE(message,*)lah_nest(1,1),nest%HLAT(1,1)
2935       CALL wrf_message(trim(message))
2936       CALL wrf_message('WRFSI LON      COMPUTED LON')
2937        WRITE(message,*)loh_nest(1,1),nest%HLON(1,1)
2938       CALL wrf_message(trim(message))
2939
2940        IF(ABS(lah_nest(1,1)-nest%HLAT(1,1)) .GE. 0.5 .OR.  &
2941           ABS(loh_nest(1,1)-nest%HLON(1,1)) .GE. 0.5)THEN
2942            CALL wrf_message('CHECK WRFSI CONFIGURATION AND INPUT HIGH RESOLUTION TOPOGRAPHY AND/OR GRID RATIO')
2943            CALL wrf_error_fatal('LATLON MISMATCH: ERROR READING static FILE FOR THE NEST')
2944        ENDIF
2945
2946        call smdhld(im,jm,avc_nest,1.0-lnd_nest,12,12)
2947
2948!-------------4-point averaging of mountains along inner boundary-------
2949
2950        do i=1,im-1
2951            avc_nest(i,2)=0.25*(avc_nest(i,1)+avc_nest(i+1,1)+ &
2952           &                    avc_nest(i,3)+avc_nest(i+1,3))
2953        enddo
2954
2955        do i=1,im-1
2956            avc_nest(i,jm-1)=0.25*(avc_nest(i,jm-2)+avc_nest(i+1,jm-2)+ &
2957           &                       avc_nest(i,jm)+avc_nest(i+1,jm))
2958        enddo
2959
2960        do j=4,jm-3,2
2961            avc_nest(1,j)=0.25*(avc_nest(1,j-1)+avc_nest(2,j-1)+ &
2962           &                    avc_nest(1,j+1)+avc_nest(2,j+1))
2963        enddo
2964
2965        do j=4,jm-3,2
2966            avc_nest(im-1,j)=0.25*(avc_nest(im-1,j-1)+avc_nest(im,j-1)+ &
2967           &                     avc_nest(im-1,j+1)+avc_nest(im,j+1))
2968        enddo
2969
2970        DO J = 1,NNYP
2971         DO I = 1,NNXP
2972            DATA1(I,J) = 9.81*avc_nest(I,J)
2973         ENDDO
2974        ENDDO
2975
2976        DEALLOCATE (avc_big,lnd_big)
2977        DEALLOCATE (avc_nest,lnd_nest)
2978!
2979      ENDIF
2980
2981      CALL wrf_dm_bcast_bytes (DATA1,NNXP*NNYP*RWORDSIZE)
2982
2983      DO J=JDS,JDE
2984       DO I =IDS,IDE
2985        IF(I.GE.ITS .AND. I .LE. MIN(ide-1,ite) .AND. J.GE.JTS  .AND. J .LE. MIN(jde-1,jte))THEN
2986          nest%hres_fis(I,J)=DATA1(I,J)
2987        ENDIF
2988       ENDDO
2989      ENDDO
2990
2991      DEALLOCATE(DATA1)
2992      CALL wrf_message('end of NEST_TERRAIN')
2993
2994END SUBROUTINE NEST_TERRAIN
2995!===========================================================================================
2996
2997
2998SUBROUTINE med_init_domain_constants_nmm ( parent, nest)   !, config_flags)
2999  ! Driver layer
3000   USE module_domain
3001   USE module_configure
3002   USE module_timing
3003   IMPLICIT NONE
3004   TYPE(domain) , POINTER                     :: parent, nest, grid
3005!
3006!
3007   INTERFACE
3008     SUBROUTINE med_initialize_nest_nmm ( grid  &   
3009!
3010# include <dummy_new_args.inc>
3011!
3012                           )
3013        USE module_domain
3014        USE module_configure
3015        USE module_timing
3016        IMPLICIT NONE
3017        TYPE(domain) , POINTER                  :: grid
3018#include <dummy_new_decl.inc>
3019     END SUBROUTINE med_initialize_nest_nmm
3020   END INTERFACE
3021
3022!------------------------------------------------------------------------------
3023!  PURPOSE:
3024!         - initialize some data, mainly 2D & 3D nmm arrays  very similar to
3025!           those done in ./dyn_nmm/module_initialize_real.f
3026!-----------------------------------------------------------------------------
3027!
3028
3029   grid => nest
3030
3031   CALL med_initialize_nest_nmm( grid &   
3032!
3033# include <actual_new_args.inc>
3034!
3035                       )
3036
3037END SUBROUTINE med_init_domain_constants_nmm
3038
3039SUBROUTINE med_initialize_nest_nmm( grid &
3040!
3041# include <dummy_new_args.inc>
3042!
3043                           )
3044
3045 USE module_domain
3046 USE module_configure
3047 USE module_timing
3048 IMPLICIT NONE
3049
3050! Local domain indices and counters.
3051
3052 INTEGER :: ids, ide, jds, jde, kds, kde, &
3053            ims, ime, jms, jme, kms, kme, &
3054            its, ite, jts, jte, kts, kte, &
3055            i, j, k, nnxp, nnyp
3056
3057 TYPE(domain) , POINTER                     :: grid
3058
3059! Local data
3060
3061 LOGICAL, EXTERNAL                   :: wrf_dm_on_monitor
3062 INTEGER                             :: KHH,KVH,JAM,JA,IHL, IHH, L
3063 INTEGER                             :: II,JJ,ISRCH,ISUM
3064 INTEGER, ALLOCATABLE, DIMENSION(:)  :: KHL2,KVL2,KHH2,KVH2,KHLA,KHHA,KVLA,KVHA
3065 INTEGER,PARAMETER                   :: KNUM=SELECTED_REAL_KIND(13)
3066!
3067 REAL(KIND=KNUM)                     :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
3068 REAL(KIND=KNUM)                     :: STPH,CTPH,TDLM,TDPH,FP,TPH,TLM,TLM0
3069 REAL                                :: TPH0D,TLM0D,ANBI,TSPH,DTAD,DTCF,DT
3070 REAL                                :: ACDT,CDDAMP,DXP
3071 REAL                                :: WBD,SBD,WBI,SBI,EBI
3072 REAL                                :: DY_NMM0
3073 REAL                                :: RSNOW,SNOFAC
3074 REAL, ALLOCATABLE, DIMENSION(:)     :: DXJ,WPDARJ,CPGFUJ,CURVJ,   &
3075                                        FCPJ,FDIVJ,EMJ,EMTJ,FADJ,  &
3076                                        HDACJ,DDMPUJ,DDMPVJ
3077!
3078 REAL, PARAMETER:: SALP=2.60
3079 REAL, PARAMETER:: SNUP=0.040
3080 REAL, PARAMETER:: W_NMM=0.08
3081 REAL, PARAMETER:: COAC=0.75
3082 REAL, PARAMETER:: CODAMP=6.4
3083 REAL, PARAMETER:: TWOM=.00014584
3084 REAL, PARAMETER:: CP=1004.6
3085 REAL, PARAMETER:: DFC=1.0   
3086 REAL, PARAMETER:: DDFC=1.0 
3087 REAL, PARAMETER:: ROI=916.6
3088 REAL, PARAMETER:: R=287.04
3089 REAL, PARAMETER:: CI=2060.0
3090 REAL, PARAMETER:: ROS=1500.
3091 REAL, PARAMETER:: CS=1339.2
3092 REAL, PARAMETER:: DS=0.050
3093 REAL, PARAMETER:: AKS=.0000005
3094 REAL, PARAMETER:: DZG=2.85
3095 REAL, PARAMETER:: DI=.1000
3096 REAL, PARAMETER:: AKI=0.000001075
3097 REAL, PARAMETER:: DZI=2.0
3098 REAL, PARAMETER:: THL=210.
3099 REAL, PARAMETER:: PLQ=70000.
3100 REAL, PARAMETER:: ERAD=6371200.
3101 REAL, PARAMETER:: DTR=0.01745329
3102
3103 CHARACTER(LEN=255) :: message
3104
3105   !  Definitions of dummy arguments to solve
3106#include <dummy_new_decl.inc>
3107
3108!#define COPY_IN
3109!#include <scalar_derefs.inc>
3110#ifdef DM_PARALLEL
3111#      include <data_calls.inc>
3112#endif
3113
3114   CALL get_ijk_from_grid (  grid ,                           &
3115                             ids, ide, jds, jde, kds, kde,    &
3116                             ims, ime, jms, jme, kms, kme,    &
3117                             its, ite, jts, jte, kts, kte     )
3118
3119
3120!=================================================================================
3121!
3122!
3123
3124    DT=grid%dt         !float(TIME_STEP)/parent_time_step_ratio
3125    NNXP=min(ITE,IDE-1)
3126    NNYP=min(JTE,JDE-1)
3127    JAM=6+2*((JDE-1)-10)         ! this should be the fix instead of JAM=6+2*(NNYP-10)
3128
3129    WRITE(message,*)'TIME STEP ON DOMAIN',grid%id,'==',dt
3130    CALL wrf_message(trim(message))
3131
3132    WRITE(message,*)'IDS,IDE ON DOMAIN',grid%id,'==',ids,ide
3133    CALL wrf_message(trim(message))
3134!
3135    ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
3136    ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
3137    ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),FADJ(JTS:NNYP))
3138    ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
3139    ALLOCATE(KHLA(JAM),KHHA(JAM))
3140    ALLOCATE(KVLA(JAM),KVHA(JAM))
3141
3142!   INITIALIZE SOME LAND/WATER SURFACE DATA ON THE BASIS OF INPUTS: SM, XICE, WEASD,
3143!   INTERPOLATED FROM MOTHER (WRFSI) DOMAIN. THIS PART OF THE CODE HAS TO BE REVISITED
3144!   LATER ON
3145
3146!   Since SM has been changed on parent domain to be 0 over sea ice it can not be used here
3147!   to find where sea ice is. That's why alogirthm here is slightly different than the
3148!   one used in module_initalize_real.f
3149
3150#ifdef HWRF
3151!zhang's doing: added to AVOID THIS COMPUTATION IF THE NEST IS STARTED USING ANALYSIS FILE
3152   IF(.not. grid%analysis)THEN
3153#endif
3154    DO J = JTS, MIN(JTE,JDE-1)
3155     DO I = ITS, MIN(ITE,IDE-1)
3156
3157      IF (grid%sm(I,J).GT.0.9) THEN                              ! OVER WATER SURFACE
3158         grid%epsr(I,J)= 0.97
3159         grid%embck(I,J)= 0.97
3160         grid%gffc(I,J)= 0.
3161         grid%albedo(I,J)=.06
3162         grid%albase(I,J)=.06
3163      ENDIF
3164
3165      IF (grid%sice(I,J).GT.0.9) THEN                            ! OVER SEA-ICE
3166         grid%sm(I,J)=0.
3167         grid%si(I,J)=0.
3168         grid%gffc(I,J)=0.
3169         grid%albedo(I,J)=.60
3170         grid%albase(I,J)=.60
3171      ENDIF
3172
3173      IF (grid%sm(I,J).LT.0.5.AND.grid%sice(I,J).LT.0.5) THEN         ! OVER LAND SURFACE
3174           grid%si(I,J)=5.0*grid%weasd(I,J)/1000. ! SNOW WATER EQ (mm) OBTAINED FROM PARENT (grid%si) IS INTERPOLATED
3175           grid%epsr(I,J)=1.0                ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
3176           grid%embck(I,J)=1.0               ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
3177           grid%gffc(I,J)=0.0                ! just leave zero as irrelevant
3178           grid%sno(I,J)=grid%si(I,J)*.20         ! LAND-SNOW COVER
3179      ENDIF
3180
3181     ENDDO
3182    ENDDO
3183
3184#if 0
3185    DO J = JTS, MIN(JTE,JDE-1)
3186     DO I = ITS, MIN(ITE,IDE-1)
3187      IF(grid%sm(I,J).GT.0.9) THEN           ! OVER WATER SURFACE
3188!
3189           IF (XICE(I,J) .gt. 0)THEN    ! XICE: SI INPUT ON PARENT, INTERPOLATED ONTO NEST
3190             grid%si(I,J)=1.0                ! INITIALIZE SI BASED ON XICE FROM INTERPOLATED INPUT
3191           ENDIF
3192!
3193           grid%epsr(I,J)= 0.97              ! VALID OVER SEA SURFACE
3194           grid%embck(I,J)= 0.97              ! VALID OVER SEA SURFACE
3195           grid%gffc(I,J)= 0.
3196           grid%albedo(I,J)=.06
3197           grid%albase(I,J)=.06
3198!
3199              IF(grid%si (I,J) .GT. 0.)THEN  ! VALID OVER SEA-ICE
3200                 grid%sm(I,J)=0.
3201                 grid%si(I,J)=0.             !
3202                 grid%sice(I,J)=1.
3203                 grid%gffc(I,J)=0.           ! just leave zero as irrelevant
3204                 grid%albedo(I,J)=.60        ! DEFINE grid%albedo
3205                 grid%albase(I,J)=.60
3206              ENDIF
3207!
3208      ELSE                              ! OVER LAND SURFACE
3209!
3210           grid%si(I,J)=5.0*grid%weasd(I,J)/1000. ! SNOW WATER EQ (mm) OBTAINED FROM PARENT (grid%si) IS INTERPOLATED
3211           grid%epsr(I,J)=1.0                ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
3212           grid%embck(I,J)=1.0                ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
3213           grid%gffc(I,J)=0.0                ! just leave zero as irrelevant
3214           grid%sice(I,J)=0.                 ! SEA ICE
3215           grid%sno(I,J)=grid%si(I,J)*.20         ! LAND-SNOW COVER
3216!
3217      ENDIF
3218!
3219     ENDDO
3220    ENDDO
3221#endif
3222
3223!   This may just be a fix and may need some Registry related changes, later on
3224
3225    DO J = JTS, MIN(JTE,JDE-1)
3226     DO I = ITS, MIN(ITE,IDE-1)
3227         grid%vegfra(I,J)=grid%vegfrc(I,J)
3228     ENDDO
3229    ENDDO
3230
3231!   DETERMINE ALBEDO OVER LAND ON THE BASIS OF INPUTS: SM, ALBASE, MXSNAL & VEGFRA
3232!   INTERPOLATED FROM MOTHER (WRFSI) DOMAIN
3233
3234 
3235    DO J = JTS, MIN(JTE,JDE-1)
3236     DO I = ITS, MIN(ITE,IDE-1)
3237
3238          IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN
3239!
3240            IF ( (grid%sno(I,J) .EQ. 0.0) .OR. &            ! SNOWFREE ALBEDO
3241                           (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN
3242              grid%albedo(I,J) = grid%albase(I,J)
3243            ELSE
3244              IF (grid%sno(I,J) .LT. SNUP) THEN             ! MODIFY ALBEDO IF SNOWCOVER:
3245                  RSNOW = grid%sno(I,J)/SNUP                ! BELOW SNOWDEPTH THRESHOLD
3246                  SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
3247              ELSE
3248                  SNOFAC = 1.0                         ! ABOVE SNOWDEPTH THRESHOLD
3249              ENDIF
3250              grid%albedo(I,J) = grid%albase(I,J) &
3251                          + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J))
3252            ENDIF
3253!
3254          END IF
3255
3256          grid%si(I,J)=5.0*grid%weasd(I,J)
3257          grid%sno(I,J)=grid%weasd(I,J)
3258! this block probably superfluous.  Meant to guarantee land/sea agreement
3259
3260        IF (grid%sm(I,J) .gt. 0.5)THEN
3261           grid%landmask(I,J)=0.0
3262        ELSE
3263           grid%landmask(I,J)=1.0
3264        ENDIF
3265
3266        IF (grid%sice(I,J) .eq. 1.0) then !!!! change vegtyp and sltyp to fit seaice (desireable??)
3267         grid%isltyp(I,J)=16
3268         grid%ivgtyp(I,J)=24
3269        ENDIF
3270
3271     ENDDO
3272    ENDDO
3273
3274!  Check land water interface
3275
3276    DO J = JTS, MIN(JTE,JDE-1)
3277     DO I = ITS,MIN(ITE,IDE-1)
3278      IF(grid%sm(I,J).GT.0.9 .AND. grid%vegfra(I,J) .NE. 0) THEN
3279        WRITE(message,*)'PROBLEM AT THE LAND-WATER INTERFACE (VEGFRA):', &
3280             I,J,grid%sm(I-1,J),grid%vegfra(I-1,j),grid%sm(I,J),grid%vegfra(I,J)
3281        CALL wrf_message(trim(message))
3282      ENDIF
3283!
3284#if defined(HWRF)
3285      ! HWRF should not perform the check below because the nmm_tsk is
3286      ! update with the correct skin temperature every timestep (on both
3287      ! land and sea points).
3288#else
3289      IF(grid%sm(I,J).GT.0.9 .AND. grid%nmm_tsk(I,J) .NE. 0) THEN
3290        WRITE(message,*)'PROBLEM AT THE LAND-WATER INTERFACE (NMM_TSK):', &
3291             I,J,grid%sm(I-1,J),grid%nmm_tsk(I-1,J),grid%sm(I,J),grid%nmm_tsk(I,J)
3292        CALL wrf_message(trim(message))
3293      ENDIF
3294#endif
3295     ENDDO
3296    ENDDO
3297
3298
3299!   hardwire root depth for time being
3300
3301        grid%rtdpth=0.
3302        grid%rtdpth(1)=0.1
3303        grid%rtdpth(2)=0.3
3304        grid%rtdpth(3)=0.6
3305
3306!   hardwire soil depth for time being
3307
3308        grid%sldpth=0.
3309        grid%sldpth(1)=0.1
3310        grid%sldpth(2)=0.3
3311        grid%sldpth(3)=0.6
3312        grid%sldpth(4)=1.0
3313
3314#ifdef HWRF
3315!zhang's doing: added to AVOID THIS COMPUTATION IF THE NEST IS STARTED USING ANALYSIS FILE
3316   ENDIF ! <------ for analysis set to false
3317#endif
3318!-----------  END OF LAND SURFACE INITIALIZATION -------------------------------------
3319!
3320    DO J = JTS, MIN(JTE,JDE-1)
3321     DO I = ITS, MIN(ITE,IDE-1)
3322       grid%res(I,J)=1.
3323     ENDDO
3324    ENDDO
3325
3326!   INITIALIZE 2D BOUNDARY MASKS
3327
3328!! grid%hbm2:
3329 
3330    grid%hbm2=0.
3331    DO J = JTS, MIN(JTE,JDE-1)
3332      DO I = ITS, MIN(ITE,IDE-1)
3333       IF((J .GE. 3 .and. J .LE. (JDE-1)-2) .AND.   &
3334          (I .GE. 2 .and. I .LE. (IDE-1)-2+MOD(J,2))) THEN
3335          grid%hbm2(I,J)=1.
3336        ENDIF
3337      ENDDO
3338    ENDDO   
3339
3340!! grid%hbm3:
3341
3342    grid%hbm3=0.
3343    DO J=JTS,MIN(JTE,JDE-1)
3344     grid%ihwg(J)=mod(J+1,2)-1
3345      IF (J .ge. 4 .and. J .le. (JDE-1)-3) THEN
3346        IHL=(IDS+1)-grid%ihwg(J)
3347        IHH=(IDE-1)-2
3348        DO I=ITS,MIN(ITE,IDE-1)
3349           IF (I .ge. IHL  .and. I .le. IHH) grid%hbm3(I,J)=1.
3350        ENDDO
3351      ENDIF
3352    ENDDO
3353   
3354!! grid%vbm2
3355
3356    grid%vbm2=0.
3357    DO J=JTS,MIN(JTE,JDE-1)
3358     DO I=ITS,MIN(ITE,IDE-1)
3359       IF((J .ge. 3 .and. J .le. (JDE-1)-2)  .AND.  &
3360          (I .ge. 2 .and. I .le. (IDE-1)-1-MOD(J,2))) THEN
3361           grid%vbm2(I,J)=1.
3362       ENDIF
3363     ENDDO
3364    ENDDO
3365
3366!! grid%vbm3
3367
3368    grid%vbm3=0.
3369    DO J=JTS,MIN(JTE,JDE-1)
3370      DO I=ITS,MIN(ITE,IDE-1)
3371        IF((J .ge. 4 .and. J .le. (JDE-1)-3)  .AND.  &
3372           (I .ge. 3-MOD(J,2) .and. I .le. (IDE-1)-2)) THEN
3373           grid%vbm3(I,J)=1.
3374        ENDIF
3375      ENDDO
3376    ENDDO
3377
3378    TPH0D  = grid%CEN_LAT
3379    TLM0D  = grid%CEN_LON
3380    TPH0   = TPH0D*DTR
3381    WBD    = grid%WBD0   ! gopal's doing: may use Registry WBD0 now
3382    WB     = WBD*DTR
3383    SBD    = grid%SBD0   ! gopal's doing: may use Registry SBD0 now
3384    SB     = SBD*DTR
3385    DLM    = grid%dlmd*DTR  ! input now from med_nest_egrid_configure
3386    DPH    = grid%dphd*DTR  ! input now from med_nest_egrid_configure
3387    TDLM   = DLM+DLM
3388    TDPH   = DPH+DPH
3389    WBI    = WB+TDLM
3390    SBI    = SB+TDPH
3391    EBI    = WB+((ide-1)-2)*TDLM  ! gopal's doing: check this for nested domain
3392    ANBI   = SB+((jde-1)-3)*DPH   ! gopal's doing: check this for nested domain
3393    STPH0  = SIN(TPH0)
3394    CTPH0  = COS(TPH0)
3395    TSPH   = 3600./grid%DT
3396    DTAD   = 1.0
3397    DTCF   = 4.0   
3398    DY_NMM0= grid%dy_nmm ! ERAD*DPH; input now from med_nest_egrid_configure
3399
3400!   CORIOLIS PARAMETER  (There appears to be some roundoff in computing TLM & STPH and other terms,
3401!   in the nested domain. The problem needs to be revisited
3402
3403    DO J=JTS,MIN(JTE,JDE-1)
3404      TLM0=WB-TDLM+MOD(J,2)*DLM           ! remember this is a wind point
3405      TPH =SB+float(J-1)*DPH
3406      STPH=SIN(TPH)
3407      CTPH=COS(TPH)
3408      DO I=ITS,MIN(ITE,IDE-1)
3409         TLM=TLM0 + I*TDLM
3410         FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
3411         grid%f(I,J)=0.5*grid%DT*FP
3412      ENDDO
3413    ENDDO
3414
3415
3416    DO J=JTS,MIN(JTE,JDE-1)
3417      KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
3418      KVL2(J)=(IDE-1)*(J-1)-J/2+2
3419      KHH2(J)=(IDE-1)*J-J/2-1
3420      KVH2(J)=(IDE-1)*J-(J+1)/2-1
3421    ENDDO
3422
3423
3424    TPH=SB-DPH
3425    DO J=JTS,MIN(JTE,JDE-1)
3426      TPH=SB+float(J-1)*DPH
3427      DXP=ERAD*DLM*COS(TPH)
3428      DXJ(J)=DXP
3429      WPDARJ(J)=-W_NMM*((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)/  &
3430                (grid%DT*32.*DXP*DY_NMM0)
3431      CPGFUJ(J)=-grid%DT/(48.*DXP)
3432      CURVJ(J)=.5*grid%DT*TAN(TPH)/ERAD
3433      FCPJ(J)=grid%DT/(CP*192.*DXP*DY_NMM0)
3434      FDIVJ(J)=1./(12.*DXP*DY_NMM0)
3435      FADJ(J)=-grid%DT/(48.*DXP*DY_NMM0)*DTAD
3436      ACDT=grid%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)
3437      CDDAMP=CODAMP*ACDT
3438      HDACJ(J)=COAC*ACDT/(4.*DXP*DY_NMM0)
3439      DDMPUJ(J)=CDDAMP/DXP
3440      DDMPVJ(J)=CDDAMP/DY_NMM0
3441    ENDDO
3442
3443! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
3444
3445     WRITE(message,*)'NEW CHANGE',grid%f4d,grid%ef4t,grid%f4q
3446     CALL wrf_message(trim(message))
3447
3448      DO L=KDS,KDE-1
3449        grid%rdeta(L)=1./grid%deta(L)
3450        grid%f4q2(L)=-.25*grid%DT*DTAD/grid%deta(L)
3451      ENDDO
3452
3453       DO J=JTS,MIN(JTE,JDE-1)
3454        DO I=ITS,MIN(ITE,IDE-1)
3455          grid%dx_nmm(I,J)=DXJ(J)
3456          grid%wpdar(I,J)=WPDARJ(J)*grid%hbm2(I,J)
3457          grid%cpgfu(I,J)=CPGFUJ(J)*grid%vbm2(I,J)
3458          grid%curv(I,J)=CURVJ(J)*grid%vbm2(I,J)
3459          grid%fcp(I,J)=FCPJ(J)*grid%hbm2(I,J)
3460          grid%fdiv(I,J)=FDIVJ(J)*grid%hbm2(I,J)
3461          grid%fad(I,J)=FADJ(J)
3462          grid%hdacv(I,J)=HDACJ(J)*grid%vbm2(I,J)
3463          grid%hdac(I,J)=HDACJ(J)*1.25*grid%hbm2(I,J)
3464        ENDDO
3465       ENDDO
3466
3467       DO J=JTS, MIN(JTE,JDE-1)
3468        IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
3469          KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
3470          DO I=ITS,MIN(ITE,IDE-1)
3471             IF (I .ge. 2 .and. I .le. KHH) THEN
3472               grid%hdac(I,J)=grid%hdac(I,J)* DFC
3473             ENDIF
3474          ENDDO
3475        ELSE
3476          KHH=2+MOD(J,2)
3477          DO I=ITS,MIN(ITE,IDE-1)
3478             IF (I .ge. 2 .and. I .le. KHH) THEN
3479               grid%hdac(I,J)=grid%hdac(I,J)* DFC
3480             ENDIF
3481          ENDDO
3482          KHH=(IDE-1)-2+MOD(J,2)
3483
3484          DO I=ITS,MIN(ITE,IDE-1)
3485             IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
3486               grid%hdac(I,J)=grid%hdac(I,J)* DFC
3487             ENDIF
3488          ENDDO
3489        ENDIF
3490      ENDDO
3491
3492      DO J=JTS,min(JTE,JDE-1)
3493      DO I=ITS,min(ITE,IDE-1)
3494        grid%ddmpu(I,J)=DDMPUJ(J)*grid%vbm2(I,J)
3495        grid%ddmpv(I,J)=DDMPVJ(J)*grid%vbm2(I,J)
3496        grid%hdacv(I,J)=grid%hdacv(I,J)*grid%vbm2(I,J)
3497      ENDDO
3498      ENDDO
3499
3500! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
3501
3502        DO J=JTS,MIN(JTE,JDE-1)
3503        IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
3504          KVH=(IDE-1)-1-MOD(J,2)
3505          DO I=ITS,MIN(ITE,IDE-1)
3506            IF (I .ge. 2 .and.  I .le. KVH) THEN
3507             grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
3508             grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
3509             grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
3510            ENDIF
3511          ENDDO
3512        ELSE
3513          KVH=3-MOD(J,2)
3514          DO I=ITS,MIN(ITE,IDE-1)
3515           IF (I .ge. 2 .and.  I .le. KVH) THEN
3516            grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
3517            grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
3518            grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
3519           ENDIF
3520          ENDDO
3521          KVH=(IDE-1)-1-MOD(J,2)
3522          DO I=ITS,MIN(ITE,IDE-1)
3523           IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
3524            grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
3525            grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
3526            grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
3527           ENDIF
3528          ENDDO
3529        ENDIF
3530      ENDDO
3531
3532! This one was left over for nested domain
3533 
3534     DO J = JTS, MIN(JTE,JDE-1)
3535       DO I = ITS, MIN(ITE,IDE-1)
3536          grid%GLAT(I,J)=grid%HLAT(I,J)*DTR
3537          grid%GLON(I,J)=grid%HLON(I,J)*DTR
3538       ENDDO
3539     ENDDO
3540
3541!!   compute EMT, EM on global domain, and only on task 0.
3542
3543!    IF (wrf_dm_on_monitor()) THEN   !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
3544     
3545     ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
3546     DO J=JDS,JDE-1
3547       TPH=SB+float(J-1)*DPH
3548       DXP=ERAD*DLM*COS(TPH)
3549       EMJ(J)= grid%DT/( 4.*DXP)*DTAD
3550       EMTJ(J)=grid%DT/(16.*DXP)*DTAD
3551     ENDDO
3552
3553          JA=0
3554          DO 161 J=3,5
3555          JA=JA+1
3556          KHLA(JA)=2
3557          KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
3558 161      grid%emt(JA)=EMTJ(J)
3559          DO 162 J=(JDE-1)-4,(JDE-1)-2
3560          JA=JA+1
3561          KHLA(JA)=2
3562          KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
3563 162      grid%emt(JA)=EMTJ(J)
3564          DO 163 J=6,(JDE-1)-5
3565          JA=JA+1
3566          KHLA(JA)=2
3567          KHHA(JA)=2+MOD(J,2)
3568 163      grid%emt(JA)=EMTJ(J)
3569          DO 164 J=6,(JDE-1)-5
3570          JA=JA+1
3571          KHLA(JA)=(IDE-1)-2
3572          KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
3573 164      grid%emt(JA)=EMTJ(J)
3574
3575! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
3576
3577          JA=0
3578          DO 171 J=3,5
3579          JA=JA+1
3580          KVLA(JA)=2
3581          KVHA(JA)=(IDE-1)-1-MOD(J,2)
3582 171      grid%em(JA)=EMJ(J)
3583          DO 172 J=(JDE-1)-4,(JDE-2)-2
3584          JA=JA+1
3585          KVLA(JA)=2
3586          KVHA(JA)=(IDE-1)-1-MOD(J,2)
3587 172      grid%em(JA)=EMJ(J)
3588          DO 173 J=6,(JDE-1)-5
3589          JA=JA+1
3590          KVLA(JA)=2
3591          KVHA(JA)=2+MOD(J+1,2)
3592 173      grid%em(JA)=EMJ(J)
3593          DO 174 J=6,(JDE-1)-5
3594          JA=JA+1
3595          KVLA(JA)=(IDE-1)-2
3596          KVHA(JA)=(IDE-1)-1-MOD(J,2)
3597 174      grid%em(JA)=EMJ(J)
3598
3599!        ENDIF ! wrf_dm_on_monitor
3600
3601!! must be a better place to put this, but will eliminate "phantom"
3602!! wind points here (no wind point on eastern boundary of odd numbered rows)
3603!!
3604                                                                !   phantom
3605        IF (ABS(IDE-1-ITE) .eq. 1 ) THEN                        !      |
3606         CALL wrf_message('zero phantom winds')                 !  H  [x]    H    V
3607         DO K=KDS,KDE-1                                         ! 
3608          DO J=JDS,JDE-1,2                                      !  V  [H]    V    H
3609           IF (J .ge. JTS .and. J .le. JTE) THEN                !
3610             grid%u(IDE-1,J,K)=0.                                    !  H  [x]    H    V
3611             grid%v(IDE-1,J,K)=0.                                    !  ------    ------ 
3612           ENDIF                                                !   ide-1      ide
3613          ENDDO                                                 !   NMM/si     WRF
3614         ENDDO                                                  !   domain    domain
3615        ENDIF                                                   !             (dummy)   
3616
3617
3618! just a test for gravity waves
3619
3620!  PD=62000.
3621!   grid%u=0.0
3622!   grid%v=0.0
3623!   T=300.
3624!   Q=0.0
3625!   Q2=0.0
3626!   CWM=0.0
3627!   FIS=0.0
3628
3629! testx
3630!  DO J = JTS, MIN(JTE,JDE-1)
3631!   DO K = KTS,KTE
3632!    DO I = ITS, MIN(ITE,IDE-1)
3633!      grid%sm(I,J)=I
3634!       grid%u(I,K,J)=J
3635!    ENDDO
3636!   ENDDO
3637!  ENDDO
3638!
3639
3640!   deallocs
3641
3642    DEALLOCATE(KHL2,KVL2,KHH2,KVH2)
3643    DEALLOCATE(DXJ,WPDARJ,CPGFUJ,CURVJ)
3644    DEALLOCATE(FCPJ,FDIVJ,FADJ)
3645    DEALLOCATE(HDACJ,DDMPUJ,DDMPVJ)
3646    DEALLOCATE(KHLA,KHHA)
3647    DEALLOCATE(KVLA,KVHA)
3648
3649
3650END SUBROUTINE med_initialize_nest_nmm
3651!======================================================================
3652
3653      subroutine smdhld(ime,jme,h,s,lines,nsmud)
3654      character(len=255) :: message
3655      dimension ihw(jme),ihe(jme)
3656      dimension h(ime,jme),s(ime,jme) &
3657     &         ,hbms(ime,jme),hne(ime,jme),hse(ime,jme)
3658!-----------------------------------------------------------------------
3659          do j=1,jme
3660      ihw(j)=-mod(j,2)
3661      ihe(j)=ihw(j)+1
3662          enddo
3663!-----------------------------------------------------------------------
3664
3665              do j=1,jme
3666          do i=1,ime
3667      hbms(i,j)=1.-s(i,j)
3668          enddo
3669              enddo
3670!
3671      jmelin=jme-lines+1
3672      ibas=lines/2
3673      m2l=mod(lines,2)
3674!
3675              do j=lines,jmelin
3676          ihl=ibas+mod(j,2)+m2l*mod(j+1,2)
3677          ihh=ime-ibas-m2l*mod(j+1,2)
3678
3679!
3680          do i=ihl,ihh
3681      hbms(i,j)=0.
3682          enddo
3683              enddo
3684!-----------------------------------------------------------------------
3685                  do ks=1,nsmud
3686
3687                write(message,*) 'H(1,1): ', h(1,1)
3688                CALL wrf_message(trim(message))
3689                write(message,*) 'H(3,1): ', h(1,1)
3690                CALL wrf_message(trim(message))
3691!-----------------------------------------------------------------------
3692              do j=1,jme-1
3693          do i=1,ime-1
3694      hne(i,j)=h(i+ihe(j),j+1)-h(i,j)
3695          enddo
3696              enddo
3697              do j=2,jme
3698          do i=1,ime-1
3699      hse(i,j)=h(i+ihe(j),j-1)-h(i,j)
3700          enddo
3701              enddo
3702!
3703              do j=2,jme-1
3704          do i=1+mod(j,2),ime-1
3705      h(i,j)=(hne(i,j)-hne(i+ihw(j),j-1) &
3706     &       +hse(i,j)-hse(i+ihw(j),j+1))*hbms(i,j)*0.125+h(i,j)
3707          enddo
3708              enddo
3709!-----------------------------------------------------------------------
3710
3711!!!         smooth around boundary somehow?
3712
3713!       special treatment for four corners
3714
3715        if (hbms(1,1) .eq. 1) then
3716        h(1,1)=0.75*h(1,1)+0.125*h(1+ihe(1),2)+ &
3717     &                            0.0625*(h(2,1)+h(1,3))
3718        endif
3719
3720        if (hbms(ime,1) .eq. 1) then
3721        h(ime,1)=0.75*h(ime,1)+0.125*h(ime+ihw(1),2)+ &
3722     &                            0.0625*(h(ime-1,1)+h(ime,3))
3723        endif
3724
3725        if (hbms(1,jme) .eq. 1) then
3726        h(1,jme)=0.75*h(1,jme)+0.125*h(1+ihe(jme),jme-1)+ &
3727     &                            0.0625*(h(2,jme)+h(1,jme-2))
3728        endif
3729
3730        if (hbms(ime,jme) .eq. 1) then
3731        h(ime,jme)=0.75*h(ime,jme)+0.125*h(ime+ihw(jme),jme-1)+ &
3732     &                            0.0625*(h(ime-1,jme)+h(ime,jme-2))
3733        endif
3734
3735
3736!       S bound
3737
3738        J=1
3739        do I=2,ime-1
3740        if (hbms(I,J) .eq. 1) then
3741        h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
3742        endif
3743        enddo
3744
3745!       N bound
3746
3747        J=JME
3748        do I=2,ime-1
3749        if (hbms(I,J) .eq. 1) then
3750        h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
3751        endif
3752        enddo
3753
3754!       W bound
3755
3756        I=1
3757        do J=3,jme-2
3758        if (hbms(I,J) .eq. 1) then
3759        h(I,J)=0.75*h(I,J)+0.125*(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
3760        endif
3761        enddo
3762
3763!       E bound
3764
3765        I=IME
3766        do J=3,jme-2
3767        if (hbms(I,J) .eq. 1) then
3768        h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
3769        endif
3770        enddo
3771
3772
3773                      enddo   ! end ks loop
3774
3775
3776!-----------------------------------------------------------------------
3777      return
3778      end subroutine smdhld
3779
3780!--------------------------------------------------------------------------------------
3781#if 0
3782SUBROUTINE initial_nest_pivot ( parent , nest, iloc, jloc )
3783
3784!==========================================================================================
3785!
3786! This program produces i_start and j_start for the nested domain depending on the
3787! central lat-lon of the storm.
3788!
3789!==========================================================================================
3790
3791 USE module_domain
3792 USE module_configure
3793 USE module_timing
3794 USE module_dm
3795
3796 IMPLICIT NONE
3797 TYPE(domain) , POINTER              :: parent , nest
3798 INTEGER, INTENT(OUT)                :: ILOC,JLOC
3799 INTEGER                             :: IMS,IME,JMS,JME,KMS,KME
3800 INTEGER                             :: IDS,IDE,JDS,JDE,KDS,KDE
3801 INTEGER                             :: IMS,IME,JMS,JME,KMS,KME
3802 INTEGER                             :: ITS,ITE,JTS,JTE,KTS,KTE
3803 INTEGER                             :: NIDE,NJDE              ! nest dimension
3804 INTEGER                             :: I,J,ITER,IDUM,JDUM
3805 REAL                                :: ALAT,ALON,DIFF1,DIFF2,ERR
3806 REAL                                :: parent_CLAT,parent_CLON,parent_SLAT,parent_SLON
3807 REAL                                :: parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
3808!========================================================================================
3809
3810!   First obtain central latitude and longitude for the parent domain
3811
3812    CALL nl_get_cen_lat (parent%ID, parent_CLAT)
3813    CALL nl_get_cen_lon (parent%ID, parent_CLON)
3814!    CALL nl_get_storm_lat (parent%ID, parent_SLAT)
3815!    CALL nl_get_storm_lon (parent%ID, parent_SLON)
3816
3817!   Parent grid configuration, including, western and southern boundary
3818
3819    IDS = parent%sd31
3820    IDE = parent%ed31
3821    JDS = parent%sd32
3822    JDE = parent%ed32
3823    KDS = parent%sd33
3824    KDE = parent%ed33
3825
3826    IMS = parent%sm31
3827    IME = parent%em31
3828    JMS = parent%sm32
3829    JME = parent%em32
3830    KMS = parent%sm33
3831    KME = parent%em33
3832
3833    ITS  = parent%sp31
3834    ITE  = parent%ep31
3835    JTS  = parent%sp32
3836    JTE  = parent%ep32
3837    KTS  = parent%sp33
3838    KTE  = parent%ep33
3839
3840    NIDE = nest%ed31
3841    NJDE = nest%ed32
3842
3843    parent_DLMD = parent%dx          ! DLMD: dlamda in degrees
3844    parent_DPHD = parent%dy          ! DPHD: dphi in degrees
3845    parent_WBD  = -(IDE-2)*parent%dx ! WBD0: in deg;factor 2 takes care of dummy last column
3846    parent_SBD  = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd
3847    ALAT  = parent_SLAT - 0.5*(NJDE-2)*parent_DPHD/nest%parent_grid_ratio
3848    ALON  = parent_SLON - 1.0*(NIDE-2)*parent_DLMD/nest%parent_grid_ratio
3849
3850    CALL EARTH_LATLON ( parent%HLAT,parent%HLON,parent%VLAT,parent%VLON, & !output
3851                        parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,                   & !inputs
3852                        parent_CLAT,parent_CLON,                                         &
3853                        IDS,IDE,JDS,JDE,KDS,KDE,                                         &
3854                        IMS,IME,JMS,JME,KMS,KME,                                         &
3855                        ITS,ITE,JTS,JTE,KTS,KTE                          )
3856
3857!   start iteration
3858
3859      ILOC=-99
3860      JLOC=-99
3861      ERR=0.1
3862      ITER=1
3863100   CONTINUE
3864
3865     DO J = JTS,min(JTE,JDE-1)
3866      DO I = ITS,min(ITE,IDE-1)
3867        DIFF1 = ABS(ALAT - parent%HLAT(I,J))
3868        DIFF2 = ABS(ALON - parent%HLON(I,J))
3869        IF(DIFF1 .LE. ERR .AND. DIFF2 .LE. ERR)THEN
3870         ILOC=I
3871         JLOC=J
3872        ENDIF
3873      ENDDO
3874     ENDDO
3875
3876     CALL wrf_dm_maxval_integer ( ILOC, idum, jdum )
3877     CALL wrf_dm_maxval_integer ( JLOC, idum, jdum )
3878
3879     IF(ILOC .EQ. -99 .AND. JLOC .EQ. -99)THEN
3880       ERR=ERR+0.1
3881       ITER=ITER+1
3882       IF(ITER .LE. 100)GO TO 100
3883     ENDIF
3884
3885     IF(ILOC .NE. -99 .AND. JLOC .NE. -99)THEN
3886       WRITE(message,*)'NOTE: I_PARENT_START AND J_PARENT_START FOUND FOR THE NESTED DOMAIN CONFIGURATION AT ITER=',ITER
3887       CALL wrf_message(trim(message))
3888       WRITE(message,*)'istart=',ILOC
3889       CALL wrf_message(trim(message))
3890       WRITE(message,*)'jstart=',JLOC
3891       CALL wrf_message(trim(message))
3892     ELSE
3893       ILOC=IDE/3
3894       JLOC=JDE/3
3895!
3896       WRITE(message,*)'WARNING: COULD NOT LOCATE I_PARENT_START AND J_PARENT_START FROM INPUT STORM INFO'
3897       CALL wrf_message(trim(message))
3898       WRITE(message,*)'ISTART=',IDE/3
3899       CALL wrf_message(trim(message))
3900       WRITE(message,*)'JSTART=',JDE/3
3901       CALL wrf_message(trim(message))
3902     ENDIF
3903
3904     RETURN
3905END SUBROUTINE initial_nest_pivot
3906
3907!============================================================================================
3908#endif
3909
3910LOGICAL FUNCTION cd_feedback_mask_orig( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
3911   INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
3912   LOGICAL, INTENT(IN) :: xstag, ystag
3913
3914   INTEGER ioff, joff
3915
3916   ioff = 0 ; joff = 0
3917   IF ( xstag  ) ioff = 1
3918   IF ( ystag  ) joff = 1
3919
3920   cd_feedback_mask_orig = ( pig .ge. ips_save+1        .and.      &
3921                            pjg .ge. jps_save+1        .and.      &
3922                            pig .le. ipe_save-1  +ioff .and.      &
3923                            pjg .le. jpe_save-1  +joff           )
3924
3925END FUNCTION cd_feedback_mask_orig
3926
3927LOGICAL FUNCTION cd_feedback_mask( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
3928   INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
3929   LOGICAL, INTENT(IN) :: xstag, ystag
3930
3931   INTEGER ioff, joff
3932
3933   ioff = 0 ; joff = 0
3934   IF ( xstag  ) ioff = 1
3935   IF ( ystag  ) joff = 1
3936
3937   cd_feedback_mask = ( pig .ge. ips_save+1 .and. &
3938                            pjg .ge. jps_save+2 .and. &
3939                            pig .le. ipe_save-1-mod(pjg-jps_save,2) .and. &
3940                            pjg .le. jpe_save-2 )
3941
3942END FUNCTION cd_feedback_mask
3943
3944LOGICAL FUNCTION cd_feedback_mask_v( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
3945   INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
3946   LOGICAL, INTENT(IN) :: xstag, ystag
3947
3948   INTEGER ioff, joff
3949
3950   ioff = 0 ; joff = 0
3951   IF ( xstag  ) ioff = 1
3952   IF ( ystag  ) joff = 1
3953   
3954   cd_feedback_mask_v = ( pig .ge. ips_save+1 .and. &
3955                            pjg .ge. jps_save+2 .and. &
3956                            pig .le. ipe_save-1-mod(pjg-jps_save+1,2) .and. &
3957                            pjg .le. jpe_save-2 )
3958
3959END FUNCTION cd_feedback_mask_v
3960
3961
3962!----------------------------------------------------------------------------
3963#else
3964SUBROUTINE stub_nmm_nest_stub
3965END SUBROUTINE stub_nmm_nest_stub
3966#endif
3967
3968RECURSIVE SUBROUTINE find_ijstart_level ( grid, i_start, j_start, level )
3969
3970! Dusan Jovic
3971
3972   USE module_domain
3973
3974   IMPLICIT NONE
3975
3976   !  Input data.
3977
3978   TYPE(domain) :: grid
3979   INTEGER, INTENT (OUT) :: i_start, j_start, level
3980   INTEGER :: iadd
3981
3982      if (grid%parent_id == 0 ) then
3983         i_start = 1
3984         j_start = 1
3985         level = 0
3986      else
3987         call find_ijstart_level ( grid%parents(1)%ptr, i_start, j_start, level )
3988         if (level > 0) then
3989             iadd = (i_start-1)*3
3990             if ( mod(j_start,2).ne.0 .and. mod(grid%j_parent_start,2).ne.0 ) iadd = iadd - 1
3991             if ( mod(j_start,2).eq.0 .and. mod(grid%j_parent_start,2).eq.0 ) iadd = iadd + 2
3992         else
3993             iadd = -mod(grid%j_parent_start,2)
3994         end if
3995         i_start = iadd + grid%i_parent_start*3 - 1
3996         j_start = ( (j_start-1) + (grid%j_parent_start-1) ) * 3 + 1
3997         level = level + 1
3998      end if
3999
4000END SUBROUTINE find_ijstart_level
Note: See TracBrowser for help on using the repository browser.