sflux_subs5.f90

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

! 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