!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! MODULE LLXY_MODULE ! ! This module handles transformations between model grid coordinates and ! latitude-longitude coordinates. The actual transformations are done through ! the map_utils module. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module llxy_module use gridinfo_module use list_module use map_utils use module_debug use misc_definitions_module ! Parameters integer, parameter :: MAX_SOURCE_LEVELS = 20 ! Variables integer :: current_nest_number integer :: SOURCE_PROJ = 0 ! The following arrays hold values for all available domains ! NOTE: The entries in the arrays for "domain 0" are used for projection ! information of user-specified source data type (proj_info), dimension(-MAX_SOURCE_LEVELS:MAX_DOMAINS) :: proj_stack, proj_stack_nestlev ! The projection and domain that we have computed constants for integer :: computed_proj = INVALID integer :: computed_domain = INVALID contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: push_source_projection ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine push_source_projection(iprojection, user_stand_lon, user_truelat1, user_truelat2, & user_dxkm, user_dykm, user_dlat, user_dlon, user_known_x, & user_known_y, user_known_lat, user_known_lon, earth_radius) implicit none ! Arguments integer, intent(in) :: iprojection real, intent(in) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, & user_dlat, user_dlon, & user_known_x, user_known_y, user_known_lat, user_known_lon real, intent(in), optional :: earth_radius SOURCE_PROJ = SOURCE_PROJ-1 if (SOURCE_PROJ < -MAX_SOURCE_LEVELS) then call mprintf(.true.,ERROR,'In push_user_projection(), too many levels of user projections.') end if call map_init(proj_stack(SOURCE_PROJ)) if (iprojection == PROJ_LATLON) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & latinc=user_dlat, & loninc=user_dlon, & r_earth=earth_radius) else if (iprojection == PROJ_MERC) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & truelat1=user_truelat1, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_CYL) then call mprintf(.true.,ERROR,'Should not have PROJ_CYL as projection for ' & //'source data in push_source_projection()') else if (iprojection == PROJ_CASSINI) then call mprintf(.true.,ERROR,'Should not have PROJ_CASSINI as projection for ' & //'source data in push_source_projection()') else if (iprojection == PROJ_LC) then print*, 'call map_set from push_source_projection' print*, 'user_known_x, user_known_y, user_known_lat, user_known_lon: ', user_known_x, user_known_y, user_known_lat, user_known_lon call map_set(iprojection, proj_stack(SOURCE_PROJ), & truelat1=user_truelat1, & truelat2=user_truelat2, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_ALBERS_NAD83) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & truelat1=user_truelat1, & truelat2=user_truelat2, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_PS) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & truelat1=user_truelat1, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_PS_WGS84) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & truelat1=user_truelat1, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_GAUSS) then call map_set(iprojection, proj_stack(SOURCE_PROJ), & lat1=user_known_lat, & lon1=user_known_lon, & nlat=nint(user_dlat), & loninc=user_dlon, & r_earth=earth_radius) else if (iprojection == PROJ_ROTLL) then ! BUG: Implement this projection. end if end subroutine push_source_projection !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: pop_source_projection ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine pop_source_projection() implicit none SOURCE_PROJ = SOURCE_PROJ+1 call mprintf((SOURCE_PROJ > 0), ERROR, & 'In pop_user_projection(), projection stack has overflowed.') end subroutine pop_source_projection #ifdef _METGRID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: set_domain_projection ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine set_domain_projection(iprojection, user_stand_lon, user_truelat1, user_truelat2, & user_dxkm, user_dykm, user_dlat, user_dlon, & user_xdim, user_ydim, user_known_x, & user_known_y, user_known_lat, user_known_lon, & user_pole_lat, user_pole_lon, earth_radius) implicit none ! Arguments integer, intent(in) :: iprojection integer, intent(in) :: user_xdim, user_ydim real, intent(in) :: user_stand_lon, user_truelat1, user_truelat2, & user_dxkm, user_dykm, user_dlat, user_dlon, & user_known_x, user_known_y, user_known_lat, user_known_lon, & user_pole_lat, user_pole_lon real, intent(in), optional :: earth_radius current_nest_number = 1 call map_init(proj_stack(current_nest_number)) if (iprojection == PROJ_LATLON) then call map_set(iprojection, proj_stack(current_nest_number), & lat1=user_known_lat, & lon1=user_known_lon, & latinc=user_dlat, & loninc=user_dlon, & r_earth=earth_radius) else if (iprojection == PROJ_MERC) then call map_set(iprojection, proj_stack(current_nest_number), & truelat1=user_truelat1, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_CYL) then call map_set(iprojection, proj_stack(current_nest_number), & latinc=user_dlat, & loninc=user_dlon, & stdlon=user_stand_lon, & r_earth=earth_radius) else if (iprojection == PROJ_CASSINI) then call map_set(iprojection, proj_stack(current_nest_number), & latinc=user_dlat, & loninc=user_dlon, & dx=user_dxkm, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & lat0=user_pole_lat, & lon0=user_pole_lon, & knowni=user_known_x, & knownj=user_known_y, & r_earth=earth_radius) else if (iprojection == PROJ_LC) then print*, 'call map_set from set_domain_projection' print*, 'user_known_x, user_known_y, user_known_lat, user_known_lon: ', & user_known_x, user_known_y, user_known_lat, user_known_lon call map_set(iprojection, proj_stack(current_nest_number), & truelat1=user_truelat1, & truelat2=user_truelat2, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_ALBERS_NAD83) then call map_set(iprojection, proj_stack(current_nest_number), & truelat1=user_truelat1, & truelat2=user_truelat2, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_PS) then call map_set(iprojection, proj_stack(current_nest_number), & truelat1=user_truelat1, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_PS_WGS84) then call map_set(iprojection, proj_stack(current_nest_number), & truelat1=user_truelat1, & stdlon=user_stand_lon, & lat1=user_known_lat, & lon1=user_known_lon, & knowni=user_known_x, & knownj=user_known_y, & dx=user_dxkm) else if (iprojection == PROJ_GAUSS) then call map_set(iprojection, proj_stack(current_nest_number), & lat1=user_known_lat, & lon1=user_known_lon, & nlat=nint(user_dlat), & loninc=user_dlon, & r_earth=earth_radius) else if (iprojection == PROJ_ROTLL) then call map_set(iprojection, proj_stack(current_nest_number), & ixdim=user_xdim, & jydim=user_ydim, & phi=user_dlat, & lambda=user_dlon, & lat1=user_known_lat, & lon1=user_known_lon, & stagger=HH, & latinc=user_dykm, & loninc=user_dxkm, & r_earth=earth_radius) else if (iprojection == PROJ_ROTLLB) then PRINT*, 'call map_set (a)' call map_set(iprojection, proj_stack(current_nest_number), & ixdim=user_xdim, & jydim=user_ydim, & phi=user_dlat, & lambda=user_dlon, & lat1=user_known_lat, & lon1=user_known_lon, & stagger=HH, & !these latinc and lonic werent in previous b-grid generation latinc=user_dykm, & loninc=user_dxkm, & !these latinc and lonic werent in previous b-grid generation ! domcenlat_loc=proj_stack(current_nest_number)%dom_cen_lat, & ! domcenlon_loc=proj_stack(current_nest_number)%dom_cen_lon, & r_earth=earth_radius) end if end subroutine set_domain_projection #endif #ifdef _GEOGRID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: compute_nest_locations ! ! Purpose: This routine computes the variables necessary in determining the ! location of all nests without reference to the parent or coarse domains. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine compute_nest_locations() implicit none ! Local variables integer :: i, i_parent_start, j_parent_start, i_parent_end, j_parent_end integer :: ixdim_proj, jydim_proj, ii, jj, nest_level real :: temp_known_x, temp_known_y, temp_known_lat, temp_known_lon, & temp_dxkm, temp_dykm, temp_dlat, temp_dlon, val real :: xcenter, ycenter, wbd, sbd, wbd_guess, sbd_guess real :: tst_lat, tst_lon, cen_lat, cen_lon, rmin write(0,*) ' ' write(0,*) ' ' write(0,*) ' compute_nest_locations ' write(0,*) ' ' ! Set location of coarse/mother domain call map_init(proj_stack(1)) if (iproj_type == PROJ_LATLON) then call map_set(iproj_type, proj_stack(1), & lat1=known_lat, & lon1=known_lon, & latinc=dykm, & loninc=dxkm) else if (iproj_type == PROJ_MERC) then call map_set(iproj_type, proj_stack(1), & truelat1=truelat1, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dxkm) else if (iproj_type == PROJ_CYL) then call map_set(iproj_type, proj_stack(1), & latinc=dlatdeg, & loninc=dlondeg, & stdlon=stand_lon) else if (iproj_type == PROJ_CASSINI) then call map_set(iproj_type, proj_stack(1), & latinc=dlatdeg, & loninc=dlondeg, & dx=dxkm, & stdlon=stand_lon, & knowni=known_x, & knownj=known_y, & lat0=pole_lat, & lon0=pole_lon, & lat1=known_lat, & lon1=known_lon) else if (iproj_type == PROJ_LC) then print*, 'call map_set from compute_nest_locations' print*, 'user_known_x, user_known_y, user_known_lat, user_known_lon: ', & known_x, known_y, known_lat, known_lon call map_set(iproj_type, proj_stack(1), & truelat1=truelat1, & truelat2=truelat2, & stdlon=stand_lon, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dxkm) else if (iproj_type == PROJ_ALBERS_NAD83) then call map_set(iproj_type, proj_stack(1), & truelat1=truelat1, & truelat2=truelat2, & stdlon=stand_lon, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dxkm) else if (iproj_type == PROJ_PS) then call map_set(iproj_type, proj_stack(1), & truelat1=truelat1, & stdlon=stand_lon, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dxkm) else if (iproj_type == PROJ_PS_WGS84) then call map_set(iproj_type, proj_stack(1), & truelat1=truelat1, & stdlon=stand_lon, & lat1=known_lat, & lon1=known_lon, & knowni=known_x, & knownj=known_y, & dx=dxkm) else if (iproj_type == PROJ_GAUSS) then call map_set(iproj_type, proj_stack(current_nest_number), & lat1=known_lat, & lon1=known_lon, & nlat=nint(dykm), & loninc=dxkm) else if (iproj_type == PROJ_ROTLL) then PRINT*, 'call map_set (b)' call map_set(iproj_type, proj_stack(1), & ixdim=ixdim(1), & jydim=jydim(1), & phi=phi, & lambda=lambda, & lat1=known_lat, & lon1=known_lon, & latinc=dykm, & loninc=dxkm, & stagger=HH) else if (iproj_type == PROJ_ROTLLB) then PRINT*, 'call map_set for domain 1 for PROJ_ROTLLB' write(0,*)' phi lambda ', phi ,lambda call map_set(iproj_type, proj_stack(1), & ixdim=ixdim(1), & jydim=jydim(1), & phi=phi, & lambda=lambda, & wbd=-lambda, & sbd=-phi, & lat1=known_lat, & lon1=known_lon, & latinc=dykm, & loninc=dxkm, & domcenlat_loc=dom_cen_lat(1), & domcenlon_loc=dom_cen_lon(1), & stagger=HH) end if ! Now we can compute lat/lon <-> x/y for coarse domain call select_domain(1) ! Call a recursive procedure to find the lat/lon of the centerpoint for ! each domain do i=2,n_domains temp_known_x = real(ixdim(i))/2. temp_known_y = real(jydim(i))/2. call find_known_latlon(i, temp_known_x, temp_known_y, & temp_known_lat, temp_known_lon, & temp_dxkm, temp_dykm, temp_dlat, temp_dlon) if (iproj_type == PROJ_LATLON) then call map_set(iproj_type, proj_stack(i), & lat1=temp_known_lat, & lon1=temp_known_lon, & latinc=temp_dlat, & loninc=temp_dlon) else if (iproj_type == PROJ_MERC) then call map_set(iproj_type, proj_stack(i), & truelat1=truelat1, & lat1=temp_known_lat, & lon1=temp_known_lon, & knowni=temp_known_x, & knownj=temp_known_y, & dx=temp_dxkm) else if (iproj_type == PROJ_CYL) then call mprintf(.true.,ERROR,'Don''t know how to do nesting with PROJ_CYL ' & //'in compute_nest_locations()') else if (iproj_type == PROJ_CASSINI) then call map_set(iproj_type, proj_stack(i), & latinc=temp_dlat, & loninc=temp_dlon, & dx=temp_dxkm, & stdlon=stand_lon, & knowni=temp_known_x, & knownj=temp_known_y, & lat0=pole_lat, & lon0=pole_lon, & lat1=temp_known_lat, & lon1=temp_known_lon) else if (iproj_type == PROJ_LC) then print*, 'dx: ', temp_dxkm print*, 'call map_set from lower in compute_nest_locations' print*, 'user_known_x, user_known_y, user_known_lat, user_known_lon: ', temp_known_x, temp_known_y, temp_known_lat, temp_known_lon call map_set(iproj_type, proj_stack(i), & truelat1=truelat1, & truelat2=truelat2, & stdlon=stand_lon, & lat1=temp_known_lat, & lon1=temp_known_lon, & knowni=temp_known_x, & knownj=temp_known_y, & dx=temp_dxkm) else if (iproj_type == PROJ_ALBERS_NAD83) then call map_set(iproj_type, proj_stack(i), & truelat1=truelat1, & truelat2=truelat2, & stdlon=stand_lon, & lat1=temp_known_lat, & lon1=temp_known_lon, & knowni=temp_known_x, & knownj=temp_known_y, & dx=temp_dxkm) else if (iproj_type == PROJ_PS) then call map_set(iproj_type, proj_stack(i), & truelat1=truelat1, & stdlon=stand_lon, & lat1=temp_known_lat, & lon1=temp_known_lon, & knowni=temp_known_x, & knownj=temp_known_y, & dx=temp_dxkm) else if (iproj_type == PROJ_PS_WGS84) then call map_set(iproj_type, proj_stack(i), & truelat1=truelat1, & stdlon=stand_lon, & lat1=temp_known_lat, & lon1=temp_known_lon, & knowni=temp_known_x, & knownj=temp_known_y, & dx=temp_dxkm) else if (iproj_type == PROJ_GAUSS) then call map_set(iproj_type, proj_stack(current_nest_number), & lat1=temp_known_lat, & lon1=temp_known_lon, & nlat=nint(temp_dykm), & loninc=temp_dxkm) else if (iproj_type == PROJ_ROTLL .or. iproj_type == PROJ_ROTLLB) then nest_level=get_nest_level(i) ixdim_proj = ixdim(1)*(3**nest_level)-(3**nest_level-1) jydim_proj = jydim(1)*(3**nest_level)-(3**nest_level-1) rmin=5*temp_dxkm do jj = 1, jydim_proj do ii = 1, ixdim_proj call ijll_rotlatlon_noproj(real(ii), real(jj), ixdim_proj, jydim_proj, & temp_dxkm, temp_dykm, known_lat, known_lon, & tst_lat, tst_lon ) val=( ( tst_lat-dom_cen_lat(i) )**2. + & ( tst_lon-dom_cen_lon(i) )**2. ) ** 0.5 if (val .lt. rmin .and. mod((II-1),3) .eq. 0 .and. mod((JJ-1),3) .eq. 0) then rmin=val print*, ' II, JJ, tst_lat, tst_lon: ', II, JJ, tst_lat, tst_lon cen_lat=tst_lat cen_lon=tst_lon endif enddo enddo write(0,*) 'cen_lat, dom_cen_lat(i): ', cen_lat, dom_cen_lat(i) write(0,*) 'cen_lon, dom_cen_lon(i): ', cen_lon, dom_cen_lon(i) ! call tll(dom_cen_lon(i),dom_cen_lat(i),xcenter, ycenter, known_lat, known_lon) call tll(cen_lon,cen_lat,xcenter, ycenter, known_lat, known_lon) write(0,*) ' xcenter, ycenter', xcenter, ycenter wbd_guess = xcenter - (ixdim(i)-1)*0.5*temp_dxkm sbd_guess = ycenter - (jydim(i)-1)*0.5*temp_dykm i_parent_start = nint ( 1 + ( wbd_guess - proj_stack(parent_id(i))%wbd ) /proj_stack(parent_id(i))%loninc ) j_parent_start = nint ( 1 + ( sbd_guess - proj_stack(parent_id(i))%sbd ) /proj_stack(parent_id(i))%latinc ) i_parent_end = i_parent_start + (ixdim(i)-1)/parent_grid_ratio(i) j_parent_end = j_parent_start + (jydim(i)-1)/parent_grid_ratio(i) wbd = proj_stack(parent_id(i))%wbd + (i_parent_start-1)*proj_stack(parent_id(i))%loninc sbd = proj_stack(parent_id(i))%sbd + (j_parent_start-1)*proj_stack(parent_id(i))%latinc call map_set(iproj_type, proj_stack(i), & ixdim=ixdim(i), & jydim=jydim(i), & phi=phi, & lambda=lambda, & wbd=wbd, & sbd=sbd, & lat1=known_lat, & lon1=known_lon, & latinc=temp_dykm, & loninc=temp_dxkm, & domcenlat_loc=dom_cen_lat(i), & domcenlon_loc=dom_cen_lon(i), & i_parent_start=i_parent_start, & j_parent_start=j_parent_start, & i_parent_end=i_parent_end, & j_parent_end=j_parent_end, & stagger=HH) write(0,*) 'ixdim=',ixdim(i) write(0,*) 'jydim=',jydim(i) write(0,*) 'phi=',phi write(0,*) 'lambda=',lambda write(0,*) 'wbd=',wbd write(0,*) 'sbd=',sbd write(0,*) 'lat1=',known_lat write(0,*) 'lon1=',known_lon write(0,*) 'latinc=',temp_dykm write(0,*) 'loninc=',temp_dxkm write(0,*) 'domcenlat_loc=',dom_cen_lat(i) write(0,*) 'domcenlon_loc=',dom_cen_lon(i) write(0,*) 'i_parent_start=',i_parent_start write(0,*) 'i_parent_end =',i_parent_end write(0,*) 'j_parent_start=',j_parent_start write(0,*) 'j_parent_end =',j_parent_end end if end do end subroutine compute_nest_locations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: find_known_latlon ! ! Purpose: This recursive routine computes the latitude and longitude for a ! specified x/y location in the given nest number, and also computes the ! grid spacing ! ! NOTE: This routine assumes that xytoll will work correctly for the ! coarse domain. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! recursive subroutine find_known_latlon(n, rx, ry, rlat, rlon, dx, dy, dlat, dlon) implicit none ! Arguments integer, intent(in) :: n real, intent(in) :: rx, ry real, intent(out) :: rlat, rlon, dx, dy, dlat, dlon ! Local variables real :: x_in_parent, y_in_parent if (n == 1) then ! Stopping case for the recursion dx = dxkm dy = dykm dlat = dlatdeg dlon = dlondeg call ij_to_latlon(PROJ=proj_stack(current_nest_number), i=rx, j=ry, lat=rlat, lon=rlon) return else ! Recursive case x_in_parent = (rx - ((parent_grid_ratio(n)+1.)/2.)) & / parent_grid_ratio(n) + proj_stack(n)%i_parent_start y_in_parent = (ry - ((parent_grid_ratio(n)+1.)/2.)) & / parent_grid_ratio(n) + proj_stack(n)%j_parent_start call find_known_latlon(parent_id(n), x_in_parent, y_in_parent, rlat, rlon, dx, dy, dlat, dlon) dx = dx / parent_grid_ratio(n) dy = dy / parent_grid_ratio(n) dlat = dlat / parent_grid_ratio(n) dlon = dlon / parent_grid_ratio(n) end if end subroutine find_known_latlon !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: compute_nest_level_info ! ! Purpose: This routine computes the parameters describing a nesting level for ! NMM grids. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine compute_nest_level_info() implicit none ! Local variables integer :: i, nest_level, temp type (list) :: level_list call list_init(level_list) ! Set location of coarse/mother domain call map_init(proj_stack(1)) call map_set(PROJ_ROTLL, proj_stack(1), & ixdim=ixdim(1), & jydim=jydim(1), & phi=phi, & lambda=lambda, & lat1=known_lat, & lon1=known_lon, & latinc=dykm, & loninc=dxkm, & stagger=HH) do i=2,n_domains nest_level = get_nest_level(i) if (.not. list_search(level_list, ikey=nest_level, ivalue=temp)) then call list_insert(level_list, ikey=nest_level, ivalue=nest_level) ixdim(nest_level) = ixdim(1)*(3**(nest_level-1))-(3**(nest_level-1)-1) jydim(nest_level) = jydim(1)*(3**(nest_level-1))-(3**(nest_level-1)-1) call map_set(PROJ_ROTLL, proj_stack(nest_level), & ixdim = ixdim(nest_level), & jydim = jydim(nest_level), & phi = phi, & lambda = lambda, & lat1=known_lat, & lon1=known_lon, & latinc=(dykm/real((3**(nest_level-1)))), & loninc=(dxkm/real((3**(nest_level-1)))), & stagger=HH) end if end do call list_destroy(level_list) end subroutine compute_nest_level_info !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_domain_resolution ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_domain_resolution(dom_dx, dom_dy) implicit none ! Arguments real, intent(out) :: dom_dx, dom_dy ! The proj_info structure only stores dx, so set both dom_dx and dom_dy to dx dom_dx = proj_stack(current_nest_number)%dx dom_dy = proj_stack(current_nest_number)%dx end subroutine get_domain_resolution !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_nest_level ! ! Purpose: This function returns, given a grid ID number, the nesting level of ! that domain; the coarse domain is taken to have nesting level 1. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function get_nest_level(i) implicit none ! Arguments integer, intent(in) :: i ! Local variables integer :: j ! Return value integer :: get_nest_level ! If argument is the coarse domain, return if (i == 1) then get_nest_level = 1 return end if if (i > MAX_DOMAINS) then call mprintf(.true., ERROR, & 'get_nest_level() called with invalid grid ID of %i.',i1=i) end if ! If not the coarse domain, then nesting level is at least 2 ! Yes, this looks silly. But we do not have a grid_id array, so ! we must check on parent_id get_nest_level = 2 j = i do while (parent_id(j) /= 1) j = parent_id(j) get_nest_level = get_nest_level + 1 ! Sanity check if (get_nest_level > MAX_DOMAINS) then call mprintf(.true., ERROR, & 'Spooky nesting setup encountered in get_nest_level().') end if end do end function get_nest_level #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: select_domain ! ! Purpose: This routine is used to select which nest x/y <-> lat/lon ! conversions will be with respect to. For example, selecting domain 2 will ! cause the llxy routine to compute x/y locations with respect to domain 2 ! given a lat/lon. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine select_domain(domain_num) implicit none ! Arguments integer, intent(in) :: domain_num #ifdef _GEOGRID if (domain_num > n_domains) then call mprintf(.true.,ERROR,'In select_domain(), selected domain is greater than n_domains.') end if #endif #ifdef _METGRID if (domain_num > 1) then call mprintf(.true.,ERROR,'In select_domain(), selected domain is greater than 1.') end if #endif current_nest_number = domain_num end subroutine select_domain !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: iget_selected_domain ! ! Purpose: This function returns the number of the currently selected nest. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function iget_selected_domain() implicit none ! Return value integer :: iget_selected_domain iget_selected_domain = current_nest_number end function iget_selected_domain !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: lltoxy ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine lltoxy(xlat, xlon, x, y, stagger) implicit none ! Arguments integer, intent(in) :: stagger real, intent(in) :: xlat, xlon real, intent(out) :: x, y ! Account for grid staggering if (stagger == HH) then proj_stack(current_nest_number)%stagger = HH else if (stagger == VV) then proj_stack(current_nest_number)%stagger = VV end if call latlon_to_ij(proj_stack(current_nest_number), xlat, xlon, x, y) ! Account for grid staggering if (stagger == U) then x = x + 0.5 else if (stagger == V) then y = y + 0.5 end if end subroutine lltoxy !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: lltoxy ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine xytoll(x, y, xlat, xlon, stagger) implicit none ! Arguments integer, intent(in) :: stagger real, intent(in) :: x, y real, intent(out) :: xlat, xlon ! Local variables real :: rx, ry ! Account for grid staggering; we cannot modify x and y, so modify local ! copies of them if (stagger == U) then rx = x - 0.5 ry = y else if (stagger == V) then rx = x ry = y - 0.5 else if (stagger == HH) then proj_stack(current_nest_number)%stagger = HH rx = x ry = y else if (stagger == VV) then proj_stack(current_nest_number)%stagger = VV rx = x ry = y else rx = x ry = y end if call ij_to_latlon(PROJ=proj_stack(current_nest_number), i=rx, j=ry, lat=xlat, lon=xlon) end subroutine xytoll subroutine tll(almd,aphd,tlmd,tphd,tph0d,tlm0d) !------------------------------------------------------------------------------- real, intent(in) :: almd, aphd real, intent(out) :: tlmd, tphd real, intent(in) :: tph0d, tlm0d !------------------------------------------------------------------------------- real, parameter :: pi=3.141592654 real, parameter :: dtr=pi/180.0 ! real :: tph0, ctph0, stph0, relm, srlm, crlm real :: aph, sph, cph, cc, anum, denom !------------------------------------------------------------------------------- ! if (tlm0d==0.0.and.tph0d==0.0) then tlmd=almd tphd=aphd else tph0=tph0d*dtr ctph0=cos(tph0) stph0=sin(tph0) ! relm=(almd-tlm0d)*dtr srlm=sin(relm) crlm=cos(relm) aph=aphd*dtr sph=sin(aph) cph=cos(aph) cc=cph*crlm anum=cph*srlm denom=ctph0*cc+stph0*sph ! tlmd=atan2(anum,denom)/dtr tphd=asin(ctph0*sph-stph0*cc)/dtr end if ! return ! end subroutine tll end module llxy_module