sflux_subs5.f90

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                                       !
!                           Heat exchange sub-model of ELCIRC				!
!                       	Version 5 (Sept. 05, 2003)                              !
!                                                                                       !
!                 Center for Coastal and Land-Margin Research                           !
!             Department of Environmental Science and Engineering                       !
!                   OGI School of Science and Engineering,                              !
!                     Oregon Health & Science University                                !
!                       Beaverton, Oregon 97006, USA                                    !
!                                                                                       !
!                   Scientific direction: Antonio Baptista                              !
!                   Code development: Mike A. Zulauf                  			!
!                                                                                       !
!               Copyright 1999-2003 Oregon Health and Science University                !
!                              All Rights Reserved                                      !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  

!-----------------------------------------------------------------------
! Note: the following global variables (from module global) are used in
!       this code.  This list does not include variables passed in as
!       arguments. . .
!
!       mnp
!       np
!       x
!       y
!       kfp
!       uu2
!       vv2
!       tnd
!       snd
!-----------------------------------------------------------------------
      subroutine get_wind (time, u_air_node, v_air_node, p_air_node, &
     &                     t_air_node, q_air_node)

!       implicit none
        use global
        implicit real*8(a-h,o-z),integer(i-n)

! define some new names for things in header file
        integer max_nodes
        parameter (max_nodes = mnp)

! input/output variables
        real*8 u_air_node(max_nodes), v_air_node(max_nodes)
        real*8 t_air_node(max_nodes), q_air_node(max_nodes)
        real*8 p_air_node(max_nodes), time

! local variables
        integer max_ni, max_nj, max_times, max_files
        integer max_elems_in, max_nodes_in
        parameter (max_ni = 1024)
        parameter (max_nj = 1024)
        parameter (max_elems_in = (max_ni-1) * (max_nj-1) * 2)
        parameter (max_nodes_in = max_ni * max_nj)
        parameter (max_times = 1000)
        parameter (max_files = 999)
        integer num_nodes
        integer in_elem_to_out_node_1(max_nodes)
        integer in_elem_to_out_node_2(max_nodes)
        integer ni_1, nj_1, num_times_1, num_nodes_in_1, num_elems_in_1
        integer node_i_1(max_nodes_in), node_j_1(max_nodes_in)
        integer node_num_in_1(max_nodes_in)
        integer elem_nodes_in_1(max_elems_in,3)
        integer ni_2, nj_2, num_times_2, num_nodes_in_2, num_elems_in_2
        integer node_i_2(max_nodes_in), node_j_2(max_nodes_in)
        integer node_num_in_2(max_nodes_in)
        integer elem_nodes_in_2(max_elems_in,3)
        integer num_files_1, num_files_2
        integer max_rank, rank
        parameter (max_rank = 3)
        integer dim_sizes(max_rank)
        real*8 weight_wind_node_1(max_nodes,3)
        real*8 weight_wind_node_2(max_nodes,3)
        real*8 x_in_1(max_ni,max_nj), y_in_1(max_ni,max_nj)
        real*8 x_in_2(max_ni,max_nj), y_in_2(max_ni,max_nj)
        real*8 start_day_1
        real*8 relative_weight_1, relative_weight_2
        real*8 max_window_1, max_window_2
        parameter (relative_weight_1 = 1.0)
        parameter (relative_weight_2 = 2.0)
        parameter (max_window_1 = 24.0)
        parameter (max_window_2 = 2.0)
        real*8 frac_day, secs_per_day, utc_start
        parameter (secs_per_day = 86400.0)
        parameter (utc_start = 8.0)
        real*8 temp_arr_1(max_ni,max_nj)
        real*4 temp_arr_2(max_ni,max_nj), temp_arr_3(max_ni,max_nj)
        real*4 temp_arr_4(max_ni,max_nj), temp_arr_5(max_ni,max_nj)
        real*8 temp_arr_6(max_elems_in)
        real*8 temp_arr_8(max_nodes), temp_arr_9(max_nodes)
        real*4 temp_sca
        real*4 wind_times_1(max_times), wind_times_2(max_times)
        character wind_set_1*50, wind_set_2*50, start_day_file*50
        character wind_time_files_1(max_times)*50
        character wind_time_files_2(max_times)*50
        parameter (wind_set_1 = 'hdf/wind_file_1')
        parameter (wind_set_2 = 'hdf/wind_file_2')
        parameter (start_day_file = 'hdf/start_day.txt')
        logical first_call, have_wind_2, have_start_day_file
        data first_call/.true./

