sflux_subs5.f90

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

        if (have_wind_2) then
          write(39,*) 'num_files_2 = ', num_files_2
          write(39,*) 'num_times_2 = ', num_times_2
        endif
        write(39,*)
        write(39,*) 'start_day = ', start_day_1
        write(39,*) 'time      = ', time
        write(39,*) 'frac_day  = ', frac_day

! now retrieve data from the wind files

        call combine_set (real(frac_day), 'u                   ', &
     &                    wind_set_1, wind_set_2, &
     &                    temp_arr_2, temp_arr_3, &
     &                    temp_arr_4, temp_arr_5, temp_arr_1, &
     &                    wind_times_1, num_times_1, ni_1, nj_1, &
     &                    wind_times_2, num_times_2, ni_2, nj_2, &
     &                    num_files_1, num_files_2, max_files, &
     &                    temp_arr_8, temp_arr_9, u_air_node, &
     &                    in_elem_to_out_node_1, weight_wind_node_1, &
     &                    in_elem_to_out_node_2, weight_wind_node_2, &
     &                    elem_nodes_in_1, node_i_1, node_j_1, &
     &                    elem_nodes_in_2, node_i_2, node_j_2, &
     &                    max_ni, max_nj, max_nodes_in, &
     &                    num_nodes, max_nodes, &
     &                    num_nodes_in_1, num_elems_in_1, &
     &                    num_nodes_in_2, num_elems_in_2, &
     &                    max_elems_in, have_wind_2, &
     &                    relative_weight_1, relative_weight_2, &
     &                    max_window_1, max_window_2, &
     &                    wind_time_files_1, wind_time_files_2)

        call combine_set (real(frac_day), 'v                   ', &
     &                    wind_set_1, wind_set_2, &
     &                    temp_arr_2, temp_arr_3, &
     &                    temp_arr_4, temp_arr_5, temp_arr_1, &
     &                    wind_times_1, num_times_1, ni_1, nj_1, &
     &                    wind_times_2, num_times_2, ni_2, nj_2, &
     &                    num_files_1, num_files_2, max_files, &
     &                    temp_arr_8, temp_arr_9, v_air_node, &
     &                    in_elem_to_out_node_1, weight_wind_node_1, &
     &                    in_elem_to_out_node_2, weight_wind_node_2, &
     &                    elem_nodes_in_1, node_i_1, node_j_1, &
     &                    elem_nodes_in_2, node_i_2, node_j_2, &
     &                    max_ni, max_nj, max_nodes_in, &
     &                    num_nodes, max_nodes, &
     &                    num_nodes_in_1, num_elems_in_1, &
     &                    num_nodes_in_2, num_elems_in_2, &
     &                    max_elems_in, have_wind_2, &
     &                    relative_weight_1, relative_weight_2, &
     &                    max_window_1, max_window_2, &
     &                    wind_time_files_1, wind_time_files_2)

        call combine_set (real(frac_day), 'p_msl               ', &
     &                    wind_set_1, wind_set_2, &
     &                    temp_arr_2, temp_arr_3, &
     &                    temp_arr_4, temp_arr_5, temp_arr_1, &
     &                    wind_times_1, num_times_1, ni_1, nj_1, &
     &                    wind_times_2, num_times_2, ni_2, nj_2, &
     &                    num_files_1, num_files_2, max_files, &
     &                    temp_arr_8, temp_arr_9, p_air_node, &
     &                    in_elem_to_out_node_1, weight_wind_node_1, &
     &                    in_elem_to_out_node_2, weight_wind_node_2, &
     &                    elem_nodes_in_1, node_i_1, node_j_1, &
     &                    elem_nodes_in_2, node_i_2, node_j_2, &
     &                    max_ni, max_nj, max_nodes_in, &
     &                    num_nodes, max_nodes, &
     &                    num_nodes_in_1, num_elems_in_1, &
     &                    num_nodes_in_2, num_elems_in_2, &
     &                    max_elems_in, have_wind_2, &
     &                    relative_weight_1, relative_weight_2, &
     &                    max_window_1, max_window_2, &
     &                    wind_time_files_1, wind_time_files_2)

        call combine_set (real(frac_day), 'T                   ', &
     &                    wind_set_1, wind_set_2, &
     &                    temp_arr_2, temp_arr_3, &
     &                    temp_arr_4, temp_arr_5, temp_arr_1, &
     &                    wind_times_1, num_times_1, ni_1, nj_1, &
     &                    wind_times_2, num_times_2, ni_2, nj_2, &
     &                    num_files_1, num_files_2, max_files, &
     &                    temp_arr_8, temp_arr_9, t_air_node, &
     &                    in_elem_to_out_node_1, weight_wind_node_1, &
     &                    in_elem_to_out_node_2, weight_wind_node_2, &
     &                    elem_nodes_in_1, node_i_1, node_j_1, &
     &                    elem_nodes_in_2, node_i_2, node_j_2, &
     &                    max_ni, max_nj, max_nodes_in, &
     &                    num_nodes, max_nodes, &
     &                    num_nodes_in_1, num_elems_in_1, &
     &                    num_nodes_in_2, num_elems_in_2, &
     &                    max_elems_in, have_wind_2, &
     &                    relative_weight_1, relative_weight_2, &
     &                    max_window_1, max_window_2, &
     &                    wind_time_files_1, wind_time_files_2)

        call combine_set (real(frac_day), 'qv                  ', &
     &                    wind_set_1, wind_set_2, &
     &                    temp_arr_2, temp_arr_3, &
     &                    temp_arr_4, temp_arr_5, temp_arr_1, &
     &                    wind_times_1, num_times_1, ni_1, nj_1, &
     &                    wind_times_2, num_times_2, ni_2, nj_2, &
     &                    num_files_1, num_files_2, max_files, &
     &                    temp_arr_8, temp_arr_9, q_air_node, &
     &                    in_elem_to_out_node_1, weight_wind_node_1, &
     &                    in_elem_to_out_node_2, weight_wind_node_2, &
     &                    elem_nodes_in_1, node_i_1, node_j_1, &
     &                    elem_nodes_in_2, node_i_2, node_j_2, &
     &                    max_ni, max_nj, max_nodes_in, &
     &                    num_nodes, max_nodes, &
     &                    num_nodes_in_1, num_elems_in_1, &
     &                    num_nodes_in_2, num_elems_in_2, &
     &                    max_elems_in, have_wind_2, &
     &                    relative_weight_1, relative_weight_2, &
     &                    max_window_1, max_window_2, &
     &                    wind_time_files_1, wind_time_files_2)

