/wrfv2_fire/dyn_nmm/NMM_NEST_UTILS1.F

Large files files are truncated, but you can click here to view the full file

#if (NMM_NEST == 1)
!===========================================================================
!
! E-GRID NESTING UTILITIES: This is gopal's doing
!
!===========================================================================

SUBROUTINE med_nest_egrid_configure ( parent , nest )
 USE module_domain
 USE module_configure
 USE module_timing

 IMPLICIT NONE
 TYPE(domain) , POINTER             :: parent , nest
 REAL, PARAMETER                    :: ERAD=6371200.
 REAL, PARAMETER                    :: DTR=0.01745329
 REAL, PARAMETER                    :: DTAD=1.0
 REAL, PARAMETER                    :: CP=1004.6
 INTEGER                            :: IDS,IDE,JDS,JDE,KDS,KDE
 INTEGER                            :: IMS,IME,JMS,JME,KMS,KME
 INTEGER                            :: ITS,ITE,JTS,JTE,KTS,KTE
 CHARACTER(LEN=255)                 :: message

!----------------------------------------------------------------------------
!  PURPOSE: 
!         - Initialize nested domain configurations including setting up 
!           wbd,sbd and some other variables and 1D arrays. 
!         - Note that in order to obtain coincident grid points, which  
!           is a basic requirement for RSL, WRF infrastructure, we use 
!           western and southern boundaries of nested domain (nest%wbd0
!           and nest%sbd0 derived from the parent domain. In this case
!           the nested domain may be considered as a part of the parent 
!           domain with a higher resolution (telescoping ?). 
!         - Also note that in this case, the central lat/lons for nested 
!           domain should coincide with the central lat/lons of the parent,
!           although the nested domain NEED NOT be located at the center of
!           the domain.
!----------------------------------------------------------------------------
!
!   BASIC TEST FOR PARENT DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
!   IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT

    IF(MOD(parent%ed32,2) .NE. 0)THEN
     CALL wrf_error_fatal("PARENT DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
    ENDIF

!   BASIC TEST FOR NESTED DOMAIN: CHECK IF JMAX IS ODD. SINCE JDE IN THE NAMELIST
!   IS JMAX + 1, WE NEED TO CHECK IF JDE IS EVEN IN WRF CONTEXT

    IF(MOD(nest%ed32,2) .NE. 0)THEN
     CALL wrf_error_fatal("NESTED DOMAIN: JMAX IS EVEN, INCREASE e_sn IN THE namelist.input BY 1")
    ENDIF

!   Parent grid configuration, including, western and southern boundary

    IDS = parent%sd31
    IDE = parent%ed31
    JDS = parent%sd32
    JDE = parent%ed32
    KDS = parent%sd33
    KDE = parent%ed33

    IMS = parent%sm31
    IME = parent%em31
    JMS = parent%sm32
    JME = parent%em32
    KMS = parent%sm33
    KME = parent%em33

    ITS = parent%sp31
    ITE = parent%ep31
    JTS = parent%sp32
    JTE = parent%ep32
    KTS = parent%sp33
    KTE = parent%ep33

!   grid configuration

    ! calculate wbd0 and sbd0 only for MOAD i.e. grid with parent_id == 0
    if (parent%parent_id == 0 ) then        ! Dusan's doing
       parent%wbd0 = -(IDE-2)*parent%dx     ! WBD0: in degrees;factor 2 takes care of dummy last column
       parent%sbd0 = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd  
    end if
    nest%wbd0   = parent%wbd0 + (nest%i_parent_start-1)*2.*parent%dx + mod(nest%j_parent_start+1,2)*parent%dx
    nest%sbd0   = parent%sbd0 + (nest%j_parent_start-1)*parent%dy
    nest%dx     = parent%dx/nest%parent_grid_ratio
    nest%dy     = parent%dy/nest%parent_grid_ratio

    write(message,*)" - i_parent_start = ",nest%i_parent_start
    CALL wrf_message(trim(message))
    write(message,*)" - j_parent_start = ",nest%j_parent_start
    CALL wrf_message(trim(message))
    write(message,*)" - parent%wbd0    = ",parent%wbd0
    CALL wrf_message(trim(message))
    write(message,*)" - parent%sbd0    = ",parent%sbd0
    CALL wrf_message(trim(message))
    write(message,*)" - nest%wbd0      = ",nest%wbd0
    CALL wrf_message(trim(message))
    write(message,*)" - nest%sbd0      = ",nest%sbd0
    CALL wrf_message(trim(message))
    write(message,*)" - nest%dx        = ",nest%dx
    CALL wrf_message(trim(message))
    write(message,*)" - nest%dy        = ",nest%dy
    CALL wrf_message(trim(message))
!
    CALL nl_set_dx (nest%id , nest%dx)   ! for output purpose
    CALL nl_set_dy (nest%id , nest%dy)   ! for output purpose

!   set lat-lons; parent set to nested domain

    CALL nl_get_cen_lat (parent%id, parent%cen_lat) ! cen_lat of parent set to nested domain
    CALL nl_get_cen_lon (parent%id, parent%cen_lon) ! cen_lon of parent set to nested domain

    nest%cen_lat=parent%cen_lat
    nest%cen_lon=parent%cen_lon
!
    CALL nl_set_cen_lat ( nest%id , nest%cen_lat)  ! for output purpose
    CALL nl_set_cen_lon ( nest%id , nest%cen_lon)  ! for output purpose

    write(message,*)" - nest%cen_lat   = ",nest%cen_lat
    CALL wrf_message(trim(message))
    write(message,*)" - nest%cen_lon   = ",nest%cen_lon
    CALL wrf_message(trim(message))


!   soil configuration

#ifdef HWRF
!zhang 
    if ( .not.nest%analysis ) then
#endif
    nest%sldpth  = parent%sldpth
    nest%dzsoil  = parent%dzsoil
    nest%rtdpth  = parent%rtdpth
#ifdef HWRF 
    endif
#endif

!   numerical set up

    nest%deta        = parent%deta
    nest%aeta        = parent%aeta
    nest%etax        = parent%etax
    nest%dfl         = parent%dfl
    nest%deta1       = parent%deta1
    nest%aeta1       = parent%aeta1
    nest%eta1        = parent%eta1
    nest%deta2       = parent%deta2
    nest%aeta2       = parent%aeta2
    nest%eta2        = parent%eta2
    nest%pdtop       = parent%pdtop
    nest%pt          = parent%pt
    nest%dfrlg       = parent%dfrlg
    nest%num_soil_layers = parent%num_soil_layers
    nest%num_moves       = parent%num_moves

! Unfortunately, some of the single value constants in used in module_initialize have 
! to be defiend here instead of the usual spot in med_initialize_nest_nmm. There
! appears to be a problem in Registry and related code in this area.
!
! state  logical upstrm   -      dyn_nmm     -      -      -


    nest%dlmd   = nest%dx
    nest%dphd   = nest%dy
    nest%dy_nmm = erad*(nest%dphd*dtr)
    nest%cpgfv  = -nest%dt/(48.*nest%dy_nmm)
    nest%en     = nest%dt/( 4.*nest%dy_nmm)*dtad
    nest%ent    = nest%dt/(16.*nest%dy_nmm)*dtad
    nest%f4d    = -.5*nest%dt*dtad
    nest%f4q    = -nest%dt*dtad
    nest%ef4t   = .5*nest%dt/cp

!  Other output configurations that will make grads happy 

   CALL nl_get_truelat1 (parent%id, parent%truelat1 )
   CALL nl_get_truelat2 (parent%id, parent%truelat2 )
#ifdef HWRF
! bao : to make the restart output identical at the restart initial time for stand_lon
   CALL nl_get_stand_lon (parent%id, parent%stand_lon )
#endif
   CALL nl_get_map_proj (parent%id, parent%map_proj )
   CALL nl_get_iswater (parent%id, parent%iswater )

   nest%truelat1=parent%truelat1
   nest%truelat2=parent%truelat2
!bao
   nest%stand_lon=parent%stand_lon
!bao
   nest%map_proj=parent%map_proj
   nest%iswater=parent%iswater

   CALL nl_set_truelat1(nest%id, nest%truelat1)
   CALL nl_set_truelat2(nest%id, nest%truelat2)
!bao
   CALL nl_set_stand_lon(nest%id, nest%stand_lon)
!bao
   CALL nl_set_map_proj(nest%id, nest%map_proj)
   CALL nl_set_iswater(nest%id, nest%iswater)

!   physics and other configurations
!   CALL nl_get_iswater (parent%id, nest%iswater) ! iswater is just based on parents
!   CALL nl_get_bl_surface_physics (nest%id,  nest%bl_surface_physics )
!   CALL nl_get_num_soil_layers( parent%num_soil_layers )
!   CALL nl_get_real_data_init_type (parent%real_data_init_type)

END SUBROUTINE med_nest_egrid_configure

SUBROUTINE med_construct_egrid_weights ( parent , nest )
 USE module_domain
 USE module_configure
 USE module_timing

 IMPLICIT NONE
 TYPE(domain) , POINTER             :: parent , nest
 LOGICAL, EXTERNAL                  :: wrf_dm_on_monitor
 INTEGER                            :: IDS,IDE,JDS,JDE,KDS,KDE
 INTEGER                            :: IMS,IME,JMS,JME,KMS,KME
 INTEGER                            :: ITS,ITE,JTS,JTE,KTS,KTE
 INTEGER                            :: I,J,II,JJ,NII,NJJ
 REAL                               :: parent_CLAT,parent_CLON,parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
 REAL                               :: nest_WBD,nest_SBD,nest_DLMD,nest_DPHD
 REAL                               :: SW_LATD,SW_LOND
 REAL                               :: ADDSUM1,ADDSUM2
 REAL                               :: xr,zr,xc
!-----------------------------------------------------------------------------------------------------------
!   PURPOSE: 
!           - Initialize lat-lons and determine weights 
!
!----------------------------------------------------------------------------------------------------------

!   First obtain central latitude and longitude for the parent domain

    CALL nl_get_cen_lat (parent%ID, parent_CLAT)
    CALL nl_get_cen_lon (parent%ID, parent_CLON)

!   Parent grid configuration, including, western and southern boundary

    IDS = parent%sd31
    IDE = parent%ed31
    JDS = parent%sd32
    JDE = parent%ed32
    KDS = parent%sd33
    KDE = parent%ed33

    IMS = parent%sm31
    IME = parent%em31
    JMS = parent%sm32
    JME = parent%em32
    KMS = parent%sm33
    KME = parent%em33

    ITS  = parent%sp31
    ITE  = parent%ep31
    JTS  = parent%sp32
    JTE  = parent%ep32
    KTS  = parent%sp33
    KTE  = parent%ep33
!
    parent_DLMD = parent%dx                ! DLMD: dlamda in degrees 
    parent_DPHD = parent%dy                ! DPHD: dphi in degrees 
    parent_WBD  = parent%wbd0
    parent_SBD  = parent%sbd0

!   Now compute Geodetic lat/lon (Positive East) of parent grid in degrees

    CALL EARTH_LATLON ( parent%HLAT,parent%HLON,parent%VLAT,parent%VLON,  & !output
                        parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,                    & !inputs 
                        parent_CLAT,parent_CLON,                                          &
                        IDS,IDE,JDS,JDE,KDS,KDE,                                          &
                        IMS,IME,JMS,JME,KMS,KME,                                          &
                        ITS,ITE,JTS,JTE,KTS,KTE                                           )

!   Nested grid configuration, including, western and southern boundary

    IDS = nest%sd31
    IDE = nest%ed31
    JDS = nest%sd32
    JDE = nest%ed32
    KDS = nest%sd33
    KDE = nest%ed33

    IMS = nest%sm31
    IME = nest%em31
    JMS = nest%sm32
    JME = nest%em32
    KMS = nest%sm33
    KME = nest%em33

    ITS  = nest%sp31
    ITE  = nest%ep31
    JTS  = nest%sp32
    JTE  = nest%ep32
    KTS  = nest%sp33
    KTE  = nest%ep33
!
    nest_DLMD = nest%dx
    nest_DPHD = nest%dy
    nest_WBD  = nest%wbd0
    nest_SBD  = nest%sbd0

!
!   Now compute Geodetic lat/lon (Positive East) of nest in degrees, with the same central lat-lon
!   as the parent grid
!

    CALL EARTH_LATLON ( nest%HLAT,nest%HLON,nest%VLAT,nest%VLON, & ! output
                        nest_DLMD,nest_DPHD,nest_WBD,nest_SBD,                   & ! nest inputs
                        parent_CLAT,parent_CLON,                                 & ! parent central lat/lon
                        IDS,IDE,JDS,JDE,KDS,KDE,                                 & ! nested domain dimension
                        IMS,IME,JMS,JME,KMS,KME,                                 &
                        ITS,ITE,JTS,JTE,KTS,KTE                                  )

!   Determine the weights of nested grid h points nearest to H points of parent domain 

  if(nest%vortex_tracker /= 1) then
    CALL G2T2H_new(    nest%IIH,nest%JJH,                            & ! output grid index in parent grid
                       nest%HBWGT1,nest%HBWGT2,                      & ! output weights in terms of parent grid
                       nest%HBWGT3,nest%HBWGT4,                      &
                       nest%I_PARENT_START,nest%J_PARENT_START,      & ! nest start I, J in parent domain
                       3,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
                       IMS,IME,JMS,JME,KMS,KME,            &
                       ITS,ITE,JTS,JTE,KTS,KTE      )
  else
    CALL G2T2H( nest%IIH,nest%JJH,                       & ! output grid index on nested grid
                nest%HBWGT1,nest%HBWGT2,                 & ! output weights on the nested grid 
                nest%HBWGT3,nest%HBWGT4,                 &
                nest%HLAT,nest%HLON,                     & ! target (nest) input lat lon in degrees
                parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,   & ! parent res, western and south boundaries
                parent_CLAT,parent_CLON,                         & ! parent central lat,lon, all in degrees
                parent%ed31,parent%ed32,                         & ! parent imax and jmax
                IDS,IDE,JDS,JDE,KDS,KDE,                         & ! 
                IMS,IME,JMS,JME,KMS,KME,                         & ! nested grid configuration
                ITS,ITE,JTS,JTE,KTS,KTE                          ) !
  endif


!   Determine the weights of nested grid v points nearest to V points of parent domain

  if(nest%vortex_tracker /= 1) then
    CALL G2T2V_new(    nest%IIV,nest%JJV,                            & ! output grid index in parent grid
                       nest%VBWGT1,nest%VBWGT2,                      & ! output weights in terms of parent grid
                       nest%VBWGT3,nest%VBWGT4,                      &
                       nest%I_PARENT_START,nest%J_PARENT_START,      & ! nest start I, J in parent domain
                       3,                              & ! Ratio of parent and child grid ( always = 3 for NMM)
                       IDS,IDE,JDS,JDE,KDS,KDE,            & ! target (nest) dimensions
                       IMS,IME,JMS,JME,KMS,KME,            &
                       ITS,ITE,JTS,JTE,KTS,KTE      )
  else
    CALL G2T2V( nest%IIV,nest%JJV,                       & ! output grid index on nested grid
                nest%VBWGT1,nest%VBWGT2,                 & ! output weights on the nested grid
                nest%VBWGT3,nest%VBWGT4,                 &
                nest%VLAT,nest%VLON,                     & ! target (nest) input lat lon in degrees
                parent_DLMD,parent_DPHD,parent_WBD,parent_SBD,   & ! parent res, western and south boundaries
                parent_CLAT,parent_CLON,                         & ! parent central lat,lon, all in degrees
                parent%ed31,parent%ed32,                         & ! parent imax and jmax
                IDS,IDE,JDS,JDE,KDS,KDE,                         & !
                IMS,IME,JMS,JME,KMS,KME,                         & ! nested grid configuration
                ITS,ITE,JTS,JTE,KTS,KTE                          ) !
  endif

!*** CHECK WEIGHTS AT MASS AND VELOCITY POINTS

     CALL WEIGTS_CHECK(nest%HBWGT1,nest%HBWGT2,          & ! output weights on the nested grid
                       nest%HBWGT3,nest%HBWGT4,          &
                       nest%VBWGT1,nest%VBWGT2,          & ! output weights on the nested grid
                       nest%VBWGT3,nest%VBWGT4,          &
                       IDS,IDE,JDS,JDE,KDS,KDE,                  & !
                       IMS,IME,JMS,JME,KMS,KME,                  & ! nested grid configuration
                       ITS,ITE,JTS,JTE,KTS,KTE                   )

!*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
!
    CALL BOUNDS_CHECK( nest%IIH,nest%JJH,nest%IIV,nest%JJV,   &
                       nest%i_parent_start,nest%j_parent_start,nest%shw,      &
                       IDS,IDE,JDS,JDE,KDS,KDE,                               & !
                       IMS,IME,JMS,JME,KMS,KME,                               & ! nested grid configuration
                       ITS,ITE,JTS,JTE,KTS,KTE                                )

!------------------------------------------------------------------------------------------

END SUBROUTINE med_construct_egrid_weights 

!======================================================================================
!
! compute earth lat-lons for parent and the nest before interpolations
!------------------------------------------------------------------------------

SUBROUTINE EARTH_LATLON ( HLAT,HLON,VLAT,VLON,     & !Earth lat,lon at H and V points
                          DLMD1,DPHD1,WBD1,SBD1,   & !input res,west & south boundaries,
                          CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees   
                          IDS,IDE,JDS,JDE,KDS,KDE, &  
                          IMS,IME,JMS,JME,KMS,KME, &
                          ITS,ITE,JTS,JTE,KTS,KTE  )
!============================================================================
!
 IMPLICIT NONE
 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME 
 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE  
 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT)        :: HLAT,HLON,VLAT,VLON

! local

 INTEGER,PARAMETER                           :: KNUM=SELECTED_REAL_KIND(13) 
 INTEGER                                     :: I,J
 REAL(KIND=KNUM)                             :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
 REAL(KIND=KNUM)                             :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR
 REAL(KIND=KNUM)                             :: STPH,CTPH,STPV,CTPV,PI_2
 REAL(KIND=KNUM)                             :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV
 REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV
!-------------------------------------------------------------------------

!
      PI_2 = ACOS(0.)
      DTR  = PI_2/90.
      WB   = WBD1 * DTR                 ! WB:   western boundary in radians
      SB   = SBD1 * DTR                 ! SB:   southern boundary in radians
      DLM  = DLMD1 * DTR                ! DLM:  dlamda in radians 
      DPH  = DPHD1 * DTR                ! DPH:  dphi   in radians
      TDLM = DLM + DLM                  ! TDLM: 2.0*dlamda 
      TDPH = DPH + DPH                  ! TDPH: 2.0*DPH 

!     For earth lat lon only

      TPH0  = CENTRAL_LAT*DTR                ! TPH0: central lat in radians 
      STPH0 = SIN(TPH0)
      CTPH0 = COS(TPH0)

                                                !    .H
      DO J = JTS,MIN(JTE,JDE-1)                 ! H./    This loop takes care of zig-zag 
!                                               !   \.H  starting points along j  
         TLMH0 = WB - TDLM + MOD(J+1,2) * DLM   !  ./    TLMH (rotated lats at H points)
         TLMV0 = WB - TDLM + MOD(J,2) * DLM     !  H     (//ly for V points) 
         TPHH = SB + (J-1)*DPH                  !   TPHH (rotated lons at H points) are simple trans.
         TPHV = TPHH                            !   TPHV (rotated lons at V points) are simple trans.
         STPH = SIN(TPHH)
         CTPH = COS(TPHH)
         STPV = SIN(TPHV)
         CTPV = COS(TPHV)

                                                              !   .H
         DO I = ITS,MIN(ITE,IDE-1)                            !  / 
           TLMH = TLMH0 + I*TDLM                              !  \.H   .U   .H 
!                                                             !H./ ----><----
           SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH)     !     DLM + DLM
           GLATH(I,J)=ASIN(SPHH)                              ! GLATH: Earth Lat in radians 
           CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) &
                - TAN(GLATH(I,J))*TAN(TPH0)
           IF(CLMH .GT. 1.) CLMH = 1.0
           IF(CLMH .LT. -1.) CLMH = -1.0
           FACTH = 1.
           IF(TLMH .GT. 0.) FACTH = -1.
           GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH)

         ENDDO                                    

         DO I = ITS,MIN(ITE,IDE-1)
           TLMV = TLMV0 + I*TDLM
           SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV)
           GLATV(I,J) = ASIN(SPHV)
           CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) &
                - TAN(GLATV(I,J))*TAN(TPH0)
           IF(CLMV .GT. 1.) CLMV = 1.
           IF(CLMV .LT. -1.) CLMV = -1.
           FACTV = 1.
           IF(TLMV .GT. 0.) FACTV = -1.
           GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV)

         ENDDO

      ENDDO

