sflux_subs5.f90

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

     &                  uu2(i_node, kfp(i_node)), &
     &                  vv2(i_node, kfp(i_node)), &
     &                  tnd(i_node, kfp(i_node))
            write(38,*) 'u, v, p, T, q (air) = ', u_air(i_node), &
     &                  v_air(i_node), p_air(i_node), t_air(i_node), &
     &                  q_air(i_node)
          endif

        enddo

! calculate the turbulent fluxes at the nodes
        write(38,*) 'above turb_fluxes'
        call turb_fluxes (num_nodes, &
     &                    u_air, v_air, p_air, t_air, q_air, &
     &                    sen_flux, lat_flux, tau_xz, tau_yz)
        write(38,*) 'below turb_fluxes'

! now we must retrieve the radiative fluxes, and copy from the filled
! real*4 arrays of reduced size to full-size real*8 arrays

! first shortwave
        call get_set (real(frac_day), 'shortwave_flux      ', &
     &                flux_set, time_exst, &
     &                temp_arr_4, temp_arr_2, temp_arr_3, &
     &                flux_times, num_times, ni, nj, &
     &                num_files, max_files, window, &
     &                flux_time_files)
        call copy_arr(temp_arr_4, ni, nj, shortwave_in, max_ni, max_nj)

! next longwave
        call get_set (real(frac_day), 'longwave_flux       ', &
     &                flux_set, time_exst, &
     &                temp_arr_4, temp_arr_2, temp_arr_3, &
     &                flux_times, num_times, ni, nj, &
     &                num_files, max_files, window, &
     &                flux_time_files)
        call copy_arr(temp_arr_4, ni, nj, longwave_in, max_ni, max_nj)
!       write(*,*) 'time_exst, window = ', time_exst, window

! if the desired time is outside of the times supplied by the dataset,
! or if the data interpolation window is greater than the maximum, then
! quit with an error message
        if ( (.not. time_exst) .or. (window .gt. max_window) ) then
          write(*,*)
          write(*,*) 'Radiative flux data: time not present in file!'
          write(11,*)
          write(11,*) 'Radiative flux data: time not present in file!'
          stop
        endif

        write(38,*) 'rad fluxes obtained'

! now the input radiative fluxes must be put on a triangular grid

! use the weightings to interpolate from the input grid to the node
! points
        write(38,*) 'interpolating longwave'
        call interp_data (longwave_d, longwave_in, weight_rad_node, &
     &                    in_elem_to_out_node, elem_nodes_in, &
     &                    node_i, node_j, &
     &                    num_nodes_in, num_nodes, num_elems_in, &
     &                    max_ni, max_nj, max_nodes)

        write(38,*) 'interpolating shortwave'
        call interp_data (shortwave_d, shortwave_in, weight_rad_node, &
     &                    in_elem_to_out_node, elem_nodes_in, &
     &                    node_i, node_j, &
     &                    num_nodes_in, num_nodes, num_elems_in, &
     &                    max_ni, max_nj, max_nodes)

! now calculate upwards longwave flux at the surface, using black-body
! equation
        write(38,*) 'calculating longwave_u'
        do i_node = 1, num_nodes
          longwave_u(i_node) = emissivity * stefan * &
     &                      ( t_freeze + tnd(i_node, kfp(i_node)) ) ** 4
        enddo

! get the albedo (store in temp_arr_8)
        write(38,*) 'calculating albedo'
        call get_albedo (temp_arr_8, num_nodes, max_nodes)

! reduce the downwards shortwave flux at the surface by the albedo
! (and check for possible negative values from interpolation)
        write(38,*) 'reducing shortwave'
        do i_node = 1, num_nodes
          shortwave_d(i_node) = max( (1.0 - temp_arr_8(i_node)) * &
     &                               shortwave_d(i_node), &
     &                               0.0 )
        enddo

        do i_node = 1, num_nodes
          if (mod(i_node-1,printit) .eq. 0) then
            if (kfp(i_node) .ne. 0) then
              write(38,*)
              write(38,*) 'i_node = ', i_node
              write(38,*) 'net_sfc_flux_d = ', &
     &                     - sen_flux(i_node) - lat_flux(i_node) &
     &                     - ( longwave_u(i_node) - longwave_d(i_node) )
              write(38,*) 'shortwave_d = ', shortwave_d(i_node)
              write(38,*) 'longwave_d, longwave_u = ', &
     &                     longwave_d(i_node), longwave_u(i_node)
              write(38,*) 'sen_flux, lat_flux = ', &
     &                     sen_flux(i_node), lat_flux(i_node)
            else
              write(38,*)
              write(38,*) 'i_node = ', i_node
              write(38,*) 'dry!'
            endif
          endif
        enddo

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

        close(38)

      return
      end