! set first_call to false, so subsequent calls will know that they're
! not the first call
        first_call = .false.

        close(39)

      return
      end
!-----------------------------------------------------------------------
!
! Note that net downward heat flux (except for solar) is given by:
!
! net_sfc_flux_d = - (sen_flux + lat_flux + (longwave_u - longwave_d))
!
      subroutine surf_fluxes (time, u_air, v_air, p_air, &
     &                   t_air, q_air, shortwave_d, &
     &                   sen_flux, lat_flux, longwave_u, longwave_d, &
     &                   tau_xz, tau_yz)

!       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(max_nodes), v_air(max_nodes)
        real*8 t_air(max_nodes), q_air(max_nodes)
        real*8 p_air(max_nodes), time
        real*8 sen_flux(max_nodes), lat_flux(max_nodes)
        real*8 longwave_d(max_nodes), longwave_u(max_nodes)
        real*8 shortwave_d(max_nodes)
        real*8 tau_xz(max_nodes), tau_yz(max_nodes)

! local variables
        integer num_nodes, i_node
        real*4 temp_sca, window
        real*8 weight_rad_node(max_nodes,3)
        real*8 start_day, one_third, frac_day, secs_per_day, utc_start
        real*8 max_window
        parameter (max_window = 6.0)
        parameter (secs_per_day = 86400.0)
        parameter (utc_start = 8.0)
        integer in_elem_to_out_node(max_nodes)
        character flux_set*50, start_day_file*50
        parameter (flux_set = 'hdf/flux_file_1')
        parameter (start_day_file = 'hdf/start_day.txt')
        logical first_call, time_exst, have_start_day_file
        data first_call/.true./
        integer max_ni, max_nj, ni, nj, max_times, num_times
        integer num_files, max_files
        integer num_nodes_in, num_elems_in
        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)
        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)
        real*4 flux_times(max_times)
        real*8 shortwave_in(max_ni,max_nj), longwave_in(max_ni,max_nj)
        real*8 x_in(max_ni,max_nj), y_in(max_ni,max_nj)
        real*8 stefan, t_freeze, emissivity
        character flux_time_files(max_times)*50
        integer node_i(max_nodes_in), node_j(max_nodes_in)
        integer max_rank, rank
        integer node_num_in(max_nodes_in), elem_nodes_in(max_elems_in,3)
        parameter (max_rank = 3)
        integer dim_sizes(max_rank)
        parameter (stefan = 5.67e-8)
        parameter (t_freeze = 273.15)
        parameter (emissivity = 1.0)
        integer printit

