sflux_subs5.f90

河口模型 使用模拟盐水入侵、热量扩散等等 河口模型 使用模拟盐水入侵、热量扩散

            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) t
10      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_label
20      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