i_elem = i_elem + 1            elem_nodes_in(i_elem,1) = node_num_in(i,j)            elem_nodes_in(i_elem,2) = node_num_in(i+1,j)            elem_nodes_in(i_elem,3) = node_num_in(i+1,j+1)          enddo        enddo      return      end!-----------------------------------------------------------------------! note: this subroutine has been modified from the original version! to account for full-sized data arrays (and still reduced size integer! arrays)      subroutine interp_data (data_out, data_in, weight, &     &                        i_elem_ae_min, elem_nodes_in, &     &                        node_i, node_j, &     &                        n_nodes_in, n_nodes_out, n_elems_in, &     &                        max_ni, max_nj, max_nodes_out)        implicit none        integer n_nodes_in, n_nodes_out, n_elems_in        integer max_ni, max_nj, max_nodes_out        integer i_elem_ae_min(n_nodes_out)        integer elem_nodes_in(n_elems_in,3)        integer node_i(n_nodes_in), node_j(n_nodes_in)        real*8 data_in(max_ni,max_nj), data_out(max_nodes_out)        real*8 weight(max_nodes_out,3)        integer i_node, i_elem, i1, j1, i2, j2, i3, j3! loop over the output nodes        do i_node = 1, n_nodes_out! get the locations of the nodes for the surrounding element on the! input grid          i_elem = i_elem_ae_min(i_node)          i1 = node_i(elem_nodes_in(i_elem,1))          j1 = node_j(elem_nodes_in(i_elem,1))          i2 = node_i(elem_nodes_in(i_elem,2))          j2 = node_j(elem_nodes_in(i_elem,2))          i3 = node_i(elem_nodes_in(i_elem,3))          j3 = node_j(elem_nodes_in(i_elem,3))! the data on the output grid is simply the weighted data from the input! grid          data_out(i_node) = data_in(i1,j1) * weight(i_node,1) &     &                     + data_in(i2,j2) * weight(i_node,2) &     &                     + data_in(i3,j3) * weight(i_node,3)        enddo      return      end!-----------------------------------------------------------------------      subroutine get_albedo (albedo, num_nodes, max_nodes)        implicit none        integer max_nodes, num_nodes, i_node        real*8 albedo(max_nodes)        do i_node = 1, num_nodes          albedo(i_node) = 0.06        enddo      return      end!-----------------------------------------------------------------------      real*8 function psi_m(zeta)        implicit none        real*8 zeta, chi, half_pi        half_pi = 2.0 * atan(1.0)        chi = (1.0 - 16.0 * zeta)**0.25        psi_m = 2.0 * log( 0.5 * (1.0 + chi) ) + &     &          log( 0.5 * (1.0 + chi*chi) ) - &     &          2.0 * atan(chi) + half_pi      return      end!-----------------------------------------------------------------------      real*8 function psi_h(zeta)        implicit none        real*8 zeta, chi        chi = (1.0 - 16.0 * zeta)**0.25        psi_h = 2.0 * log( 0.5 * (1.0 + chi*chi) )      return      end!-----------------------------------------------------------------------      subroutine read_2d_arr(in_file, data, data_label, t, ni, nj)        implicit none        integer ni, nj, edges(2)        real*4 data(ni,nj), t        integer ret, read_only, sfstart, sfend, sd_id, sds_id        integer start(2), stride(2), sds_index, sfn2index        integer sfrdata, sfendacc, sfselect        parameter (read_only = 1)        character data_label*20, in_file*50, dat_time_label*33! open in_file in read only mode        sd_id = sfstart(in_file, read_only)! create the data-time label, which we'll use as the name        call get_label (dat_time_label, data_label, t)! find index for this data set        if (t .lt. 0.0) then          sds_index = sfn2index(sd_id, data_label)        else          sds_index = sfn2index(sd_id, dat_time_label)        endif! find the id for this index        sds_id = sfselect(sd_id, sds_index)! set up start, stride, and edges to read in the entire dataset        start(1) = 0        start(2) = 0        stride(1) = 1        stride(2) = 1        edges(1) = ni        edges(2) = nj! read in the dataset        ret = sfrdata(sds_id, start, stride, edges, data)        call checkret(ret)! close access to the dataset        ret = sfendacc(sds_id)        call checkret(ret)! close access to the file        ret = sfend(sd_id)        call checkret(ret)      return      end!-----------------------------------------------------------------------      subroutine read_vec_int(in_file, data, data_label, t, ni)        implicit none        integer ni        real*4 t        integer data(ni)        integer ret, read_only, sfstart, sfend, sd_id, sds_id        integer sds_index, sfn2index, sfrdata, sfendacc, sfselect        parameter (read_only = 1)        character data_label*20, in_file*50, dat_time_label*33        ! open in_file in read only mode        sd_id = sfstart(in_file, read_only)! create the data-time label, which we'll use as the name        call get_label (dat_time_label, data_label, t)        ! find index for this data set        if (t .lt. 0.0) then          sds_index = sfn2index(sd_id, data_label)        else          sds_index = sfn2index(sd_id, dat_time_label)        endif! find the id for this index        sds_id = sfselect(sd_id, sds_index)! read in the dataset        ret = sfrdata(sds_id, 0, 1, ni, data)        call checkret(ret)! close access to the dataset        ret = sfendacc(sds_id)        call checkret(ret)! close access to the file        ret = sfend(sd_id)        call checkret(ret)      return      end!-----------------------------------------------------------------------      subroutine read_scalar(in_file, data, data_label, t)        implicit none        real*4 data, t        integer ret, read_only, sfstart, sfend, sd_id, sds_id        integer sds_index, sfn2index, sfrdata, sfendacc, sfselect        parameter (read_only = 1)        character data_label*20, in_file*50, dat_time_label*33        ! open in_file in read only mode        sd_id = sfstart(in_file, read_only)! create the data-time label, which we'll use as the name        call get_label (dat_time_label, data_label, t)        ! find index for this data set        if (t .lt. 0.0) then          sds_index = sfn2index(sd_id, data_label)        else          sds_index = sfn2index(sd_id, dat_time_label)        endif! find the id for this index        sds_id = sfselect(sd_id, sds_index)! read in the dataset        ret = sfrdata(sds_id, 0, 1, 1, data)        call checkret(ret)! close access to the dataset        ret = sfendacc(sds_id)        call checkret(ret)! close access to the file        ret = sfend(sd_id)        call checkret(ret)      return      end!-----------------------------------------------------------------------      subroutine get_label (dat_time_label, data_label, t)        implicit none        character data_label*20, dat_time_label*33, time_label*12        character zero_short*12        real*4 t        integer i        logical nonzero        parameter (zero_short = '      0.    ')! create the time_label        write(time_label,10) t10      format(g12.6)! since different platforms use different formats for zero, we'll use a! standard form        nonzero = .false.        do i = 1, 12          if (time_label(i:i) .ge. '1' .and. &     &        time_label(i:i) .le. '9') nonzero = .true.        enddo                if (.not. nonzero) time_label = zero_short! create the data-time label        write(dat_time_label,20) data_label, ' ', time_label20      format(a20,a1,a12)      return      end!-----------------------------------------------------------------------      subroutine checkret (ret)        implicit none        integer ret                if (ret .ne. 0) then          write(*,*) 'nonzero HDF return code!'          write(*,*) 'ret = ', ret          write(*,*)          write(11,*) 'nonzero HDF return code!'          write(11,*) 'ret = ', ret          write(11,*)          stop!       else!         write(*,*) 'HDF file access ok. . .'        endif              return      end!-----------------------------------------------------------------------      subroutine get_weight (x_in, y_in, x_out, y_out, &     &                       elem_nodes_in, node_i, node_j, &     &                       max_ni, max_nj, &     &                       n_elems_in, n_nodes_in, &     &                       n_nodes_out, &     &                       max_nodes_out, &     &                       i_elem_ae_min, &     &                       area_in, weight)        implicit none        integer max_ni, max_nj, n_elems_in, n_nodes_in        integer n_nodes_out, max_nodes_out        integer node_i(n_nodes_in), node_j(n_nodes_in)        integer elem_nodes_in(n_elems_in,3)        real*8 x_in(max_ni,max_nj), y_in(max_ni,max_nj)        real*8 x_out(n_nodes_out), y_out(n_nodes_out)        real*8 area_in(n_elems_in)        real*8 weight(max_nodes_out,3)        integer i_elem, i_node, i_elem_ae_min(n_nodes_out)        integer i1, j1, i2, j2, i3, j3        real*8 x1, y1, x2, y2, x3, y3, x4, y4, a1, a2, a3, aa, ae        real*8 ae_min! calculate and store the areas of the input grid elements        do i_elem = 1, n_elems_in          i1 = node_i(elem_nodes_in(i_elem,1))          j1 = node_j(elem_nodes_in(i_elem,1))          x1 = x_in(i1,j1)          y1 = y_in(i1,j1)          i2 = node_i(elem_nodes_in(i_elem,2))          j2 = node_j(elem_nodes_in(i_elem,2))          x2 = x_in(i2,j2)          y2 = y_in(i2,j2)          i3 = node_i(elem_nodes_in(i_elem,3))          j3 = node_j(elem_nodes_in(i_elem,3))          x3 = x_in(i3,j3)          y3 = y_in(i3,j3)          area_in(i_elem) = 0.5 * &     &                      ( (x1-x3)*(y2-y3) + (x3-x2)*(y1-y3) )        enddo! now loop over the nodes of the output grid, searching for the! surrounding elements on the input grid        do i_node = 1, n_nodes_out          ae_min = 1.0e25          i_elem_ae_min(i_node) = 0          do i_elem = 1, n_elems_in! get the locations of the nodes for this element on the input grid            i1 = node_i(elem_nodes_in(i_elem,1))            j1 = node_j(elem_nodes_in(i_elem,1))            x1 = x_in(i1,j1)            y1 = y_in(i1,j1)            i2 = node_i(elem_nodes_in(i_elem,2))            j2 = node_j(elem_nodes_in(i_elem,2))            x2 = x_in(i2,j2)            y2 = y_in(i2,j2)            i3 = node_i(elem_nodes_in(i_elem,3))            j3 = node_j(elem_nodes_in(i_elem,3))            x3 = x_in(i3,j3)            y3 = y_in(i3,j3)! get the locations o