From 1d40660e42cf0d6ee4a177503c3731a75f865e0f Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 10 Mar 2017 12:22:44 -0700 Subject: [PATCH] Enable metgrid to process native MPAS output files (#11) This merge introduces changes to the metgrid program that enable the WPS to read and interpolate native MPAS output files, thereby allowing MPAS to provide initial and boundary conditions for a WRF simulation. From a user perspective, two namelist variables in the &metgrid namelist record must be set: 1. constants_name must be set to the name of the MPAS "static" file, prefixed with "mpas:", e.g., constants_name = "mpas:x1.40962.static.nc". 2. fg_name must provide the prefix of MPAS netCDF output files, prefixed with "mpas:", e.g., fg_name = "mpas:FOO"; the corresponding MPAS netCDF files should be named "FOO.YYYY-MM-DD_hh.nc", where YYYY, MM, DD, and hh are the valid time year, month, day, and hour of the data contained in the file. Note that MPAS files should contain only one time per file. When writing output from MPAS for use as initial and boundary conditions in WRF, the following MPAS stream definition can serve as a template: --- geogrid/src/constants_module.F | 4 + metgrid/METGRID.TBL.ARW | 30 + metgrid/src/Makefile | 12 +- metgrid/src/interp_option_module.F | 101 ++- metgrid/src/mpas_mesh.F | 310 ++++++++ metgrid/src/process_domain_module.F | 1505 ++++++++++++++++++++++++++--------- metgrid/src/remapper.F | 1341 +++++++++++++++++++++++++++++++ metgrid/src/scan_input.F | 604 ++++++++++++++ metgrid/src/target_mesh.F | 214 +++++ 9 files changed, 3724 insertions(+), 397 deletions(-) create mode 100644 metgrid/src/mpas_mesh.F create mode 100644 metgrid/src/remapper.F create mode 100644 metgrid/src/scan_input.F create mode 100644 metgrid/src/target_mesh.F diff --git a/geogrid/src/constants_module.F b/geogrid/src/constants_module.F index a06122d..d46abe3 100644 --- a/geogrid/src/constants_module.F +++ b/geogrid/src/constants_module.F @@ -25,4 +25,8 @@ module constants_module real, parameter :: EARTH_RADIUS_M = 6370000. ! same as MM5 system real, parameter :: EARTH_CIRC_M = 2.*PI*EARTH_RADIUS_M + real, parameter :: P0 = 1.0e5 ! Reference surface pressure, Pa + real, parameter :: RD = 287.0 ! Gas constant for dry air + real, parameter :: CP = 1004.0 ! Heat capacity for dry air at const. pressure + end module constants_module diff --git a/metgrid/METGRID.TBL.ARW b/metgrid/METGRID.TBL.ARW index 9cfc1da..397eb31 100644 --- a/metgrid/METGRID.TBL.ARW +++ b/metgrid/METGRID.TBL.ARW @@ -133,6 +133,7 @@ name=SNOW flag_in_output=FLAG_SNOW ======================================== name=SKINTEMP +mpas_name=skintemp interp_option=sixteen_pt+four_pt+wt_average_4pt+wt_average_16pt+search masked=both interp_land_mask = LANDSEA(1) @@ -140,6 +141,7 @@ name=SKINTEMP fill_missing=0. ======================================== name=PSFC +mpas_name=surface_pressure interp_option=four_pt+average_4pt fill_lev=200100:const(200100.) flag_in_output=FLAG_PSFC @@ -560,6 +562,7 @@ name=SOILT050 flag_in_output=FLAG_SOILT050 ======================================== name=PMSL +mpas_name=mslp interp_option=sixteen_pt+four_pt+average_4pt flag_in_output=FLAG_SLP ======================================== @@ -599,6 +602,7 @@ name=T ; output_name=TT # If we get T, use entry from TT and # write the field out as TT ======================================== name=TT +mpas_name=theta mandatory=yes # MUST HAVE THIS FIELD interp_option=sixteen_pt+four_pt+average_4pt fill_missing=0. @@ -608,6 +612,7 @@ name=U ; output_name=UU # If we get U, use entry from UU and # write the field out as UU ======================================== name=UU +mpas_name=uReconstructZonal mandatory=yes # MUST HAVE THIS FIELD interp_option=sixteen_pt+four_pt+average_4pt is_u_field=yes @@ -619,6 +624,7 @@ name=V ; output_name=VV # If we get V, use entry from VV and # write the field out as VV ======================================== name=VV +mpas_name=uReconstructMeridional mandatory=yes # MUST HAVE THIS FIELD interp_option=sixteen_pt+four_pt+average_4pt is_v_field=yes @@ -633,42 +639,49 @@ name=SST flag_in_output=FLAG_SST ======================================== name=QV +mpas_name=qv interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QV ======================================== name=QR +mpas_name=qr interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QR ======================================== name=QC +mpas_name=qc interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QC ======================================== name=QI +mpas_name=qi interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QI ======================================== name=QS +mpas_name=qs interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QS ======================================== name=QG +mpas_name=qg interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) flag_in_output=FLAG_QG ======================================== name=QNI +mpas_name=ni interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) @@ -682,6 +695,7 @@ name=QNC flag_in_output=FLAG_QNC ======================================== name=QNR +mpas_name=nr interp_option=four_pt+average_4pt fill_missing=0. fill_lev=200100:const(0.) @@ -693,6 +707,7 @@ name=VPTMP fill_lev=200100:const(0.) ======================================== name=PRESSURE +mpas_name=pressure interp_option=sixteen_pt+four_pt+average_4pt fill_missing=0. fill_lev=200100:PSFC(200100) @@ -704,6 +719,7 @@ name=PRHO flag_in_output=FLAG_PRHO ======================================== name=GHT +mpas_name=zgrid interp_option=sixteen_pt+four_pt+average_4pt fill_missing=0. fill_lev=200100:SOILHGT(200100) @@ -1002,3 +1018,17 @@ name=QNIFA z_dim_name=num_qnwfa_levels interp_option=four_pt+average_4pt ======================================== +# This entry is only used for MPAS data +name=smois +mpas_name=smois + output=no + masked=water + fill_missing=1. +======================================== +# This entry is only used for MPAS data +name=tslb +mpas_name=tslb + output=no + masked=water + fill_missing=1. +======================================== diff --git a/metgrid/src/Makefile b/metgrid/src/Makefile index f78efdb..e99bed3 100644 --- a/metgrid/src/Makefile +++ b/metgrid/src/Makefile @@ -1,6 +1,6 @@ include ../../configure.wps -OBJS = cio.o wrf_debug.o bitarray_module.o constants_module.o datatype_module.o module_stringutil.o gridinfo_module.o metgrid.o input_module.o interp_module.o interp_option_module.o list_module.o llxy_module.o met_data_module.o minheap_module.o misc_definitions_module.o module_date_pack.o module_debug.o module_map_utils.o module_mergesort.o output_module.o parallel_module.o process_domain_module.o queue_module.o read_met_module.o rotate_winds_module.o storage_module.o write_met_module.o +OBJS = cio.o wrf_debug.o bitarray_module.o constants_module.o datatype_module.o module_stringutil.o gridinfo_module.o metgrid.o input_module.o interp_module.o interp_option_module.o list_module.o llxy_module.o met_data_module.o minheap_module.o misc_definitions_module.o module_date_pack.o module_debug.o module_map_utils.o module_mergesort.o output_module.o parallel_module.o process_domain_module.o queue_module.o read_met_module.o rotate_winds_module.o storage_module.o write_met_module.o scan_input.o mpas_mesh.o target_mesh.o remapper.o all: clear ; @@ -60,7 +60,7 @@ output_module.o: gridinfo_module.o misc_definitions_module.o module_debug.o para parallel_module.o: -process_domain_module.o: module_date_pack.o bitarray_module.o gridinfo_module.o input_module.o interp_module.o interp_option_module.o list_module.o llxy_module.o misc_definitions_module.o module_debug.o module_mergesort.o output_module.o parallel_module.o read_met_module.o rotate_winds_module.o storage_module.o +process_domain_module.o: module_date_pack.o bitarray_module.o gridinfo_module.o input_module.o interp_module.o interp_option_module.o list_module.o llxy_module.o misc_definitions_module.o module_debug.o module_mergesort.o output_module.o parallel_module.o read_met_module.o rotate_winds_module.o storage_module.o scan_input.o mpas_mesh.o target_mesh.o remapper.o queue_module.o: module_debug.o @@ -74,6 +74,14 @@ wrf_debug.o: gridinfo_module.o cio.o write_met_module.o: misc_definitions_module.o module_debug.o met_data_module.o +scan_input.o: + +mpas_mesh.o: scan_input.o + +target_mesh.o: + +remapper.o: scan_input.o mpas_mesh.o target_mesh.o + clean: $(RM) $(OBJS) *.f90 *.mod $(RM) metgrid.exe diff --git a/metgrid/src/interp_option_module.F b/metgrid/src/interp_option_module.F index 0b44b25..b300e57 100644 --- a/metgrid/src/interp_option_module.F +++ b/metgrid/src/interp_option_module.F @@ -15,7 +15,8 @@ module interp_option_module logical, pointer, dimension(:) :: output_this_field, is_u_field, is_v_field, is_derived_field, is_mandatory character (len=128), pointer, dimension(:) :: fieldname, interp_method, v_interp_method, & interp_mask, interp_land_mask, interp_water_mask, & - flag_in_output, output_name, from_input, z_dim_name, level_template + flag_in_output, output_name, from_input, z_dim_name, level_template, & + mpas_name character (len=1), pointer, dimension(:) :: interp_mask_relational, interp_land_mask_relational, interp_water_mask_relational type (list), pointer, dimension(:) :: fill_lev_list type (list) :: flag_in_output_list @@ -79,6 +80,7 @@ module interp_option_module num_entries = num_entries + 1 allocate(fieldname(num_entries)) + allocate(mpas_name(num_entries)) allocate(interp_method(num_entries)) allocate(v_interp_method(num_entries)) allocate(masked(num_entries)) @@ -111,6 +113,7 @@ module interp_option_module ! do i=1,num_entries fieldname(i) = ' ' + mpas_name(i) = ' ' flag_in_output(i) = ' ' output_name(i) = ' ' from_input(i) = '*' @@ -200,6 +203,15 @@ module interp_option_module fieldname(i) = ' ' fieldname(i)(1:ispace-idx) = buffer(idx+1:ispace-1) + else if (index('mpas_name',trim(buffer(1:idx-1))) /= 0 .and. & + len_trim('mpas_name') == len_trim(buffer(1:idx-1))) then + ispace = idx+1 + do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) + ispace = ispace + 1 + end do + mpas_name(i) = ' ' + mpas_name(i)(1:ispace-idx) = buffer(idx+1:ispace-1) + else if (index('from_input',trim(buffer(1:idx-1))) /= 0 .and. & len_trim('from_input') == len_trim(buffer(1:idx-1))) then ispace = idx+1 @@ -586,6 +598,93 @@ module interp_option_module !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: mpas_name_to_idx + ! + ! Pupose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + function mpas_name_to_idx(mpasname) result(idx) + + implicit none + + ! Arguments + character (len=*), intent(in) :: mpasname + + ! Return value + integer :: idx + + ! Local variables + integer :: i + + idx = 0 + do i=1,num_entries + if (trim(mpasname) == trim(mpas_name(i))) then + idx = i + exit + end if + end do + + end function mpas_name_to_idx + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: mpas_to_intermediate_name + ! + ! Pupose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + function mpas_to_intermediate_name(mpasname) result(intermediate_name) + + implicit none + + ! Arguments + character (len=*), intent(in) :: mpasname + + ! Return value + character (len=128) :: intermediate_name + + ! Local variables + integer :: i + + intermediate_name = fieldname(num_entries) + do i=1,num_entries + if (trim(mpasname) == trim(mpas_name(i))) then + intermediate_name = fieldname(i) + exit + end if + end do + + end function mpas_to_intermediate_name + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: mpas_output_stagger + ! + ! Pupose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + function mpas_output_stagger(mpasname) result(istagger) + + implicit none + + ! Arguments + character (len=*), intent(in) :: mpasname + + ! Return value + integer :: istagger + + ! Local variables + integer :: i + + istagger = M + do i=1,num_entries + if (trim(mpasname) == trim(mpas_name(i))) then + istagger = output_stagger(i) + exit + end if + end do + + end function mpas_output_stagger + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_gcell_threshold ! ! Pupose: diff --git a/metgrid/src/mpas_mesh.F b/metgrid/src/mpas_mesh.F new file mode 100644 index 0000000..b21302d --- /dev/null +++ b/metgrid/src/mpas_mesh.F @@ -0,0 +1,310 @@ +module mpas_mesh + + use scan_input + + type mpas_mesh_type + logical :: valid = .false. + integer :: nCells = 0 + integer :: nVertices = 0 + integer :: nEdges = 0 + integer :: maxEdges = 0 + integer, dimension(:), pointer :: landmask => null() + integer, dimension(:), pointer :: nEdgesOnCell => null() + integer, dimension(:,:), pointer :: cellsOnCell => null() + integer, dimension(:,:), pointer :: verticesOnCell => null() + integer, dimension(:,:), pointer :: cellsOnVertex => null() + integer, dimension(:,:), pointer :: edgesOnCell => null() + integer, dimension(:,:), pointer :: cellsOnEdge => null() + real, dimension(:), pointer :: latCell => null() + real, dimension(:), pointer :: lonCell => null() + real, dimension(:), pointer :: latVertex => null() + real, dimension(:), pointer :: lonVertex => null() + real, dimension(:), pointer :: latEdge => null() + real, dimension(:), pointer :: lonEdge => null() + end type mpas_mesh_type + + + contains + + + integer function mpas_mesh_setup(mesh_filename, mesh) result(stat) + + implicit none + + character (len=*), intent(in) :: mesh_filename + type (mpas_mesh_type), intent(out) :: mesh + + type (input_handle_type) :: handle + type (input_field_type) :: field + + stat = 0 + + if (scan_input_open(mesh_filename, handle) /= 0) then + stat = 1 + return + end if + + ! + ! nEdgesOnCell + ! + if (scan_input_for_field(handle, 'nEdgesOnCell', field) /= 0) then + stat = 1 + return + end if + + mesh % nCells = field % dimlens(1) + mesh % nVertices = 2 * (mesh % nCells - 2) + mesh % nEdges = 3 * (mesh % nCells - 2) + + stat = scan_input_read_field(field) + allocate(mesh % nEdgesOnCell(mesh % nCells)) + mesh % nEdgesOnCell(:) = field % array1i(:) + stat = scan_input_free_field(field) + + ! + ! landmask + ! + if (scan_input_for_field(handle, 'landmask', field) == 0) then + stat = scan_input_read_field(field) + allocate(mesh % landmask(mesh % nCells)) + mesh % landmask(:) = field % array1i(:) + stat = scan_input_free_field(field) + else ! no landmask available + allocate(mesh % landmask(mesh % nCells)) + mesh % landmask(:) = 1 + end if + + ! + ! cellsOnCell + ! + if (scan_input_for_field(handle, 'cellsOnCell', field) /= 0) then + stat = 1 + return + end if + + mesh % maxEdges = field % dimlens(1) + + stat = scan_input_read_field(field) + allocate(mesh % cellsOnCell(mesh % maxEdges, mesh % nCells)) + mesh % cellsOnCell(:,:) = field % array2i(:,:) + stat = scan_input_free_field(field) + + ! + ! verticesOnCell + ! + if (scan_input_for_field(handle, 'verticesOnCell', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % verticesOnCell(mesh % maxEdges, mesh % nCells)) + mesh % verticesOnCell(:,:) = field % array2i(:,:) + stat = scan_input_free_field(field) + + ! + ! cellsOnVertex + ! + if (scan_input_for_field(handle, 'cellsOnVertex', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % cellsOnVertex(3, mesh % nVertices)) + mesh % cellsOnVertex(:,:) = field % array2i(:,:) + stat = scan_input_free_field(field) + + ! + ! edgesOnCell + ! + if (scan_input_for_field(handle, 'edgesOnCell', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % edgesOnCell(mesh % maxEdges, mesh % nCells)) + mesh % edgesOnCell(:,:) = field % array2i(:,:) + stat = scan_input_free_field(field) + + ! + ! cellsOnEdge + ! + if (scan_input_for_field(handle, 'cellsOnEdge', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % cellsOnEdge(2, mesh % nEdges)) + mesh % cellsOnEdge(:,:) = field % array2i(:,:) + stat = scan_input_free_field(field) + + ! + ! latCell + ! + if (scan_input_for_field(handle, 'latCell', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % latCell(mesh % nCells)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % latCell(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % latCell(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + ! + ! lonCell + ! + if (scan_input_for_field(handle, 'lonCell', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % lonCell(mesh % nCells)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % lonCell(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % lonCell(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + ! + ! latVertex + ! + if (scan_input_for_field(handle, 'latVertex', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % latVertex(mesh % nVertices)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % latVertex(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % latVertex(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + ! + ! lonVertex + ! + if (scan_input_for_field(handle, 'lonVertex', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % lonVertex(mesh % nVertices)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % lonVertex(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % lonVertex(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + ! + ! latEdge + ! + if (scan_input_for_field(handle, 'latEdge', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % latEdge(mesh % nEdges)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % latEdge(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % latEdge(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + ! + ! lonEdge + ! + if (scan_input_for_field(handle, 'lonEdge', field) /= 0) then + stat = 1 + return + end if + + stat = scan_input_read_field(field) + allocate(mesh % lonEdge(mesh % nEdges)) + if (field % xtype == FIELD_TYPE_REAL) then + mesh % lonEdge(:) = field % array1r(:) + else if (field % xtype == FIELD_TYPE_DOUBLE) then + mesh % lonEdge(:) = real(field % array1d(:)) + end if + stat = scan_input_free_field(field) + + stat = scan_input_close(handle) + + mesh % valid = .true. + + end function mpas_mesh_setup + + + integer function mpas_mesh_free(mesh) result(stat) + + implicit none + + type (mpas_mesh_type), intent(inout) :: mesh + + + stat = 0 + + mesh % valid = .false. + mesh % nCells = 0 + mesh % nVertices = 0 + mesh % nEdges = 0 + mesh % maxEdges = 0 + + if (associated(mesh % landmask)) then + deallocate(mesh % landmask) + end if + if (associated(mesh % nEdgesOnCell)) then + deallocate(mesh % nEdgesOnCell) + end if + if (associated(mesh % cellsOnCell)) then + deallocate(mesh % cellsOnCell) + end if + if (associated(mesh % verticesOnCell)) then + deallocate(mesh % verticesOnCell) + end if + if (associated(mesh % cellsOnVertex)) then + deallocate(mesh % cellsOnVertex) + end if + if (associated(mesh % edgesOnCell)) then + deallocate(mesh % edgesOnCell) + end if + if (associated(mesh % cellsOnEdge)) then + deallocate(mesh % cellsOnEdge) + end if + if (associated(mesh % latCell)) then + deallocate(mesh % latCell) + end if + if (associated(mesh % lonCell)) then + deallocate(mesh % lonCell) + end if + if (associated(mesh % latVertex)) then + deallocate(mesh % latVertex) + end if + if (associated(mesh % lonVertex)) then + deallocate(mesh % lonVertex) + end if + if (associated(mesh % latEdge)) then + deallocate(mesh % latEdge) + end if + if (associated(mesh % lonEdge)) then + deallocate(mesh % lonEdge) + end if + + end function mpas_mesh_free + +end module mpas_mesh diff --git a/metgrid/src/process_domain_module.F b/metgrid/src/process_domain_module.F index b1cb370..bdc5727 100644 --- a/metgrid/src/process_domain_module.F +++ b/metgrid/src/process_domain_module.F @@ -1,5 +1,14 @@ module process_domain_module + use mpas_mesh + use target_mesh + use remapper + + type (mpas_mesh_type) :: mpas_source_mesh + type (target_mesh_type) :: wrf_target_mesh_m, wrf_target_mesh_u, wrf_target_mesh_v + type (remap_info_type) :: remap_info_m, remap_info_u, remap_info_v + real, dimension(:,:), allocatable, target :: xlat_rad, xlon_rad, xlat_u_rad, xlon_u_rad, xlat_v_rad, xlon_v_rad + contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -23,6 +32,7 @@ module process_domain_module logical, intent(in) :: extra_row, extra_col ! Local variables + integer :: istatus integer :: i, t, dyn_opt, & we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, & we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & @@ -182,6 +192,22 @@ module process_domain_module if (associated(geogrid_flags)) deallocate(geogrid_flags) call storage_delete_all() + + istatus = mpas_mesh_free(mpas_source_mesh) + + if (allocated(xlat_rad)) deallocate(xlat_rad) + if (allocated(xlon_rad)) deallocate(xlon_rad) + if (allocated(xlat_u_rad)) deallocate(xlat_u_rad) + if (allocated(xlon_u_rad)) deallocate(xlon_u_rad) + if (allocated(xlat_v_rad)) deallocate(xlat_v_rad) + if (allocated(xlon_v_rad)) deallocate(xlon_v_rad) + istatus = target_mesh_free(wrf_target_mesh_m) + istatus = target_mesh_free(wrf_target_mesh_u) + istatus = target_mesh_free(wrf_target_mesh_v) + + istatus = remap_info_free(remap_info_m) + istatus = remap_info_free(remap_info_u) + istatus = remap_info_free(remap_info_v) end subroutine process_domain @@ -697,19 +723,16 @@ module process_domain_module integer, parameter :: BDR_WIDTH = 3 ! Local variables - integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, & + integer :: istatus, fg_idx, idx, idxt, i, j, bottom_top_dim, & sm1, em1, sm2, em2, sm3, em3, & sp1, ep1, sp2, ep2, sp3, ep3, & sd1, ed1, sd2, ed2, sd3, ed3, & - u_idx, bdr_wdth + u_idx integer :: nmet_flags integer :: num_metgrid_soil_levs integer, pointer, dimension(:) :: soil_levels real :: rx, ry - real :: threshold - logical :: do_gcell_interp integer, pointer, dimension(:) :: u_levels, v_levels - real, pointer, dimension(:,:) :: halo_slab real, pointer, dimension(:,:,:) :: real_array character (len=19) :: output_date character (len=128) :: cname, title @@ -722,7 +745,6 @@ integer, parameter :: BDR_WIDTH = 3 nullify(soil_levels) nullify(u_levels) nullify(v_levels) - nullify(halo_slab) nullify(real_array) @@ -739,395 +761,30 @@ integer, parameter :: BDR_WIDTH = 3 call mprintf(.true.,STDOUT, ' %s', s1=input_name) call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name) - ! Do a first pass through this fg source to get all mask fields used - ! during interpolation - call get_interp_masks(trim(input_name), do_const_processing, temp_date) - - istatus = 0 - - ! Initialize the module for reading in the met fields - call read_met_init(trim(input_name), do_const_processing, temp_date, istatus) - - if (istatus == 0) then - - ! Process all fields and levels from the current file; read_next_met_field() - ! will return a non-zero status when there are no more fields to be read. - do while (istatus == 0) - - call read_next_met_field(fg_data, istatus) - - if (istatus == 0) then - - ! Find index into fieldname, interp_method, masked, and fill_missing - ! of the current field - idxt = num_entries + 1 - do idx=1,num_entries - if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. & - (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then - - got_this_field(idx) = .true. - - if (index(input_name,trim(from_input(idx))) /= 0 .or. & - (from_input(idx) == '*' .and. idxt == num_entries + 1)) then - idxt = idx - end if - - end if - end do - idx = idxt - if (idx > num_entries) idx = num_entries ! The last entry is a default - - ! Do we need to rename this field? - if (output_name(idx) /= ' ') then - fg_data%field = output_name(idx)(1:9) - - idxt = num_entries + 1 - do idx=1,num_entries - if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. & - (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then - - got_this_field(idx) = .true. - - if (index(input_name,trim(from_input(idx))) /= 0 .or. & - (from_input(idx) == '*' .and. idxt == num_entries + 1)) then - idxt = idx - end if - - end if - end do - idx = idxt - if (idx > num_entries) idx = num_entries ! The last entry is a default - end if - - ! Do a simple check to see whether this is a global dataset - ! Note that we do not currently support regional Gaussian grids - if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) & - .or. (fg_data%iproj == PROJ_GAUSS)) then - bdr_wdth = BDR_WIDTH - allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny)) - - halo_slab(1:fg_data%nx, 1:fg_data%ny) = & - fg_data%slab(1:fg_data%nx, 1:fg_data%ny) - - halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = & - fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny) - - halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = & - fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny) - - deallocate(fg_data%slab) - else - bdr_wdth = 0 - halo_slab => fg_data%slab - nullify(fg_data%slab) - end if - - call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl) - - call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, & - fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, & - fg_data%deltalon, fg_data%starti, fg_data%startj, & - fg_data%startlat, fg_data%startlon, & - fg_data%pole_lat, fg_data%pole_lon, & - fg_data%centerlat, fg_data%centerlon, & - real(fg_data%nx+1)/2., real(fg_data%ny+1)/2., & - earth_radius=fg_data%earth_radius*1000.) - - ! Initialize fg_input structure to store the field - field%header%version = 1 - field%header%date = fg_data%hdate//' ' - if (do_const_processing) then - field%header%time_dependent = .false. - field%header%constant_field = .true. - else - field%header%time_dependent = .true. - field%header%constant_field = .false. - end if - field%header%forecast_hour = fg_data%xfcst - field%header%fg_source = 'FG' - field%header%field = ' ' - field%header%field(1:9) = fg_data%field - field%header%units = ' ' - field%header%units(1:25) = fg_data%units - field%header%description = ' ' - field%header%description(1:46) = fg_data%desc - call get_z_dim_name(fg_data%field,field%header%vertical_coord) - field%header%vertical_level = nint(fg_data%xlvl) - field%header%sr_x = 1 - field%header%sr_y = 1 - field%header%array_order = 'XY ' - field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel - field%header%array_has_missing_values = .false. - nullify(field%r_arr) - nullify(field%valid_mask) - nullify(field%modified_mask) - - if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then - output_flags(idx) = flag_in_output(idx) - end if - - ! If we should not output this field, just list it as a mask field - if (output_this_field(idx)) then - field%header%mask_field = .false. - else - field%header%mask_field = .true. - end if - - ! - ! Before actually doing any interpolation to the model grid, we must check - ! whether we will be using the average_gcell interpolator that averages all - ! source points in each model grid cell - ! - do_gcell_interp = .false. - if (index(interp_method(idx),'average_gcell') /= 0) then - - call get_gcell_threshold(interp_method(idx), threshold, istatus) - if (istatus == 0) then - if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. & - fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then - fg_data%dx = abs(fg_data%deltalon) - fg_data%dy = abs(fg_data%deltalat) - else -! BUG: Need to more correctly handle dx/dy in meters. - fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees - fg_data%dy = fg_data%dy / 111000. - end if - if (gridtype == 'C') then - if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) & - do_gcell_interp = .true. - else if (gridtype == 'E') then - if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) & - do_gcell_interp = .true. - end if - end if - end if - - ! Interpolate to U staggering - if (output_stagger(idx) == U) then - - call storage_query_field(field, iqstatus) - if (iqstatus == 0) then - call storage_get_field(field, iqstatus) - call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & - ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) - if (associated(field%modified_mask)) then - call bitarray_destroy(field%modified_mask) - nullify(field%modified_mask) - end if - else - allocate(field%valid_mask) - call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1) - end if - - ! Save a copy of the fg_input structure for the U field so that we can find it later - if (is_u_field(idx)) call dup(field, u_field) - - allocate(field%modified_mask) - call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1) - - if (do_const_processing .or. field%header%time_dependent) then - call interp_met_field(input_name, fg_data%field, U, M, & - field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, & - we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & - halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & - field%modified_mask, process_bdy_width) - else - call mprintf(.true.,INFORM,' - already processed this field from constant file.') - end if - - ! Interpolate to V staggering - else if (output_stagger(idx) == V) then - - call storage_query_field(field, iqstatus) - if (iqstatus == 0) then - call storage_get_field(field, iqstatus) - call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & - ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) - if (associated(field%modified_mask)) then - call bitarray_destroy(field%modified_mask) - nullify(field%modified_mask) - end if - else - allocate(field%valid_mask) - call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1) - end if - - ! Save a copy of the fg_input structure for the V field so that we can find it later - if (is_v_field(idx)) call dup(field, v_field) - - allocate(field%modified_mask) - call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1) - - if (do_const_processing .or. field%header%time_dependent) then - call interp_met_field(input_name, fg_data%field, V, M, & - field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, & - we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & - halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & - field%modified_mask, process_bdy_width) - else - call mprintf(.true.,INFORM,' - already processed this field from constant file.') - end if - - ! Interpolate to VV staggering - else if (output_stagger(idx) == VV) then - - call storage_query_field(field, iqstatus) - if (iqstatus == 0) then - call storage_get_field(field, iqstatus) - call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & - ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) - if (associated(field%modified_mask)) then - call bitarray_destroy(field%modified_mask) - nullify(field%modified_mask) - end if - else - allocate(field%valid_mask) - call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) - end if - - ! Save a copy of the fg_input structure for the U field so that we can find it later - if (is_u_field(idx)) call dup(field, u_field) - - ! Save a copy of the fg_input structure for the V field so that we can find it later - if (is_v_field(idx)) call dup(field, v_field) - - allocate(field%modified_mask) - call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) - - if (do_const_processing .or. field%header%time_dependent) then - call interp_met_field(input_name, fg_data%field, VV, M, & - field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & + if (index(input_name, 'mpas:') == 1) then + call process_mpas_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, & + landmask, process_bdy_width, & + u_field, v_field, & + dom_dx, dom_dy, & + xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, & + output_flags, west_east_dim, south_north_dim, & we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & - halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & - field%modified_mask, process_bdy_width) - else - call mprintf(.true.,INFORM,' - already processed this field from constant file.') - end if - - ! All other fields interpolated to M staggering for C grid, H staggering for E grid - else - - call storage_query_field(field, iqstatus) - if (iqstatus == 0) then - call storage_get_field(field, iqstatus) - call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & - ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) - if (associated(field%modified_mask)) then - call bitarray_destroy(field%modified_mask) - nullify(field%modified_mask) - end if - else - allocate(field%valid_mask) - call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) - end if - - allocate(field%modified_mask) - call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) - - if (do_const_processing .or. field%header%time_dependent) then - if (gridtype == 'C') then - call interp_met_field(input_name, fg_data%field, M, M, & - field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & - we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & - halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & - field%modified_mask, process_bdy_width, landmask) - - else if (gridtype == 'E') then - call interp_met_field(input_name, fg_data%field, HH, M, & - field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & - we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & - halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & - field%modified_mask, process_bdy_width, landmask) - end if - else - call mprintf(.true.,INFORM,' - already processed this field from constant file.') - end if - - end if - - call bitarray_merge(field%valid_mask, field%modified_mask) - - deallocate(halo_slab) - - ! Store the interpolated field - call storage_put_field(field) - - call pop_source_projection() - - end if - end do - - call read_met_close() - - call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, & - fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, & - fg_data%deltalon, fg_data%starti, fg_data%startj, & - fg_data%startlat, fg_data%startlon, & - fg_data%pole_lat, fg_data%pole_lon, & - fg_data%centerlat, fg_data%centerlon, & - real(fg_data%nx+1)/2., real(fg_data%ny+1)/2., & - earth_radius=fg_data%earth_radius*1000.) - - - ! - ! If necessary, rotate winds to earth-relative for this fg source - ! - - call storage_get_levels(u_field, u_levels) - call storage_get_levels(v_field, v_levels) - - if (associated(u_levels) .and. associated(v_levels)) then - u_idx = 1 - do u_idx = 1, size(u_levels) - u_field%header%vertical_level = u_levels(u_idx) - call storage_get_field(u_field, istatus) - v_field%header%vertical_level = v_levels(u_idx) - call storage_get_field(v_field, istatus) - - if (associated(u_field%modified_mask) .and. & - associated(v_field%modified_mask)) then - - if (u_field%header%is_wind_grid_rel) then - if (gridtype == 'C') then - call map_to_met(u_field%r_arr, u_field%modified_mask, & - v_field%r_arr, v_field%modified_mask, & - we_mem_stag_s, sn_mem_s, & - we_mem_stag_e, sn_mem_e, & - we_mem_s, sn_mem_stag_s, & - we_mem_e, sn_mem_stag_e, & - xlon_u, xlon_v, xlat_u, xlat_v) - else if (gridtype == 'E') then - call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, & - v_field%r_arr, v_field%modified_mask, & - we_mem_s, sn_mem_s, & - we_mem_e, sn_mem_e, & - xlat_v, xlon_v) - end if - end if - - call bitarray_destroy(u_field%modified_mask) - call bitarray_destroy(v_field%modified_mask) - nullify(u_field%modified_mask) - nullify(v_field%modified_mask) - call storage_put_field(u_field) - call storage_put_field(v_field) - end if - - end do - - deallocate(u_levels) - deallocate(v_levels) - - end if - - call pop_source_projection() - + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e ) else - if (do_const_processing) then - call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name) - else - call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date)) - end if + call process_intermediate_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, & + landmask, process_bdy_width, & + u_field, v_field, & + dom_dx, dom_dy, & + xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, & + output_flags, west_east_dim, south_north_dim, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e ) end if fg_idx = fg_idx + 1 @@ -1137,7 +794,7 @@ integer, parameter :: BDR_WIDTH = 3 input_name = fg_name(fg_idx) end if end do ! while (input_name /= '*') - + ! ! Rotate winds from earth-relative to grid-relative ! @@ -1192,6 +849,13 @@ integer, parameter :: BDR_WIDTH = 3 got_this_field, output_flags) ! + ! Derive some MPAS fields from others, e.g., TT from theta and pressure + ! + if (.not. do_const_processing) then + call derive_mpas_fields() + end if + + ! ! Check that every mandatory field was found in input data ! do i=1,num_entries @@ -1330,6 +994,1059 @@ integer, parameter :: BDR_WIDTH = 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: process_intermediate_fields + ! + ! Purpose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine process_intermediate_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, & + landmask, process_bdy_width, & + u_field, v_field, & + dom_dx, dom_dy, & + xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, & + output_flags, west_east_dim, south_north_dim, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e ) + + use bitarray_module + use gridinfo_module + use interp_module + use interp_option_module + use llxy_module + use misc_definitions_module + use module_debug + use output_module + use parallel_module + use read_met_module + use rotate_winds_module + use storage_module + + implicit none + +! BUG: Move this constant to misc_definitions_module? +integer, parameter :: BDR_WIDTH = 3 + + character (len=*), intent(inout) :: input_name + logical, intent(in) :: do_const_processing + character (len=*), intent(in) :: temp_date + type (met_data), intent(inout) :: fg_data + logical, dimension(:), intent(inout) :: got_this_field + real, pointer, dimension(:,:) :: landmask + integer, intent(in) :: process_bdy_width + type (fg_input), intent(inout) :: u_field, v_field + character (len=128), dimension(:), intent(inout) :: output_flags + real, intent(in) :: dom_dx, dom_dy + real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v + integer, intent(in) :: west_east_dim, south_north_dim, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e + + integer :: istatus + integer :: idx, idxt, u_idx + integer :: iqstatus + real :: threshold + real, pointer, dimension(:,:) :: halo_slab => null() + integer, pointer, dimension(:) :: u_levels => null() + integer, pointer, dimension(:) :: v_levels => null() + integer :: bdr_wdth + logical :: do_gcell_interp + type (fg_input) :: field + + + ! Do a first pass through this fg source to get all mask fields used + ! during interpolation + call get_interp_masks(trim(input_name), do_const_processing, temp_date) + + istatus = 0 + + ! Initialize the module for reading in the met fields + call read_met_init(trim(input_name), do_const_processing, temp_date, istatus) + + if (istatus == 0) then + + ! Process all fields and levels from the current file; read_next_met_field() + ! will return a non-zero status when there are no more fields to be read. + do while (istatus == 0) + + call read_next_met_field(fg_data, istatus) + + if (istatus == 0) then + + ! Find index into fieldname, interp_method, masked, and fill_missing + ! of the current field + idxt = num_entries + 1 + do idx=1,num_entries + if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. & + (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then + + got_this_field(idx) = .true. + + if (index(input_name,trim(from_input(idx))) /= 0 .or. & + (from_input(idx) == '*' .and. idxt == num_entries + 1)) then + idxt = idx + end if + + end if + end do + idx = idxt + if (idx > num_entries) idx = num_entries ! The last entry is a default + + ! Do we need to rename this field? + if (output_name(idx) /= ' ') then + fg_data%field = output_name(idx)(1:9) + + idxt = num_entries + 1 + do idx=1,num_entries + if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. & + (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then + + got_this_field(idx) = .true. + + if (index(input_name,trim(from_input(idx))) /= 0 .or. & + (from_input(idx) == '*' .and. idxt == num_entries + 1)) then + idxt = idx + end if + + end if + end do + idx = idxt + if (idx > num_entries) idx = num_entries ! The last entry is a default + end if + + ! Do a simple check to see whether this is a global dataset + ! Note that we do not currently support regional Gaussian grids + if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) & + .or. (fg_data%iproj == PROJ_GAUSS)) then + bdr_wdth = BDR_WIDTH + allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny)) + + halo_slab(1:fg_data%nx, 1:fg_data%ny) = & + fg_data%slab(1:fg_data%nx, 1:fg_data%ny) + + halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = & + fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny) + + halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = & + fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny) + + deallocate(fg_data%slab) + else + bdr_wdth = 0 + halo_slab => fg_data%slab + nullify(fg_data%slab) + end if + + call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl) + + call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, & + fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, & + fg_data%deltalon, fg_data%starti, fg_data%startj, & + fg_data%startlat, fg_data%startlon, & + fg_data%pole_lat, fg_data%pole_lon, & + fg_data%centerlat, fg_data%centerlon, & + real(fg_data%nx+1)/2., real(fg_data%ny+1)/2., & + earth_radius=fg_data%earth_radius*1000.) + + ! Initialize fg_input structure to store the field + field%header%version = 1 + field%header%date = fg_data%hdate//' ' + if (do_const_processing) then + field%header%time_dependent = .false. + field%header%constant_field = .true. + else + field%header%time_dependent = .true. + field%header%constant_field = .false. + end if + field%header%forecast_hour = fg_data%xfcst + field%header%fg_source = 'FG' + field%header%field = ' ' + field%header%field(1:9) = fg_data%field + field%header%units = ' ' + field%header%units(1:25) = fg_data%units + field%header%description = ' ' + field%header%description(1:46) = fg_data%desc + call get_z_dim_name(fg_data%field,field%header%vertical_coord) + field%header%vertical_level = nint(fg_data%xlvl) + field%header%sr_x = 1 + field%header%sr_y = 1 + field%header%array_order = 'XY ' + field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel + field%header%array_has_missing_values = .false. + nullify(field%r_arr) + nullify(field%valid_mask) + nullify(field%modified_mask) + + if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then + output_flags(idx) = flag_in_output(idx) + end if + + ! If we should not output this field, just list it as a mask field + if (output_this_field(idx)) then + field%header%mask_field = .false. + else + field%header%mask_field = .true. + end if + + ! + ! Before actually doing any interpolation to the model grid, we must check + ! whether we will be using the average_gcell interpolator that averages all + ! source points in each model grid cell + ! + do_gcell_interp = .false. + if (index(interp_method(idx),'average_gcell') /= 0) then + + call get_gcell_threshold(interp_method(idx), threshold, istatus) + if (istatus == 0) then + if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. & + fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then + fg_data%dx = abs(fg_data%deltalon) + fg_data%dy = abs(fg_data%deltalat) + else +! BUG: Need to more correctly handle dx/dy in meters. + fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees + fg_data%dy = fg_data%dy / 111000. + end if + if (gridtype == 'C') then + if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) & + do_gcell_interp = .true. + else if (gridtype == 'E') then + if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) & + do_gcell_interp = .true. + end if + end if + end if + + ! Interpolate to U staggering + if (output_stagger(idx) == U) then + + call storage_query_field(field, iqstatus) + if (iqstatus == 0) then + call storage_get_field(field, iqstatus) + call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & + ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) + if (associated(field%modified_mask)) then + call bitarray_destroy(field%modified_mask) + nullify(field%modified_mask) + end if + else + allocate(field%valid_mask) + call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1) + end if + + ! Save a copy of the fg_input structure for the U field so that we can find it later + if (is_u_field(idx)) call dup(field, u_field) + + allocate(field%modified_mask) + call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1) + + if (do_const_processing .or. field%header%time_dependent) then + call interp_met_field(input_name, fg_data%field, U, M, & + field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & + field%modified_mask, process_bdy_width) + else + call mprintf(.true.,INFORM,' - already processed this field from constant file.') + end if + + ! Interpolate to V staggering + else if (output_stagger(idx) == V) then + + call storage_query_field(field, iqstatus) + if (iqstatus == 0) then + call storage_get_field(field, iqstatus) + call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & + ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) + if (associated(field%modified_mask)) then + call bitarray_destroy(field%modified_mask) + nullify(field%modified_mask) + end if + else + allocate(field%valid_mask) + call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1) + end if + + ! Save a copy of the fg_input structure for the V field so that we can find it later + if (is_v_field(idx)) call dup(field, v_field) + + allocate(field%modified_mask) + call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1) + + if (do_const_processing .or. field%header%time_dependent) then + call interp_met_field(input_name, fg_data%field, V, M, & + field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & + field%modified_mask, process_bdy_width) + else + call mprintf(.true.,INFORM,' - already processed this field from constant file.') + end if + + ! Interpolate to VV staggering + else if (output_stagger(idx) == VV) then + + call storage_query_field(field, iqstatus) + if (iqstatus == 0) then + call storage_get_field(field, iqstatus) + call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & + ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) + if (associated(field%modified_mask)) then + call bitarray_destroy(field%modified_mask) + nullify(field%modified_mask) + end if + else + allocate(field%valid_mask) + call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) + end if + + ! Save a copy of the fg_input structure for the U field so that we can find it later + if (is_u_field(idx)) call dup(field, u_field) + + ! Save a copy of the fg_input structure for the V field so that we can find it later + if (is_v_field(idx)) call dup(field, v_field) + + allocate(field%modified_mask) + call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) + + if (do_const_processing .or. field%header%time_dependent) then + call interp_met_field(input_name, fg_data%field, VV, M, & + field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & + field%modified_mask, process_bdy_width) + else + call mprintf(.true.,INFORM,' - already processed this field from constant file.') + end if + + ! All other fields interpolated to M staggering for C grid, H staggering for E grid + else + + call storage_query_field(field, iqstatus) + if (iqstatus == 0) then + call storage_get_field(field, iqstatus) + call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// & + ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl)) + if (associated(field%modified_mask)) then + call bitarray_destroy(field%modified_mask) + nullify(field%modified_mask) + end if + else + allocate(field%valid_mask) + call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) + end if + + allocate(field%modified_mask) + call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1) + + if (do_const_processing .or. field%header%time_dependent) then + if (gridtype == 'C') then + call interp_met_field(input_name, fg_data%field, M, M, & + field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & + field%modified_mask, process_bdy_width, landmask) + + else if (gridtype == 'E') then + call interp_met_field(input_name, fg_data%field, HH, M, & + field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, & + field%modified_mask, process_bdy_width, landmask) + end if + else + call mprintf(.true.,INFORM,' - already processed this field from constant file.') + end if + + end if + + call bitarray_merge(field%valid_mask, field%modified_mask) + + deallocate(halo_slab) + + ! Store the interpolated field + call storage_put_field(field) + + call pop_source_projection() + + end if + end do + + call read_met_close() + + call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, & + fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, & + fg_data%deltalon, fg_data%starti, fg_data%startj, & + fg_data%startlat, fg_data%startlon, & + fg_data%pole_lat, fg_data%pole_lon, & + fg_data%centerlat, fg_data%centerlon, & + real(fg_data%nx+1)/2., real(fg_data%ny+1)/2., & + earth_radius=fg_data%earth_radius*1000.) + + ! + ! If necessary, rotate winds to earth-relative for this fg source + ! + + call storage_get_levels(u_field, u_levels) + call storage_get_levels(v_field, v_levels) + + if (associated(u_levels) .and. associated(v_levels)) then + u_idx = 1 + do u_idx = 1, size(u_levels) + u_field%header%vertical_level = u_levels(u_idx) + call storage_get_field(u_field, istatus) + v_field%header%vertical_level = v_levels(u_idx) + call storage_get_field(v_field, istatus) + + if (associated(u_field%modified_mask) .and. & + associated(v_field%modified_mask)) then + + if (u_field%header%is_wind_grid_rel) then + if (gridtype == 'C') then + call map_to_met(u_field%r_arr, u_field%modified_mask, & + v_field%r_arr, v_field%modified_mask, & + we_mem_stag_s, sn_mem_s, & + we_mem_stag_e, sn_mem_e, & + we_mem_s, sn_mem_stag_s, & + we_mem_e, sn_mem_stag_e, & + xlon_u, xlon_v, xlat_u, xlat_v) + else if (gridtype == 'E') then + call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, & + v_field%r_arr, v_field%modified_mask, & + we_mem_s, sn_mem_s, & + we_mem_e, sn_mem_e, & + xlat_v, xlon_v) + end if + end if + + call bitarray_destroy(u_field%modified_mask) + call bitarray_destroy(v_field%modified_mask) + nullify(u_field%modified_mask) + nullify(v_field%modified_mask) + call storage_put_field(u_field) + call storage_put_field(v_field) + end if + + end do + + deallocate(u_levels) + deallocate(v_levels) + + end if + + call pop_source_projection() + + else + if (do_const_processing) then + call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name) + else + call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date)) + end if + end if + + end subroutine process_intermediate_fields + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: process_mpas_fields + ! + ! Purpose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine process_mpas_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, & + landmask, process_bdy_width, & + u_field, v_field, & + dom_dx, dom_dy, & + xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, & + output_flags, west_east_dim, south_north_dim, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e ) + + use bitarray_module + use gridinfo_module + use interp_module + use interp_option_module + use read_met_module + use llxy_module + use misc_definitions_module + use module_debug + use output_module + use parallel_module + use rotate_winds_module + use storage_module + use scan_input + use mpas_mesh + + implicit none + +! BUG: Move this constant to misc_definitions_module? +integer, parameter :: BDR_WIDTH = 3 + + character (len=*), intent(inout) :: input_name + logical, intent(in) :: do_const_processing + character (len=*), intent(in) :: temp_date + type (met_data), intent(inout) :: fg_data + logical, dimension(:), intent(inout) :: got_this_field + real, pointer, dimension(:,:) :: landmask + integer, intent(in) :: process_bdy_width + type (fg_input), intent(inout) :: u_field, v_field + character (len=128), dimension(:), intent(inout) :: output_flags + real, intent(in) :: dom_dx, dom_dy + real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v + integer, intent(in) :: west_east_dim, south_north_dim, & + we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, & + we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, & + sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, & + we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, & + sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e + + real, parameter :: deg2rad = asin(1.0) / 90.0 + + integer :: i, j, k + integer :: idx + integer :: istat + integer :: strlen + character (len=MAX_FILENAME_LEN) :: mpas_filename + integer :: nRecords + type (input_handle_type) :: mpas_handle + type (input_field_type) :: mpas_field + type (target_field_type) :: wrf_field + type (fg_input) :: field_to_store + + strlen = len_trim(input_name) + if (do_const_processing) then + write(mpas_filename,'(a)') input_name(6:strlen) + else + write(mpas_filename,'(a)') input_name(6:strlen)//'.'//trim(temp_date)//'.nc' + end if + call mprintf(.true.,LOGFILE,'Processing MPAS fields from file %s',s1=mpas_filename) + + ! + ! If we do not already have mesh information, get that now... + ! + if (.not. mpas_source_mesh % valid) then + if (mpas_mesh_setup(mpas_filename, mpas_source_mesh) /= 0) then + call mprintf(.true.,ERROR, 'Error setting up MPAS mesh %s with scan_input_open', s1=mpas_filename) + end if + end if + + ! + ! If we have not already defined the WRF grid, do that now... + ! + if (.not. wrf_target_mesh_m % valid) then + allocate(xlat_rad(size(xlat,1), size(xlat,2))) + allocate(xlon_rad(size(xlat,1), size(xlat,2))) + xlat_rad(:,:) = xlat(:,:) * deg2rad + xlon_rad(:,:) = xlon(:,:) * deg2rad + call mprintf(.true.,LOGFILE,'Need to set up WRF target mass-grid') + if (target_mesh_setup(wrf_target_mesh_m, lat2d=xlat_rad, lon2d=xlon_rad) /= 0) then + call mprintf(.true.,ERROR, 'Error setting up WRF target grid') + end if + + call mprintf(.true.,LOGFILE,'Also computing remapping weights...') + if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_m, remap_info_m) /= 0) then + call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid') + end if + else + call mprintf(.true.,LOGFILE,'Already set up WRF target mass-grid') + end if + + if (.not. wrf_target_mesh_u % valid) then + allocate(xlat_u_rad(size(xlat_u,1), size(xlat_u,2))) + allocate(xlon_u_rad(size(xlat_u,1), size(xlat_u,2))) + xlat_u_rad(:,:) = xlat_u(:,:) * deg2rad + xlon_u_rad(:,:) = xlon_u(:,:) * deg2rad + call mprintf(.true.,LOGFILE,'Need to set up WRF target U-grid') + if (target_mesh_setup(wrf_target_mesh_u, lat2d=xlat_u_rad, lon2d=xlon_u_rad) /= 0) then + call mprintf(.true.,ERROR, 'Error setting up WRF target grid') + end if + + call mprintf(.true.,LOGFILE,'Also computing remapping weights...') + if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_u, remap_info_u) /= 0) then + call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid') + end if + else + call mprintf(.true.,LOGFILE,'Already set up WRF target U-grid') + end if + + if (.not. wrf_target_mesh_v % valid) then + allocate(xlat_v_rad(size(xlat_v,1), size(xlat_v,2))) + allocate(xlon_v_rad(size(xlat_v,1), size(xlat_v,2))) + xlat_v_rad(:,:) = xlat_v(:,:) * deg2rad + xlon_v_rad(:,:) = xlon_v(:,:) * deg2rad + call mprintf(.true.,LOGFILE,'Need to set up WRF target V-grid') + if (target_mesh_setup(wrf_target_mesh_v, lat2d=xlat_v_rad, lon2d=xlon_v_rad) /= 0) then + call mprintf(.true.,ERROR, 'Error setting up WRF target grid') + end if + + call mprintf(.true.,LOGFILE,'Also computing remapping weights...') + if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_v, remap_info_v) /= 0) then + call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid') + end if + else + call mprintf(.true.,LOGFILE,'Already set up WRF target V-grid') + end if + + + if (scan_input_open(mpas_filename, mpas_handle, nRecords) /= 0) then + call mprintf(.true.,ERROR, 'Error opening %s with scan_input_open', s1=mpas_filename) + end if + + + ! Initialize fg_input structure to store the field + field_to_store%header%version = 1 + field_to_store%header%date = '?' + if (do_const_processing) then + field_to_store%header%time_dependent = .false. + field_to_store%header%constant_field = .true. + else + field_to_store%header%time_dependent = .true. + field_to_store%header%constant_field = .false. + end if + field_to_store%header%forecast_hour = 0.0 + field_to_store%header%fg_source = 'FG' + field_to_store%header%field = ' ' + field_to_store%header%field(1:9) = '?' + field_to_store%header%units = ' ' + field_to_store%header%units(1:25) = '?' + field_to_store%header%description = ' ' + field_to_store%header%description(1:46) = '?' + field_to_store%header%vertical_coord = 'z_dim_name' + field_to_store%header%vertical_level = 0 + field_to_store%header%sr_x = 1 + field_to_store%header%sr_y = 1 + field_to_store%header%array_order = 'XY ' + field_to_store%header%is_wind_grid_rel = .false. + field_to_store%header%array_has_missing_values = .false. + nullify(field_to_store%r_arr) + nullify(field_to_store%valid_mask) + nullify(field_to_store%modified_mask) + + ! If we should not output this field, just list it as a mask field +!??? if (output_this_field(idx)) then + field_to_store%header%mask_field = .false. +!??? else +!??? field%header%mask_field = .true. +!??? end if + + + do while (scan_input_next_field(mpas_handle, mpas_field) == 0) + + if (can_remap_field(mpas_field)) then + + ! Here, rename a few MPAS fields that would be difficult to treat + ! with METGRID.TBL options; principally, these are surface fields + ! that have different names from their upper-air counterparts. + if (trim(mpas_field % name) == 'u10') then + mpas_field % name = 'uReconstructZonal' + else if (trim(mpas_field % name) == 'v10') then + mpas_field % name = 'uReconstructMeridional' + else if (trim(mpas_field % name) == 'q2') then + mpas_field % name = 'qv' + else if (trim(mpas_field % name) == 't2m') then + mpas_field % name = 'theta' + end if + + ! Mark this MPAS field as "gotten" for any later checks + ! on mandatory fields + idx = mpas_name_to_idx(trim(mpas_field % name)) + if (idx > 0) then + got_this_field(idx) = .true. + if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then + output_flags(idx) = flag_in_output(idx) + end if + end if + + istat = scan_input_read_field(mpas_field, frame=1) + + field_to_store%map%stagger = mpas_output_stagger(mpas_field % name) + if (field_to_store%map%stagger == M) then + field_to_store%header%dim1(1) = we_mem_s + field_to_store%header%dim1(2) = we_mem_e + field_to_store%header%dim2(1) = sn_mem_s + field_to_store%header%dim2(2) = sn_mem_e + if (idx > 0) then + if (masked(idx) == MASKED_WATER) then + istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.true.) + else + istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.) + end if + else + istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.) + end if + else if (field_to_store%map%stagger == U) then + field_to_store%header%dim1(1) = we_mem_stag_s + field_to_store%header%dim1(2) = we_mem_stag_e + field_to_store%header%dim2(1) = sn_mem_s + field_to_store%header%dim2(2) = sn_mem_e + istat = remap_field(remap_info_u, mpas_field, wrf_field) + else if (field_to_store%map%stagger == V) then + field_to_store%header%dim1(1) = we_mem_s + field_to_store%header%dim1(2) = we_mem_e + field_to_store%header%dim2(1) = sn_mem_stag_s + field_to_store%header%dim2(2) = sn_mem_stag_e + istat = remap_field(remap_info_v, mpas_field, wrf_field) + else + call mprintf(.true.,ERROR, 'Cannot handle requested output stagger %i for MPAS field %s ...', & + i1=field_to_store%map%stagger, s1=trim(mpas_field % name)) + end if + + if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevels') then ! 3-d MPAS atmosphere field + field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name) + + ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name + if (len_trim(field_to_store % header % field) == 0) then + field_to_store % header % field = trim(mpas_field % name) + end if + + field_to_store % header % vertical_coord = 'num_mpas_levels' + do k=1,wrf_field % dimlens(1) + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = k + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:) + end if + + ! The u_field and v_field are used later by calling code to + ! determine which fields represent the x- and y-components of + ! horizonal wind velocity for the purposes of wind rotation + if (trim(mpas_field % name) == 'uReconstructZonal') then + call dup(field_to_store, u_field) + end if + if (trim(mpas_field % name) == 'uReconstructMeridional') then + call dup(field_to_store, v_field) + end if + + call storage_put_field(field_to_store) + end do + + else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevelsP1') then ! 3-d MPAS atmosphere field + field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name) + + ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name + if (len_trim(field_to_store % header % field) == 0) then + field_to_store % header % field = trim(mpas_field % name) + end if + + ! Handle surface level + field_to_store % header % vertical_coord = 'none' + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = 200100.0 + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:) + end if + + call storage_put_field(field_to_store) + + ! Handle all other layers + field_to_store % header % vertical_coord = 'num_mpas_levels' + do k=1,wrf_field % dimlens(1) - 1 + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = k + + ! Average to layer midpoint + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3r(k,:,:) + wrf_field % array3r(k+1,:,:)) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3d(k,:,:) + wrf_field % array3d(k+1,:,:)) + end if + + call storage_put_field(field_to_store) + end do + + ! Special case: zgrid field also provides SOILHGT + if (trim(mpas_field % name) == 'zgrid') then + field_to_store % header % field = 'SOILHGT' + + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = 200100.0 + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:) + end if + + call storage_put_field(field_to_store) + + do idx=1,num_entries + if (trim(fieldname(idx)) == 'SOILHGT') then + got_this_field(idx) = .true. + if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then + output_flags(idx) = flag_in_output(idx) + end if + exit + end if + end do + end if + + else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nSoilLevels') then ! 3-d MPAS soil field + + field_to_store % header % vertical_coord = 'none' + if (trim(mpas_field % name) == 'tslb') then + do k=1,wrf_field % dimlens(1) + if (k == 1) then + field_to_store % header % field = 'ST000010' + else if (k == 2) then + field_to_store % header % field = 'ST010040' + else if (k == 3) then + field_to_store % header % field = 'ST040100' + else if (k == 4) then + field_to_store % header % field = 'ST100200' + else + call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name)) + end if + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = 200100.0 + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:) + end if + + if (idx > 0) then + if (masked(idx) == MASKED_WATER) then + where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx) + end if + end if + + call storage_put_field(field_to_store) + end do + else if (trim(mpas_field % name) == 'smois') then + do k=1,wrf_field % dimlens(1) + if (k == 1) then + field_to_store % header % field = 'SM000010' + else if (k == 2) then + field_to_store % header % field = 'SM010040' + else if (k == 3) then + field_to_store % header % field = 'SM040100' + else if (k == 4) then + field_to_store % header % field = 'SM100200' + else + call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name)) + end if + allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3)) + do j=1,wrf_field % dimlens(3) + do i=1,wrf_field % dimlens(2) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = 200100.0 + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:) + end if + + if (idx > 0) then + if (masked(idx) == MASKED_WATER) then + where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx) + end if + end if + + call storage_put_field(field_to_store) + end do + else + call mprintf(.true.,WARN, 'Skipping unknown MPAS soil field %s ...', s1=trim(mpas_field % name)) + end if + + else if (wrf_field % ndims == 2) then ! 2-d MPAS field + field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name) + + ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name + if (len_trim(field_to_store % header % field) == 0) then + field_to_store % header % field = trim(mpas_field % name) + end if + + field_to_store % header % vertical_coord = 'none' + allocate(field_to_store % r_arr(wrf_field % dimlens(1), wrf_field % dimlens(2))) + allocate(field_to_store % valid_mask) + call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(1), wrf_field % dimlens(2)) + do j=1,wrf_field % dimlens(2) + do i=1,wrf_field % dimlens(1) + call bitarray_set(field_to_store % valid_mask, i, j) + end do + end do + field_to_store % header % vertical_level = 200100.0 + + if (wrf_field % xtype == FIELD_TYPE_REAL) then + field_to_store % r_arr(:,:) = wrf_field % array2r(:,:) + else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then + field_to_store % r_arr(:,:) = wrf_field % array2d(:,:) + end if + + if (idx > 0) then + if (masked(idx) == MASKED_WATER) then + where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx) + end if + end if + + call storage_put_field(field_to_store) + end if + + istat = free_target_field(wrf_field) + end if + + istat = scan_input_free_field(mpas_field) + end do + + if (scan_input_close(mpas_handle) /= 0) then + call mprintf(.true.,ERROR, 'Error closing %s with scan_input_close', s1=mpas_filename) + end if + + end subroutine process_mpas_fields + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: derive_mpas_fields + ! + ! Purpose: + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine derive_mpas_fields() + + use bitarray_module + use gridinfo_module + use interp_module + use interp_option_module + use read_met_module + use llxy_module + use misc_definitions_module + use module_debug + use output_module + use parallel_module + use rotate_winds_module + use storage_module + use scan_input + use mpas_mesh + use constants_module + + implicit none + + integer :: k + integer :: istatus + integer, pointer, dimension(:) :: theta_levels => null() + integer, pointer, dimension(:) :: pressure_levels => null() + real, pointer, dimension(:,:) :: exner => null() + type (fg_input) :: theta_field + type (fg_input) :: pressure_field + + ! + ! Derive TT from theta and pressure + ! + theta_field%header%time_dependent = .true. + theta_field%header%constant_field = .false. + theta_field%header%field = 'TT' + theta_field%header%vertical_level = 0 + nullify(theta_field%r_arr) + nullify(theta_field%valid_mask) + nullify(theta_field%modified_mask) + + pressure_field%header%time_dependent = .true. + pressure_field%header%constant_field = .false. + pressure_field%header%field = 'PRESSURE' + pressure_field%header%vertical_level = 0 + nullify(pressure_field%r_arr) + nullify(pressure_field%valid_mask) + nullify(pressure_field%modified_mask) + + call storage_get_levels(theta_field, theta_levels) + call storage_get_levels(pressure_field, pressure_levels) + + if (associated(theta_levels) .and. associated(pressure_levels)) then + call mprintf(.true.,LOGFILE, 'Computing MPAS TT field from theta and pressure...') + + if (size(theta_levels) == size(pressure_levels)) then + do k = 1, size(theta_levels) + call mprintf(.true.,LOGFILE, 'Computing TT at level %i for MPAS dataset', i1=theta_levels(k)) + theta_field % header % vertical_level = theta_levels(k) + call storage_get_field(theta_field, istatus) + if (istatus /= 0) then + call mprintf(.true.,ERROR, 'Could not get MPAS theta field at level %i', i1=theta_levels(k)) + return + end if + + pressure_field % header % vertical_level = pressure_levels(k) + call storage_get_field(pressure_field, istatus) + if (istatus /= 0) then + call mprintf(.true.,ERROR, 'Could not get MPAS pressure field at level %i', i1=theta_levels(k)) + return + end if + + ! Compute temperature + if (.not. associated(exner)) then + allocate(exner(size(theta_field % r_arr, 1), size(theta_field % r_arr, 2))) + end if + exner(:,:) = (pressure_field % r_arr(:,:) / P0)**(RD/CP) + theta_field % r_arr(:,:) = theta_field % r_arr(:,:) * exner(:,:) + + call storage_put_field(theta_field) + end do + else + call mprintf(.true.,ERROR, 'The MPAS theta and pressure fields do not have the same number of levels!') + end if + + deallocate(theta_levels) + deallocate(pressure_levels) + if (associated(exner)) then + deallocate(exner) + end if + end if + + end subroutine derive_mpas_fields + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_interp_masks ! ! Purpose: diff --git a/metgrid/src/remapper.F b/metgrid/src/remapper.F new file mode 100644 index 0000000..b18a2a0 --- /dev/null +++ b/metgrid/src/remapper.F @@ -0,0 +1,1341 @@ +module remapper + + use mpas_mesh + use scan_input, only : input_field_type, FIELD_TYPE_REAL, FIELD_TYPE_DOUBLE, FIELD_TYPE_INTEGER + use target_mesh + + integer, parameter :: max_queue_length = 2700 ! should be at least (earth circumference / minimum grid distance) + integer, parameter :: max_dictionary_size = 82000 ! should be at least (nCells/32) + + integer :: queue_head = 0 + integer :: queue_tail = 0 + integer :: queue_size = 0 + integer, dimension(0:max_queue_length-1) :: queue_array + + integer :: int_size = 32 + integer, dimension(max_dictionary_size) :: dictionary_array + + + type remap_info_type + integer :: method = -1 + type (mpas_mesh_type), pointer :: src_mesh + type (target_mesh_type), pointer :: dst_mesh + + ! For nearest-neighbor + integer, dimension(:,:), pointer :: nearestCell => null() + integer, dimension(:,:), pointer :: nearestVertex => null() + integer, dimension(:,:), pointer :: nearestEdge => null() + + ! For Wachspress interpolation + real, dimension(:,:,:), pointer :: cellWeights => null() + real, dimension(:,:,:), pointer :: vertexWeights => null() + real, dimension(:,:,:), pointer :: edgeWeights => null() + integer, dimension(:,:,:), pointer :: sourceCells => null() + integer, dimension(:,:,:), pointer :: sourceVertices => null() + integer, dimension(:,:,:), pointer :: sourceEdges => null() + + ! For masked interpolation + real, dimension(:,:,:), pointer :: cellMaskedWeights => null() + integer, dimension(:,:,:), pointer :: sourceMaskedCells => null() + end type remap_info_type + + + type target_field_type + character (len=64) :: name + integer :: ndims = -1 + integer :: xtype = -1 + logical :: isTimeDependent = .false. + integer, dimension(:), pointer :: dimlens => null() + character (len=64), dimension(:), pointer :: dimnames => null() + + ! Members to store field data + real :: array0r + real, dimension(:), pointer :: array1r => null() + real, dimension(:,:), pointer :: array2r => null() + real, dimension(:,:,:), pointer :: array3r => null() + real, dimension(:,:,:,:), pointer :: array4r => null() + double precision :: array0d + double precision, dimension(:), pointer :: array1d => null() + double precision, dimension(:,:), pointer :: array2d => null() + double precision, dimension(:,:,:), pointer :: array3d => null() + double precision, dimension(:,:,:,:), pointer :: array4d => null() + integer :: array0i + integer, dimension(:), pointer :: array1i => null() + integer, dimension(:,:), pointer :: array2i => null() + integer, dimension(:,:,:), pointer :: array3i => null() + integer, dimension(:,:,:,:), pointer :: array4i => null() + end type target_field_type + + + private :: nearest_cell, & + nearest_vertex, & + sphere_distance, & + mpas_arc_length, & + mpas_triangle_signed_area_sphere, & + mpas_wachspress_coordinates, & + convert_lx, & + index2d + + + contains + + + integer function remap_info_setup(src_mesh, dst_mesh, remap_info) result(stat) + + implicit none + + type (mpas_mesh_type), intent(in), target :: src_mesh + type (target_mesh_type), intent(in), target :: dst_mesh + type (remap_info_type), intent(out) :: remap_info + + integer :: idx + integer :: j + integer :: nn + integer :: ix, iy + integer :: irank + real :: sumWeights + real, dimension(:,:), allocatable :: vertCoords + real, dimension(3) :: pointInterp + + stat = 0 + + remap_info % method = 1 + remap_info % src_mesh => src_mesh + remap_info % dst_mesh => dst_mesh + + ! + ! For nearest neighbor + ! + allocate(remap_info % nearestCell(dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % nearestVertex(dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % nearestEdge(dst_mesh % nlon, dst_mesh % nlat)) + + irank = dst_mesh % irank + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % maxEdges, & + src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) + remap_info % nearestCell(ix, iy) = idx + end do + end do + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & + src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & + src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & + src_mesh % latVertex, src_mesh % lonVertex ) + remap_info % nearestVertex(ix, iy) = idx + end do + end do + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % nEdges, src_mesh % maxEdges, 2, & + src_mesh % nEdgesOnCell, src_mesh % edgesOnCell, & + src_mesh % cellsOnEdge, src_mesh % latCell, src_mesh % lonCell, & + src_mesh % latEdge, src_mesh % lonEdge ) + remap_info % nearestEdge(ix, iy) = idx + end do + end do + + + ! + ! For Wachspress interpolation + ! + allocate(vertCoords(3,3)) + + allocate(remap_info % cellWeights(3, dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % sourceCells(3, dst_mesh % nlon, dst_mesh % nlat)) + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & + src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & + src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & + src_mesh % latVertex, src_mesh % lonVertex ) + remap_info % sourceCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx) + pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) + do j=1,3 + vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), & + src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), & + 6371229.0) + end do + remap_info % cellWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp) + end do + end do + + deallocate(vertCoords) + + + allocate(vertCoords(3,src_mesh % maxEdges)) + + allocate(remap_info % vertexWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % sourceVertices(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % maxEdges, & + src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) + nn = src_mesh % nEdgesOnCell(idx) + remap_info % sourceVertices(:,ix,iy) = 1 + remap_info % sourceVertices(1:nn,ix,iy) = src_mesh % verticesOnCell(1:nn,idx) + pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) + do j=1,nn + vertCoords(:,j) = convert_lx(src_mesh % latVertex(src_mesh % verticesOnCell(j,idx)), & + src_mesh % lonVertex(src_mesh % verticesOnCell(j,idx)), & + 6371229.0) + end do + remap_info % vertexWeights(:,ix,iy) = 0.0 + remap_info % vertexWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp) + end do + end do + + deallocate(vertCoords) + + + allocate(vertCoords(3,src_mesh % maxEdges)) + + allocate(remap_info % edgeWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % sourceEdges(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % maxEdges, & + src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) + nn = src_mesh % nEdgesOnCell(idx) + remap_info % sourceEdges(:,ix,iy) = 1 + remap_info % sourceEdges(1:nn,ix,iy) = src_mesh % edgesOnCell(1:nn,idx) + pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) + do j=1,nn + vertCoords(:,j) = convert_lx(src_mesh % latEdge(src_mesh % edgesOnCell(j,idx)), & + src_mesh % lonEdge(src_mesh % edgesOnCell(j,idx)), & + 6371229.0) + end do + remap_info % edgeWeights(:,ix,iy) = 0.0 + remap_info % edgeWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp) + end do + end do + + deallocate(vertCoords) + + + ! + ! For masked interpolation + ! + allocate(vertCoords(3,3)) + + allocate(remap_info % cellMaskedWeights(3, dst_mesh % nlon, dst_mesh % nlat)) + allocate(remap_info % sourceMaskedCells(3, dst_mesh % nlon, dst_mesh % nlat)) + + idx = 1 + do iy=1,dst_mesh % nlat + do ix=1,dst_mesh % nlon + idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & + src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & + src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & + src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & + src_mesh % latVertex, src_mesh % lonVertex ) + remap_info % sourceMaskedCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx) + pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) + do j=1,3 + vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), & + src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), & + 6371229.0) + end do + remap_info % cellMaskedWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp) + sumWeights = 0.0 + do j=1,3 + if (src_mesh % landmask(remap_info % sourceMaskedCells(j,ix,iy)) == 1) then + sumWeights = sumWeights + remap_info % cellMaskedWeights(j,ix,iy) + else + remap_info % cellMaskedWeights(j,ix,iy) = 0.0 + end if + end do + if (sumWeights > 0.0) then + remap_info % cellMaskedWeights(:,ix,iy) = remap_info % cellMaskedWeights(:,ix,iy) / sumWeights + else + idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), & + remap_info % sourceMaskedCells(1,ix,iy), & + src_mesh % nCells, src_mesh % maxEdges, & + src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) + call search_for_cells(src_mesh % nCells, src_mesh % maxEdges, src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, & + src_mesh % landmask, dst_mesh % lats(index2d(irank,ix),iy), & + dst_mesh % lons(ix,index2d(irank,iy)), & + src_mesh % latCell, src_mesh % lonCell, idx, remap_info % sourceMaskedCells(:,ix,iy), & + remap_info % cellMaskedWeights(:,ix,iy)) + end if + end do + end do + + deallocate(vertCoords) + + end function remap_info_setup + + + integer function remap_info_free(remap_info) result(stat) + + implicit none + + type (remap_info_type), intent(inout) :: remap_info + + + stat = 0 + + remap_info % method = -1 + nullify(remap_info % src_mesh) + nullify(remap_info % dst_mesh) + + if (associated(remap_info % nearestCell)) then + deallocate(remap_info % nearestCell) + end if + if (associated(remap_info % nearestVertex)) then + deallocate(remap_info % nearestVertex) + end if + if (associated(remap_info % nearestEdge)) then + deallocate(remap_info % nearestEdge) + end if + if (associated(remap_info % cellWeights)) then + deallocate(remap_info % cellWeights) + end if + if (associated(remap_info % vertexWeights)) then + deallocate(remap_info % vertexWeights) + end if + if (associated(remap_info % edgeWeights)) then + deallocate(remap_info % edgeWeights) + end if + if (associated(remap_info % sourceCells)) then + deallocate(remap_info % sourceCells) + end if + if (associated(remap_info % sourceVertices)) then + deallocate(remap_info % sourceVertices) + end if + if (associated(remap_info % sourceEdges)) then + deallocate(remap_info % sourceEdges) + end if + if (associated(remap_info % cellMaskedWeights)) then + deallocate(remap_info % cellMaskedWeights) + end if + if (associated(remap_info % sourceMaskedCells)) then + deallocate(remap_info % sourceMaskedCells) + end if + + end function remap_info_free + + + logical function can_remap_field(field) + + implicit none + + type (input_field_type), intent(in) :: field + + integer :: decompDim + + + can_remap_field = .true. + + if (field % xtype /= FIELD_TYPE_INTEGER .and. & + field % xtype /= FIELD_TYPE_REAL .and. & + field % xtype /= FIELD_TYPE_DOUBLE) then + can_remap_field = .false. + return + end if + + if (field % ndims == 0 .or. & + (field % ndims == 1 .and. field % isTimeDependent)) then + can_remap_field = .false. + return + end if + + decompDim = field % ndims + if (field % isTimeDependent) then + decompDim = decompDim - 1 + end if + + if (trim(field % dimnames(decompDim)) /= 'nCells' .and. & + trim(field % dimnames(decompDim)) /= 'nVertices' .and. & + trim(field % dimnames(decompDim)) /= 'nEdges') then + + can_remap_field = .false. + return + end if + + end function can_remap_field + + + integer function remap_field_dryrun(remap_info, src_field, dst_field) result(stat) + + implicit none + + type (remap_info_type), intent(in) :: remap_info + type (input_field_type), intent(in) :: src_field + type (target_field_type), intent(out) :: dst_field + + integer :: idim + + stat = 0 + + dst_field % name = src_field % name + dst_field % xtype = src_field % xtype + if (src_field % isTimeDependent) then + ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon, + ! but the time dimension is not counted in the target field + dst_field % ndims = src_field % ndims + else + ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon + dst_field % ndims = src_field % ndims + 1 + end if + allocate(dst_field % dimnames(dst_field % ndims)) + allocate(dst_field % dimlens(dst_field % ndims)) + dst_field % isTimeDependent = src_field % isTimeDependent + do idim=1,dst_field % ndims-2 + dst_field % dimlens(idim) = src_field % dimlens(idim) + dst_field % dimnames(idim) = src_field % dimnames(idim) + end do + + dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon + dst_field % dimnames(dst_field % ndims-1) = 'nLons' + dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat + dst_field % dimnames(dst_field % ndims) = 'nLats' + + end function remap_field_dryrun + + + integer function remap_field(remap_info, src_field, dst_field, masked) result(stat) + + implicit none + + type (remap_info_type), intent(in) :: remap_info + type (input_field_type), intent(in) :: src_field + type (target_field_type), intent(out) :: dst_field + logical, intent(in), optional :: masked + + integer :: decompDim + integer :: idim + integer :: j + integer :: iy, ix + logical :: local_masked + integer, dimension(:,:), pointer :: nearestIndex + integer, dimension(:,:,:), pointer :: sourceNodes + real, dimension(:,:,:), pointer :: nodeWeights + + stat = 0 + + decompDim = src_field % ndims + if (src_field % isTimeDependent) then + decompDim = decompDim - 1 + end if + + local_masked = .false. + if (present(masked)) then + local_masked = masked + end if + + if (trim(src_field % dimnames(decompDim)) == 'nCells') then + nearestIndex => remap_info % nearestCell + if (local_masked) then + sourceNodes => remap_info % sourceMaskedCells + nodeWeights => remap_info % cellMaskedWeights + else + sourceNodes => remap_info % sourceCells + nodeWeights => remap_info % cellWeights + end if + else if (trim(src_field % dimnames(decompDim)) == 'nVertices') then + nearestIndex => remap_info % nearestVertex + sourceNodes => remap_info % sourceVertices + nodeWeights => remap_info % vertexWeights + else if (trim(src_field % dimnames(decompDim)) == 'nEdges') then + nearestIndex => remap_info % nearestEdge + sourceNodes => remap_info % sourceEdges + nodeWeights => remap_info % edgeWeights + else + write(0,*) 'Remap exception: unhandled decomposed dim' + stat = 1 + return + end if + + dst_field % name = src_field % name + dst_field % xtype = src_field % xtype + if (src_field % isTimeDependent) then + ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon, + ! but the time dimension is not counted in the target field + dst_field % ndims = src_field % ndims + else + ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon + dst_field % ndims = src_field % ndims + 1 + end if + allocate(dst_field % dimnames(dst_field % ndims)) + allocate(dst_field % dimlens(dst_field % ndims)) + dst_field % isTimeDependent = src_field % isTimeDependent + do idim=1,dst_field % ndims-2 + dst_field % dimlens(idim) = src_field % dimlens(idim) + dst_field % dimnames(idim) = src_field % dimnames(idim) + end do + dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon + dst_field % dimnames(dst_field % ndims-1) = 'nLons' + dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat + dst_field % dimnames(dst_field % ndims) = 'nLats' + + if (src_field % xtype == FIELD_TYPE_REAL) then + if (dst_field % ndims == 2) then + allocate(dst_field % array2r(dst_field % dimlens(1), & + dst_field % dimlens(2))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array2r(ix,iy) = src_field % array1r(nearestIndex(ix,iy)) + end do + end do +#else + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array2r(ix,iy) = 0.0 + do j=1,size(sourceNodes, dim=1) + dst_field % array2r(ix,iy) = dst_field % array2r(ix,iy) + & + nodeWeights(j,ix,iy) & + * src_field % array1r(sourceNodes(j,ix,iy)) + end do + end do + end do +#endif + else if (dst_field % ndims == 3) then + allocate(dst_field % array3r(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array3r(:,ix,iy) = src_field % array2r(:,nearestIndex(ix,iy)) + end do + end do +#else + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array3r(:,ix,iy) = 0.0 + do j=1,3 + dst_field % array3r(:,ix,iy) = dst_field % array3r(:,ix,iy) + & + nodeWeights(j,ix,iy) & + * src_field % array2r(:,sourceNodes(j,ix,iy)) + end do + end do + end do +#endif + else if (dst_field % ndims == 4) then + allocate(dst_field % array4r(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3), & + dst_field % dimlens(4))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array4r(:,:,ix,iy) = src_field % array3r(:,:,nearestIndex(ix,iy)) + end do + end do +#endif + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array4r(:,:,ix,iy) = 0.0 + do j=1,3 + dst_field % array4r(:,:,ix,iy) = dst_field % array4r(:,:,ix,iy) + & + remap_info % cellWeights(j,ix,iy) & + * src_field % array3r(:,:,remap_info % sourceCells(j,ix,iy)) + end do + end do + end do + else + write(0,*) 'Remap exception: unhandled dimension for real ', dst_field % ndims + end if + else if (src_field % xtype == FIELD_TYPE_DOUBLE) then + if (dst_field % ndims == 2) then + allocate(dst_field % array2d(dst_field % dimlens(1), & + dst_field % dimlens(2))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array2d(ix,iy) = src_field % array1d(nearestIndex(ix,iy)) + end do + end do +#endif + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array2d(ix,iy) = 0.0 + do j=1,size(sourceNodes, dim=1) + dst_field % array2d(ix,iy) = dst_field % array2d(ix,iy) + & + nodeWeights(j,ix,iy) & + * src_field % array1d(sourceNodes(j,ix,iy)) + end do + end do + end do + else if (dst_field % ndims == 3) then + allocate(dst_field % array3d(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array3d(:,ix,iy) = src_field % array2d(:,nearestIndex(ix,iy)) + end do + end do +#endif + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array3d(:,ix,iy) = 0.0 + do j=1,3 + dst_field % array3d(:,ix,iy) = dst_field % array3d(:,ix,iy) + & + remap_info % cellWeights(j,ix,iy) & + * src_field % array2d(:,remap_info % sourceCells(j,ix,iy)) + end do + end do + end do + else if (dst_field % ndims == 4) then + allocate(dst_field % array4d(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3), & + dst_field % dimlens(4))) +#ifdef NEAREST_NEIGHBOR + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array4d(:,:,ix,iy) = src_field % array3d(:,:,nearestIndex(ix,iy)) + end do + end do +#endif + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array4d(:,:,ix,iy) = 0.0 + do j=1,3 + dst_field % array4d(:,:,ix,iy) = dst_field % array4d(:,:,ix,iy) + & + remap_info % cellWeights(j,ix,iy) & + * src_field % array3d(:,:,remap_info % sourceCells(j,ix,iy)) + end do + end do + end do + else + write(0,*) 'Remap exception: unhandled dimension for dbl ', dst_field % ndims + end if + else if (src_field % xtype == FIELD_TYPE_INTEGER) then + if (dst_field % ndims == 2) then + allocate(dst_field % array2i(dst_field % dimlens(1), & + dst_field % dimlens(2))) + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array2i(ix,iy) = src_field % array1i(nearestIndex(ix,iy)) + end do + end do + else if (dst_field % ndims == 3) then + allocate(dst_field % array3i(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3))) + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array3i(:,ix,iy) = src_field % array2i(:,nearestIndex(ix,iy)) + end do + end do + else if (dst_field % ndims == 4) then + allocate(dst_field % array4i(dst_field % dimlens(1), & + dst_field % dimlens(2), & + dst_field % dimlens(3), & + dst_field % dimlens(4))) + do iy=1,size(nearestIndex, dim=2) + do ix=1,size(nearestIndex, dim=1) + dst_field % array4i(:,:,ix,iy) = src_field % array3i(:,:,nearestIndex(ix,iy)) + end do + end do + else + write(0,*) 'Remap exception: unhandled dimension for int ', dst_field % ndims + end if + else + write(0,*) 'Remap exception: unhandled type' + end if + + end function remap_field + + + integer function remap_get_target_latitudes(remap_info, lat_field) result(stat) + + implicit none + + type (remap_info_type), intent(in) :: remap_info + type (target_field_type), intent(out) :: lat_field + + real, parameter :: rad2deg = 90.0 / asin(1.0) + + stat = 0 + + + lat_field % name = 'lat' + lat_field % xtype = FIELD_TYPE_REAL + lat_field % ndims = 1 + lat_field % isTimeDependent = .false. + + allocate(lat_field % dimnames(lat_field % ndims)) + allocate(lat_field % dimlens(lat_field % ndims)) + + lat_field % dimnames(1) = 'nLats' + lat_field % dimlens(1) = remap_info % dst_mesh % nlat + + allocate(lat_field % array1r(lat_field % dimlens(1))) + lat_field % array1r(:) = remap_info % dst_mesh % lats(1,:) * rad2deg + + end function remap_get_target_latitudes + + + integer function remap_get_target_longitudes(remap_info, lon_field) result(stat) + + implicit none + + type (remap_info_type), intent(in) :: remap_info + type (target_field_type), intent(out) :: lon_field + + real, parameter :: rad2deg = 90.0 / asin(1.0) + + stat = 0 + + + lon_field % name = 'lon' + lon_field % xtype = FIELD_TYPE_REAL + lon_field % ndims = 1 + lon_field % isTimeDependent = .false. + + allocate(lon_field % dimnames(lon_field % ndims)) + allocate(lon_field % dimlens(lon_field % ndims)) + + lon_field % dimnames(1) = 'nLons' + lon_field % dimlens(1) = remap_info % dst_mesh % nlon + + allocate(lon_field % array1r(lon_field % dimlens(1))) + lon_field % array1r(:) = remap_info % dst_mesh % lons(:,1) * rad2deg + + end function remap_get_target_longitudes + + + integer function free_target_field(field) result(stat) + + implicit none + + type (target_field_type), intent(inout) :: field + + stat = 0 + + if (associated(field % dimlens)) then + deallocate(field % dimlens) + end if + if (associated(field % dimnames)) then + deallocate(field % dimnames) + end if + + if (associated(field % array1r)) then + deallocate(field % array1r) + end if + if (associated(field % array2r)) then + deallocate(field % array2r) + end if + if (associated(field % array3r)) then + deallocate(field % array3r) + end if + if (associated(field % array4r)) then + deallocate(field % array4r) + end if + + if (associated(field % array1d)) then + deallocate(field % array1d) + end if + if (associated(field % array2d)) then + deallocate(field % array2d) + end if + if (associated(field % array3d)) then + deallocate(field % array3d) + end if + if (associated(field % array4d)) then + deallocate(field % array4d) + end if + + if (associated(field % array1i)) then + deallocate(field % array1i) + end if + if (associated(field % array2i)) then + deallocate(field % array2i) + end if + if (associated(field % array3i)) then + deallocate(field % array3i) + end if + if (associated(field % array4i)) then + deallocate(field % array4i) + end if + + end function free_target_field + + + integer function nearest_cell(target_lat, target_lon, start_cell, nCells, maxEdges, & + nEdgesOnCell, cellsOnCell, latCell, lonCell) + + implicit none + + real, intent(in) :: target_lat, target_lon + integer, intent(in) :: start_cell + integer, intent(in) :: nCells, maxEdges + integer, dimension(nCells), intent(in) :: nEdgesOnCell + integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell + real, dimension(nCells), intent(in) :: latCell, lonCell + + integer :: i + integer :: iCell + integer :: current_cell + real :: current_distance, d + real :: nearest_distance + + nearest_cell = start_cell + current_cell = -1 + + do while (nearest_cell /= current_cell) + current_cell = nearest_cell + current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, & + target_lon, 1.0) + nearest_cell = current_cell + nearest_distance = current_distance + do i = 1, nEdgesOnCell(current_cell) + iCell = cellsOnCell(i,current_cell) + if (iCell <= nCells) then + d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0) + if (d < nearest_distance) then + nearest_cell = iCell + nearest_distance = d + end if + end if + end do + end do + + end function nearest_cell + + + integer function nearest_vertex( target_lat, target_lon, & + start_vertex, & + nCells, nVertices, maxEdges, vertexDegree, & + nEdgesOnCell, verticesOnCell, & + cellsOnVertex, latCell, lonCell, & + latVertex, lonVertex ) + + implicit none + + real, intent(in) :: target_lat, target_lon + integer, intent(in) :: start_vertex + integer, intent(in) :: nCells, nVertices, maxEdges, vertexDegree + integer, dimension(nCells), intent(in) :: nEdgesOnCell + integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell + integer, dimension(vertexDegree,nVertices), intent(in) :: cellsOnVertex + real, dimension(nCells), intent(in) :: latCell, lonCell + real, dimension(nVertices), intent(in) :: latVertex, lonVertex + + + integer :: i, cell1, cell2, cell3, iCell + integer :: iVtx + integer :: current_vertex + real :: cell1_dist, cell2_dist, cell3_dist + real :: current_distance, d + real :: nearest_distance + + nearest_vertex = start_vertex + current_vertex = -1 + + do while (nearest_vertex /= current_vertex) + current_vertex = nearest_vertex + current_distance = sphere_distance(latVertex(current_vertex), lonVertex(current_vertex), & + target_lat, target_lon, 1.0) + nearest_vertex = current_vertex + nearest_distance = current_distance + cell1 = cellsOnVertex(1,current_vertex) + cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0) + cell2 = cellsOnVertex(2,current_vertex) + cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0) + if (vertexDegree == 3) then + cell3 = cellsOnVertex(3,current_vertex) + cell3_dist = sphere_distance(latCell(cell3), lonCell(cell3), target_lat, target_lon, 1.0) + end if + if (vertexDegree == 3) then + if (cell1_dist < cell2_dist) then + if (cell1_dist < cell3_dist) then + iCell = cell1 + else + iCell = cell3 + end if + else + if (cell2_dist < cell3_dist) then + iCell = cell2 + else + iCell = cell3 + end if + end if + else + if (cell1_dist < cell2_dist) then + iCell = cell1 + else + iCell = cell2 + end if + end if + do i = 1, nEdgesOnCell(iCell) + iVtx = verticesOnCell(i,iCell) + d = sphere_distance(latVertex(iVtx), lonVertex(iVtx), target_lat, target_lon, 1.0) + if (d < nearest_distance) then + nearest_vertex = iVtx + nearest_distance = d + end if + end do + end do + + end function nearest_vertex + + + real function sphere_distance(lat1, lon1, lat2, lon2, radius) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a + ! sphere with given radius. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + real, intent(in) :: lat1, lon1, lat2, lon2, radius + + real :: arg1 + + arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + & + cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 ) + sphere_distance = 2.*radius*asin(arg1) + + end function sphere_distance + + + !*********************************************************************** + ! + ! function mpas_wachspress_coordinates + ! + !> \brief Compute the barycentric Wachspress coordinates for a polygon + !> \author Phillip Wolfram + !> \date 01/26/2015 + !> \details + !> Computes the barycentric Wachspress coordinates for a polygon with nVertices + !> points in R3, vertCoords for a particular pointInterp with normalized radius. + !> Follows Gillette, A., Rand, A., Bajaj, C., 2011. + !> Error estimates for generalized barycentric interpolation. + !> Advances in computational mathematics 37 (3), 417–439. + !> Optimized version of mpas_wachspress_coordinates uses optional cached B_i areas + !------------------------------------------------------------------------ + function mpas_wachspress_coordinates(nVertices, vertCoords, pointInterp, areaBin) + + implicit none + + ! input points + integer, intent(in) :: nVertices + real, dimension(3, nVertices), intent(in) :: vertCoords + real, dimension(3), intent(in) :: pointInterp + real, dimension(nVertices), optional, intent(in) :: areaBin + + ! output + real, dimension(nVertices) :: mpas_wachspress_coordinates + + ! computational intermediates + real, dimension(nVertices) :: wach ! The wachpress area-product + real :: wach_total ! The wachpress total weight + integer :: i, j ! Loop indices + integer :: im1, i0, ip1 ! im1 = (i-1), i0 = i, ip1 = (i+1) + + ! triangle areas to compute wachspress coordinate + real, dimension(nVertices) :: areaA + real, dimension(nVertices) :: areaB + + real :: radiusLocal + + radiusLocal = sqrt(sum(vertCoords(:,1)**2)) + + if (.not. present(areaBin)) then + ! compute areas + do i = 1, nVertices + ! compute first area B_i + ! get vertex indices + im1 = mod(nVertices + i - 2, nVertices) + 1 + i0 = mod(nVertices + i - 1, nVertices) + 1 + ip1 = mod(nVertices + i , nVertices) + 1 + + ! precompute B_i areas + ! always the same because B_i independent of xp,yp,zp + ! (COULD CACHE AND USE RESULT FROM ARRAY FOR FURTHER OPTIMIZATION) + areaB(i) = mpas_triangle_signed_area_sphere(vertCoords(:, im1), vertCoords(:, i0), vertCoords(:, ip1), radiusLocal) + end do + else + ! assign areas + do i = 1, nVertices + areaB(i) = areaBin(i) + end do + end if + + ! compute areas + do i = 1, nVertices + ! compute first area B_i + ! get vertex indices + im1 = mod(nVertices + i - 2, nVertices) + 1 + i0 = mod(nVertices + i - 1, nVertices) + 1 + ip1 = mod(nVertices + i , nVertices) + 1 + + ! compute A_ij areas + ! must be computed each time + areaA(i0) = mpas_triangle_signed_area_sphere(pointInterp, vertCoords(:, i0), vertCoords(:, ip1), radiusLocal) + + ! precomputed B_i areas, cached + end do + + + ! for each vertex compute wachpress coordinate + do i = 1, nVertices + wach(i) = areaB(i) + do j = (i + 1), (i + nVertices - 2) + i0 = mod(nVertices + j - 1, nVertices) + 1 + ! accumulate products for A_ij subareas + wach(i) = wach(i) * areaA(i0) + end do + end do + + ! get summed weights for normalization + wach_total = 0 + do i = 1, nVertices + wach_total = wach_total + wach(i) + end do + + ! compute lambda + mpas_wachspress_coordinates= 0.0 + do i = 1, nVertices + mpas_wachspress_coordinates(i) = wach(i)/wach_total + end do + + end function mpas_wachspress_coordinates + + + !*********************************************************************** + ! + ! routine mpas_triangle_signed_area_sphere + ! + !> \brief Calculates area of a triangle on a sphere + !> \author Matthew Hoffman + !> \date 13 January 2015 + !> \details + !> This routine calculates the area of a triangle on the surface of a sphere. + !> Uses the spherical analog of Heron's formula. + !> Copied from mesh generator. A CCW winding angle is positive. + !----------------------------------------------------------------------- + real function mpas_triangle_signed_area_sphere(a, b, c, radius) + + implicit none + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + real, dimension(3), intent(in) :: a, b, c !< Input: 3d (x,y,z) points forming the triangle in which to calculate the bary weights + real, intent(in) :: radius !< sphere radius + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + real :: ab, bc, ca, semiperim, tanqe + real, dimension(3) :: ablen, aclen, Dlen + + ab = mpas_arc_length(a(1), a(2), a(3), b(1), b(2), b(3))/radius + bc = mpas_arc_length(b(1), b(2), b(3), c(1), c(2), c(3))/radius + ca = mpas_arc_length(c(1), c(2), c(3), a(1), a(2), a(3))/radius + semiperim = 0.5 * (ab + bc + ca) + + tanqe = sqrt(max(0.0,tan(0.5 * semiperim) * tan(0.5 * (semiperim - ab)) & + * tan(0.5 * (semiperim - bc)) * tan(0.5 * (semiperim - ca)))) + + mpas_triangle_signed_area_sphere = 4.0 * radius * radius * atan(tanqe) + + ! computing correct signs (in similar fashion to mpas_sphere_angle) + ablen(1) = b(1) - a(1) + ablen(2) = b(2) - a(2) + ablen(3) = b(3) - a(3) + + aclen(1) = c(1) - a(1) + aclen(2) = c(2) - a(2) + aclen(3) = c(3) - a(3) + + dlen(1) = (ablen(2) * aclen(3)) - (ablen(3) * aclen(2)) + dlen(2) = -((ablen(1) * aclen(3)) - (ablen(3) * aclen(1))) + dlen(3) = (ablen(1) * aclen(2)) - (ablen(2) * aclen(1)) + + if ((Dlen(1)*a(1) + Dlen(2)*a(2) + Dlen(3)*a(3)) < 0.0) then + mpas_triangle_signed_area_sphere = -mpas_triangle_signed_area_sphere + end if + + end function mpas_triangle_signed_area_sphere + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! FUNCTION MPAS_ARC_LENGTH + ! + ! Returns the length of the great circle arc from A=(ax, ay, az) to + ! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the + ! same sphere centered at the origin. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + real function mpas_arc_length(ax, ay, az, bx, by, bz) + + implicit none + + real, intent(in) :: ax, ay, az, bx, by, bz + + real :: r, c + real :: cx, cy, cz + + cx = bx - ax + cy = by - ay + cz = bz - az + + r = sqrt(ax*ax + ay*ay + az*az) + c = sqrt(cx*cx + cy*cy + cz*cz) + + mpas_arc_length = r * 2.0 * asin(c/(2.0*r)) + + end function mpas_arc_length + + + function convert_lx(lat, lon, radius) result(vec) + + implicit none + + real, intent(in) :: lat, lon, radius + + real, dimension(3) :: vec + + vec(1) = radius * cos(lon) * cos(lat) + vec(2) = radius * sin(lon) * cos(lat) + vec(3) = radius * sin(lat) + + end function convert_lx + + + subroutine search_for_cells(nCells, maxEdges, nEdgesOnCell, cellsOnCell, cellMask, & + targetLat, targetLon, latCell, lonCell, & + startIdx, sourceMaskedCells, cellMaskedWeights) + + implicit none + + integer, intent(in) :: nCells + integer, intent(in) :: maxEdges + integer, dimension(nCells), intent(in) :: nEdgesOnCell + integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell + integer, dimension(nCells), intent(in) :: cellMask + real, intent(in) :: targetLat + real, intent(in) :: targetLon + real, dimension(nCells), intent(in) :: latCell + real, dimension(nCells), intent(in) :: lonCell + integer, intent(in) :: startIdx + integer, dimension(3), intent(inout) :: sourceMaskedCells + real, dimension(3), intent(inout) :: cellMaskedWeights + + integer :: i + integer :: scan_cell + integer :: neighbor + integer :: unscanned + logical :: no_more_queueing + real :: d, d_min + + ! + ! Reset data structures + ! + call queue_reset() + call dictionary_reset() + + no_more_queueing = .false. + unscanned = 0 + + d_min = 1000.0 + + sourceMaskedCells(:) = 1 + cellMaskedWeights(:) = 0.0 + + ! + ! Insert the origin cell into the queue for processing + ! + call queue_insert(startIdx) + call dictionary_insert(startIdx) + + do while (queue_size > 0) + scan_cell = queue_remove() + + ! + ! Each cell index removed from the queue represents a unique grid cell + ! that falls within the specified radius of the origin point + ! Here, we can do any processing we like for these cells + ! + if (cellMask(scan_cell) == 1) then + d = sphere_distance(targetLat, targetLon, latCell(scan_cell), lonCell(scan_cell), 1.0) + if (d < d_min) then + d_min = d + sourceMaskedCells(1) = scan_cell + cellMaskedWeights(1) = 1.0 + if (.not. no_more_queueing) then + unscanned = queue_size + end if + no_more_queueing = .true. + end if + end if + + ! + ! Add any neighbors of scan_cell within specified radius to the queue for processing + ! + if (.not. no_more_queueing .or. unscanned > 0) then + do i=1,nEdgesOnCell(scan_cell) + neighbor = cellsOnCell(i,scan_cell) + if (.not. dictionary_search(neighbor)) then + call queue_insert(neighbor) + call dictionary_insert(neighbor) + end if + end do + end if + + if (no_more_queueing) then + unscanned = unscanned - 1 + end if + + end do + + end subroutine search_for_cells + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Insert a new integer into the queue + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine queue_insert(i) + + implicit none + + integer, intent(in) :: i + + if (queue_size == max_queue_length) then + write(0,*) 'Error: queue overrun' + return + end if + queue_size = queue_size + 1 + queue_array(queue_head) = i + queue_head = mod(queue_head + 1, max_queue_length) + + end subroutine queue_insert + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Remove the oldest integer from the queue + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + integer function queue_remove() + + implicit none + + if (queue_size <= 0) then + write(0,*) 'Error: queue underrun' + queue_remove = -1 + return + end if + queue_size = queue_size - 1 + queue_remove = queue_array(queue_tail) + queue_tail = mod(queue_tail + 1, max_queue_length) + + end function queue_remove + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Reset the queue to an empty state + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine queue_reset() + + implicit none + + queue_head = 0 + queue_tail = 0 + queue_size = 0 + + end subroutine queue_reset + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Insert an integer into the dictionary + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine dictionary_insert(i) + + implicit none + + integer, intent(in) :: i + + integer :: n_integer + integer :: n_bit + + n_integer = ((i-1) / int_size) + 1 + n_bit = mod((i-1), int_size) + + if (n_integer > max_dictionary_size) then + write(0,*) 'Error: dictionary insert out of bounds' + return + end if + + dictionary_array(n_integer) = ibset(dictionary_array(n_integer), n_bit) + + end subroutine dictionary_insert + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Search for an integer in the dictionary + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + logical function dictionary_search(i) + + implicit none + + integer, intent(in) :: i + + integer :: n_integer + integer :: n_bit + + n_integer = ((i-1) / int_size) + 1 + n_bit = mod((i-1), int_size) + + if (n_integer > max_dictionary_size) then + write(0,*) 'Error: dictionary search out of bounds' + dictionary_search = .false. + return + end if + + dictionary_search = btest(dictionary_array(n_integer), n_bit) + + end function dictionary_search + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Reset the dictionary to an empty state + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine dictionary_reset() + + implicit none + + dictionary_array(:) = 0 + + end subroutine dictionary_reset + + + function index2d(irank, idx) result(i) + + implicit none + + integer, intent(in) :: irank, idx + + integer :: i + + i = irank * (idx - 1) + 1 + + end function index2d + +end module remapper diff --git a/metgrid/src/scan_input.F b/metgrid/src/scan_input.F new file mode 100644 index 0000000..030e86a --- /dev/null +++ b/metgrid/src/scan_input.F @@ -0,0 +1,604 @@ +module scan_input + + use netcdf + + type input_handle_type + integer :: ncid + integer :: num_vars = 0 + integer :: current_var = 0 + integer, dimension(:), pointer :: varids => null() + integer :: unlimited_dimid + end type input_handle_type + + type input_field_type + character (len=64) :: name + logical :: isTimeDependent = .false. + integer :: varid = -1 + integer :: xtype = -1 + integer :: ndims = -1 + character (len=64), dimension(:), pointer :: dimnames + integer, dimension(:), pointer :: dimlens + integer, dimension(:), pointer :: dimids + type (input_handle_type), pointer :: file_handle + + ! Members to store field data + real :: array0r + real, dimension(:), pointer :: array1r => null() + real, dimension(:,:), pointer :: array2r => null() + real, dimension(:,:,:), pointer :: array3r => null() + real, dimension(:,:,:,:), pointer :: array4r => null() + double precision :: array0d + double precision, dimension(:), pointer :: array1d => null() + double precision, dimension(:,:), pointer :: array2d => null() + double precision, dimension(:,:,:), pointer :: array3d => null() + double precision, dimension(:,:,:,:), pointer :: array4d => null() + integer :: array0i + integer, dimension(:), pointer :: array1i => null() + integer, dimension(:,:), pointer :: array2i => null() + integer, dimension(:,:,:), pointer :: array3i => null() + end type input_field_type + + integer, parameter :: FIELD_TYPE_UNSUPPORTED = -1, & + FIELD_TYPE_REAL = 1, & + FIELD_TYPE_DOUBLE = 2, & + FIELD_TYPE_INTEGER = 3, & + FIELD_TYPE_CHARACTER = 4 + + + contains + + + integer function scan_input_open(filename, handle, nRecords) result(stat) + + implicit none + + character (len=*), intent(in) :: filename + type (input_handle_type), intent(out) :: handle + integer, intent(out), optional :: nRecords + + integer :: i + + stat = 0 + + stat = nf90_open(trim(filename), NF90_NOWRITE, handle % ncid) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + stat = nf90_inquire(handle % ncid, nVariables=handle % num_vars) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + allocate(handle % varids(handle % num_vars)) + +#ifdef HAVE_NF90_INQ_VARIDS + stat = nf90_inq_varids(handle % ncid, handle % num_vars, handle % varids) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if +#else + ! Newer versions of the netCDF4 library (perhaps newer than 4.2.0?) + ! provide a function to return a list of all variable IDs in a file; if + ! we are using an older version of the netCDF library, we can apparently + ! assume that the variable IDs are numbered 1 through nVars. + ! See http://www.unidata.ucar.edu/software/netcdf/docs/tutorial_ncids.html + do i=1,handle % num_vars + handle % varids(i) = i + end do +#endif + + stat = nf90_inquire(handle % ncid, unlimitedDimId=handle % unlimited_dimid) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + if (present(nRecords)) then + stat = nf90_inquire_dimension(handle % ncid, handle % unlimited_dimid, len=nRecords) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + ! In case we have an input file that no time-varying records but + ! does have time-invariant fields, set nRecords = 1 so that we can + ! at least extract these fields + if ((nRecords == 0) .and. (handle % num_vars > 0)) then + nRecords = 1 + end if + end if + + handle % current_var = 1 + + end function scan_input_open + + + integer function scan_input_close(handle) result(stat) + + implicit none + + type (input_handle_type), intent(inout) :: handle + + + stat = 0 + + stat = nf90_close(handle % ncid) + if (stat /= NF90_NOERR) then + stat = 1 + end if + + if (associated(handle % varids)) then + deallocate(handle % varids) + end if + handle % current_var = 0 + + end function scan_input_close + + + integer function scan_input_rewind(handle) result(stat) + + implicit none + + type (input_handle_type), intent(inout) :: handle + + + stat = 0 + + handle % current_var = 1 + + end function scan_input_rewind + + + integer function scan_input_next_field(handle, field) result(stat) + + implicit none + + type (input_handle_type), intent(inout), target :: handle + type (input_field_type), intent(out) :: field + + integer :: idim + + + stat = 0 + + if (handle % current_var < 1 .or. handle % current_var > handle % num_vars) then + stat = 1 + return + end if + + field % varid = handle % varids(handle % current_var) + stat = nf90_inquire_variable(handle % ncid, field % varid, & + name=field % name, & + xtype=field % xtype, & + ndims=field % ndims) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + if (field % xtype == NF90_FLOAT) then + field % xtype = FIELD_TYPE_REAL + else if (field % xtype == NF90_DOUBLE) then + field % xtype = FIELD_TYPE_DOUBLE + else if (field % xtype == NF90_INT) then + field % xtype = FIELD_TYPE_INTEGER + else if (field % xtype == NF90_CHAR) then + field % xtype = FIELD_TYPE_CHARACTER + else + field % xtype = FIELD_TYPE_UNSUPPORTED + end if + + allocate(field % dimids(field % ndims)) + + stat = nf90_inquire_variable(handle % ncid, field % varid, & + dimids=field % dimids) + if (stat /= NF90_NOERR) then + stat = 1 + deallocate(field % dimids) + return + end if + + allocate(field % dimlens(field % ndims)) + allocate(field % dimnames(field % ndims)) + + do idim=1,field % ndims + stat = nf90_inquire_dimension(handle % ncid, field % dimids(idim), & + name=field % dimnames(idim), & + len=field % dimlens(idim)) + if (field % dimids(idim) == handle % unlimited_dimid) then + field % isTimeDependent = .true. + end if + end do + + field % file_handle => handle + + handle % current_var = handle % current_var + 1 + + end function scan_input_next_field + + + integer function scan_input_for_field(handle, fieldname, field) result(stat) + + implicit none + + type (input_handle_type), intent(inout), target :: handle + character (len=*), intent(in) :: fieldname + type (input_field_type), intent(out) :: field + + integer :: idim + + stat = 0 + + stat = nf90_inq_varid(handle % ncid, trim(fieldname), field % varid) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + stat = nf90_inquire_variable(handle % ncid, field % varid, & + name=field % name, & + xtype=field % xtype, & + ndims=field % ndims) + if (stat /= NF90_NOERR) then + stat = 1 + return + end if + + if (field % xtype == NF90_FLOAT) then + field % xtype = FIELD_TYPE_REAL + else if (field % xtype == NF90_DOUBLE) then + field % xtype = FIELD_TYPE_DOUBLE + else if (field % xtype == NF90_INT) then + field % xtype = FIELD_TYPE_INTEGER + else if (field % xtype == NF90_CHAR) then + field % xtype = FIELD_TYPE_CHARACTER + else + field % xtype = FIELD_TYPE_UNSUPPORTED + end if + + allocate(field % dimids(field % ndims)) + + stat = nf90_inquire_variable(handle % ncid, field % varid, & + dimids=field % dimids) + if (stat /= NF90_NOERR) then + stat = 1 + deallocate(field % dimids) + return + end if + + allocate(field % dimlens(field % ndims)) + allocate(field % dimnames(field % ndims)) + + do idim=1,field % ndims + stat = nf90_inquire_dimension(handle % ncid, field % dimids(idim), & + name=field % dimnames(idim), & + len=field % dimlens(idim)) + if (field % dimids(idim) == handle % unlimited_dimid) then + field % isTimeDependent = .true. + end if + end do + + field % file_handle => handle + + end function scan_input_for_field + + + integer function scan_input_read_field(field, frame) result(stat) + + implicit none + + type (input_field_type), intent(inout) :: field + integer, intent(in), optional :: frame + + integer :: local_frame + integer, dimension(5) :: start, count + real, dimension(1) :: temp1r + double precision, dimension(1) :: temp1d + integer, dimension(1) :: temp1i + + + stat = 0 + + local_frame = 1 + if (present(frame)) then + local_frame = frame + end if + + if (field % xtype == FIELD_TYPE_REAL) then + if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = local_frame + count(1) = 1 + stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1r, & + start=start(1:1), count=count(1:1)) + field % array0r = temp1r(1) + else + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0r) + end if + else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = local_frame + count(2) = 1 + allocate(field % array1r(count(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1r, & + start=start(1:2), count=count(1:2)) + else + allocate(field % array1r(field%dimlens(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1r) + end if + else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = local_frame + count(3) = 1 + allocate(field % array2r(count(1),count(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2r, & + start=start(1:3), count=count(1:3)) + else + allocate(field % array2r(field%dimlens(1),field%dimlens(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2r) + end if + else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = 1 + count(3) = field % dimlens(3) + start(4) = local_frame + count(4) = 1 + allocate(field % array3r(count(1),count(2),count(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3r, & + start=start(1:4), count=count(1:4)) + else + allocate(field % array3r(field%dimlens(1),field%dimlens(2),field%dimlens(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3r) + end if + else if (field % ndims == 4 .or. (field % ndims == 5 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = 1 + count(3) = field % dimlens(3) + start(4) = 1 + count(4) = field % dimlens(4) + start(5) = local_frame + count(5) = 1 + allocate(field % array4r(count(1),count(2),count(3),count(4))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4r, & + start=start(1:5), count=count(1:5)) + else + allocate(field % array4r(field%dimlens(1),field%dimlens(2),field%dimlens(3),field%dimlens(4))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4r) + end if + end if + else if (field % xtype == FIELD_TYPE_DOUBLE) then + if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = local_frame + count(1) = 1 + stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1d, & + start=start(1:1), count=count(1:1)) + field % array0d = temp1d(1) + else + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0d) + end if + else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = local_frame + count(2) = 1 + allocate(field % array1d(count(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1d, & + start=start(1:2), count=count(1:2)) + else + allocate(field % array1d(field%dimlens(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1d) + end if + else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = local_frame + count(3) = 1 + allocate(field % array2d(count(1),count(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2d, & + start=start(1:3), count=count(1:3)) + else + allocate(field % array2d(field%dimlens(1),field%dimlens(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2d) + end if + else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = 1 + count(3) = field % dimlens(3) + start(4) = local_frame + count(4) = 1 + allocate(field % array3d(count(1),count(2),count(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3d, & + start=start(1:4), count=count(1:4)) + else + allocate(field % array3d(field%dimlens(1),field%dimlens(2),field%dimlens(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3d) + end if + else if (field % ndims == 4 .or. (field % ndims == 5 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = 1 + count(3) = field % dimlens(3) + start(4) = 1 + count(4) = field % dimlens(4) + start(5) = local_frame + count(5) = 1 + allocate(field % array4d(count(1),count(2),count(3),count(4))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4d, & + start=start(1:5), count=count(1:5)) + else + allocate(field % array4d(field%dimlens(1),field%dimlens(2),field%dimlens(3),field%dimlens(4))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4d) + end if + end if + + else if (field % xtype == FIELD_TYPE_INTEGER) then + if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = local_frame + count(1) = 1 + stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1i, & + start=start(1:1), count=count(1:1)) + field % array0i = temp1i(1) + else + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0i) + end if + else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = local_frame + count(2) = 1 + allocate(field % array1i(count(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1i, & + start=start(1:2), count=count(1:2)) + else + allocate(field % array1i(field%dimlens(1))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1i) + end if + else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = local_frame + count(3) = 1 + allocate(field % array2i(count(1),count(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2i, & + start=start(1:3), count=count(1:3)) + else + allocate(field % array2i(field%dimlens(1),field%dimlens(2))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2i) + end if + else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then + if (field % isTimeDependent) then + start(1) = 1 + count(1) = field % dimlens(1) + start(2) = 1 + count(2) = field % dimlens(2) + start(3) = 1 + count(3) = field % dimlens(3) + start(4) = local_frame + count(4) = 1 + allocate(field % array3i(count(1),count(2),count(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3i, & + start=start(1:4), count=count(1:4)) + else + allocate(field % array3i(field%dimlens(1),field%dimlens(2),field%dimlens(3))) + stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3i) + end if + end if + + else if (field % xtype == FIELD_TYPE_CHARACTER) then + write(0,*) ' ' + write(0,*) 'Processing of character fields is not supported; skipping read of field '//trim(field % name) + write(0,*) ' ' + + else + write(0,*) ' ' + write(0,*) 'Unsupported type; skipping read of field '//trim(field % name) + write(0,*) ' ' + end if + + if (stat /= NF90_NOERR) then + write(0,*) ' ' + write(0,*) 'NetCDF error: reading '//trim(field % name)//' returned ', stat + write(0,*) ' ' + stat = 1 + else + stat = 0 + end if + + end function scan_input_read_field + + + integer function scan_input_free_field(field) result(stat) + + implicit none + + type (input_field_type), intent(inout) :: field + + + stat = 0 + + if (associated(field % dimids)) then + deallocate(field % dimids) + end if + if (associated(field % dimlens)) then + deallocate(field % dimlens) + end if + if (associated(field % dimnames)) then + deallocate(field % dimnames) + end if + + if (associated(field % array1r)) then + deallocate(field % array1r) + end if + if (associated(field % array2r)) then + deallocate(field % array2r) + end if + if (associated(field % array3r)) then + deallocate(field % array3r) + end if + if (associated(field % array4r)) then + deallocate(field % array4r) + end if + + if (associated(field % array1d)) then + deallocate(field % array1d) + end if + if (associated(field % array2d)) then + deallocate(field % array2d) + end if + if (associated(field % array3d)) then + deallocate(field % array3d) + end if + if (associated(field % array4d)) then + deallocate(field % array4d) + end if + + if (associated(field % array1i)) then + deallocate(field % array1i) + end if + if (associated(field % array2i)) then + deallocate(field % array2i) + end if + if (associated(field % array3i)) then + deallocate(field % array3i) + end if + + nullify(field % file_handle) + + end function scan_input_free_field + +end module scan_input diff --git a/metgrid/src/target_mesh.F b/metgrid/src/target_mesh.F new file mode 100644 index 0000000..91023c4 --- /dev/null +++ b/metgrid/src/target_mesh.F @@ -0,0 +1,214 @@ +module target_mesh + + type target_mesh_type + logical :: valid = .false. + integer :: irank = 0 + integer :: nLat = 0 + integer :: nLon = 0 + real :: startLat = 0.0 + real :: endLat = 0.0 + real :: startLon = 0.0 + real :: endLon = 0.0 + real, dimension(:,:), pointer :: lats => null() + real, dimension(:,:), pointer :: lons => null() + end type target_mesh_type + + + contains + + + integer function target_mesh_setup(mesh, lat2d, lon2d) result(stat) + + implicit none + + type (target_mesh_type), intent(out) :: mesh + real, dimension(:,:), target, optional :: lat2d + real, dimension(:,:), target, optional :: lon2d + + integer :: i, j + integer :: iostatus + integer :: eqIdx + real :: delta + logical :: exists + character (len=64) :: spec + real, parameter :: pi_const = 2.0 * asin(1.0) + + stat = 0 + + ! + ! If 2-d arrays of latitude and longitude are provided, we can just + ! point to those arrays rather than generate lat/lon values based on + ! a specified target domain + ! + if (present(lat2d) .and. present(lon2d)) then + + mesh % irank = 1 + mesh % nLat = size(lat2d,2) + mesh % nLon = size(lon2d,1) + mesh % lats => lat2d + mesh % lons => lon2d + mesh % valid = .true. + + return + end if + + + ! + ! Try to parse nLat, nLon from target_domain file + ! + inquire(file='target_domain', exist=exists) + if (exists) then + write(0,*) ' ' + write(0,*) 'Reading target domain specification from file ''target_domain''' + + mesh % startLat = -90.0 + mesh % endLat = 90.0 + mesh % startLon = -180.0 + mesh % endLon = 180.0 + mesh % nLat = 360 + mesh % nLon = 720 + + open(22, file='target_domain', form='formatted') + read(22,fmt='(a)',iostat=iostatus) spec + j = 1 + do while (iostatus >= 0) + call despace(spec) + eqIdx = index(spec, '=') + if (eqIdx /= 0) then + if (spec(1:eqIdx-1) == 'nlat') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % nLat + write(0,*) 'Setting nlat = ', mesh % nLat + else if (spec(1:eqIdx-1) == 'nlon') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % nLon + write(0,*) 'Setting nlon = ', mesh % nLon + else if (spec(1:eqIdx-1) == 'startlat') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % startLat + write(0,*) 'Setting startlat = ', mesh % startLat + else if (spec(1:eqIdx-1) == 'endlat') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % endLat + write(0,*) 'Setting endlat = ', mesh % endLat + else if (spec(1:eqIdx-1) == 'startlon') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % startLon + write(0,*) 'Setting startlon = ', mesh % startLon + else if (spec(1:eqIdx-1) == 'endlon') then + read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % endLon + write(0,*) 'Setting endlon = ', mesh % endLon + else + write(0,*) 'Unrecognized keyword on line ', j, ' of file ''target_domain'': '//spec(1:eqIdx-1) + stat = 1 + close(22) + return + end if + else + write(0,*) 'Syntax error on line ', j, ' of file ''target_domain'': ''='' not found' + stat = 1 + close(22) + return + end if + read(22,fmt='(a)',iostat=iostatus) spec + j = j + 1 + end do + close(22) + else + write(0,*) ' ' + write(0,*) 'Target domain specification file ''target_domain'' not found.' + write(0,*) 'Default 0.5-degree global target domain will be used.' + write(0,*) ' ' + + mesh % startLat = -90.0 + mesh % endLat = 90.0 + mesh % startLon = -180.0 + mesh % endLon = 180.0 + mesh % nLat = 360 + mesh % nLon = 720 + end if + + + allocate(mesh % lats(1, mesh % nLat)) + allocate(mesh % lons(mesh % nLon, 1)) + + delta = (mesh % endLat - mesh % startLat) / real(mesh % nLat) + do i=0,mesh % nLat-1 + mesh % lats(1,i+1) = mesh % startLat + 0.5 * delta + real(i) * delta + mesh % lats(1,i+1) = mesh % lats(1,i+1) * pi_const / 180.0 + end do + + delta = (mesh % endLon - mesh % startLon) / real(mesh % nLon) + do i=0,mesh % nLon-1 + mesh % lons(i+1,1) = mesh % startLon + 0.5 * delta + real(i) * delta + mesh % lons(i+1,1) = mesh % lons(i+1,1) * pi_const / 180.0 + end do + + mesh % valid = .true. + + end function target_mesh_setup + + + integer function target_mesh_free(mesh) result(stat) + + implicit none + + type (target_mesh_type), intent(inout) :: mesh + + + stat = 0 + + mesh % valid = .false. + mesh % nLat = 0 + mesh % nLon = 0 + mesh % startLat = 0.0 + mesh % endLat = 0.0 + mesh % startLon = 0.0 + mesh % endLon = 0.0 + + ! + ! When irank == 0, we allocated the lats and lons arrays + ! internally and should therefore deallocate them + ! + if (mesh % irank == 0) then + if (associated(mesh % lats)) deallocate(mesh % lats) + if (associated(mesh % lons)) deallocate(mesh % lons) + end if + + end function target_mesh_free + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: despace + ! + ! Purpose: Remove all space and tab characters from a string, thus + ! compressing the string to the left. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine despace(string) + + implicit none + + ! Arguments + character (len=*), intent(inout) :: string + + ! Local variables + integer :: i, j, length, iquoted + + length = len(string) + + iquoted = 0 + j = 1 + do i=1,length + ! Check for a quote mark + if (string(i:i) == '"' .or. string(i:i) == '''') iquoted = mod(iquoted+1,2) + + ! Check for non-space, non-tab character, or if we are inside quoted + ! text + if ((string(i:i) /= ' ' .and. string(i:i) /= achar(9)) .or. iquoted == 1) then + string(j:j) = string(i:i) + j = j + 1 + end if + end do + + do i=j,length + string(i:i) = ' ' + end do + + end subroutine despace + +end module target_mesh -- 2.11.4.GIT