!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module monin_obukhov_mod !======================================================================= ! ! MONIN-OBUKHOV MODULE ! ! Routines for computing surface drag coefficients ! from data at the lowest model level ! and for computing the profile of fields ! between the lowest model level and the ground ! using Monin-Obukhov scaling ! !======================================================================= use constants_mod, only : grav, vonkarm use fms_mod, only: error_mesg, FATAL, file_exist, & check_nml_error, open_namelist_file, & mpp_pe, mpp_root_pe, close_file, stdlog, & write_version_number implicit none private !======================================================================= public monin_obukhov_init public monin_obukhov_end public mo_drag public mo_profile public mo_diff public stable_mix !======================================================================= interface mo_drag module procedure mo_drag_0d, mo_drag_1d, mo_drag_2d end interface interface mo_profile module procedure mo_profile_0d, mo_profile_1d, mo_profile_2d, & mo_profile_0d_n, mo_profile_1d_n, mo_profile_2d_n end interface interface mo_diff module procedure mo_diff_0d_n, mo_diff_0d_1, & mo_diff_1d_n, mo_diff_1d_1, & mo_diff_2d_n, mo_diff_2d_1 end interface interface stable_mix module procedure stable_mix_0d, stable_mix_1d, & stable_mix_2d, stable_mix_3d end interface !--------------------- version number --------------------------------- character(len=128) :: version = '$Id: monin_obukhov.F90,v 17.0 2009/07/21 02:55:38 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !======================================================================= ! DEFAULT VALUES OF NAMELIST PARAMETERS: real :: rich_crit = 2.0 real :: drag_min = 1.e-05 logical :: neutral = .false. integer :: stable_option = 1 real :: zeta_trans = 0.5 namelist /monin_obukhov_nml/ rich_crit, neutral, drag_min, & stable_option, zeta_trans !======================================================================= ! MODULE VARIABLES real, parameter :: small = 1.e-04 real :: b_stab, r_crit, sqrt_drag_min, lambda, rich_trans logical :: module_is_initialized = .false. contains !======================================================================= subroutine monin_obukhov_init integer :: unit, ierr, io, logunit !------------------- read namelist input ------------------------------- if (file_exist('input.nml')) then unit = open_namelist_file () ierr=1; do while (ierr /= 0) read (unit, nml=monin_obukhov_nml, iostat=io, end=10) ierr = check_nml_error(io,'monin_obukhov_nml') enddo 10 call close_file (unit) endif !---------- output namelist to log------------------------------------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) logunit = stdlog() write (logunit, nml=monin_obukhov_nml) endif !---------------------------------------------------------------------- if(rich_crit.le.0.25) call error_mesg( & 'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', & 'rich_crit in monin_obukhov_mod must be > 0.25', FATAL) if(drag_min.le.0.0) call error_mesg( & 'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', & 'drag_min in monin_obukhov_mod must be >= 0.0', FATAL) if(stable_option < 1 .or. stable_option > 2) call error_mesg( & 'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', & 'the only allowable values of stable_option are 1 and 2', FATAL) if(stable_option == 2 .and. zeta_trans < 0) call error_mesg( & 'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', & 'zeta_trans must be positive', FATAL) b_stab = 1.0/rich_crit r_crit = 0.95*rich_crit ! convergence can get slow if one is ! close to rich_crit sqrt_drag_min = 0.0 if(drag_min.ne.0.0) sqrt_drag_min = sqrt(drag_min) lambda = 1.0 + (5.0 - b_stab)*zeta_trans ! used only if stable_option = 2 rich_trans = zeta_trans/(1.0 + 5.0*zeta_trans) ! used only if stable_option = 2 module_is_initialized = .true. return end subroutine monin_obukhov_init !======================================================================= subroutine monin_obukhov_end module_is_initialized = .false. end subroutine monin_obukhov_end !======================================================================= subroutine mo_drag_1d & (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, & u_star, b_star, avail) real, intent(in) , dimension(:) :: pt, pt0, z, z0, zt, zq, speed real, intent(inout), dimension(:) :: drag_m, drag_t, drag_q, u_star, b_star logical, intent(in), optional, dimension(:) :: avail logical :: lavail, avail_dummy(1) integer :: n, ier integer, parameter :: max_iter = 20 real , parameter :: error=1.e-04, zeta_min=1.e-06, small=1.e-04 ! #include "monin_obukhov_interfaces.h" if(.not.module_is_initialized) call monin_obukhov_init n = size(pt) lavail = .false. if(present(avail)) lavail = .true. if(lavail) then if (count(avail) .eq. 0) return call monin_obukhov_drag_1d(grav, vonkarm, & & error, zeta_min, max_iter, small, & & neutral, stable_option, rich_crit, zeta_trans, drag_min, & & n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, & & drag_q, u_star, b_star, lavail, avail, ier) else call monin_obukhov_drag_1d(grav, vonkarm, & & error, zeta_min, max_iter, small, & & neutral, stable_option, rich_crit, zeta_trans, drag_min, & & n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, & & drag_q, u_star, b_star, lavail, avail_dummy, ier) endif end subroutine mo_drag_1d !======================================================================= subroutine mo_profile_1d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_t, del_q, avail) real, intent(in) :: zref, zref_t real, intent(in) , dimension(:) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out), dimension(:) :: del_m, del_t, del_q logical, intent(in) , optional, dimension(:) :: avail logical :: dummy_avail(1) integer :: n, ier ! #include "monin_obukhov_interfaces.h" if(.not. module_is_initialized) call monin_obukhov_init n = size(z) if(present(avail)) then if (count(avail) .eq. 0) return call monin_obukhov_profile_1d(vonkarm, & & neutral, stable_option, rich_crit, zeta_trans, & & n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, & & del_m, del_t, del_q, .true., avail, ier) else call monin_obukhov_profile_1d(vonkarm, & & neutral, stable_option, rich_crit, zeta_trans, & & n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, & & del_m, del_t, del_q, .false., dummy_avail, ier) endif end subroutine mo_profile_1d !======================================================================= subroutine stable_mix_3d(rich, mix) real, intent(in) , dimension(:,:,:) :: rich real, intent(out), dimension(:,:,:) :: mix integer :: n, ier if(.not. module_is_initialized) call monin_obukhov_init n = size(rich,1)*size(rich,2)*size(rich,3) call monin_obukhov_stable_mix(stable_option, rich_crit, zeta_trans, & & n, rich, mix, ier) end subroutine stable_mix_3d !======================================================================= subroutine mo_diff_2d_n(z, u_star, b_star, k_m, k_h) real, intent(in), dimension(:,:,:) :: z real, intent(in), dimension(:,:) :: u_star, b_star real, intent(out), dimension(:,:,:) :: k_m, k_h integer :: ni, nj, nk, ier real, parameter :: ustar_min = 1.e-10 if(.not.module_is_initialized) call monin_obukhov_init ni = size(z, 1); nj = size(z, 2); nk = size(z, 3) call monin_obukhov_diff(vonkarm, & & ustar_min, & & neutral, stable_option, rich_crit, zeta_trans, & & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier) end subroutine mo_diff_2d_n !======================================================================= ! The following routines are used by the public interfaces above !======================================================================= subroutine solve_zeta(rich, z, z0, zt, zq, f_m, f_t, f_q, mask) real , intent(in) , dimension(:) :: rich, z, z0, zt, zq logical, intent(in) , dimension(:) :: mask real , intent(out), dimension(:) :: f_m, f_t, f_q real, parameter :: error = 1.e-04 real, parameter :: zeta_min = 1.e-06 integer, parameter :: max_iter = 20 real :: max_cor integer :: iter real, dimension(size(rich(:))) :: & d_rich, rich_1, correction, corr, z_z0, z_zt, z_zq, & ln_z_z0, ln_z_zt, ln_z_zq, zeta, & phi_m, phi_m_0, phi_t, phi_t_0, rzeta, & zeta_0, zeta_t, zeta_q, df_m, df_t logical, dimension(size(rich(:))) :: mask_1 z_z0 = z/z0 z_zt = z/zt z_zq = z/zq ln_z_z0 = log(z_z0) ln_z_zt = log(z_zt) ln_z_zq = log(z_zq) corr = 0.0 mask_1 = mask ! initial guess where(mask_1) zeta = rich*ln_z_z0*ln_z_z0/ln_z_zt elsewhere zeta = 0.0 end where where (mask_1 .and. rich >= 0.0) zeta = zeta/(1.0 - rich/rich_crit) end where iter_loop: do iter = 1, max_iter where (mask_1 .and. abs(zeta).lt.zeta_min) zeta = 0.0 f_m = ln_z_z0 f_t = ln_z_zt f_q = ln_z_zq mask_1 = .false. ! don't do any more calculations at these pts end where where (mask_1) rzeta = 1.0/zeta zeta_0 = zeta/z_z0 zeta_t = zeta/z_zt zeta_q = zeta/z_zq elsewhere zeta_0 = 0.0 zeta_t = 0.0 zeta_q = 0.0 end where call mo_derivative_m(phi_m , zeta , mask_1) call mo_derivative_m(phi_m_0, zeta_0, mask_1) call mo_derivative_t(phi_t , zeta , mask_1) call mo_derivative_t(phi_t_0, zeta_t, mask_1) call mo_integral_m(f_m, zeta, zeta_0, ln_z_z0, mask_1) call mo_integral_tq(f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask_1) where (mask_1) df_m = (phi_m - phi_m_0)*rzeta df_t = (phi_t - phi_t_0)*rzeta rich_1 = zeta*f_t/(f_m*f_m) d_rich = rich_1*( rzeta + df_t/f_t - 2.0 *df_m/f_m) correction = (rich - rich_1)/d_rich corr = min(abs(correction),abs(correction/zeta)) ! the criterion corr < error seems to work ok, but is a bit arbitrary ! when zeta is small the tolerance is reduced end where max_cor= maxval(corr) if(max_cor > error) then mask_1 = mask_1 .and. (corr > error) ! change the mask so computation proceeds only on non-converged points where(mask_1) zeta = zeta + correction end where cycle iter_loop else return end if end do iter_loop call error_mesg ('solve_zeta in monin_obukhov_mod', & 'surface drag iteration did not converge', FATAL) end subroutine solve_zeta !======================================================================= subroutine mo_derivative_m(phi_m, zeta, mask) ! the differential similarity function for momentum real , intent(out), dimension(:) :: phi_m real , intent(in), dimension(:) :: zeta logical , intent(in), dimension(:) :: mask logical, dimension(size(zeta(:))) :: stable, unstable real , dimension(size(zeta(:))) :: x stable = mask .and. zeta >= 0.0 unstable = mask .and. zeta < 0.0 where (unstable) x = (1 - 16.0*zeta )**(-0.5) phi_m = sqrt(x) ! phi_m = (1 - 16.0*zeta)**(-0.25) end where if(stable_option == 1) then where (stable) phi_m = 1.0 + zeta *(5.0 + b_stab*zeta)/(1.0 + zeta) end where else if(stable_option == 2) then where (stable .and. zeta < zeta_trans) phi_m = 1 + 5.0*zeta end where where (stable .and. zeta >= zeta_trans) phi_m = lambda + b_stab*zeta end where endif return end subroutine mo_derivative_m !======================================================================= subroutine mo_derivative_t(phi_t, zeta, mask) ! the differential similarity function for buoyancy and tracers real , intent(out), dimension(:) :: phi_t real , intent(in), dimension(:) :: zeta logical , intent(in), dimension(:) :: mask logical, dimension(size(zeta(:))) :: stable, unstable stable = mask .and. zeta >= 0.0 unstable = mask .and. zeta < 0.0 where (unstable) phi_t = (1 - 16.0*zeta)**(-0.5) end where if(stable_option == 1) then where (stable) phi_t = 1.0 + zeta*(5.0 + b_stab*zeta)/(1.0 + zeta) end where else if(stable_option == 2) then where (stable .and. zeta < zeta_trans) phi_t = 1 + 5.0*zeta end where where (stable .and. zeta >= zeta_trans) phi_t = lambda + b_stab*zeta end where endif return end subroutine mo_derivative_t !======================================================================= subroutine mo_integral_tq (psi_t, psi_q, zeta, zeta_t, zeta_q, & ln_z_zt, ln_z_zq, mask) ! the integral similarity function for moisture and tracers real , intent(out), dimension(:) :: psi_t, psi_q real , intent(in), dimension(:) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq logical , intent(in), dimension(:) :: mask real, dimension(size(zeta(:))) :: x, x_t, x_q logical, dimension(size(zeta(:))) :: stable, unstable, & weakly_stable, strongly_stable stable = mask .and. zeta >= 0.0 unstable = mask .and. zeta < 0.0 where(unstable) x = sqrt(1 - 16.0*zeta) x_t = sqrt(1 - 16.0*zeta_t) x_q = sqrt(1 - 16.0*zeta_q) psi_t = ln_z_zt - 2.0*log( (1.0 + x)/(1.0 + x_t) ) psi_q = ln_z_zq - 2.0*log( (1.0 + x)/(1.0 + x_q) ) end where if( stable_option == 1) then where (stable) psi_t = ln_z_zt + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_t)) & + b_stab*(zeta - zeta_t) psi_q = ln_z_zq + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_q)) & + b_stab*(zeta - zeta_q) end where else if (stable_option == 2) then weakly_stable = stable .and. zeta <= zeta_trans strongly_stable = stable .and. zeta > zeta_trans where (weakly_stable) psi_t = ln_z_zt + 5.0*(zeta - zeta_t) psi_q = ln_z_zq + 5.0*(zeta - zeta_q) end where where(strongly_stable) x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans) endwhere where (strongly_stable .and. zeta_t <= zeta_trans) psi_t = ln_z_zt + x + 5.0*(zeta_trans - zeta_t) end where where (strongly_stable .and. zeta_t > zeta_trans) psi_t = lambda*ln_z_zt + b_stab*(zeta - zeta_t) endwhere where (strongly_stable .and. zeta_q <= zeta_trans) psi_q = ln_z_zq + x + 5.0*(zeta_trans - zeta_q) end where where (strongly_stable .and. zeta_q > zeta_trans) psi_q = lambda*ln_z_zq + b_stab*(zeta - zeta_q) endwhere end if return end subroutine mo_integral_tq !======================================================================= subroutine mo_integral_m (psi_m, zeta, zeta_0, ln_z_z0, mask) ! the integral similarity function for momentum real , intent(out), dimension(:) :: psi_m real , intent(in), dimension(:) :: zeta, zeta_0, ln_z_z0 logical , intent(in), dimension(:) :: mask real, dimension(size(zeta(:))) :: x, x_0, x1, x1_0, num, denom, y logical, dimension(size(zeta(:))) :: stable, unstable, & weakly_stable, strongly_stable stable = mask .and. zeta >= 0.0 unstable = mask .and. zeta < 0.0 where(unstable) x = sqrt(1 - 16.0*zeta) x_0 = sqrt(1 - 16.0*zeta_0) x = sqrt(x) x_0 = sqrt(x_0) x1 = 1.0 + x x1_0 = 1.0 + x_0 num = x1*x1*(1.0 + x*x) denom = x1_0*x1_0*(1.0 + x_0*x_0) y = atan(x) - atan(x_0) psi_m = ln_z_z0 - log(num/denom) + 2*y end where if( stable_option == 1) then where (stable) psi_m = ln_z_z0 + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_0)) & + b_stab*(zeta - zeta_0) end where else if (stable_option == 2) then weakly_stable = stable .and. zeta <= zeta_trans strongly_stable = stable .and. zeta > zeta_trans where (weakly_stable) psi_m = ln_z_z0 + 5.0*(zeta - zeta_0) end where where(strongly_stable) x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans) endwhere where (strongly_stable .and. zeta_0 <= zeta_trans) psi_m = ln_z_z0 + x + 5.0*(zeta_trans - zeta_0) end where where (strongly_stable .and. zeta_0 > zeta_trans) psi_m = lambda*ln_z_z0 + b_stab*(zeta - zeta_0) endwhere end if return end subroutine mo_integral_m !======================================================================= ! The following routines allow the public interfaces to be used ! with different dimensions of the input and output ! !======================================================================= subroutine mo_drag_2d & (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star) real, intent(in) , dimension(:,:) :: z, speed, pt, pt0, z0, zt, zq real, intent(out) , dimension(:,:) :: drag_m, drag_t, drag_q real, intent(inout), dimension(:,:) :: u_star, b_star integer :: j do j = 1, size(pt,2) call mo_drag_1d (pt(:,j), pt0(:,j), z(:,j), z0(:,j), zt(:,j), zq(:,j), & speed(:,j), drag_m(:,j), drag_t(:,j), drag_q(:,j), & u_star(:,j), b_star(:,j)) end do return end subroutine mo_drag_2d !======================================================================= subroutine mo_drag_0d & (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star) real, intent(in) :: z, speed, pt, pt0, z0, zt, zq real, intent(out) :: drag_m, drag_t, drag_q, u_star, b_star real, dimension(1) :: pt_1, pt0_1, z_1, z0_1, zt_1, zq_1, speed_1, & drag_m_1, drag_t_1, drag_q_1, u_star_1, b_star_1 pt_1 (1) = pt pt0_1 (1) = pt0 z_1 (1) = z z0_1 (1) = z0 zt_1 (1) = zt zq_1 (1) = zq speed_1(1) = speed call mo_drag_1d (pt_1, pt0_1, z_1, z0_1, zt_1, zq_1, speed_1, & drag_m_1, drag_t_1, drag_q_1, u_star_1, b_star_1) drag_m = drag_m_1(1) drag_t = drag_t_1(1) drag_q = drag_q_1(1) u_star = u_star_1(1) b_star = b_star_1(1) return end subroutine mo_drag_0d !======================================================================= subroutine mo_profile_2d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_h, del_q) real, intent(in) :: zref, zref_t real, intent(in) , dimension(:,:) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out), dimension(:,:) :: del_m, del_h, del_q integer :: j do j = 1, size(z,2) call mo_profile_1d (zref, zref_t, z(:,j), z0(:,j), zt(:,j), & zq(:,j), u_star(:,j), b_star(:,j), q_star(:,j), & del_m(:,j), del_h (:,j), del_q (:,j)) enddo return end subroutine mo_profile_2d !======================================================================= subroutine mo_profile_0d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_h, del_q) real, intent(in) :: zref, zref_t real, intent(in) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out) :: del_m, del_h, del_q real, dimension(1) :: z_1, z0_1, zt_1, zq_1, u_star_1, b_star_1, q_star_1, & del_m_1, del_h_1, del_q_1 z_1 (1) = z z0_1 (1) = z0 zt_1 (1) = zt zq_1 (1) = zq u_star_1(1) = u_star b_star_1(1) = b_star q_star_1(1) = q_star call mo_profile_1d (zref, zref_t, z_1, z0_1, zt_1, zq_1, & u_star_1, b_star_1, q_star_1, & del_m_1, del_h_1, del_q_1) del_m = del_m_1(1) del_h = del_h_1(1) del_q = del_q_1(1) return end subroutine mo_profile_0d !======================================================================= subroutine mo_profile_1d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_t, del_q, avail) real, intent(in), dimension(:) :: zref real, intent(in) , dimension(:) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out), dimension(:,:) :: del_m, del_t, del_q logical, intent(in) , optional, dimension(:) :: avail integer :: k do k = 1, size(zref(:)) if(present(avail)) then call mo_profile_1d (zref(k), zref(k), z, z0, zt, zq, & u_star, b_star, q_star, del_m(:,k), del_t(:,k), del_q(:,k), avail) else call mo_profile_1d (zref(k), zref(k), z, z0, zt, zq, & u_star, b_star, q_star, del_m(:,k), del_t(:,k), del_q(:,k)) endif enddo return end subroutine mo_profile_1d_n !======================================================================= subroutine mo_profile_0d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_t, del_q) real, intent(in), dimension(:) :: zref real, intent(in) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out), dimension(:) :: del_m, del_t, del_q integer :: k do k = 1, size(zref(:)) call mo_profile_0d (zref(k), zref(k), z, z0, zt, zq, & u_star, b_star, q_star, del_m(k), del_t(k), del_q(k)) enddo return end subroutine mo_profile_0d_n !======================================================================= subroutine mo_profile_2d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, & del_m, del_t, del_q) real, intent(in), dimension(:) :: zref real, intent(in), dimension(:,:) :: z, z0, zt, zq, u_star, b_star, q_star real, intent(out), dimension(:,:,:) :: del_m, del_t, del_q integer :: k do k = 1, size(zref(:)) call mo_profile_2d (zref(k), zref(k), z, z0, zt, zq, & u_star, b_star, q_star, del_m(:,:,k), del_t(:,:,k), del_q(:,:,k)) enddo return end subroutine mo_profile_2d_n !======================================================================= subroutine mo_diff_2d_1(z, u_star, b_star, k_m, k_h) real, intent(in), dimension(:,:) :: z, u_star, b_star real, intent(out), dimension(:,:) :: k_m, k_h real , dimension(size(z,1),size(z,2),1) :: z_n, k_m_n, k_h_n z_n(:,:,1) = z call mo_diff_2d_n(z_n, u_star, b_star, k_m_n, k_h_n) k_m = k_m_n(:,:,1) k_h = k_h_n(:,:,1) return end subroutine mo_diff_2d_1 !======================================================================= subroutine mo_diff_1d_1(z, u_star, b_star, k_m, k_h) real, intent(in), dimension(:) :: z, u_star, b_star real, intent(out), dimension(:) :: k_m, k_h real, dimension(size(z),1,1) :: z_n, k_m_n, k_h_n real, dimension(size(z),1) :: u_star_n, b_star_n z_n (:,1,1) = z u_star_n(:,1) = u_star b_star_n(:,1) = b_star call mo_diff_2d_n(z_n, u_star_n, b_star_n, k_m_n, k_h_n) k_m = k_m_n(:,1,1) k_h = k_h_n(:,1,1) return end subroutine mo_diff_1d_1 !======================================================================= subroutine mo_diff_1d_n(z, u_star, b_star, k_m, k_h) real, intent(in), dimension(:,:) :: z real, intent(in), dimension(:) :: u_star, b_star real, intent(out), dimension(:,:) :: k_m, k_h real, dimension(size(z,1),1) :: u_star2, b_star2 real, dimension(size(z,1),1, size(z,2)) :: z2, k_m2, k_h2 integer :: n do n = 1, size(z,2) z2 (:,1,n) = z(:,n) enddo u_star2(:,1) = u_star b_star2(:,1) = b_star call mo_diff_2d_n(z2, u_star2, b_star2, k_m2, k_h2) do n = 1, size(z,2) k_m(:,n) = k_m2(:,1,n) k_h(:,n) = k_h2(:,1,n) enddo return end subroutine mo_diff_1d_n !======================================================================= subroutine mo_diff_0d_1(z, u_star, b_star, k_m, k_h) real, intent(in) :: z, u_star, b_star real, intent(out) :: k_m, k_h integer :: ni, nj, nk, ier real, parameter :: ustar_min = 1.e-10 if(.not.module_is_initialized) call monin_obukhov_init ni = 1; nj = 1; nk = 1 call monin_obukhov_diff(vonkarm, & & ustar_min, & & neutral, stable_option, rich_crit, zeta_trans, & & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier) end subroutine mo_diff_0d_1 !======================================================================= subroutine mo_diff_0d_n(z, u_star, b_star, k_m, k_h) real, intent(in), dimension(:) :: z real, intent(in) :: u_star, b_star real, intent(out), dimension(:) :: k_m, k_h integer :: ni, nj, nk, ier real, parameter :: ustar_min = 1.e-10 if(.not.module_is_initialized) call monin_obukhov_init ni = 1; nj = 1; nk = size(z(:)) call monin_obukhov_diff(vonkarm, & & ustar_min, & & neutral, stable_option, rich_crit, zeta_trans, & & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier) end subroutine mo_diff_0d_n !======================================================================= subroutine stable_mix_2d(rich, mix) real, intent(in) , dimension(:,:) :: rich real, intent(out), dimension(:,:) :: mix real, dimension(size(rich,1),size(rich,2),1) :: rich_3d, mix_3d rich_3d(:,:,1) = rich call stable_mix_3d(rich_3d, mix_3d) mix = mix_3d(:,:,1) return end subroutine stable_mix_2d !======================================================================= subroutine stable_mix_1d(rich, mix) real, intent(in) , dimension(:) :: rich real, intent(out), dimension(:) :: mix real, dimension(size(rich),1,1) :: rich_3d, mix_3d rich_3d(:,1,1) = rich call stable_mix_3d(rich_3d, mix_3d) mix = mix_3d(:,1,1) return end subroutine stable_mix_1d !======================================================================= subroutine stable_mix_0d(rich, mix) real, intent(in) :: rich real, intent(out) :: mix real, dimension(1,1,1) :: rich_3d, mix_3d rich_3d(1,1,1) = rich call stable_mix_3d(rich_3d, mix_3d) mix = mix_3d(1,1,1) return end subroutine stable_mix_0d !======================================================================= end module monin_obukhov_mod