! retain the values of some local variables between calls
        save first_call, start_day, in_elem_to_out_node, &
     &       weight_rad_node, num_nodes, &
     &       flux_times, num_times, num_files, ni, nj, &
     &       num_nodes_in, num_elems_in, node_i, node_j, node_num_in, &
     &       elem_nodes_in, flux_time_files

        printit = 100

        open(38,file='fort.38')
        rewind(38)
        write(38,*)
        write(38,*) 'enter surf_fluxes'
        write(38,*) 'first_call = ', first_call

! define some constants, etc
        one_third = 1.0 / 3.0

! 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 flux_set
          if (have_start_day_file) then
            open (unit=77, file=start_day_file, status='old')
            read(77,*) temp_sca
            close (unit=77)
            start_day = temp_sca
          else
            flux_time_files(1) = 'hdf/flux_file_1.001.hdf'
            call read_scalar(flux_time_files(1), temp_sca, &
     &                       'start_day           ', 0.0)
            start_day = temp_sca
          endif

! get the times of the fluxes available in the flux_set
          call get_times(flux_times, flux_set, &
     &                   'shortwave_flux      ', &
     &                   flux_time_files, num_times, &
     &                   num_files, max_times, max_files)

! get the dimensions of the flux datasets (use first dataset)
          call get_dims(flux_time_files(1), 'shortwave_flux      ', &
     &                  flux_times(1), rank, dim_sizes)
          ni = dim_sizes(1)
          nj = dim_sizes(2)

! check the dimensions of the flux datasets, to ensure they don't exceed
! the maximums
          if (ni .gt. max_ni .or. nj .gt. max_nj) then
            write(*,*)
            write(*,*) 'Radiative flux data: max dimensions exceeded!'
            write(11,*)
            write(11,*) 'Radiative flux data: max dimensions exceeded!'
            stop
          endif

! calculate the total number of nodes and elements for the radiative
! fluxes input grid
          num_nodes_in = ni * nj
          num_elems_in = (ni-1) * (nj-1) * 2

! check the elems/nodes of the input dataset, to ensure they don't
! exceed the maximums
          if (num_elems_in .gt. max_elems_in .or. &
     &        num_nodes_in .gt. max_nodes_in) then
            write(*,*)
            write(*,*) 'Radiative flux data: max elems/nodes exceeded!'
            write(11,*)
            write(11,*) 'Radiative flux data: max elems/nodes exceeded!'
            stop
          endif

! create list of all nodes for input grid
          call list_nodes (node_i, node_j, node_num_in, &
     &                     num_nodes_in, ni, nj)

! now create the list of all the elements (and the nodes defining them)
          call list_elems (elem_nodes_in, node_num_in, &
     &                     ni, nj, num_elems_in)

! read in the x and y values for the input grid, and copy to full size
! real*8 arrays
          call read_2d_arr(flux_time_files(1), temp_arr_4, &
     &                     'x                   ', 0.0, &
     &                     ni, nj)
          call read_2d_arr(flux_time_files(1), temp_arr_5, &
     &                     'y                   ', 0.0, &
     &                     ni, nj)
          call copy_arr(temp_arr_4, ni, nj, x_in, max_ni, max_nj)
          call copy_arr(temp_arr_5, ni, nj, y_in, max_ni, max_nj)

! calculate the weightings from rad grid to elcirc nodes (this is slow)
          write(*,*)
          write(*,*) 'calculating grid weightings for rad fluxes'
          write(16,*)
          write(16,*) 'calculating grid weightings for rad fluxes'
          call get_weight (x_in, y_in, x, y, &
     &                     elem_nodes_in, node_i, node_j, &
     &                     max_ni, max_nj, &
     &                     num_elems_in, num_nodes_in, &
     &                     num_nodes, &
     &                     max_nodes, &
     &                     in_elem_to_out_node, &
     &                     temp_arr_6, weight_rad_node)

! output starting date and time
          write(*,*)
          write(*,*) 'rad fluxes file starting Julian date: ', &
     &               start_day
          write(*,*) 'rad fluxes file assumed UTC starting time: ', &
     &               utc_start
          write(16,*)
          write(16,*) 'rad fluxes file starting Julian date: ', &
     &                start_day
          write(16,*) 'rad fluxes 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 + time/secs_per_day + utc_start/24.0

! output info to debug file
        write(38,*) 'num_nodes = ', num_nodes
        write(38,*) 'num_files = ', num_files
        write(38,*) 'num_times = ', num_times
        write(38,*)
        write(38,*) 'start_day = ', start_day
        write(38,*) 'time      = ', time
        write(38,*) 'frac_day  = ', frac_day

! convert air temperatures to Celcius
        do i_node = 1, num_nodes
          t_air(i_node) = t_air(i_node) - t_freeze

          if (mod(i_node-1,printit) .eq. 0) then
            write(38,*)
            write(38,*) 'i_node, sfc u, v, T = ', i_node, &