!     Conversion to degrees (may not be required, eventually)

      DO J = JTS, MIN(JTE,JDE-1)
       DO I = ITS, MIN(ITE,IDE-1)
          HLAT(I,J) = GLATH(I,J) / DTR
          HLON(I,J)= -GLONH(I,J)/DTR
          IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J)  - 360.
          IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360.
!
          VLAT(I,J) = GLATV(I,J) / DTR
          VLON(I,J) = -GLONV(I,J) / DTR
          IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J)  - 360.
          IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360.

       ENDDO
      ENDDO

END SUBROUTINE EARTH_LATLON

!-----------------------------------------------------------------------------  

  SUBROUTINE G2T2H( IIH,JJH,                     & ! output grid index and weights 
                    HBWGT1,HBWGT2,               & ! output weights in terms of parent grid
                    HBWGT3,HBWGT4,               & 
                    HLAT,HLON,                   & ! target (nest) input lat lon in degrees
                    DLMD1,DPHD1,WBD1,SBD1,       & ! parent res, west and south boundaries
                    CENTRAL_LAT,CENTRAL_LON,     & ! parent central lat,lon, all in degrees
                    P_IDE,P_JDE,                 & ! parent imax and jmax
                    IDS,IDE,JDS,JDE,KDS,KDE,     & ! target (nest) dIMEnsions
                    IMS,IME,JMS,JME,KMS,KME,     &
                    ITS,ITE,JTS,JTE,KTS,KTE      )