! retain the values of some local variables between calls
        save first_call, start_day_1, &
     &    in_elem_to_out_node_1, weight_wind_node_1, &
     &    in_elem_to_out_node_2, weight_wind_node_2, &
     &    num_nodes, have_wind_2, &
     &    wind_times_1, num_times_1, num_files_1, ni_1, nj_1, &
     &    num_nodes_in_1, num_elems_in_1, node_i_1, node_j_1, &
     &    node_num_in_1, elem_nodes_in_1, &
     &    wind_times_2, num_times_2, num_files_2, ni_2, nj_2, &
     &    num_nodes_in_2, num_elems_in_2, node_i_2, node_j_2, &
     &    node_num_in_2, elem_nodes_in_2, &
     &    wind_time_files_1, wind_time_files_2

        open(39,file='fort.39')
        rewind(39)
        write(39,*)
        write(39,*) 'enter get_wind'
        write(39,*) 'first_call = ', first_call

! if this is the first call to this routine then get some things ready
        if (first_call) then

! define the local variables num_nodes
          num_nodes = np

! check to see if start_day_file exists
          call file_exst (start_day_file, have_start_day_file, .false.)

! if the start day file does exist, get start_day from it, otherwise
! use the first start_day in wind_set_1
          if (have_start_day_file) then
            open (unit=77, file=start_day_file, status='old')
            read(77,*) temp_sca
            close (unit=77)
            start_day_1 = temp_sca
          else
            wind_time_files_1(1) = 'hdf/wind_file_1.001.hdf'
            call read_scalar(wind_time_files_1(1), temp_sca, &
     &                       'start_day           ', 0.0)
            start_day_1 = temp_sca
          endif

! check to see if _any_ wind_file_2 exists (use first possible name)
          wind_time_files_2(1) = 'hdf/wind_file_2.001.hdf'
          call file_exst (wind_time_files_2(1), have_wind_2, .false.)
          if (.not. have_wind_2) then
            write(39,*)
            write(39,*) wind_time_files_2(1), ' not exist. . .'
          endif

! get the times of the data available in wind_set_1
          call get_times(wind_times_1, wind_set_1, &
     &                   'u                   ', &
     &                   wind_time_files_1, num_times_1, &
     &                   num_files_1, max_times, max_files)

! get the dimensions of the datasets in wind_set_1 (use first dataset)
          call get_dims(wind_time_files_1(1), 'u                   ', &
     &                  wind_times_1(1), rank, dim_sizes)
          ni_1 = dim_sizes(1)
          nj_1 = dim_sizes(2)

! check the dimensions of wind_set_1, to ensure they don't exceed the
! maximums
          if (ni_1 .gt. max_ni .or. nj_1 .gt. max_nj) then
            write(*,*)
            write(*,*) 'wind_file_1: max dimensions exceeded!'
            write(11,*)
            write(11,*) 'wind_file_1: max dimensions exceeded!'
            stop
          endif

! calculate the total number of nodes and elements for wind_set_1
          num_nodes_in_1 = ni_1 * nj_1
          num_elems_in_1 = (ni_1-1) * (nj_1-1) * 2

! check the elems/nodes of wind_set_1, to ensure they don't exceed the
! maximums
          if (num_elems_in_1 .gt. max_elems_in .or. &
     &        num_nodes_in_1 .gt. max_nodes_in) then
            write(*,*)
            write(*,*) 'wind_file_1: max elems/nodes exceeded!'
            write(11,*)
            write(11,*) 'wind_file_1: max elems/nodes exceeded!'
            stop
          endif

! create list of all nodes for wind_set_1
          call list_nodes (node_i_1, node_j_1, node_num_in_1, &
     &                     num_nodes_in_1, ni_1, nj_1)

! now create the list of all the elements (and the nodes defining them)
! for wind_set_1
          call list_elems (elem_nodes_in_1, node_num_in_1, &
     &                     ni_1, nj_1, num_elems_in_1)

! do the same as above for wind_set_2 (if it exists)
          if (have_wind_2) then

! get the times of the data available in wind_set_2
            call get_times(wind_times_2, wind_set_2, &
     &                     'u                   ', &
     &                     wind_time_files_2, num_times_2, &
     &                     num_files_2, max_times, max_files)

! get the dimensions of the datasets in wind_set_2 (use first dataset)
            call get_dims(wind_time_files_2(1), 'u                   ', &
     &                    wind_times_2(1), rank, dim_sizes)
            ni_2 = dim_sizes(1)
            nj_2 = dim_sizes(2)

! check the dimensions of wind_set_2, to ensure they don't exceed the
! maximums
            if (ni_2 .gt. max_ni .or. nj_2 .gt. max_nj) then
              write(*,*)
              write(*,*) 'wind_file_2: max dimensions exceeded!'
              write(11,*)
              write(11,*) 'wind_file_2: max dimensions exceeded!'
              stop
            endif

! calculate the total number of nodes and elements for wind_set_2
            num_nodes_in_2 = ni_2 * nj_2
            num_elems_in_2 = (ni_2-1) * (nj_2-1) * 2

