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