! 
!***  Tom Black   - Initial Version
!***  Gopal       - Revised Version for WRF (includes coincident grid points) 
!*** 
!***  GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
!***  AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE 
!***  INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE
!***  h POINTS OF THE NESTED DOMAIN   
!
!============================================================================
!
 IMPLICIT NONE
 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
 INTEGER,    INTENT(IN   )                            :: P_IDE,P_JDE
 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(IN)   :: HLAT,HLON
 REAL,    DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
 INTEGER, DIMENSION(IMS:IME,JMS:JME),    INTENT(OUT)  :: IIH,JJH

! local

 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
 INTEGER           :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
 INTEGER           :: NROW,NCOL,KROWS
 REAL(KIND=KNUM)   :: X,Y,Z,TLAT,TLON
 REAL(KIND=KNUM)   :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB                
 REAL(KIND=KNUM)   :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE
 REAL(KIND=KNUM)   :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
 REAL(KIND=KNUM)   :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3                    ! Q
 REAL(KIND=KNUM)   :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO 
 REAL(KIND=KNUM)   :: DTEMP
 REAL   , DIMENSION(IMS:IME,JMS:JME)    :: TLATHX,TLONHX
 INTEGER, DIMENSION(IMS:IME,JMS:JME)    :: KOUTB