! check the elems/nodes of wind_set_2, to ensure they don't exceed the
! maximums
            if (num_elems_in_2 .gt. max_elems_in .or. &
     &          num_nodes_in_2 .gt. max_nodes_in) then
              write(*,*)
              write(*,*) 'wind_file_2: max elems/nodes exceeded!'
              write(11,*)
              write(11,*) 'wind_file_2: max elems/nodes exceeded!'
              stop
            endif

! create list of all nodes for wind_set_2
            call list_nodes (node_i_2, node_j_2, node_num_in_2, &
     &                       num_nodes_in_2, ni_2, nj_2)

! now create the list of all the elements (and the nodes defining them)
! for wind_set_2
            call list_elems (elem_nodes_in_2, node_num_in_2, &
     &                       ni_2, nj_2, num_elems_in_2)

          endif ! end of have_wind_2 block

! read in the x and y values for wind_set_1, and copy to full size
! real*8 arrays
          call read_2d_arr(wind_time_files_1(1), temp_arr_2, &
     &                     'x                   ', 0.0, &
     &                     ni_1, nj_1)
          call read_2d_arr(wind_time_files_1(1), temp_arr_3, &
     &                     'y                   ', 0.0, &
     &                     ni_1, nj_1)
          call copy_arr(temp_arr_2, ni_1, nj_1, x_in_1, &
     &                    max_ni, max_nj)
          call copy_arr(temp_arr_3, ni_1, nj_1, y_in_1, &
     &                    max_ni, max_nj)

! calculate the weightings from wind_set_1 to elcirc nodes
! (this is slow)
          write(*,*)
          write(*,*) &
     &      'begin calculating grid weightings for wind_file_1'
          write(16,*)
          write(16,*) &
     &      'begin calculating grid weightings for wind_file_1'
          call get_weight (x_in_1, y_in_1, x, y, &
     &                     elem_nodes_in_1, node_i_1, node_j_1, &
     &                     max_ni, max_nj, &
     &                     num_elems_in_1, num_nodes_in_1, &
     &                     num_nodes, &
     &                     max_nodes, &
     &                     in_elem_to_out_node_1, &
     &                     temp_arr_6, weight_wind_node_1)
          write(*,*) &
     &      'done calculating grid weightings for wind_file_1'
          write(16,*) &
     &      'done calculating grid weightings for wind_file_1'

! do the same but for wind_set_2 (if it exists)
          if (have_wind_2) then

! read in the x and y values for wind_set_2, and copy to full size
! real*8 arrays
            call read_2d_arr(wind_time_files_2(1), temp_arr_2, &
     &                       'x                   ', 0.0, &
     &                       ni_2, nj_2)
            call read_2d_arr(wind_time_files_2(1), temp_arr_3, &
     &                       'y                   ', 0.0, &
     &                       ni_2, nj_2)
            call copy_arr(temp_arr_2, ni_2, nj_2, x_in_2, &
     &                    max_ni, max_nj)
            call copy_arr(temp_arr_3, ni_2, nj_2, y_in_2, &
     &                    max_ni, max_nj)

! calculate the weightings from wind_set_2 to elcirc nodes
! (this is slow)
            write(*,*)
            write(*,*) &
     &        'begin calculating grid weightings for wind_file_2'
            write(16,*)
            write(16,*) &
     &        'begin calculating grid weightings for wind_file_2'
            call get_weight (x_in_2, y_in_2, x, y, &
     &                       elem_nodes_in_2, node_i_2, node_j_2, &
     &                       max_ni, max_nj, &
     &                       num_elems_in_2, num_nodes_in_2, &
     &                       num_nodes, &
     &                       max_nodes, &
     &                       in_elem_to_out_node_2, &
     &                       temp_arr_6, weight_wind_node_2)
            write(*,*) &
     &        'done calculating grid weightings for wind_file_2'
            write(16,*) &
     &        'done calculating grid weightings for wind_file_2'

          endif

! output starting date and time
          write(*,*)
          write(*,*) 'wind file starting Julian date: ', start_day_1
          write(*,*) 'wind file assumed UTC starting time: ', utc_start
          write(16,*)
          write(16,*) 'wind file starting Julian date: ', start_day_1
          write(16,*) 'wind file assumed UTC starting time: ', utc_start

        endif ! (end of first_call block)

! define frac_day - the fractional Julian date
! include offset for starting time in UTC (in hours)
        frac_day = start_day_1 + time/secs_per_day + utc_start/24.0

! output info to debug file
        write(39,*) 'num_nodes = ', num_nodes
        write(39,*) 'num_files_1 = ', num_files_1
        write(39,*) 'num_times_1 = ', num_times_1