! caulculate theta_star and q_star, depending on zeta            if (zeta_t .lt. zeta_h) then                ! very unstable              theta_star = karman * delta_theta &     &                   / ( log(zeta_h*monin/z_0_t) &     &                       - psi_h(zeta_h) &! extra term?     &                       + psi_h(z_0_t/monin) &     &                       + 0.8 * ((-zeta_h)**(-one_third) - &     &                                (-zeta_t)**(-one_third)) )              q_star = karman * delta_q &     &                   / ( log(zeta_h*monin/z_0_t) &     &                       - psi_h(zeta_h) &! extra term?     &                       + psi_h(z_0_t/monin) &     &                       + 0.8 * ((-zeta_h)**(-one_third) - &     &                                (-zeta_t)**(-one_third)) )            else if (zeta_t .lt. 0.0) then              ! unstable              theta_star = karman * delta_theta &     &                   / ( log(z_t/z_0_t) &     &                       - psi_h(zeta_t) &! extra term?     &                       + psi_h(z_0_t/monin) &     &                     )              q_star = karman * delta_q &     &                   / ( log(z_t/z_0_t) &     &                       - psi_h(zeta_t) &! extra term?     &                       + psi_h(z_0_t/monin) &     &                     )            else if (zeta_t .lt. 1.0) then              ! neutral/stable              theta_star = karman * delta_theta &     &                   / ( log(z_t/z_0_t) &     &                       + 5.0*zeta_t &! extra term?     &                       - 5.0*z_0_t/monin &     &                     )              q_star = karman * delta_q &     &                   / ( log(z_t/z_0_t) &     &                       + 5.0*zeta_t &! extra term?     &                       - 5.0*z_0_t/monin &     &                     )            else                                        ! very stable              theta_star = karman * delta_theta &     &                   / ( log(monin/z_0_t) + 5.0 &     &                       + 5.0*log(zeta_t) &! extra term?     &                       - 5.0*z_0_t/monin &     &                       + zeta_t - 1.0 )              q_star = karman * delta_q &     &                   / ( log(monin/z_0_t) + 5.0 &     &                       + 5.0*log(zeta_t) &! extra term?     &                       - 5.0*z_0_t/monin &     &                       + zeta_t - 1.0 )            endif! calculate theta_v_star and monin            theta_v_star = theta_star * (1.0 + 0.608 * mix_ratio) &     &                   + 0.608 * theta_air * q_star            monin = theta_v_air * u_star * u_star &     &            / (karman * g * theta_v_star)! depending on surface layer stability, calculate the effective! near-surface wind speed! (ie relative to the flowing water surface)            if (delta_theta_v .ge. 0.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! calculate the convective velocity scale              w_star = (-g * theta_v_star * u_star * z_i / theta_v_air) &     &                 ** one_third              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            if (mod(i_node-1,printit) .eq. 0) then              write(38,*) 'iter, u_star, q_star, theta_star = ', &     &                     iter, u_star, q_star, theta_star              write(38,*) 'iter, theta_v_star, monin, speed = ', &     &                     iter, theta_v_star, monin, speed              write(38,*) 'iter, zeta_u, zeta_t = ', &     &                     iter, zeta_u, zeta_t            endif! bottom of main iteration loop          if (.not. converged .and. iter .lt. max_iter) goto 100! calculate fluxes          sen_flux(i_node) = - rho_air * c_p_air * u_star * theta_star          lat_flux(i_node) = - rho_air * latent * u_star * q_star! calculate wind stresses          speed_res = &     &          sqrt( (u_air(i_node) - uu2(i_node, kfp(i_node)))**2 + &     &                (v_air(i_node) - vv2(i_node, kfp(i_node)))**2 )          if (speed_res .gt. 0.0) then            tau = rho_air * u_star * u_star * speed_res / speed            tau_xz(i_node) = - tau &     &                     * (u_air(i_node) - uu2(i_node, kfp(i_node))) &     &                     / speed_res            tau_yz(i_node) = - tau &     &                     * (v_air(i_node) - vv2(i_node, kfp(i_node))) &     &                     / speed_res          else            tau_xz(i_node) = 0.0            tau_yz(i_node) = 0.0          endif          if (mod(i_node-1,printit) .eq. 0) then            write(38,*) 'sen_flux, lat_flux = ', &     &                   sen_flux(i_node), lat_flux(i_node)            write(38,*) 'tau_xz, tau_yz = ', &     &                   tau_xz(i_node), tau_yz(i_node)          endif!         write(38,*) tnd(i_node, kfp(i_node)), sen_flux(i_node),!    +                lat_flux(i_node)! end of wet/dry block        endif! end of loop over points        enddo        write(38,*) 'exit turb_fluxes'      return      end!-----------------------------------------------------------------------!! Calculate saturation vapor pressure using the eighth order relative! error norm method of Flatau et al., J Applied Meteo, v31, p 1507,! Dec 1992.!      real*4 function esat_flat_r(t)        implicit none        real*4 t        real*4 c0, c1, c2, c3, c4, c5, c6, c7, c8, t_eff        parameter ( &     &        c0= 6.11583699e+02,  c1= 0.444606896e+02, &     &        c2= 0.143177157e+01, c3= 0.264224321e-01, &     &        c4= 0.299291081e-03, c5= 0.203154182e-05, &     &        c6= 0.702620698e-08, c7= 0.379534310e-11, &     &        c8=-0.321582393e-13)! t     : temperature in K! t_eff : effective temperature in C        t_eff = max(-85.,t-273.16)        esat_flat_r = c0+t_eff*(c1+t_eff*(c2+t_eff*(c3+t_eff*(c4+t_eff* &     &                         (c5+t_eff*(c6+t_eff*(c7+t_eff*c8)))))))      return      end!-----------------------------------------------------------------------      subroutine copy_arr(array_1, ni_1, nj_1, array_2, ni_2, nj_2)        implicit none        integer ni_1, nj_1, ni_2, nj_2, i, j        real*4 array_1(ni_1,nj_1)        real*8 array_2(ni_2,nj_2)        do j = 1, nj_1          do i = 1, ni_1            array_2(i,j) = array_1(i,j)          enddo        enddo      return      end!-----------------------------------------------------------------------      subroutine file_exst (file_name, exst, fail)        implicit none        character file_name*50        logical exst, fail        inquire(file=file_name, exist=exst)        if (exst .and. fail) then          write(*,*)          write(*,*) file_name, ' already exists!'          write(11,*)          write(11,*) file_name, ' already exists!'         stop        endif      return      end!-----------------------------------------------------------------------      subroutine get_bracket (time, input_times, i_time, bracket, &     &                        num_times)        implicit none        integer num_times, i_time        real*4 time, input_times(num_times)        logical bracket        i_time = 0        bracket = .false.10      continue          i_time = i_time + 1          bracket = ( input_times(i_time) .le. time .and. &     &                time .le. input_times(i_time+1) )        if (.not. bracket .and. i_time .lt. num_times-1) goto 10      return      end!-----------------------------------------------------------------------      subroutine interp_set(in_set, data, data_label, time, &     &                      data_1, data_2, input_times, &     &                      i_time, num_times, ni, nj, &     &                      time_files, num_files, max_files)        implicit none        integer ni, nj, i, j, num_times, i_time, num_files, max_files        real*4 data_1(ni,nj), data_2(ni,nj)        real*4 data(ni,nj), time, input_times(num_times), ratio        character in_set*50, data_label*20        character time_files(num_times)*50! read in the data at the first of the bracketing times        call read_2d_arr (time_files(i_time), data_1, data_label, &     &                    input_times(i_time), ni, nj)! read in the data at the second of the bracketing times        call read_2d_arr (time_files(i_time+1), data_2, data_label, &     &                    input_times(i_time+1), ni, nj)! now interpolate these to the desired time        ratio = (time                  - input_times(i_time)) &     &        / (input_times(i_time+1) - input_times(i_time))        do j = 1, nj          do i = 1, ni            data(i,j) = data_1(i,j) &     &                + (data_2(i,j) - data_1(i,j)) * ratio          enddo        enddo      return      end!-----------------------------------------------------------------------      subroutine get_dims(in_file, data_label, t, rank, dim_sizes)        implicit none        real*4 t        integer ret, read_only, sfstart, sfend, sd_id, sds_id        integer sds_index, sfn2index, sfendacc, sfselect        integer max_rank, rank, data_type, n_attrs, sfginfo        parameter (max_rank = 3)        parameter (read_only = 1)        integer dim_sizes(max_rank)        character data_label*20, in_file*50        character dat_time_label*33, sds_name*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! get info on dataset        sds_id = sfselect(sd_id, sds_index)        ret = sfginfo(sds_id, sds_name, rank, dim_sizes, &     &                data_type, n_attrs)        call checkret(ret)! find the id for this index        sds_id = sfselect(sd_id, sds_index)! 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 list_nodes (node_i, node_j, node_num, &     &                       n_nodes_in, ni, nj)        implicit none        integer ni, nj, n_nodes_in, i, j, i_node        integer node_i(n_nodes_in), node_j(n_nodes_in)        integer node_num(ni,nj)        i_node = 0        do j = 1, nj          do i = 1, ni            i_node = i_node + 1            node_i(i_node) = i            node_j(i_node) = j            node_num(i,j) = i_node          enddo        enddo      return      end!-----------------------------------------------------------------------      subroutine list_elems (elem_nodes_in, node_num_in, &     &                       ni, nj, n_elems_in)        implicit none        integer ni, nj, n_elems_in, i, j, i_elem        integer node_num_in(ni,nj), elem_nodes_in(n_elems_in,3)        i_elem = 0        do j = 1, nj-1          do i = 1, ni-1! define the first element in this grid box            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+1)            elem_nodes_in(i_elem,3) = node_num_in(i,j+1)! define the second element in this grid box