!------------------------------------------------------------------------------- 

  IMT=2*P_IDE-2             ! parent i dIMEnsions 
  JMT=P_JDE/2               ! parent j dIMEnsions 
  PI_2=ACOS(0.)
  D2R=PI_2/90.
  R2D=1./D2R
  DPH=DPHD1*D2R
  DLM=DLMD1*D2R
  TPH0= CENTRAL_LAT*D2R
  TLM0=-CENTRAL_LON*D2R        ! NOTE THE MINUS HERE
  WB=WBD1*D2R                   ! CONVERT NESTED GRID H POINTS FROM GEODETIC
  SB=SBD1*D2R
  SLP0=DPHD1/DLMD1
  DSLP0=NINT(R2D*ATAN(SLP0))
  DS1=SQRT(DPH*DPH+DLM*DLM)    ! Q
  AN1=ASIN(DLM/DS1)
  AN2=ASIN(DPH/DS1)

  DO J =  JTS,MIN(JTE,JDE-1) 
    DO I = ITS,MIN(ITE,IDE-1) 

!***
!***  LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND
!***  DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST 
!***  CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED 
!***  COORDINATE ON THE PARENT GRID
!

      GLAT=HLAT(I,J)*D2R
      GLON= (360. - HLON(I,J))*D2R
      X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
      Y=-COS(GLAT)*SIN(GLON-TLM0)
      Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
      TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
      TLON=R2D*ATAN(Y/X)