!-----------------------------------------------------------------------
!
! Calculate bulk aerodynamic surface fluxes over water using method of
! Zeng et al., J Clim, v11, p 2628, Oct 1998.
!
! Note: this subroutine uses actual temperatures instead of potential
!       temperatures.  Since the temperature height is only 2m, the
!       difference should be negligible. . .
!
      subroutine turb_fluxes (num_nodes, &
     &                        u_air, v_air, p_air, t_air, q_air, &
     &                        sen_flux, lat_flux, 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
        integer num_nodes
        real*8 u_air(max_nodes), v_air(max_nodes), p_air(max_nodes)
        real*8 t_air(max_nodes), q_air(max_nodes)
        real*8 sen_flux(max_nodes), lat_flux(max_nodes)
        real*8 tau_xz(max_nodes), tau_yz(max_nodes)

! local variables
        integer i_node, iter, max_iter
        parameter (max_iter = 10)
        real*8 u_star, theta_star, q_star, z_0, monin
        real*8 zeta_u, zeta_t
        real*8 one_third, two_thirds, z_t, z_u
        real*8 a1, a2, b1, b2, nu, g, beta, w_star, mix_ratio
        real*8 t_v, z_i, speed, karman, zeta_m, psi_m, zeta_h, psi_h
        real*8 re, z_0_t, e_sfc, t_freeze, q_sfc, epsilon_r
        real*8 rho_air, c_p_air, latent, r_air, rb
        real*8 theta_air, theta_v_air, delta_theta, delta_q
        real*8 delta_theta_v, theta_v_star, speed_res, tau
        real*4 esat_flat_r
        parameter (z_t = 2.0)
        parameter (z_u = 10.0)
        parameter (a1 = 0.013)
        parameter (a2 = 0.11)
        parameter (b1 = 2.67)
        parameter (b2 = -2.57)
        parameter (nu = 1.46e-5)
        parameter (beta = 1.0)
        parameter (g = 9.81)
        parameter (z_i = 1000.0)
        parameter (karman = 0.4)
        parameter (zeta_m = -1.574)
        parameter (zeta_h = -0.465)
        parameter (t_freeze = 273.15)
        parameter (epsilon_r = 0.6220)
        parameter (c_p_air = 1004.0)
        parameter (latent = 2.501e6)
        parameter (r_air = 287.0)
        integer printit
        logical converged
        printit = 100

        write(38,*) 'enter turb_fluxes'

! precalculate constants
        one_third = 1.0 / 3.0
        two_thirds = 2.0 / 3.0

! now loop over all points
        do i_node = 1, num_nodes
          if (mod(i_node-1,printit) .eq. 0) then
            write(38,*)
            write(38,*) 'i_node = ', i_node
          endif

! if this point isn't dry, then calculate fluxes (if dry, then skip)
        if (kfp(i_node) .ne. 0) then

! calculate q_sfc from e_sfc
! (e_sfc reduced for salinity using eqn from Smithsonian Met Tables)
          e_sfc = (1.0 - 0.000537 * snd(i_node, kfp(i_node))) &
     &          * esat_flat_r(real(tnd(i_node, kfp(i_node)) + t_freeze))
          q_sfc = epsilon_r * e_sfc &
     &          / ( p_air(i_node) - e_sfc * (1.0 - epsilon_r) )

! calculate the water vapor mixing ratio of the air
          mix_ratio = q_air(i_node) / (1.0 - q_air(i_node))

! calculate theta_air, theta_v_air, delta_theta, delta_q,
! and delta_theta_v
          theta_air = (t_air(i_node) + t_freeze) + 0.0098*z_t
          theta_v_air = theta_air * (1.0 + 0.608 * mix_ratio)
          delta_theta = theta_air - &
     &                  (tnd(i_node, kfp(i_node)) + t_freeze)
          delta_q = q_air(i_node) - q_sfc
          delta_theta_v = delta_theta * (1.0 + 0.608 * mix_ratio) &
     &                  + 0.608 * theta_air * delta_q

! calculate the air virtual temperature and density
          t_v = (t_air(i_node) + t_freeze) * (1.0 + 0.608 * mix_ratio)
          rho_air = p_air(i_node) / (r_air * t_v)
          if (mod(i_node-1,printit) .eq. 0) then
            write(38,*) 'e_sfc, q_sfc, mix_ratio = ', &
     &                   e_sfc, q_sfc, mix_ratio
            write(38,*) 'theta_air, theta_v_air, delta_theta = ', &
     &                   theta_air, theta_v_air, delta_theta
            write(38,*) 'delta_q, delta_theta_v, t_v = ', &
     &                   delta_q, delta_theta_v, t_v
            write(38,*) 'rho_air = ', rho_air
          endif

! begin with initial values of u_star, w_star, and speed
          u_star = 0.06
          w_star = 0.5
          if (delta_theta_v .ge. 0) then                    ! stable
            speed = &
     &        max( sqrt( &
     &               (u_air(i_node) - uu2(i_node, kfp(i_node)))**2 + &
     &               (v_air(i_node) - vv2(i_node, kfp(i_node)))**2 ), &
     &             0.1)
          else                                              ! unstable
            speed = &
     &        sqrt( (u_air(i_node) - uu2(i_node, kfp(i_node)))**2 + &
     &              (v_air(i_node) - vv2(i_node, kfp(i_node)))**2 + &
     &              (beta * w_star)**2 )
          endif

! now loop to obtain good initial values for u_star and z_0
          do iter = 1, 5
            z_0 = a1 * u_star * u_star / g + a2 * nu / u_star
            u_star = karman * speed / log(z_u/z_0)
          enddo

! calculate rb (some stability parameter from Xubin's code?)
          rb = g * z_u * delta_theta_v / (theta_v_air * speed * speed)

! calculate initial values for zeta_u, monin, zeta_t
          if (rb .ge. 0) then                      ! neutral or stable
            zeta_u = rb * log(z_u/z_0) &
     &             / (1.0 - 0.5*min(rb,0.19))
          else
            zeta_u = rb * log(z_u/z_0)
          endif
          monin = z_u / zeta_u
          zeta_t = z_t / monin
          if (mod(i_node-1,printit) .eq. 0) then
            write(38,*) 'speed, z_0, u_star = ', &
     &                   speed, z_0, u_star
            write(38,*) 'zeta_u, zeta_t, monin = ', &
     &                   zeta_u, zeta_t, monin
          endif

! iterate a maximum of max_iter times
          iter = 0
          converged = .false.
100       continue
            iter = iter + 1

! Calculate the roughness lengths
            z_0 = a1 * u_star * u_star / g + a2 * nu / u_star
            re = u_star * z_0 / nu
            z_0_t = z_0 / exp(b1 * (re**0.25) + b2)
            if (mod(i_node-1,printit) .eq. 0) then
              write(38,*) 're, z_0, z_0_t = ', &
     &                     re, z_0, z_0_t
            endif

! calculate the zetas
            zeta_u = z_u / monin
            zeta_t = z_t / monin

! apply asymptotic limit to stable conditions
            if (zeta_t .gt. 2.5) then
              converged = .true.
              zeta_t = 2.5
              monin = z_t / zeta_t
              zeta_u = z_u / monin
              if (mod(i_node-1,printit) .eq. 0) then
                write(38,*) 'limiting zeta_u, zeta_t, monin!'
              endif
            endif

! caulculate u_star, depending on zeta
            if (zeta_u .lt. zeta_m) then                ! very unstable

              u_star = speed * karman &
     &               / ( log(zeta_m*monin/z_0) &
     &                   - psi_m(zeta_m) &
! extra term?
     &                   + psi_m(z_0/monin) &
     &                   + 1.14 * ((-zeta_u)**(one_third) - &
     &                           (-zeta_m)**(one_third)) )

            else if (zeta_u .lt. 0.0) then              ! unstable

              u_star = speed * karman &
     &               / ( log(z_u/z_0) &
     &                   - psi_m(zeta_u) &
! extra term?
     &                   + psi_m(z_0/monin) &
     &                 )

            else if (zeta_u .le. 1.0) then              ! neutral/stable

              u_star = speed * karman &
     &               / ( log(z_u/z_0) + 5.0*zeta_u &
! extra term?
     &                   - 5.0*z_0/monin &
     &                 )

            else                                        ! very stable

              u_star = speed * karman &
     &               / ( log(monin/z_0) + 5.0 &
     &                   + 5.0*log(zeta_u) &
! extra term?
     &                   - 5.0*z_0/monin &
     &                   + zeta_u - 1.0 )

            endif