/wrfv2_fire/dyn_nmm/NMM_NEST_UTILS1.F
FORTRAN Legacy | 3802 lines | 2481 code | 422 blank | 899 comment | 16 complexity | 9b1eb832c545956ed2b6043efc59e2ba MD5 | raw file
Possible License(s): AGPL-1.0
- #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=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
- !
- !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX
- !
- IIV(I,J)=NINT(0.5*IIV(I,J))
- VBWGT1(I,J)=MAX(VBWGT1(I,J),0.0) ! all weights must be GE zero (postive def)
- VBWGT2(I,J)=MAX(VBWGT2(I,J),0.0) ! all weights must be GE zero (postive def)
- VBWGT3(I,J)=MAX(VBWGT3(I,J),0.0) ! all weights must be GE zero (postive def)
- VBWGT4(I,J)=MAX(VBWGT4(I,J),0.0) ! all weights must be GE zero (postive def)
- ENDDO
- ENDDO
- RETURN
- END SUBROUTINE G2T2V
- !-----------------------------------------------------------------------------
- SUBROUTINE G2T2H_new( IIH,JJH, & ! output grid index and weights
- HBWGT1,HBWGT2, & ! output weights in terms of parent grid
- HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- !
- !*** XUEJIN ZHANG --- Initial version (09/08/2008)
- !*** XUEJIN ZHANG --- Modified for parallel purpose (09/10/2009)
- !*** XUEJIN ZHANG --- Modified for parallel purpose (09/29/2009)
- !Function: Bilnear interpolation weight and indexing for E-grid
- !
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: JP
- !*** ARRANGEMENT OF 4 VERTEXES FROM PARENT DOMAIN
- !***
- !*** 4
- !***
- !*** h
- !*** 1 2
- !***
- !***
- !*** 3
- !
- !*************************************************************
- !*** POINT (i,j) SPANS 9 NESTED POINTS
- !*** A VERTEX IN THE NEST DOMAIN AND INDEXING FOR PARENT DOMAIN
- !***
- !***
- !*** H
- !***
- !*** h h
- !*** / \
- !*** h h h
- !*** / \ / \
- !*** H h h H
- !*** \ / \ /
- !*** h h h
- !*** \ /
- !*** h h
- !***
- !*** H
- !***
- !***
- !***
- !*************************************************************
- !*** MOVING NEST ALWAYS STARTING FROM MASS GRID H
- !*** PLEASE REFER TO E-GRID ARRANGEMENT FIGURE FOR MASS POINTS
- DO J=JTS,MIN(JTE,JDE-1)
- DO I=ITS,MIN(ITE,IDE-1)
- JP = MOD(J,RATIO*2)
- SELECT CASE(JP)
- CASE ( 1 )
- CALL SUB1H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 2 )
- CALL SUB2H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 3 )
- CALL SUB3H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 4 )
- CALL SUB4H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 5 )
- CALL SUB5H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 0 )
- CALL SUB6H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- END SELECT
- #if defined(EXPENSIVE_HWRF_DEBUG_STUFF)
- 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), &
- HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J)
- #endif
- 105 format(a,4i4,5f7.3)
- END DO
- END DO
- RETURN
- END SUBROUTINE G2T2H_new
- SUBROUTINE SUB1H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 0.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 4.0/9.0
- HBWGT2(I,J) = 1.0/9.0
- HBWGT3(I,J) = 2.0/9.0
- HBWGT4(I,J) = 2.0/9.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0/9.0
- HBWGT2(I,J) = 4.0/9.0
- HBWGT3(I,J) = 2.0/9.0
- HBWGT4(I,J) = 2.0/9.0
- END SELECT
- RETURN
- END SUBROUTINE SUB1H
- SUBROUTINE SUB2H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 2.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 1.0/3.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 2.0/9.0
- HBWGT2(I,J) = 2.0/9.0
- HBWGT3(I,J) = 1.0/9.0
- HBWGT4(I,J) = 4.0/9.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)+1
- HBWGT1(I,J) = 1.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 2.0/3.0
- HBWGT4(I,J) = 0.0
- END SELECT
- RETURN
- END SUBROUTINE SUB2H
- SUBROUTINE SUB3H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)-1
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)+1
- HBWGT1(I,J) = 2.0/9.0
- HBWGT2(I,J) = 2.0/9.0
- HBWGT3(I,J) = 4.0/9.0
- HBWGT4(I,J) = 1.0/9.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 2.0/3.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)+1
- HBWGT1(I,J) = 2.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 1.0/3.0
- HBWGT4(I,J) = 0.0
- END SELECT
- RETURN
- END SUBROUTINE SUB3H
- SUBROUTINE SUB4H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)-1
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0/9.0
- HBWGT2(I,J) = 4.0/9.0
- HBWGT3(I,J) = 2.0/9.0
- HBWGT4(I,J) = 2.0/9.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 0.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 4.0/9.0
- HBWGT2(I,J) = 1.0/9.0
- HBWGT3(I,J) = 2.0/9.0
- HBWGT4(I,J) = 2.0/9.0
- END SELECT
- RETURN
- END SUBROUTINE SUB4H
- SUBROUTINE SUB5H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)-1
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 2.0/9.0
- HBWGT2(I,J) = 2.0/9.0
- HBWGT3(I,J) = 1.0/9.0
- HBWGT4(I,J) = 4.0/9.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)+1
- HBWGT1(I,J) = 1.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 2.0/3.0
- HBWGT4(I,J) = 0.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 2.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 1.0/3.0
- END SELECT
- RETURN
- END SUBROUTINE SUB5H
- SUBROUTINE SUB6H(I,J,IIH,JJH,HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIH,JJH
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- HBWGT1(I,J) = 2.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 1.0/3.0
- HBWGT4(I,J) = 0.0
- CASE ( 2 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- HBWGT1(I,J) = 2.0/9.0
- HBWGT2(I,J) = 2.0/9.0
- HBWGT3(I,J) = 4.0/9.0
- HBWGT4(I,J) = 1.0/9.0
- CASE ( 0 )
- IIH(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJH(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- HBWGT1(I,J) = 1.0/3.0
- HBWGT2(I,J) = 0.0
- HBWGT3(I,J) = 0.0
- HBWGT4(I,J) = 2.0/3.0
- END SELECT
- RETURN
- END SUBROUTINE SUB6H
- !-----------------------------------------------------------------------------
- SUBROUTINE G2T2V_new( IIV,JJV, & ! output grid index and weights
- VBWGT1,VBWGT2, & ! output weights in terms of parent grid
- VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- !
- !*** XUEJIN ZHANG --- Initial version (09/08/2008)
- !*** XUEJIN ZHANG --- Modified for parallel purpose (09/10/2009)
- !*** XUEJIN ZHANG --- Modified for parallel purpose (09/29/2009)
- !Function: Bilnear interpolation weight and indexing for E-grid
- !
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: JP
- !*** ARRANGEMENT OF 4 VERTEXES FROM PARENT DOMAIN
- !***
- !*** 4
- !***
- !*** v
- !*** 1 2
- !***
- !***
- !*** 3
- !
- !*************************************************************
- !*** POINT (i,j) SPANS 9 NESTED POINTS
- !*** A VERTEX IN THE NEST DOMAIN AND INDEXING FOR PARENT DOMAIN
- !***
- !***
- !*** V
- !***
- !*** v h v
- !*** / \
- !*** v h v h v
- !*** / \ / \
- !*** V H v h v h V H
- !*** \ / \ /
- !*** v h v h v
- !*** \ /
- !*** v h v
- !***
- !*** V
- !***
- !***
- !***
- !*************************************************************
- !*** MOVING NEST ALWAYS STARTING FROM MASS GRID H
- !*** PLEASE REFER TO E-GRID ARRANGEMENT FIGURE FOR WIND POINTS
- DO J=JTS,MIN(JTE,JDE-1)
- DO I=ITS,MIN(ITE,IDE-1)
- JP = MOD(J,RATIO*2)
- SELECT CASE(JP)
- CASE ( 1 )
- CALL SUB1V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 2 )
- CALL SUB2V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 3 )
- CALL SUB3V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 4 )
- CALL SUB4V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 5 )
- CALL SUB5V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- CASE ( 0 )
- CALL SUB6V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- END SELECT
- ! WRITE(0,105)'NEW V WEIGHTS:',I,J,VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J), &
- ! VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J),IIV(i,j),JJV(i,j)
- ! 105 format(a,2i4,5f7.3,2i4)
- END DO
- END DO
- RETURN
- END SUBROUTINE G2T2V_new
- SUBROUTINE SUB1V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO) - 1
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0/9.0
- VBWGT2(I,J) = 4.0/9.0
- VBWGT3(I,J) = 2.0/9.0
- VBWGT4(I,J) = 2.0/9.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 0.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 4.0/9.0
- VBWGT2(I,J) = 1.0/9.0
- VBWGT3(I,J) = 2.0/9.0
- VBWGT4(I,J) = 2.0/9.0
- END SELECT
- RETURN
- END SUBROUTINE SUB1V
- SUBROUTINE SUB2V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO) - 1
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 2.0/9.0
- VBWGT2(I,J) = 2.0/9.0
- VBWGT3(I,J) = 1.0/9.0
- VBWGT4(I,J) = 4.0/9.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 1.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 2.0/3.0
- VBWGT4(I,J) = 0.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 2.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 1.0/3.0
- END SELECT
- RETURN
- END SUBROUTINE SUB2V
- SUBROUTINE SUB3V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 2.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 1.0/3.0
- VBWGT4(I,J) = 0.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 2.0/9.0
- VBWGT2(I,J) = 2.0/9.0
- VBWGT3(I,J) = 4.0/9.0
- VBWGT4(I,J) = 1.0/9.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 2.0/3.0
- END SELECT
- RETURN
- END SUBROUTINE SUB3V
- SUBROUTINE SUB4V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 0.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 4.0/9.0
- VBWGT2(I,J) = 1.0/9.0
- VBWGT3(I,J) = 2.0/9.0
- VBWGT4(I,J) = 2.0/9.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0/9.0
- VBWGT2(I,J) = 4.0/9.0
- VBWGT3(I,J) = 2.0/9.0
- VBWGT4(I,J) = 2.0/9.0
- END SELECT
- RETURN
- END SUBROUTINE SUB4V
- SUBROUTINE SUB5V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 2.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 1.0/3.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 2.0/9.0
- VBWGT2(I,J) = 2.0/9.0
- VBWGT3(I,J) = 1.0/9.0
- VBWGT4(I,J) = 4.0/9.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 1.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 2.0/3.0
- VBWGT4(I,J) = 0.0
- END SELECT
- RETURN
- END SUBROUTINE SUB5V
- SUBROUTINE SUB6V(I,J,IIV,JJV,VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- I_PARENT_START,J_PARENT_START, & ! nest start I and J in parent domain
- RATIO, & ! 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 )
- 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 ) :: I_PARENT_START,J_PARENT_START
- INTEGER, INTENT(IN ) :: RATIO
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- INTEGER, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: IIV,JJV
- !LOCAL VARIABLES
- INTEGER :: I,J
- INTEGER :: IP
- IP = MOD(I,RATIO)
- SELECT CASE(IP)
- CASE ( 1 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO) - 1
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 2.0/9.0
- VBWGT2(I,J) = 2.0/9.0
- VBWGT3(I,J) = 4.0/9.0
- VBWGT4(I,J) = 1.0/9.0
- CASE ( 2 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO)
- VBWGT1(I,J) = 1.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 0.0
- VBWGT4(I,J) = 2.0/3.0
- CASE ( 0 )
- IIV(I,J) = I_PARENT_START + INT((I-1)/RATIO)
- JJV(I,J) = J_PARENT_START + INT((J-1)/RATIO) + 1
- VBWGT1(I,J) = 2.0/3.0
- VBWGT2(I,J) = 0.0
- VBWGT3(I,J) = 1.0/3.0
- VBWGT4(I,J) = 0.0
- END SELECT
- RETURN
- END SUBROUTINE SUB6V
- !------------------------------------------------------------------------------
- !
- SUBROUTINE WEIGTS_CHECK(HBWGT1,HBWGT2,HBWGT3,HBWGT4, &
- VBWGT1,VBWGT2,VBWGT3,VBWGT4, &
- 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, &
- IMS,IME,JMS,JME,KMS,KME, &
- ITS,ITE,JTS,JTE,KTS,KTE
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
- ! local
- REAL , PARAMETER :: EPSI=1.0E-3
- INTEGER :: I,J
- REAL :: ADDSUM
- CHARACTER(LEN=255):: message
- !-------------------------------------------------------------------------------------
- ! DUE TO THE NEED FOR HALO EXCHANGES IN PARALLEL RUNS ONE HAS TO ENSURE CONSISTENT
- ! USAGE OF NUMBER OF PROCESSORS BEFORE ANY FURTHER COMPUTATIONS. WE INTRODUCE THIS
- ! CHECK FIRST
- IF((ITE-ITS) .LE. 5 .OR. (JTE-JTS) .LE. 5)THEN
- WRITE(message,*)'ITE-ITS=',ITE-ITS,'JTE-JTS=',JTE-JTS
- CALL wrf_message(trim(message))
- CALL wrf_error_fatal ('NESTED DOMAIN:PLEASE OPTIMIZE THE NUMBER OF PROCESSES; TRY SQUARES OF NUMBERS')
- ENDIF
- !
- ! NOW CHECK WEIGHTS
- !
- ADDSUM=0.
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- ADDSUM=HBWGT1(I,J)+HBWGT2(I,J)+HBWGT3(I,J)+HBWGT4(I,J)
- IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
- WRITE(message,*)'I=',I,'J=',J,'WEIGHTS=',HBWGT1(I,J),HBWGT2(I,J),HBWGT3(I,J),HBWGT4(I,J),1-ADDSUM
- CALL wrf_message(trim(message))
- CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT MASS POINTS')
- ENDIF
- ENDDO
- ENDDO
- ADDSUM=0.
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- ADDSUM=VBWGT1(I,J)+VBWGT2(I,J)+VBWGT3(I,J)+VBWGT4(I,J)
- IF(ABS(1.0-ADDSUM) .GE. EPSI)THEN
- WRITE(message,*)'I=',I,'J=',J,'WEIGHTS=',VBWGT1(I,J),VBWGT2(I,J),VBWGT3(I,J),VBWGT4(I,J),1-ADDSUM
- CALL wrf_message(trim(message))
- CALL wrf_error_fatal ('NESTED DOMAIN:SOMETHING IS WRONG WITH WEIGHTS COMPUTATION AT VELOCITY POINTS')
- ENDIF
- ENDDO
- ENDDO
- END SUBROUTINE WEIGTS_CHECK
- !-----------------------------------------------------------------------------------
- SUBROUTINE BOUNDS_CHECK( IIH,JJH,IIV,JJV, &
- IPOS,JPOS,SHW, &
- IDS,IDE,JDS,JDE,KDS,KDE, & !
- IMS,IME,JMS,JME,KMS,KME, & ! nested grid configuration
- ITS,ITE,JTS,JTE,KTS,KTE )
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: IPOS,JPOS,SHW, &
- IDS,IDE,JDS,JDE,KDS,KDE, &
- IMS,IME,JMS,JME,KMS,KME, &
- ITS,ITE,JTS,JTE,KTS,KTE
- INTEGER, DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IIH,JJH,IIV,JJV
- ! local variables
- INTEGER :: I,J
- CHARACTER(LEN=255) :: message
- !*** Gopal - Initial version
- !***
- !*** CHECK DOMAIN BOUNDS BEFORE PROCEEDING TO INTERPOLATION
- !
- !============================================================================
- IF(IPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs X-BOUNDARY')
- IF(JPOS .LE. SHW)CALL wrf_error_fatal('NESTED DOMAIN TOO CLOSE TO PARENTs Y-BOUNDARY')
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF(IIH(I,J) .EQ. 0)CALL wrf_error_fatal ('IIH=0: SOMETHING IS WRONG')
- IF(JJH(I,J) .EQ. 0)CALL wrf_error_fatal ('JJH=0: SOMETHING IS WRONG')
- ENDDO
- ENDDO
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF(IIH(I,J) .LT. (IPOS-SHW) .OR. JJH(I,J) .LT. (JPOS-SHW) .OR. &
- IIV(I,J) .LT. (IPOS-SHW) .OR. JJV(I,J) .LT. (JPOS-SHW))THEN
- WRITE(message,*)I,J,IIH(I,J),IPOS,JJH(I,J),JPOS,SHW
- CALL wrf_message(trim(message))
- WRITE(message,*)I,J,IIV(I,J),IPOS,JJV(I,J),JPOS,SHW
- CALL wrf_message(trim(message))
- CALL wrf_error_fatal ('CHECK NESTED DOMAIN BOUNDS: TRY INCREASING STENCIL WIDTH')
- ENDIF
- ENDDO
- ENDDO
- END SUBROUTINE BOUNDS_CHECK
- !==========================================================================================
- SUBROUTINE BASE_STATE_PARENT ( Z3d,Q3d,T3d,PSTD, &
- PINT,T,Q,CWM, &
- FIS,QS,PD,PDTOP,PTOP, &
- ETA1,ETA2, &
- DETA1,DETA2, &
- IDS,IDE,JDS,JDE,KDS,KDE, &
- IMS,IME,JMS,JME,KMS,KME, &
- ITS,ITE,JTS,JTE,KTS,KTE )
- !
- USE MODULE_MODEL_CONSTANTS
- 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 ) :: PDTOP,PTOP
- REAL, DIMENSION(KMS:KME), INTENT(IN) :: ETA1,ETA2,DETA1,DETA2
- REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: FIS,PD,QS
- REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME), INTENT(IN) :: PINT,T,Q,CWM
- REAL, DIMENSION(KMS:KME), INTENT(OUT):: PSTD
- REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME), INTENT(OUT):: Z3d,Q3d,T3d
- ! local
- INTEGER,PARAMETER :: JTB=134
- INTEGER :: I,J,K,ILOC,JLOC
- REAL, PARAMETER :: LAPSR=6.5E-3, GI=1./G,D608=0.608
- REAL, PARAMETER :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
- REAL, PARAMETER :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
- REAL, PARAMETER :: P_REF=103000.
- REAL :: A,B,APELP,RTOPP,DZ,ZMID
- REAL, DIMENSION(IMS:IME,JMS:JME) :: SLP,TSFC,ZMSLP
- REAL, DIMENSION(IMS:IME,JMS:JME,KMS:KME) :: Z3d_IN
- REAL,DIMENSION(JTB) :: PIN,ZIN,Y2,PIO,ZOUT,DUM1,DUM2
- REAL,DIMENSION(JTB) :: QIN,QOUT,TIN,TOUT
- !--------------------------------------------------------------------------------------
- ! CLEAN Z3D ARRAY FIRST
- DO K=KDS,KDE
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- Z3d(I,J,K)=0.0
- T3d(I,J,K)=0.0
- Q3d(I,J,K)=0.0
- ENDDO
- ENDDO
- ENDDO
- ! DETERMINE THE HEIGHTS ON THE PARENT DOMAIN
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- Z3d_IN(I,J,1)=FIS(I,J)*GI
- ENDDO
- ENDDO
- DO K = KDS,KDE-1
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- APELP = (PINT(I,J,K+1)+PINT(I,J,K))
- ! RTOPP = TRG*T(I,J,K)*(1.0+Q(I,J,K)*P608-CWM(I,J,K))/APELP
- RTOPP = TRG*T(I,J,K)*(1.0+Q(I,J,K)*P608)/APELP
- DZ = RTOPP*(DETA1(K)*PDTOP+DETA2(K)*PD(I,J)) ! (RTv/P_TOT)*D(P_HYDRO)
- Z3d_IN(I,J,K+1) = Z3d_IN(I,J,K) + DZ
- ENDDO
- ENDDO
- ENDDO
- ! CONSTRUCT STANDARD ISOBARIC SURFACES
- DO K=KDS,KDE ! target points in model interface levels (pint)
- PSTD(K) = ETA1(K)*PDTOP + ETA2(K)*(P_REF -PDTOP - PTOP) + PTOP
- ENDDO
- ! DETERMINE THE MSLP USE THAT TO CREATE HEIGHTS AT 1000. mb LEVEL. THESE HEIGHTS
- ! MAY ONLY BE USED IN VERTICAL INTERPOLATION TO ISOBARIC SURFACES WHICH ARE LOCATED
- ! BELOW GROUND LEVEL.
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- 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
- A = LAPSR*Z3d_IN(I,J,1)/TSFC(I,J)
- SLP(I,J) = PINT(I,J,1)*(1-A)**COEF2 ! sea level pressure
- B = (PSTD(1)/SLP(I,J))**COEF3
- ZMSLP(I,J)= TSFC(I,J)*LAPSI*(1.0 - B) ! Height at 1000. mb level
- ENDDO
- ENDDO
- ! INTERPOLATE Z3d_IN TO STANDARD PRESSURE INTERFACES. FOR LEVELS BELOW
- ! GROUND USE ZMSLP(I,J)
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- !
- ! clean local array before use of spline
- PIN=0.;ZIN=0.;Y2=0;PIO=0.;ZOUT=0.;DUM1=0.;DUM2=0.
- DO K=KDS,KDE ! inputs at model interfaces
- PIN(K) = PINT(I,J,KDE-K+1)
- ZIN(K) = Z3d_IN(I,J,KDE-K+1)
- ENDDO
- IF(PINT(I,J,1) .LE. PSTD(1))THEN
- PIN(KDE) = PSTD(1)
- ZIN(KDE) = ZMSLP(I,J)
- ENDIF
- !
- Y2(1 )=0.
- Y2(KDE)=0.
- !
- DO K=KDS,KDE
- PIO(K)=PSTD(K)
- ENDDO
- !
- CALL SPLINE1(I,J,JTB,KDE,PIN,ZIN,Y2,KDE,PIO,ZOUT,DUM1,DUM2) ! interpolate
- !
- DO K=KDS,KDE ! inputs at model interfaces
- Z3d(I,J,K)=ZOUT(K)
- ENDDO
- ENDDO
- ENDDO
- !
- ! INTERPOLATE TEMPERATURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW
- ! GROUND USE A LAPSE RATE ATMOSPHERE
- !
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- !
- ! clean local array before use of spline or linear interpolation
- !
- PIN=0.;TIN=0.;Y2=0;PIO=0.;TOUT=0.;DUM1=0.;DUM2=0.
- DO K=KDS+1,KDE ! inputs at model levels
- PIN(K-1) = EXP((ALOG(PINT(I,J,KDE-K+1))+ALOG(PINT(I,J,KDE-K+2)))*0.5)
- TIN(K-1) = T(I,J,KDE-K+1)
- ENDDO
- IF(PINT(I,J,1) .LE. PSTD(1))THEN
- PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
- ZMID = 0.5*(Z3d_IN(I,J,1)+Z3d_IN(I,J,2))
- TIN(KDE-1) = T(I,J,1) + LAPSR*(ZMID-ZMSLP(I,J))
- ENDIF
- !
- Y2(1 )=0.
- Y2(KDE-1)=0.
- !
- DO K=KDS,KDE-1
- PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
- ENDDO
- CALL SPLINE1(I,J,JTB,KDE-1,PIN,TIN,Y2,KDE-1,PIO,TOUT,DUM1,DUM2) ! interpolate
- DO K=KDS,KDE-1 ! inputs at model levels
- T3d(I,J,K)=TOUT(K)
- ENDDO
- ENDDO
- ENDDO
- !
- ! INTERPOLATE MOISTURE ONTO THE STANDARD PRESSURE LEVELS. FOR LEVELS BELOW
- ! GROUND USE THE SURFACE MOISTURE
- !
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- !
- ! clean local array before use of spline or linear interpolation
- PIN=0.;QIN=0.;Y2=0;PIO=0.;QOUT=0.;DUM1=0.;DUM2=0.
- DO K=KDS+1,KDE ! inputs at model levels
- PIN(K-1) = EXP((ALOG(PINT(I,J,KDE-K+1))+ALOG(PINT(I,J,KDE-K+2)))*0.5)
- QIN(K-1) = Q(I,J,KDE-K+1)
- ENDDO
- IF(PINT(I,J,1) .LE. PSTD(1))THEN
- PIN(KDE-1) = EXP((ALOG(PSTD(1))+ALOG(PSTD(2)))*0.5)
- ! QIN(KDE-1) = QS(I,J)
- ENDIF
- Y2(1 )=0.
- Y2(KDE-1)=0.
- !
- DO K=KDS,KDE-1
- PIO(K)=EXP((ALOG(PSTD(K))+ALOG(PSTD(K+1)))*0.5)
- ENDDO
- CALL SPLINE1(I,J,JTB,KDE-1,PIN,QIN,Y2,KDE-1,PIO,QOUT,DUM1,DUM2) ! interpolate
- DO K=KDS,KDE-1 ! inputs at model levels
- Q3d(I,J,K)=QOUT(K)
- ENDDO
- ENDDO
- ENDDO
- END SUBROUTINE BASE_STATE_PARENT
- !=============================================================================
- SUBROUTINE SPLINE1(I,J,JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
- !
- ! ******************************************************************
- ! * *
- ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
- ! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
- ! * *
- ! * PROGRAMER Z. JANJIC *
- ! * *
- ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
- ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
- ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
- ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
- ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
- ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
- ! * SPECIFIED. *
- ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
- ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
- ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
- ! * AND LE XOLD(NOLD). *
- ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
- ! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
- ! * *
- ! ******************************************************************
- !---------------------------------------------------------------------
- IMPLICIT NONE
- !---------------------------------------------------------------------
- INTEGER,INTENT(IN) :: I,J,JTBX,NNEW,NOLD
- REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
- REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
- REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
- !
- INTEGER :: II,JJ,K,K1,K2,KOLD,NOLDM1
- REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
- ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
- CHARACTER(LEN=255) :: message
- !---------------------------------------------------------------------
- ! debug
- II=9999 !67 !35 !50 !4
- JJ=9999 !31 !73 !115 !192
- #if defined(EXPENSIVE_HWRF_DEBUG_STUFF)
- IF(I.eq.II.and.J.eq.JJ)THEN
- WRITE(message,*)'DEBUG in SPLINE1:HSO= ',xnew(1:nold)
- CALL wrf_debug(1,trim(message))
- DO K=1,NOLD
- WRITE(message,*)'DEBUG in SPLINE1:L,ZETAI,PINTI= ' &
- ,K,YOLD(K),XOLD(K)
- CALL wrf_debug(1,trim(message))
- ENDDO
- ENDIF
- #endif
- !
- NOLDM1=NOLD-1
- !
- DXL=XOLD(2)-XOLD(1)
- DXR=XOLD(3)-XOLD(2)
- DYDXL=(YOLD(2)-YOLD(1))/DXL
- DYDXR=(YOLD(3)-YOLD(2))/DXR
- RTDXC=0.5/(DXL+DXR)
- !
- P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
- Q(1)=-RTDXC*DXR
- !
- IF(NOLD.EQ.3)GO TO 150
- !---------------------------------------------------------------------
- K=3
- !
- 100 DXL=DXR
- DYDXL=DYDXR
- DXR=XOLD(K+1)-XOLD(K)
- DYDXR=(YOLD(K+1)-YOLD(K))/DXR
- DXC=DXL+DXR
- DEN=1./(DXL*Q(K-2)+DXC+DXC)
- !
- P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
- Q(K-1)=-DEN*DXR
- !
- K=K+1
- IF(K.LT.NOLD)GO TO 100
- !-----------------------------------------------------------------------
- 150 K=NOLDM1
- !
- 200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
- !
- K=K-1
- IF(K.GT.1)GO TO 200
- !-----------------------------------------------------------------------
- K1=1
- !
- 300 XK=XNEW(K1)
- !
- DO 400 K2=2,NOLD
- !
- IF(XOLD(K2).GT.XK)THEN
- KOLD=K2-1
- GO TO 450
- ENDIF
- !
- 400 CONTINUE
- !
- YNEW(K1)=YOLD(NOLD)
- GO TO 600
- !
- 450 IF(K1.EQ.1)GO TO 500
- IF(K.EQ.KOLD)GO TO 550
- !
- 500 K=KOLD
- !
- Y2K=Y2(K)
- Y2KP1=Y2(K+1)
- DX=XOLD(K+1)-XOLD(K)
- RDX=1./DX
- !
- AK=.1666667*RDX*(Y2KP1-Y2K)
- BK=0.5*Y2K
- CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
- !
- 550 X=XK-XOLD(K)
- XSQ=X*X
- !
- YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)
- ! debug
- #if defined(EXPENSIVE_HWRF_DEBUG_STUFF)
- if(i.eq.ii.and.j.eq.jj)then
- write(message,*) 'DEBUG:: k1,xnew(k1),ynew(k1): ', k1,xnew(k1),ynew(k1)
- CALL wrf_debug(1,trim(message))
- endif
- #endif
- !
- 600 K1=K1+1
- IF(K1.LE.NNEW)GO TO 300
- RETURN
- END SUBROUTINE SPLINE1
- !---------------------------------------------------------------------
- SUBROUTINE NEST_TERRAIN ( nest, config_flags )
- USE module_domain
- USE module_configure
- USE module_timing
- USE module_TERRAIN
- USE wrfsi_static
- USE module_SMOOTH_TERRAIN
- IMPLICIT NONE
- TYPE(domain) , POINTER :: nest
- TYPE(grid_config_rec_type) , INTENT(IN) :: config_flags
- !
- ! Local variables
- !
- 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_parent_start, j_parent_start
- INTEGER :: parent_grid_ratio
- INTEGER :: n,i,j,ii,jj,nnxp,nnyp
- INTEGER :: i_start,j_start,level
- REAL, DIMENSION(1,1), TARGET :: nothing
- REAL :: lah_nest_11, loh_nest_11
- INTEGER :: im_big, jm_big, i_add
- INTEGER :: im, jm
- CHARACTER(LEN=6) :: nestpath
- integer :: input_type
- character(len=128) :: input_fname
- character (len=32) :: cname
- integer :: ndim
- character (len=3) :: memorder
- character (len=32) :: stagger
- integer, dimension(3) :: domain_start, domain_end
- integer :: wrftype
- character (len=128), dimension(3) :: dimnames
- integer :: istatus
- integer :: handle
- integer :: comm_1, comm_2
- real, allocatable, dimension(:,:,:) :: real_domain
- character (len=10), dimension(24) :: name = (/ "XLAT_M ", &
- "XLONG_M ", &
- "XLAT_V ", &
- "XLONG_V ", &
- "E ", &
- "F ", &
- "LANDMASK ", &
- "LANDUSEF ", &
- "LU_INDEX ", &
- "HCNVX ", &
- "HSTDV ", &
- "HASYW ", &
- "HASYS ", &
- "HASYSW ", &
- "HASYNW ", &
- "HLENW ", &
- "HLENS ", &
- "HLENSW ", &
- "HLENNW ", &
- "HANIS ", &
- "HSLOP ", &
- "HANGL ", &
- "HZMAX ", &
- "HGT_M " /)
- integer, parameter :: IO_BIN=1, IO_NET=2
- integer :: io_form_input
- integer :: itsok,iteok,jtsok,jteok
- CHARACTER(LEN=512) :: message
- write(message,'("Nest d",I2," entering nest_terrain")') nest%id
- call wrf_debug(1,trim(message))
- call START_TIMING()
- write(message,*)"in NEST_TERRAIN config_flags%io_form_input = ", config_flags%io_form_input
- CALL wrf_debug(2,trim(message))
- write(message,*)"in NEST_TERRAIN config_flags%auxinput1_inname = ", config_flags%auxinput1_inname
- CALL wrf_debug(2,trim(message))
- io_form_input = config_flags%io_form_input
- if (config_flags%auxinput1_inname(1:7) == "met_nmm") then
- input_type = 2
- else
- input_type = 1
- end if
- !----------------------------------------------------------------------------------
- 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
- i_parent_start = nest%i_parent_start
- j_parent_start = nest%j_parent_start
- parent_grid_ratio = nest%parent_grid_ratio
- NNXP=IDE-1
- NNYP=JDE-1
- im = NNXP
- jm = NNYP
- ! Find nesting depth:
- call find_ijstart_level (nest,i_start,j_start,level)
- write(message,*)" nest%id =", nest%id , " i_start,j_start,level =", i_start,j_start,level
- CALL wrf_debug(2,trim(message))
- if ( level <= 0 ) then
- CALL wrf_error_fatal('this routine NEST_TERRAIN should not be called for top-level domain')
- end if
- ! Monitor process stores high-resolution topography:
- #ifdef DM_PARALLEL
- monitor_only: IF ( wrf_dm_on_monitor() ) THEN
- #endif
- call wrf_debug(1,'NEST_TERRAIN MASTER PROCESS')
- call MASTER(IDS,IDE,JDS,JDE)
- #ifdef DM_PARALLEL
- ELSE
- call wrf_debug(1,'NEST_TERRAIN SLAVE PROCESS')
- call SLAVE(IDS,IDE,JDS,JDE)
- ENDIF monitor_only
- #endif
- if(config_flags%terrain_smoothing==2) then
- call wrf_debug(1,'Call fast smoother (smooth_terrain)')
- call smooth_terrain(nest,12,12, &
- IDS,IDE,JDS,JDE,KDS,KDE, &
- IMS,IME,JMS,JME,KMS,KME, &
- ITS,ITE,JTS,JTE,KTS,KTE)
- elseif(config_flags%terrain_smoothing==1) then
- continue ! already handled this case in the call to MASTER
- elseif(config_flags%terrain_smoothing==0) then
- call wrf_debug(1,'Terrain smoothing is disabled.')
- else
- write(message,*) 'Invalid option for terrain_smoothing: ',config_flags%terrain_smoothing
- call wrf_error_fatal(message)
- endif
- DO J = jts,jte
- DO I = its,ite
- #ifdef IDEAL_NMM_TC
- nest%hres_fis(I,J)= 0.0 ! idealized
- #else
- nest%hres_fis(i,j)=9.81*nest%hres_avc(i,j)
- #endif
- ENDDO
- ENDDO
- write(message,'("Nest d",I0," nest_terrain")') nest%id
- call END_TIMING(trim(message))
- CONTAINS
- #ifdef DM_PARALLEL
- SUBROUTINE SLAVE(IDS,IDE,JDS,JDE)
- IMPLICIT NONE
- integer, intent(in) :: IDS,IDE,JDS,JDE
- REAL, DIMENSION(1,1) :: avc_nest,lnd_nest
- call wrf_debug(1,'call wrf_global_to_patch_real in nest_terrain')
- call wrf_global_to_patch_real(avc_nest,nest%hres_avc,nest%domdesc,'z','xy', &
- ids, ide-1, jds, jde-1, 1, 1, &
- ims, ime, jms, jme, 1, 1, &
- its, ite, jts, jte, 1, 1)
- call wrf_global_to_patch_real(lnd_nest,nest%hres_lnd,nest%domdesc,'z','xy', &
- ids, ide-1, jds, jde-1, 1, 1, &
- ims, ime, jms, jme, 1, 1, &
- its, ite, jts, jte, 1, 1)
- call wrf_debug(1,'back from wrf_global_to_patch_real in nest_terrain')
- END SUBROUTINE SLAVE
- #endif
- SUBROUTINE MASTER(IDS,IDE,JDS,JDE)
- IMPLICIT NONE
- integer, intent(in) :: IDS,IDE,JDS,JDE
- REAL, DIMENSION(IDS:IDE,JDS:JDE) :: avc_nest, lnd_nest
- type(nmm_terrain), pointer :: tr
- nullify(tr)
- avc_nest = 0.0
- lnd_nest = 0.0
-
- tr=>terrain_for(level,input_type,io_form_input)
- ! select subdomain from big fine grid
- i_add = mod(j_start+1,2)
- DO j=1,jde
- DO i=1,ide
- avc_nest(i,j) = tr%avc(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
- lnd_nest(i,j) = tr%lnd(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
- END DO
- END DO
- i=1 ; j=1
- lah_nest_11 = tr%lah(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
- loh_nest_11 = tr%loh(i_start+i-1 + mod(j+1,2)*i_add, j_start+j-1)
- IF(ABS(lah_nest_11-nest%HLAT(1,1)) .GE. 0.5 .OR. &
- ABS(loh_nest_11-nest%HLON(1,1)) .GE. 0.5)THEN
- WRITE(message,*)'SOME MATCHING TEST i_parent_start, j_parent_start',i_parent_start,j_parent_start
- CALL wrf_message(trim(message))
- CALL wrf_message('WRFSI LAT COMPUTED LAT')
- WRITE(message,*)lah_nest_11,nest%HLAT(1,1)
- CALL wrf_message(trim(message))
- CALL wrf_message('WRFSI LON COMPUTED LON')
- WRITE(message,*)loh_nest_11,nest%HLON(1,1)
- CALL wrf_message(trim(message))
- CALL wrf_message('CHECK WRFSI CONFIGURATION AND INPUT HIGH RESOLUTION TOPOGRAPHY AND/OR GRID RATIO')
- #ifndef IDEAL_NMM_TC
- CALL wrf_error_fatal('LATLON MISMATCH: ERROR READING static FILE FOR THE NEST')
- #endif
- ENDIF
- if(config_flags%terrain_smoothing==1) then
- call wrf_debug(1,'Call slow smoother (smdhld).')
- call smdhld(ids,ide,jds,jde,avc_nest,lnd_nest,12,12)
- endif
- #ifdef DM_PARALLEL
- call wrf_debug(1,'call wrf_global_to_patch_real in nest_terrain')
- call wrf_global_to_patch_real(avc_nest,nest%hres_avc,nest%domdesc,'z','xy', &
- ids, ide-1, jds, jde-1, 1, 1, &
- ims, ime, jms, jme, 1, 1, &
- its, ite, jts, jte, 1, 1)
- call wrf_global_to_patch_real(lnd_nest,nest%hres_lnd,nest%domdesc,'z','xy', &
- ids, ide-1, jds, jde-1, 1, 1, &
- ims, ime, jms, jme, 1, 1, &
- its, ite, jts, jte, 1, 1)
- call wrf_debug(1,'back from wrf_global_to_patch_real in nest_terrain')
- #endif
- END SUBROUTINE MASTER
-
- END SUBROUTINE NEST_TERRAIN
- !===========================================================================================
- SUBROUTINE med_init_domain_constants_nmm ( parent, nest) !, config_flags)
- ! Driver layer
- USE module_domain
- USE module_configure
- USE module_timing
- IMPLICIT NONE
- TYPE(domain) , POINTER :: parent, nest, grid
- !
- !
- INTERFACE
- SUBROUTINE med_initialize_nest_nmm ( grid &
- !
- # include <dummy_new_args.inc>
- !
- )
- USE module_domain
- USE module_configure
- USE module_timing
- IMPLICIT NONE
- TYPE(domain) , POINTER :: grid
- #include <dummy_new_decl.inc>
- END SUBROUTINE med_initialize_nest_nmm
- END INTERFACE
- !------------------------------------------------------------------------------
- ! PURPOSE:
- ! - initialize some data, mainly 2D & 3D nmm arrays very similar to
- ! those done in ./dyn_nmm/module_initialize_real.f
- !-----------------------------------------------------------------------------
- !
- grid => nest
- CALL med_initialize_nest_nmm( grid &
- !
- # include <actual_new_args.inc>
- !
- )
- END SUBROUTINE med_init_domain_constants_nmm
- SUBROUTINE med_initialize_nest_nmm( grid &
- !
- # include <dummy_new_args.inc>
- !
- )
- USE module_domain
- USE module_configure
- USE module_timing
- IMPLICIT NONE
- ! Local domain indices and counters.
- INTEGER :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte, &
- i, j, k, nnxp, nnyp
- TYPE(domain) , POINTER :: grid
- ! Local data
- LOGICAL, EXTERNAL :: wrf_dm_on_monitor
- INTEGER :: KHH,KVH,JAM,JA,IHL, IHH, L
- INTEGER :: II,JJ,ISRCH,ISUM
- INTEGER, ALLOCATABLE, DIMENSION(:) :: KHL2,KVL2,KHH2,KVH2,KHLA,KHHA,KVLA,KVHA
- INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
- !
- REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
- REAL(KIND=KNUM) :: STPH,CTPH,TDLM,TDPH,FP,TPH,TLM,TLM0
- REAL :: TPH0D,TLM0D,ANBI,TSPH,DTAD,DTCF,DT
- REAL :: ACDT,CDDAMP,DXP
- REAL :: WBD,SBD,WBI,SBI,EBI
- REAL :: DY_NMM0
- REAL :: RSNOW,SNOFAC
- REAL, ALLOCATABLE, DIMENSION(:) :: DXJ,WPDARJ,CPGFUJ,CURVJ, &
- FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
- HDACJ,DDMPUJ,DDMPVJ
- !
- REAL, PARAMETER:: SALP=2.60
- REAL, PARAMETER:: SNUP=0.040
- REAL, PARAMETER:: W_NMM=0.08
- ! REAL, PARAMETER:: COAC=0.75
- REAL, PARAMETER:: CODAMP=6.4
- REAL, PARAMETER:: TWOM=.00014584
- REAL, PARAMETER:: CP=1004.6
- REAL, PARAMETER:: DFC=1.0
- REAL, PARAMETER:: DDFC=1.0
- REAL, PARAMETER:: ROI=916.6
- REAL, PARAMETER:: R=287.04
- REAL, PARAMETER:: CI=2060.0
- REAL, PARAMETER:: ROS=1500.
- REAL, PARAMETER:: CS=1339.2
- REAL, PARAMETER:: DS=0.050
- REAL, PARAMETER:: AKS=.0000005
- REAL, PARAMETER:: DZG=2.85
- REAL, PARAMETER:: DI=.1000
- REAL, PARAMETER:: AKI=0.000001075
- REAL, PARAMETER:: DZI=2.0
- REAL, PARAMETER:: THL=210.
- REAL, PARAMETER:: PLQ=70000.
- REAL, PARAMETER:: ERAD=6371200.
- REAL, PARAMETER:: DTR=0.01745329
- REAL :: COAC
- CHARACTER(LEN=255) :: message
- ! Definitions of dummy arguments to solve
- #include <dummy_new_decl.inc>
- !#define COPY_IN
- !#include <scalar_derefs.inc>
- #ifdef DM_PARALLEL
- # include <data_calls.inc>
- #endif
- CALL get_ijk_from_grid ( grid , &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- !=================================================================================
- !
- !
- call nl_get_coac(grid%id,coac)
- DT=grid%dt !float(TIME_STEP)/parent_time_step_ratio
- NNXP=min(ITE,IDE-1)
- NNYP=min(JTE,JDE-1)
- JAM=6+2*((JDE-1)-10) ! this should be the fix instead of JAM=6+2*(NNYP-10)
- WRITE(message,*)'TIME STEP ON DOMAIN',grid%id,'==',dt
- CALL wrf_message(trim(message))
- WRITE(message,*)'IDS,IDE ON DOMAIN',grid%id,'==',ids,ide
- CALL wrf_message(trim(message))
- !
- ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
- ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
- ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),FADJ(JTS:NNYP))
- ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
- ALLOCATE(KHLA(JAM),KHHA(JAM))
- ALLOCATE(KVLA(JAM),KVHA(JAM))
- ! INITIALIZE SOME LAND/WATER SURFACE DATA ON THE BASIS OF INPUTS: SM, XICE, WEASD,
- ! INTERPOLATED FROM MOTHER (WRFSI) DOMAIN. THIS PART OF THE CODE HAS TO BE REVISITED
- ! LATER ON
- ! Since SM has been changed on parent domain to be 0 over sea ice it can not be used here
- ! to find where sea ice is. That's why alogirthm here is slightly different than the
- ! one used in module_initalize_real.f
- #ifdef HWRF
- !zhang's doing: added to AVOID THIS COMPUTATION IF THE NEST IS STARTED USING ANALYSIS FILE
- IF(.not. grid%analysis)THEN
- #endif
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF (grid%sm(I,J).GT.0.9) THEN ! OVER WATER SURFACE
- grid%epsr(I,J)= 0.97
- grid%embck(I,J)= 0.97
- grid%gffc(I,J)= 0.
- grid%albedo(I,J)=.06
- grid%albase(I,J)=.06
- ENDIF
- IF (grid%sice(I,J).GT.0.9) THEN ! OVER SEA-ICE
- grid%sm(I,J)=0.
- grid%si(I,J)=0.
- grid%gffc(I,J)=0.
- grid%albedo(I,J)=.60
- grid%albase(I,J)=.60
- ENDIF
- IF (grid%sm(I,J).LT.0.5.AND.grid%sice(I,J).LT.0.5) THEN ! OVER LAND SURFACE
- grid%si(I,J)=5.0*grid%weasd(I,J)/1000. ! SNOW WATER EQ (mm) OBTAINED FROM PARENT (grid%si) IS INTERPOLATED
- grid%epsr(I,J)=1.0 ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
- grid%embck(I,J)=1.0 ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
- grid%gffc(I,J)=0.0 ! just leave zero as irrelevant
- grid%sno(I,J)=grid%si(I,J)*.20 ! LAND-SNOW COVER
- ENDIF
- ENDDO
- ENDDO
- #if 0
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF(grid%sm(I,J).GT.0.9) THEN ! OVER WATER SURFACE
- !
- IF (XICE(I,J) .gt. 0)THEN ! XICE: SI INPUT ON PARENT, INTERPOLATED ONTO NEST
- grid%si(I,J)=1.0 ! INITIALIZE SI BASED ON XICE FROM INTERPOLATED INPUT
- ENDIF
- !
- grid%epsr(I,J)= 0.97 ! VALID OVER SEA SURFACE
- grid%embck(I,J)= 0.97 ! VALID OVER SEA SURFACE
- grid%gffc(I,J)= 0.
- grid%albedo(I,J)=.06
- grid%albase(I,J)=.06
- !
- IF(grid%si (I,J) .GT. 0.)THEN ! VALID OVER SEA-ICE
- grid%sm(I,J)=0.
- grid%si(I,J)=0. !
- grid%sice(I,J)=1.
- grid%gffc(I,J)=0. ! just leave zero as irrelevant
- grid%albedo(I,J)=.60 ! DEFINE grid%albedo
- grid%albase(I,J)=.60
- ENDIF
- !
- ELSE ! OVER LAND SURFACE
- !
- grid%si(I,J)=5.0*grid%weasd(I,J)/1000. ! SNOW WATER EQ (mm) OBTAINED FROM PARENT (grid%si) IS INTERPOLATED
- grid%epsr(I,J)=1.0 ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
- grid%embck(I,J)=1.0 ! EMISSIVITY DEFINED OVER LAND IN THE NESTED DOMAIN
- grid%gffc(I,J)=0.0 ! just leave zero as irrelevant
- grid%sice(I,J)=0. ! SEA ICE
- grid%sno(I,J)=grid%si(I,J)*.20 ! LAND-SNOW COVER
- !
- ENDIF
- !
- ENDDO
- ENDDO
- #endif
- ! This may just be a fix and may need some Registry related changes, later on
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- grid%vegfra(I,J)=grid%vegfrc(I,J)
- ENDDO
- ENDDO
- ! DETERMINE ALBEDO OVER LAND ON THE BASIS OF INPUTS: SM, ALBASE, MXSNAL & VEGFRA
- ! INTERPOLATED FROM MOTHER (WRFSI) DOMAIN
-
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN
- !
- IF ( (grid%sno(I,J) .EQ. 0.0) .OR. & ! SNOWFREE ALBEDO
- (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN
- grid%albedo(I,J) = grid%albase(I,J)
- ELSE
- IF (grid%sno(I,J) .LT. SNUP) THEN ! MODIFY ALBEDO IF SNOWCOVER:
- RSNOW = grid%sno(I,J)/SNUP ! BELOW SNOWDEPTH THRESHOLD
- SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
- ELSE
- SNOFAC = 1.0 ! ABOVE SNOWDEPTH THRESHOLD
- ENDIF
- grid%albedo(I,J) = grid%albase(I,J) &
- + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J))
- ENDIF
- !
- END IF
- grid%si(I,J)=5.0*grid%weasd(I,J)
- grid%sno(I,J)=grid%weasd(I,J)
- ! this block probably superfluous. Meant to guarantee land/sea agreement
- IF (grid%sm(I,J) .gt. 0.5)THEN
- grid%landmask(I,J)=0.0
- ELSE
- grid%landmask(I,J)=1.0
- ENDIF
- IF (grid%sice(I,J) .eq. 1.0) then !!!! change vegtyp and sltyp to fit seaice (desireable??)
- grid%isltyp(I,J)=16
- grid%ivgtyp(I,J)=24
- ENDIF
- ENDDO
- ENDDO
- ! Check land water interface
- ! The write(20,*) statements below were very slow on Jet when using
- ! the intel compiler (added hours to the runtime). The fix
- ! suggested by Chris Harrop was to send messages to wrf_message
- ! instead:
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS,MIN(ITE,IDE-1)
- IF(grid%sm(I,J).GT.0.9 .AND. grid%vegfra(I,J) .NE. 0) THEN
- #if defined(USE_SLOW_INTERFACE_CHECK)
- WRITE(20,*)'PROBLEM AT THE LAND-WATER INTERFACE (VEGFRA):', &
- I,J,grid%sm(I-1,J),grid%vegfra(I-1,j),grid%sm(I,J),grid%vegfra(I,J)
- #else
- WRITE(message,*)'PROBLEM AT THE LAND-WATER INTERFACE (VEGFRA):', &
- I,J,grid%sm(I-1,J),grid%vegfra(I-1,j),grid%sm(I,J),grid%vegfra(I,J)
- CALL wrf_message(trim(message))
- #endif
- ENDIF
- !
- #if defined(HWRF)
- ! HWRF should not perform the check below because the nmm_tsk is
- ! update with the correct skin temperature every timestep (on both
- ! land and sea points).
- #else
- IF(grid%sm(I,J).GT.0.9 .AND. grid%nmm_tsk(I,J) .NE. 0) THEN
- #if defined(USE_SLOW_INTERFACE_CHECK)
- WRITE(20,*)'PROBLEM AT THE LAND-WATER INTERFACE (NMM_TSK):', &
- I,J,grid%sm(I-1,J),grid%nmm_tsk(I-1,J),grid%sm(I,J),grid%nmm_tsk(I,J)
- #else
- WRITE(message,*)'PROBLEM AT THE LAND-WATER INTERFACE (NMM_TSK):', &
- I,J,grid%sm(I-1,J),grid%nmm_tsk(I-1,J),grid%sm(I,J),grid%nmm_tsk(I,J)
- CALL wrf_message(trim(message))
- #endif
- ENDIF
- #endif
- ENDDO
- ENDDO
- ! hardwire root depth for time being
- grid%rtdpth=0.
- grid%rtdpth(1)=0.1
- grid%rtdpth(2)=0.3
- grid%rtdpth(3)=0.6
- ! hardwire soil depth for time being
- grid%sldpth=0.
- grid%sldpth(1)=0.1
- grid%sldpth(2)=0.3
- grid%sldpth(3)=0.6
- grid%sldpth(4)=1.0
- #ifdef HWRF
- !zhang's doing: added to AVOID THIS COMPUTATION IF THE NEST IS STARTED USING ANALYSIS FILE
- ENDIF ! <------ for analysis set to false
- #endif
- !----------- END OF LAND SURFACE INITIALIZATION -------------------------------------
- !
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- grid%res(I,J)=1.
- ENDDO
- ENDDO
- ! INITIALIZE 2D BOUNDARY MASKS
- !! grid%hbm2:
-
- grid%hbm2=0.
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- IF((J .GE. 3 .and. J .LE. (JDE-1)-2) .AND. &
- (I .GE. 2 .and. I .LE. (IDE-1)-2+MOD(J,2))) THEN
- grid%hbm2(I,J)=1.
- ENDIF
- ENDDO
- ENDDO
- !! grid%hbm3:
- grid%hbm3=0.
- DO J=JTS,MIN(JTE,JDE-1)
- grid%ihwg(J)=mod(J+1,2)-1
- IF (J .ge. 4 .and. J .le. (JDE-1)-3) THEN
- IHL=(IDS+1)-grid%ihwg(J)
- IHH=(IDE-1)-2
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. IHL .and. I .le. IHH) grid%hbm3(I,J)=1.
- ENDDO
- ENDIF
- ENDDO
-
- !! grid%vbm2
- grid%vbm2=0.
- DO J=JTS,MIN(JTE,JDE-1)
- DO I=ITS,MIN(ITE,IDE-1)
- IF((J .ge. 3 .and. J .le. (JDE-1)-2) .AND. &
- (I .ge. 2 .and. I .le. (IDE-1)-1-MOD(J,2))) THEN
- grid%vbm2(I,J)=1.
- ENDIF
- ENDDO
- ENDDO
- !! grid%vbm3
- grid%vbm3=0.
- DO J=JTS,MIN(JTE,JDE-1)
- DO I=ITS,MIN(ITE,IDE-1)
- IF((J .ge. 4 .and. J .le. (JDE-1)-3) .AND. &
- (I .ge. 3-MOD(J,2) .and. I .le. (IDE-1)-2)) THEN
- grid%vbm3(I,J)=1.
- ENDIF
- ENDDO
- ENDDO
- TPH0D = grid%CEN_LAT
- TLM0D = grid%CEN_LON
- TPH0 = TPH0D*DTR
- WBD = grid%WBD0 ! gopal's doing: may use Registry WBD0 now
- WB = WBD*DTR
- SBD = grid%SBD0 ! gopal's doing: may use Registry SBD0 now
- SB = SBD*DTR
- DLM = grid%dlmd*DTR ! input now from med_nest_egrid_configure
- DPH = grid%dphd*DTR ! input now from med_nest_egrid_configure
- TDLM = DLM+DLM
- TDPH = DPH+DPH
- WBI = WB+TDLM
- SBI = SB+TDPH
- EBI = WB+((ide-1)-2)*TDLM ! gopal's doing: check this for nested domain
- ANBI = SB+((jde-1)-3)*DPH ! gopal's doing: check this for nested domain
- STPH0 = SIN(TPH0)
- CTPH0 = COS(TPH0)
- TSPH = 3600./grid%DT
- DTAD = 1.0
- DTCF = 4.0
- DY_NMM0= grid%dy_nmm ! ERAD*DPH; input now from med_nest_egrid_configure
- ! CORIOLIS PARAMETER (There appears to be some roundoff in computing TLM & STPH and other terms,
- ! in the nested domain. The problem needs to be revisited
- DO J=JTS,MIN(JTE,JDE-1)
- TLM0=WB-TDLM+MOD(J,2)*DLM ! remember this is a wind point
- TPH =SB+float(J-1)*DPH
- STPH=SIN(TPH)
- CTPH=COS(TPH)
- DO I=ITS,MIN(ITE,IDE-1)
- TLM=TLM0 + I*TDLM
- FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
- grid%f(I,J)=0.5*grid%DT*FP
- ENDDO
- ENDDO
- DO J=JTS,MIN(JTE,JDE-1)
- KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
- KVL2(J)=(IDE-1)*(J-1)-J/2+2
- KHH2(J)=(IDE-1)*J-J/2-1
- KVH2(J)=(IDE-1)*J-(J+1)/2-1
- ENDDO
- TPH=SB-DPH
- DO J=JTS,MIN(JTE,JDE-1)
- TPH=SB+float(J-1)*DPH
- DXP=ERAD*DLM*COS(TPH)
- DXJ(J)=DXP
- WPDARJ(J)=-W_NMM*((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)/ &
- (grid%DT*32.*DXP*DY_NMM0)
- CPGFUJ(J)=-grid%DT/(48.*DXP)
- CURVJ(J)=.5*grid%DT*TAN(TPH)/ERAD
- FCPJ(J)=grid%DT/(CP*192.*DXP*DY_NMM0)
- FDIVJ(J)=1./(12.*DXP*DY_NMM0)
- FADJ(J)=-grid%DT/(48.*DXP*DY_NMM0)*DTAD
- ACDT=grid%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY_NMM0**2)
- CDDAMP=CODAMP*ACDT
- HDACJ(J)=COAC*ACDT/(4.*DXP*DY_NMM0)
- DDMPUJ(J)=CDDAMP/DXP
- DDMPVJ(J)=CDDAMP/DY_NMM0
- ENDDO
- ! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
- WRITE(message,*)'NEW CHANGE',grid%f4d,grid%ef4t,grid%f4q
- CALL wrf_message(trim(message))
- DO L=KDS,KDE-1
- grid%rdeta(L)=1./grid%deta(L)
- grid%f4q2(L)=-.25*grid%DT*DTAD/grid%deta(L)
- ENDDO
- DO J=JTS,MIN(JTE,JDE-1)
- DO I=ITS,MIN(ITE,IDE-1)
- grid%dx_nmm(I,J)=DXJ(J)
- grid%wpdar(I,J)=WPDARJ(J)*grid%hbm2(I,J)
- grid%cpgfu(I,J)=CPGFUJ(J)*grid%vbm2(I,J)
- grid%curv(I,J)=CURVJ(J)*grid%vbm2(I,J)
- grid%fcp(I,J)=FCPJ(J)*grid%hbm2(I,J)
- grid%fdiv(I,J)=FDIVJ(J)*grid%hbm2(I,J)
- grid%fad(I,J)=FADJ(J)
- grid%hdacv(I,J)=HDACJ(J)*grid%vbm2(I,J)
- grid%hdac(I,J)=HDACJ(J)*1.25*grid%hbm2(I,J)
- ENDDO
- ENDDO
- DO J=JTS, MIN(JTE,JDE-1)
- IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
- KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. 2 .and. I .le. KHH) THEN
- grid%hdac(I,J)=grid%hdac(I,J)* DFC
- ENDIF
- ENDDO
- ELSE
- KHH=2+MOD(J,2)
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. 2 .and. I .le. KHH) THEN
- grid%hdac(I,J)=grid%hdac(I,J)* DFC
- ENDIF
- ENDDO
- KHH=(IDE-1)-2+MOD(J,2)
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
- grid%hdac(I,J)=grid%hdac(I,J)* DFC
- ENDIF
- ENDDO
- ENDIF
- ENDDO
- DO J=JTS,min(JTE,JDE-1)
- DO I=ITS,min(ITE,IDE-1)
- grid%ddmpu(I,J)=DDMPUJ(J)*grid%vbm2(I,J)
- grid%ddmpv(I,J)=DDMPVJ(J)*grid%vbm2(I,J)
- grid%hdacv(I,J)=grid%hdacv(I,J)*grid%vbm2(I,J)
- ENDDO
- ENDDO
- ! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
- DO J=JTS,MIN(JTE,JDE-1)
- IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
- KVH=(IDE-1)-1-MOD(J,2)
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. 2 .and. I .le. KVH) THEN
- grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
- grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
- grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
- ENDIF
- ENDDO
- ELSE
- KVH=3-MOD(J,2)
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. 2 .and. I .le. KVH) THEN
- grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
- grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
- grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
- ENDIF
- ENDDO
- KVH=(IDE-1)-1-MOD(J,2)
- DO I=ITS,MIN(ITE,IDE-1)
- IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
- grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
- grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
- grid%hdacv(I,J)=grid%hdacv(I,J)*DFC
- ENDIF
- ENDDO
- ENDIF
- ENDDO
- ! This one was left over for nested domain
-
- DO J = JTS, MIN(JTE,JDE-1)
- DO I = ITS, MIN(ITE,IDE-1)
- grid%GLAT(I,J)=grid%HLAT(I,J)*DTR
- grid%GLON(I,J)=grid%HLON(I,J)*DTR
- ENDDO
- ENDDO
- !! compute EMT, EM on global domain, and only on task 0.
- ! IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
-
- ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
- DO J=JDS,JDE-1
- TPH=SB+float(J-1)*DPH
- DXP=ERAD*DLM*COS(TPH)
- EMJ(J)= grid%DT/( 4.*DXP)*DTAD
- EMTJ(J)=grid%DT/(16.*DXP)*DTAD
- ENDDO
- JA=0
- DO 161 J=3,5
- JA=JA+1
- KHLA(JA)=2
- KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
- 161 grid%emt(JA)=EMTJ(J)
- DO 162 J=(JDE-1)-4,(JDE-1)-2
- JA=JA+1
- KHLA(JA)=2
- KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
- 162 grid%emt(JA)=EMTJ(J)
- DO 163 J=6,(JDE-1)-5
- JA=JA+1
- KHLA(JA)=2
- KHHA(JA)=2+MOD(J,2)
- 163 grid%emt(JA)=EMTJ(J)
- DO 164 J=6,(JDE-1)-5
- JA=JA+1
- KHLA(JA)=(IDE-1)-2
- KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
- 164 grid%emt(JA)=EMTJ(J)
- ! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
- JA=0
- DO 171 J=3,5
- JA=JA+1
- KVLA(JA)=2
- KVHA(JA)=(IDE-1)-1-MOD(J,2)
- 171 grid%em(JA)=EMJ(J)
- DO 172 J=(JDE-1)-4,(JDE-2)-2
- JA=JA+1
- KVLA(JA)=2
- KVHA(JA)=(IDE-1)-1-MOD(J,2)
- 172 grid%em(JA)=EMJ(J)
- DO 173 J=6,(JDE-1)-5
- JA=JA+1
- KVLA(JA)=2
- KVHA(JA)=2+MOD(J+1,2)
- 173 grid%em(JA)=EMJ(J)
- DO 174 J=6,(JDE-1)-5
- JA=JA+1
- KVLA(JA)=(IDE-1)-2
- KVHA(JA)=(IDE-1)-1-MOD(J,2)
- 174 grid%em(JA)=EMJ(J)
- ! ENDIF ! wrf_dm_on_monitor
- !! must be a better place to put this, but will eliminate "phantom"
- !! wind points here (no wind point on eastern boundary of odd numbered rows)
- !!
- ! phantom
- IF (ABS(IDE-1-ITE) .eq. 1 ) THEN ! |
- CALL wrf_message('zero phantom winds') ! H [x] H V
- DO K=KDS,KDE-1 !
- DO J=JDS,JDE-1,2 ! V [H] V H
- IF (J .ge. JTS .and. J .le. JTE) THEN !
- grid%u(IDE-1,J,K)=0. ! H [x] H V
- grid%v(IDE-1,J,K)=0. ! ------ ------
- ENDIF ! ide-1 ide
- ENDDO ! NMM/si WRF
- ENDDO ! domain domain
- ENDIF ! (dummy)
- ! just a test for gravity waves
- ! PD=62000.
- ! grid%u=0.0
- ! grid%v=0.0
- ! T=300.
- ! Q=0.0
- ! Q2=0.0
- ! CWM=0.0
- ! FIS=0.0
- ! testx
- ! DO J = JTS, MIN(JTE,JDE-1)
- ! DO K = KTS,KTE
- ! DO I = ITS, MIN(ITE,IDE-1)
- ! grid%sm(I,J)=I
- ! grid%u(I,K,J)=J
- ! ENDDO
- ! ENDDO
- ! ENDDO
- !
- ! deallocs
- DEALLOCATE(KHL2,KVL2,KHH2,KVH2)
- DEALLOCATE(DXJ,WPDARJ,CPGFUJ,CURVJ)
- DEALLOCATE(FCPJ,FDIVJ,FADJ)
- DEALLOCATE(HDACJ,DDMPUJ,DDMPVJ)
- DEALLOCATE(KHLA,KHHA)
- DEALLOCATE(KVLA,KVHA)
- END SUBROUTINE med_initialize_nest_nmm
- !======================================================================
- !--------------------------------------------------------------------------------------
- #if 0
- SUBROUTINE initial_nest_pivot ( parent , nest, iloc, jloc )
- !==========================================================================================
- !
- ! This program produces i_start and j_start for the nested domain depending on the
- ! central lat-lon of the storm.
- !
- !==========================================================================================
- USE module_domain
- USE module_configure
- USE module_timing
- USE module_dm
- IMPLICIT NONE
- TYPE(domain) , POINTER :: parent , nest
- INTEGER, INTENT(OUT) :: ILOC,JLOC
- INTEGER :: IMS,IME,JMS,JME,KMS,KME
- INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
- INTEGER :: IMS,IME,JMS,JME,KMS,KME
- INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
- INTEGER :: NIDE,NJDE ! nest dimension
- INTEGER :: I,J,ITER,IDUM,JDUM
- REAL :: ALAT,ALON,DIFF1,DIFF2,ERR
- REAL :: parent_CLAT,parent_CLON,parent_SLAT,parent_SLON
- REAL :: parent_WBD,parent_SBD,parent_DLMD,parent_DPHD
- !========================================================================================
- ! 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)
- ! CALL nl_get_storm_lat (parent%ID, parent_SLAT)
- ! CALL nl_get_storm_lon (parent%ID, parent_SLON)
- ! 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
- NIDE = nest%ed31
- NJDE = nest%ed32
- parent_DLMD = parent%dx ! DLMD: dlamda in degrees
- parent_DPHD = parent%dy ! DPHD: dphi in degrees
- parent_WBD = -(IDE-2)*parent%dx ! WBD0: in deg;factor 2 takes care of dummy last column
- parent_SBD = -((JDE-1)/2)*parent%dy ! SBD0: in degrees; note that JDE-1 should be odd
- ALAT = parent_SLAT - 0.5*(NJDE-2)*parent_DPHD/nest%parent_grid_ratio
- ALON = parent_SLON - 1.0*(NIDE-2)*parent_DLMD/nest%parent_grid_ratio
- 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 )
- ! start iteration
- ILOC=-99
- JLOC=-99
- ERR=0.1
- ITER=1
- 100 CONTINUE
- DO J = JTS,min(JTE,JDE-1)
- DO I = ITS,min(ITE,IDE-1)
- DIFF1 = ABS(ALAT - parent%HLAT(I,J))
- DIFF2 = ABS(ALON - parent%HLON(I,J))
- IF(DIFF1 .LE. ERR .AND. DIFF2 .LE. ERR)THEN
- ILOC=I
- JLOC=J
- ENDIF
- ENDDO
- ENDDO
- CALL wrf_dm_maxval_integer ( ILOC, idum, jdum )
- CALL wrf_dm_maxval_integer ( JLOC, idum, jdum )
- IF(ILOC .EQ. -99 .AND. JLOC .EQ. -99)THEN
- ERR=ERR+0.1
- ITER=ITER+1
- IF(ITER .LE. 100)GO TO 100
- ENDIF
- IF(ILOC .NE. -99 .AND. JLOC .NE. -99)THEN
- WRITE(message,*)'NOTE: I_PARENT_START AND J_PARENT_START FOUND FOR THE NESTED DOMAIN CONFIGURATION AT ITER=',ITER
- CALL wrf_message(trim(message))
- WRITE(message,*)'istart=',ILOC
- CALL wrf_message(trim(message))
- WRITE(message,*)'jstart=',JLOC
- CALL wrf_message(trim(message))
- ELSE
- ILOC=IDE/3
- JLOC=JDE/3
- !
- WRITE(message,*)'WARNING: COULD NOT LOCATE I_PARENT_START AND J_PARENT_START FROM INPUT STORM INFO'
- CALL wrf_message(trim(message))
- WRITE(message,*)'ISTART=',IDE/3
- CALL wrf_message(trim(message))
- WRITE(message,*)'JSTART=',JDE/3
- CALL wrf_message(trim(message))
- ENDIF
- RETURN
- END SUBROUTINE initial_nest_pivot
- !============================================================================================
- #endif
- LOGICAL FUNCTION cd_feedback_mask_orig( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
- INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
- LOGICAL, INTENT(IN) :: xstag, ystag
- INTEGER ioff, joff
- ioff = 0 ; joff = 0
- IF ( xstag ) ioff = 1
- IF ( ystag ) joff = 1
- cd_feedback_mask_orig = ( pig .ge. ips_save+1 .and. &
- pjg .ge. jps_save+1 .and. &
- pig .le. ipe_save-1 +ioff .and. &
- pjg .le. jpe_save-1 +joff )
- END FUNCTION cd_feedback_mask_orig
- LOGICAL FUNCTION cd_feedback_mask( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
- INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
- LOGICAL, INTENT(IN) :: xstag, ystag
- INTEGER ioff, joff
- ioff = 0 ; joff = 0
- IF ( xstag ) ioff = 1
- IF ( ystag ) joff = 1
- cd_feedback_mask = ( pig .ge. ips_save+1 .and. &
- pjg .ge. jps_save+2 .and. &
- pig .le. ipe_save-1-mod(pjg-jps_save,2) .and. &
- pjg .le. jpe_save-2 )
- END FUNCTION cd_feedback_mask
- LOGICAL FUNCTION cd_feedback_mask_v( pig, ips_save, ipe_save , pjg, jps_save, jpe_save, xstag, ystag )
- INTEGER, INTENT(IN) :: pig, ips_save, ipe_save , pjg, jps_save, jpe_save
- LOGICAL, INTENT(IN) :: xstag, ystag
- INTEGER ioff, joff
- ioff = 0 ; joff = 0
- IF ( xstag ) ioff = 1
- IF ( ystag ) joff = 1
-
- cd_feedback_mask_v = ( pig .ge. ips_save+1 .and. &
- pjg .ge. jps_save+2 .and. &
- pig .le. ipe_save-1-mod(pjg-jps_save+1,2) .and. &
- pjg .le. jpe_save-2 )
- END FUNCTION cd_feedback_mask_v
- !----------------------------------------------------------------------------
- #else
- SUBROUTINE stub_nmm_nest_stub
- END SUBROUTINE stub_nmm_nest_stub
- #endif
- RECURSIVE SUBROUTINE find_ijstart_level ( grid, i_start, j_start, level )
- ! Dusan Jovic
- USE module_domain
- IMPLICIT NONE
- ! Input data.
- TYPE(domain) :: grid
- INTEGER, INTENT (OUT) :: i_start, j_start, level
- INTEGER :: iadd
- if (grid%parent_id == 0 ) then
- i_start = 1
- j_start = 1
- level = 0
- else
- call find_ijstart_level ( grid%parents(1)%ptr, i_start, j_start, level )
- if (level > 0) then
- iadd = (i_start-1)*3
- if ( mod(j_start,2).ne.0 .and. mod(grid%j_parent_start,2).ne.0 ) iadd = iadd - 1
- if ( mod(j_start,2).eq.0 .and. mod(grid%j_parent_start,2).eq.0 ) iadd = iadd + 2
- else
- iadd = -mod(grid%j_parent_start,2)
- end if
- i_start = iadd + grid%i_parent_start*3 - 1
- j_start = ( (j_start-1) + (grid%j_parent_start-1) ) * 3 + 1
- level = level + 1
- end if
- END SUBROUTINE find_ijstart_level