!      ROW=TLAT/DPHD1+JMT         ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
!      COL=TLON/DLMD1+P_IDE-1     ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN 

      ROW=(TLAT-SBD1)/DPHD1+1     ! Dusan's doing
      COL=(TLON-WBD1)/DLMD1+1     ! Dusan's doing

      NROW=INT(ROW + 0.001)     ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
      NCOL=INT(COL + 0.001)     
      TLAT=TLAT*D2R
      TLON=TLON*D2R

!***  
!***
!***  FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT
!***
!***              V      H
!***
!***
!***                 H           
!***              H      V
!***
!***  THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
!***
      IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR.     &     
         MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN        
           TLAT1=(NROW-JMT)*DPH                           
           TLAT2=TLAT1+DPH                              
           TLON1=(NCOL-(P_IDE-1))*DLM                                   
           TLON2=TLON1+DLM                                 
           DLM1=TLON-TLON1                                 
           DLM2=TLON-TLON2
!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
           D1=ACOS(DTEMP)
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           D2=ACOS(DTEMP)
            IF(D1.GT.D2)THEN
             NROW=NROW+1                    ! FIND THE NEAREST H ROW
             NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
            ENDIF 
      ELSE
!***
!***  NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT
!***
!***              H      V
!***
!***
!***                 H 
!***              V      H
!***
!***  THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
!***
!***
           TLAT1=(NROW+1-JMT)*DPH
           TLAT2=TLAT1-DPH
           TLON1=(NCOL-(P_IDE-1))*DLM
           TLON2=TLON1+DLM
           DLM1=TLON-TLON1
           DLM2=TLON-TLON2
!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
           D1=ACOS(DTEMP)
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           D2=ACOS(DTEMP)
             IF(D1.LT.D2)THEN
              NROW=NROW+1                    ! FIND THE NEAREST H ROW
             ELSE
              NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
             ENDIF
      ENDIF

      KROWS=((NROW-1)/2)*IMT
      IF(MOD(NROW,2).EQ.1)THEN
        K=KROWS+(NCOL+1)/2
      ELSE
        K=KROWS+P_IDE-1+NCOL/2
      ENDIF

!***
!***  WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
!***  NEAREST TO THE CENTER K AS SEEN BELOW.  WE MUST FIND
!***  WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS
!***  A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION.
!***
!**
!***                  H
!***
!***
!***
!***            H     V     H
!***
!***
!***               H
!***      H     V     H     V     H
!***
!***
!***
!***            H     V     H
!***
!***
!***
!***                  H
!***
!***
!***  FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H.
!***
    N2R=K/IMT
    MK=MOD(K,IMT)
!
    IF(MK.EQ.0)THEN
      TLATHC=SB+(2*N2R-1)*DPH
    ELSE
      TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH
    ENDIF
!
    IF(MK.LE.(P_IDE-1))THEN
      TLONHC=WB+2*(MK-1)*DLM
    ELSE
      TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM
    ENDIF
  
!
!***  EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
!***  DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE
!***  CAREFUL HERE      
!

    IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON
    IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT
    DENOM=(TLON-TLONHC)
!
!***
!***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
!***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
!***
!*** COINCIDENT CONDITIONS

    IF(DENOM.EQ.0.0)THEN

      IF(TLATHC.EQ.TLAT)THEN
        KOUTB(I,J)=K
        IIH(I,J) = NCOL
        JJH(I,J) = NROW
        TLATHX(I,J)=TLATHC
        TLONHX(I,J)=TLONHC
        HBWGT1(I,J)=1.0
        HBWGT2(I,J)=0.0
        HBWGT3(I,J)=0.0
        HBWGT4(I,J)=0.0
      ELSE                                      ! SAME LONGITUDE BUT DIFFERENT LATS
!
         IF(TLATHC .GT. TLAT)THEN      ! NESTED POINT SOUTH OF PARENT
          KOUTB(I,J)=K-(P_IDE-1)
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW-1
          TLATHX(I,J)=TLATHC-DPH
          TLONHX(I,J)=TLONHC-DLM
         ELSE                                   ! NESTED POINT NORTH OF PARENT
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW+1
          TLATHX(I,J)=TLATHC+DPH
          TLONHX(I,J)=TLONHC-DLM
         ENDIF
!***
!***
!***                  4
!***
!***                  H
!***             1         2
!***
!***                  3
!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX

       TLATO=TLATHX(I,J)
       TLONO=TLONHX(I,J)
       DLM1=TLON-TLONO
       DLA1=TLAT-TLATO                                               ! Q
!      DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q 
!
       TLATO=TLATHX(I,J)
       TLONO=TLONHX(I,J)+2.*DLM
       DLM2=TLON-TLONO
       DLA2=TLAT-TLATO                                               ! Q
!      DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
! 
       TLATO=TLATHX(I,J)-DPH
       TLONO=TLONHX(I,J)+DLM
       DLM3=TLON-TLONO
       DLA3=TLAT-TLATO                                               ! Q
!      DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
 
       TLATO=TLATHX(I,J)+DPH
       TLONO=TLONHX(I,J)+DLM
       DLM4=TLON-TLONO
       DLA4=TLAT-TLATO                                               ! Q
!      DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q


!      THE BILINEAR WEIGHTS
!***
!***
       AN3=ATAN2(DLA1,DLM1)                                          ! Q
       R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
       S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
       R1=R1/DS1
       S1=S1/DS1
       DL1I=(1.-R1)*(1.-S1)
       DL2I=R1*S1
       DL3I=R1*(1.-S1)
       DL4I=(1.-R1)*S1
!
       HBWGT1(I,J)=DL1I
       HBWGT2(I,J)=DL2I
       HBWGT3(I,J)=DL3I
       HBWGT4(I,J)=DL4I
! 
      ENDIF

    ELSE
