&                  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