!
!*** NON-COINCIDENT POINTS   
!
      SLOPE=(TLAT-TLATHC)/DENOM
      DSLOPE=NINT(R2D*ATAN(SLOPE))

      IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
        IF(TLON.GT.TLONHC)THEN
          IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K
          IIH(I,J) = NCOL
          JJH(I,J) = NROW
          TLATHX(I,J)=TLATHC
          TLONHX(I,J)=TLONHC
        ELSE
          IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-1
          IIH(I,J) = NCOL-2
          JJH(I,J) = NROW
          TLATHX(I,J)=TLATHC
          TLONHX(I,J)=TLONHC -2.*DLM
        ENDIF

!
      ELSEIF(DSLOPE.GT.DSLP0)THEN
        IF(TLON.GT.TLONHC)THEN
          IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW+1
          TLATHX(I,J)=TLATHC+DPH
          TLONHX(I,J)=TLONHC-DLM
        ELSE
          IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-(P_IDE-1)
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW-1
          TLATHX(I,J)=TLATHC-DPH
          TLONHX(I,J)=TLONHC-DLM
        ENDIF

!
      ELSEIF(DSLOPE.LT.-DSLP0)THEN
        IF(TLON.GT.TLONHC)THEN
          IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-(P_IDE-1)
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW-1
          TLATHX(I,J)=TLATHC-DPH
          TLONHX(I,J)=TLONHC-DLM
        ELSE
          IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIH(I,J) = NCOL-1
          JJH(I,J) = NROW+1
          TLATHX(I,J)=TLATHC+DPH
          TLONHX(I,J)=TLONHC-DLM
        ENDIF
      ENDIF

!
!***  NOW WE WILL MOVE AS FOLLOWS:
!***
!***
!***                      4
!***
!***
!***  
!***                   H
!***             1                 2
!***
!***
!***
!***
!***                      3
!***
!***
!***
!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
      
      TLATO=TLATHX(I,J)
      TLONO=TLONHX(I,J)
      DLM1=TLON-TLONO
      DLA1=TLAT-TLATO                                               ! Q
!     DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
!
      TLATO=TLATHX(I,J)                                             ! redundant computations
      TLONO=TLONHX(I,J)+2.*DLM
      DLM2=TLON-TLONO
      DLA2=TLAT-TLATO                                               ! Q
!     DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
!
      TLATO=TLATHX(I,J)-DPH
      TLONO=TLONHX(I,J)+DLM
      DLM3=TLON-TLONO
      DLA3=TLAT-TLATO                                               ! Q
!     DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q
!
      TLATO=TLATHX(I,J)+DPH
      TLONO=TLONHX(I,J)+DLM
      DLM4=TLON-TLONO
      DLA4=TLAT-TLATO                                               ! Q
!     DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q

!     THE BILINEAR WEIGHTS
!***
      AN3=ATAN2(DLA1,DLM1)                                          ! Q
      R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
      S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
      R1=R1/DS1
      S1=S1/DS1
      DL1I=(1.-R1)*(1.-S1)
      DL2I=R1*S1
      DL3I=R1*(1.-S1)
      DL4I=(1.-R1)*S1
!
      HBWGT1(I,J)=DL1I
      HBWGT2(I,J)=DL2I
      HBWGT3(I,J)=DL3I
      HBWGT4(I,J)=DL4I
!
    ENDIF 

!
!***  FINALLY STORE IIH IN TERMS OF E-GRID INDEX
!
      IIH(I,J)=NINT(0.5*IIH(I,J))

      HBWGT1(I,J)=MAX(HBWGT1(I,J),0.0)   ! all weights must be GE zero (postive def)
      HBWGT2(I,J)=MAX(HBWGT2(I,J),0.0)   ! all weights must be GE zero (postive def)
      HBWGT3(I,J)=MAX(HBWGT3(I,J),0.0)   ! all weights must be GE zero (postive def)
      HBWGT4(I,J)=MAX(HBWGT4(I,J),0.0)   ! all weights must be GE zero (postive def)

!      write(message,105)"H WEIGHTS:",I,J,HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J), &
!                               HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J),IIH(i,j),JJH(i,j)
!      CALL wrf_message(trim(message))
! 105  format(a,2i4,5f7.3,2i4)

   ENDDO
  ENDDO


  RETURN 
  END SUBROUTINE G2T2H
!========================================================================================


  SUBROUTINE G2T2V( IIV,JJV,                     & ! output grid index and weights
                    VBWGT1,VBWGT2,               & ! output weights in terms of parent grid
                    VBWGT3,VBWGT4,               &
                    VLAT,VLON,                   & ! target (nest) input lat lon in degrees
                    DLMD1,DPHD1,WBD1,SBD1,       & ! parent res, west and south boundaries
                    CENTRAL_LAT,CENTRAL_LON,     & ! parent central lat,lon, all in degrees
                    P_IDE,P_JDE,                 & ! parent imax and jmax
                    IDS,IDE,JDS,JDE,KDS,KDE,     & ! target (nest) dIMEnsions
                    IMS,IME,JMS,JME,KMS,KME,     &
                    ITS,ITE,JTS,JTE,KTS,KTE      )

!
!***  Tom Black   - Initial Version 
!***  Gopal       - Revised Version for WRF (includes coincIDEnt grid points)
!***
!***  GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
!***  AND THE NESTED GRID LAT-LONS AT V POINTS, THIS ROUTINE FIRST LOCATES THE
!***  INDICES,IIV,JJV, OF THE PARENT DOMAIN'S V POINTS THAT LIES CLOSEST TO THE
!***  V POINTS OF THE NESTED DOMAIN
!
!============================================================================

 IMPLICIT NONE
 INTEGER,    INTENT(IN   )                            :: IDS,IDE,JDS,JDE,KDS,KDE
 INTEGER,    INTENT(IN   )                            :: IMS,IME,JMS,JME,KMS,KME
 INTEGER,    INTENT(IN   )                            :: ITS,ITE,JTS,JTE,KTS,KTE
 INTEGER,    INTENT(IN   )                            :: P_IDE,P_JDE
 REAL,       INTENT(IN   )                            :: DLMD1,DPHD1,WBD1,SBD1
 REAL,       INTENT(IN   )                            :: CENTRAL_LAT,CENTRAL_LON
 REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(IN)    :: VLAT,VLON
 REAL,    DIMENSION(IMS:IME,JMS:JME),   INTENT(OUT)   :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
 INTEGER, DIMENSION(IMS:IME,JMS:JME),   INTENT(OUT)   :: IIV,JJV

! local

 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)     ! (6) single precision
 INTEGER           :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE
 INTEGER           :: NROW,NCOL,KROWS
 REAL(KIND=KNUM)   :: X,Y,Z,TLAT,TLON
 REAL(KIND=KNUM)   :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
 REAL(KIND=KNUM)   :: ROW,COL,SLP0,TLATVC,TLONVC,DENOM,SLOPE
 REAL(KIND=KNUM)   :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
 REAL(KIND=KNUM)   :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3                    ! Q
 REAL(KIND=KNUM)   :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
 REAL(KIND=KNUM)   :: DTEMP
 REAL  , DIMENSION(IMS:IME,JMS:JME)      :: TLATVX,TLONVX
 INTEGER, DIMENSION(IMS:IME,JMS:JME)     :: KOUTB
!-------------------------------------------------------------------------------------

  IMT=2*P_IDE-2             ! parent i dIMEnsions
  JMT=P_JDE/2               ! parent j dIMEnsions
  PI_2=ACOS(0.)
  D2R=PI_2/90.
  R2D=1./D2R
  DPH=DPHD1*D2R
  DLM=DLMD1*D2R
  TPH0= CENTRAL_LAT*D2R
  TLM0=-CENTRAL_LON*D2R        ! NOTE THE MINUS HERE
  WB=WBD1*D2R                   ! DEGREES TO RADIANS
  SB=SBD1*D2R
  SLP0=DPHD1/DLMD1
  DSLP0=NINT(R2D*ATAN(SLP0))
  DS1=SQRT(DPH*DPH+DLM*DLM)    ! Q
  AN1=ASIN(DLM/DS1)
  AN2=ASIN(DPH/DS1)

  DO J =  JTS,MIN(JTE,JDE-1)
    DO I = ITS,MIN(ITE,IDE-1)
!***
!***  LOCATE TARGET V POINTS (VLAT AND VLON) ON THE PARENT DOMAIN AND
!***  DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
!***  CONVERT NESTED GRID V POINTS FROM GEODETIC TO TRANSFORMED
!***  COORDINATE ON THE PARENT GRID
!

      GLAT=VLAT(I,J)*D2R
      GLON=(360. - VLON(I,J))*D2R
      X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
      Y=-COS(GLAT)*SIN(GLON-TLM0)
      Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
      TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
      TLON=R2D*ATAN(Y/X)

!      ROW=TLAT/DPHD1+JMT         ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
!      COL=TLON/DLMD1+P_IDE-1     ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN

      ROW=(TLAT-SBD1)/DPHD1+1     ! Dusan's doing
      COL=(TLON-WBD1)/DLMD1+1     ! Dusan's doing

      NROW=INT(ROW + 0.001)     ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
      NCOL=INT(COL + 0.001)
      TLAT=TLAT*D2R
      TLON=TLON*D2R

!***
!***
!***  FIRST CONSIDER THE SITUATION WHERE THE POINT V IS AT
!***
!***              H      V
!***
!***
!***                 V
!***              V      H
!***
!***  THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
!***

      IF(MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.1.OR.     &
         MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.0)THEN
           TLAT1=(NROW-JMT)*DPH
           TLAT2=TLAT1+DPH
           TLON1=(NCOL-(P_IDE-1))*DLM
           TLON2=TLON1+DLM
           DLM1=TLON-TLON1
           DLM2=TLON-TLON2
!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
           D1=ACOS(DTEMP)
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           D2=ACOS(DTEMP)
            IF(D1.GT.D2)THEN
             NROW=NROW+1                    ! FIND THE NEAREST V ROW
             NCOL=NCOL+1                    ! FIND THE NEAREST V COLUMN
            ENDIF
      ELSE

!***
!***  NOW CONSIDER THE SITUATION WHERE THE POINT V IS AT
!***
!***              V      H
!***
!***
!***                 V
!***              H      V
!***
!*** THEN LOCATE THE NEAREST V POINT ON THE PARENT GRID
!***
           TLAT1=(NROW+1-JMT)*DPH
           TLAT2=TLAT1-DPH
           TLON1=(NCOL-(P_IDE-1))*DLM
           TLON2=TLON1+DLM
           DLM1=TLON-TLON1
           DLM2=TLON-TLON2
!           D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
!           D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
           D1=ACOS(DTEMP)
           DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
           D2=ACOS(DTEMP)
             IF(D1.LT.D2)THEN
              NROW=NROW+1                    ! FIND THE NEAREST H ROW
             ELSE
              NCOL=NCOL+1                    ! FIND THE NEAREST H COLUMN
             ENDIF

      ENDIF

      KROWS=((NROW-1)/2)*IMT
      IF(MOD(NROW,2).EQ.1)THEN
        K=KROWS+NCOL/2
      ELSE
        K=KROWS+P_IDE-2+(NCOL+1)/2     ! check this one should this not be P_IDE-2 ????
      ENDIF

!***
!***  WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
!***  NEAREST TO THE CENTER K AS SEEN BELOW.  WE MUST FIND
!***  WHICH OF THE FOUR v-BOXES (OF WHICH THIS V POINT IS
!***  A VERTEX) SURROUNDS THE INNER GRID V POINT IN QUESTION.
!***
!***
!***                  V
!***
!***
!***
!***            V     H     V
!***
!***
!***               V
!***      V     H     V     H     V
!***
!***
!***
!***            V     H     V
!***
!***
!***
!***                  V
!***
!***
!***  FIND THE SLOPE OF THE LINE CONNECTING V AND THE CENTER v.
!***
      N2R=K/IMT
      MK=MOD(K,IMT)
!
      IF(MK.EQ.0)THEN
        TLATVC=SB+(2*N2R-1)*DPH
      ELSE
        TLATVC=SB+(2*N2R+MK/(P_IDE-1))*DPH
      ENDIF
!
      IF(MK.LE.(P_IDE-1)-1)THEN
        TLONVC=WB+(2*MK-1)*DLM
      ELSE
        TLONVC=WB+2*(MK-(P_IDE-1))*DLM
      ENDIF

!
!***  EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
!***  DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE V BOXES, WE NEED TO BE
!***  CAREFUL HERE
!
       IF(ABS(TLON-TLONVC) .LE. 1.E-4)TLONVC=TLON
       IF(ABS(TLAT-TLATVC) .LE. 1.E-4)TLATVC=TLAT
       DENOM=(TLON-TLONVC)
!
!***
!***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
!***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
!***
!*** COINCIDENT CONDITIONS

     IF(DENOM.EQ.0.0)THEN

       IF(TLATVC.EQ.TLAT)THEN
         KOUTB(I,J)=K
         IIV(I,J) = NCOL
         JJV(I,J) = NROW
         TLATVX(I,J)=TLATVC
         TLONVX(I,J)=TLONVC
         VBWGT1(I,J)=1.0
         VBWGT2(I,J)=0.0
         VBWGT3(I,J)=0.0
         VBWGT4(I,J)=0.0
       ELSE                              ! SAME LONGITUDE BUT DIFFERENT LATS 
          
         IF(TLATVC .GT. TLAT)THEN      ! NESTED POINT SOUTH OF PARENT
          KOUTB(I,J)=K-(P_IDE-1)
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW-1
          TLATVX(I,J)=TLATVC-DPH
          TLONVX(I,J)=TLONVC-DLM
         ELSE                                   ! NESTED POINT NORTH OF PARENT
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW+1
          TLATVX(I,J)=TLATVC+DPH
          TLONVX(I,J)=TLONVC-DLM
         ENDIF

!***
!***
!***                  4
!***
!***                  V 
!***             1         2
!***
!***                  3
!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX

       TLATO=TLATVX(I,J)
       TLONO=TLONVX(I,J)
       DLM1=TLON-TLONO
       DLA1=TLAT-TLATO                                               ! Q
!      DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
!
       TLATO=TLATVX(I,J)
       TLONO=TLONVX(I,J)+2.*DLM
       DLM2=TLON-TLONO
       DLA2=TLAT-TLATO                                               ! Q
!      DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q

       TLATO=TLATVX(I,J)-DPH
       TLONO=TLONVX(I,J)+DLM
       DLM3=TLON-TLONO
       DLA3=TLAT-TLATO                                               ! Q
!      DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL3=SQRT(DLM3*DLM3+DLA3*DLA3)                                 ! Q

       TLATO=TLATVX(I,J)+DPH
       TLONO=TLONVX(I,J)+DLM
       DLM4=TLON-TLONO
       DLA4=TLAT-TLATO                                               ! Q
!      DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
       DL4=SQRT(DLM4*DLM4+DLA4*DLA4)                                 ! Q
                 
!      THE BILINEAR WEIGHTS
!***
       AN3=ATAN2(DLA1,DLM1)                                          ! Q
       R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
       S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
       R1=R1/DS1
       S1=S1/DS1
       DL1I=(1.-R1)*(1.-S1)
       DL2I=R1*S1
       DL3I=R1*(1.-S1)
       DL4I=(1.-R1)*S1
!
       VBWGT1(I,J)=DL1I
       VBWGT2(I,J)=DL2I
       VBWGT3(I,J)=DL3I
       VBWGT4(I,J)=DL4I      

      ENDIF

    ELSE

!
!*** NON-COINCIDENT POINTS
!
      SLOPE=(TLAT-TLATVC)/DENOM
      DSLOPE=NINT(R2D*ATAN(SLOPE))

      IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
        IF(TLON.GT.TLONVC)THEN
          IF(TLONVC.GE.-WB-DLM)CALL wrf_error_fatal("1V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K
          IIV(I,J)=NCOL
          JJV(I,J)=NROW
          TLATVX(I,J)=TLATVC
          TLONVX(I,J)=TLONVC
        ELSE
          IF(TLONVC.LE.WB+DLM)CALL wrf_error_fatal("2V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-1
          IIV(I,J) = NCOL-2
          JJV(I,J) = NROW
          TLATVX(I,J)=TLATVC
          TLONVX(I,J)=TLONVC-2.*DLM
        ENDIF
 
      ELSEIF(DSLOPE.GT.DSLP0)THEN
        IF(TLON.GT.TLONVC)THEN
          IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("3V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW+1
          TLATVX(I,J)=TLATVC+DPH
          TLONVX(I,J)=TLONVC-DLM
        ELSE
          IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("4V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-(P_IDE-1)
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW-1
          TLATVX(I,J)=TLATVC-DPH
          TLONVX(I,J)=TLONVC-DLM
        ENDIF
 
      ELSEIF(DSLOPE.LT.-DSLP0)THEN
        IF(TLON.GT.TLONVC)THEN
          IF(TLATVC.LE.SB+DPH)CALL wrf_error_fatal("5V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K-(P_IDE-1)
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW-1
          TLATVX(I,J)=TLATVC-DPH
          TLONVX(I,J)=TLONVC-DLM
        ELSE
          IF(TLATVC.GE.-SB-DPH)CALL wrf_error_fatal("6V:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
          KOUTB(I,J)=K+(P_IDE-1)-1
          IIV(I,J) = NCOL-1
          JJV(I,J) = NROW+1
          TLATVX(I,J)=TLATVC+DPH
          TLONVX(I,J)=TLONVC-DLM
        ENDIF
      ENDIF
!
!***  NOW WE WILL MOVE AS FOLLOWS:
!***
!***
!***                      4
!***
!***
!*** 
!***                   V 
!***             1                 2
!***
!***
!***
!***
!***                      3
!***
!***
!***
!***  DL 1-4 ARE THE ANGULAR DISTANCES FROM V TO EACH VERTEX

      TLATO=TLATVX(I,J)
      TLONO=TLONVX(I,J)
      DLM1=TLON-TLONO
      DLA1=TLAT-TLATO                                               ! Q
!     DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL1=SQRT(DLM1*DLM1+DLA1*DLA1)                                 ! Q
!
      TLATO=TLATVX(I,J)
      TLONO=TLONVX(I,J)+2.*DLM
      DLM2=TLON-TLONO
      DLA2=TLAT-TLATO                                               ! Q
!     DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
      DL2=SQRT(DLM2*DLM2+DLA2*DLA2)                                 ! Q
!
      TLATO…