module diffusivity_mod !======================================================================= ! ! DIFFUSIVITY MODULE ! ! Routines for computing atmospheric diffusivities in the ! planetary boundary layer and in the free atmosphere ! !======================================================================= use constants_mod, only : grav, vonkarm, cp_air, rdgas, rvgas use mpp_mod, only : input_nml_file use fms_mod, only : error_mesg, FATAL, file_exist, & check_nml_error, open_namelist_file, & mpp_pe, mpp_root_pe, close_file, & write_version_number, stdlog use monin_obukhov_mod, only : mo_diff implicit none private ! public interfaces !======================================================================= public diffusivity, pbl_depth, molecular_diff !======================================================================= ! form of iterfaces !======================================================================= ! subroutine diffusivity (t, q, u, v, p_full, p_half, z_full, z_half, ! u_star, b_star, h, k_m, k_t) ! input: ! t : real, dimension(:,:,:) -- (:,:,pressure), third index running ! from top of atmosphere to bottom ! temperature (K) ! ! q : real, dimension(:,:,:) ! water vapor specific humidity (nondimensional) ! ! u : real, dimension(:,:) ! zonal wind (m/s) ! ! v : real, dimension(:,:,:) ! meridional wind (m/s) ! ! z_full : real, dimension(:,:,: ! height of full levels (m) ! 1 = top of atmosphere; size(p_half,3) = surface ! size(z_full,3) = size(t,3) ! ! z_half : real, dimension(:,:,:) ! height of half levels (m) ! size(z_half,3) = size(t,3) +1 ! z_half(:,:,size(z_half,3)) must be height of surface! ! (if you are not using eta-model) ! ! u_star: real, dimension(:,:) ! friction velocity (m/s) ! ! b_star: real, dimension(:,:) ! buoyancy scale (m/s**2) ! (u_star and b_star can be obtained by calling ! mo_drag in monin_obukhov_mod) ! output: ! h : real, dimension(:,:,) ! depth of planetary boundary layer (m) ! ! k_m : real, dimension(:,:,:) ! diffusivity for momentum (m**2/s) ! ! defined at half-levels ! size(k_m,3) should be at least as large as size(t,3) ! only the returned values at ! levels 2 to size(t,3) are meaningful ! other values will be returned as zero ! ! k_t : real, dimension(:,:,:) ! diffusivity for temperature and scalars (m**2/s) ! ! !======================================================================= !--------------------- version number ---------------------------------- character(len=128) :: version = '$Id: diffusivity.F90,v 19.0 2012/01/06 20:05:43 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' !======================================================================= ! DEFAULT VALUES OF NAMELIST PARAMETERS: logical :: fixed_depth = .false. real :: depth_0 = 5000.0 real :: frac_inner = 0.1 real :: rich_crit_pbl = 1.0 real :: entr_ratio = 0.2 real :: parcel_buoy = 2.0 real :: znom = 1000.0 logical :: free_atm_diff = .false. logical :: free_atm_skyhi_diff = .false. logical :: pbl_mcm = .false. real :: rich_crit_diff = 0.25 real :: mix_len = 30. real :: rich_prandtl = 1.00 real :: background_m = 0.0 real :: background_t = 0.0 logical :: ampns = .false. ! include delta z factor in ! defining ri ? real :: ampns_max = 1.0E20 ! limit to reduction factor ! applied to ri due to delta z ! factor logical :: do_entrain =.true. logical :: do_simple =.false. namelist /diffusivity_nml/ fixed_depth, depth_0, frac_inner,& rich_crit_pbl, entr_ratio, parcel_buoy,& znom, free_atm_diff, free_atm_skyhi_diff,& pbl_mcm, rich_crit_diff, mix_len, rich_prandtl,& background_m, background_t, ampns, ampns_max, & do_entrain, do_simple !======================================================================= ! OTHER MODULE VARIABLES real :: small = 1.e-04 real :: gcp = grav/cp_air logical :: module_is_initialized = .false. real :: beta = 1.458e-06 real :: rbop1 = 110.4 real :: rbop2 = 1.405 real, parameter :: d608 = (rvgas-rdgas)/rdgas contains !======================================================================= subroutine diffusivity_init integer :: unit, ierr, io, logunit !------------------- read namelist input ------------------------------- if (file_exist('input.nml')) then #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=diffusivity_nml, iostat=io) ierr = check_nml_error(io,"diffusivity_nml") #else unit = open_namelist_file () ierr=1; do while (ierr /= 0) read (unit, nml=diffusivity_nml, iostat=io, end=10) ierr = check_nml_error(io,'diffusivity_nml') enddo 10 call close_file (unit) #endif !------------------- dummy checks -------------------------------------- if (frac_inner .le. 0. .or. frac_inner .ge. 1.) & call error_mesg ('diffusivity_init', & 'frac_inner must be between 0 and 1', FATAL) if (rich_crit_pbl .lt. 0.) & call error_mesg ('diffusivity_init', & 'rich_crit_pbl must be greater than or equal to zero', FATAL) if (entr_ratio .lt. 0.) & call error_mesg ('diffusivity_init', & 'entr_ratio must be greater than or equal to zero', FATAL) if (znom .le. 0.) & call error_mesg ('diffusivity_init', & 'znom must be greater than zero', FATAL) if (.not.free_atm_diff .and. free_atm_skyhi_diff)& call error_mesg ('diffusivity_init', & 'free_atm_diff must be set to true if '//& 'free_atm_skyhi_diff = .true.', FATAL) if (rich_crit_diff .le. 0.) & call error_mesg ('diffusivity_init', & 'rich_crit_diff must be greater than zero', FATAL) if (mix_len .lt. 0.) & call error_mesg ('diffusivity_init', & 'mix_len must be greater than or equal to zero', FATAL) if (rich_prandtl .lt. 0.) & call error_mesg ('diffusivity_init', & 'rich_prandtl must be greater than or equal to zero', FATAL) if (background_m .lt. 0.) & call error_mesg ('diffusivity_init', & 'background_m must be greater than or equal to zero', FATAL) if (background_t .lt. 0.) & call error_mesg ('diffusivity_init', & 'background_t must be greater than or equal to zero', FATAL) if (ampns_max .lt. 1.) & call error_mesg ('diffusivity_init', & 'ampns_max must be greater than or equal to one', FATAL) if (ampns .and. .not. free_atm_skyhi_diff) & call error_mesg ('diffusivity_init', & 'ampns is only valid when free_atm_skyhi_diff is & & also true', FATAL) endif !end of reading input.nml !---------- output namelist to log------------------------------------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) logunit = stdlog() write (logunit, nml=diffusivity_nml) endif module_is_initialized = .true. return end subroutine diffusivity_init !======================================================================= subroutine diffusivity_end module_is_initialized = .false. end subroutine diffusivity_end !======================================================================= subroutine diffusivity(t, q, u, v, p_full, p_half, z_full, z_half, & u_star, b_star, h, k_m, k_t, kbot) real, intent(in), dimension(:,:,:) :: t, q, u, v real, intent(in), dimension(:,:,:) :: p_full, p_half real, intent(in), dimension(:,:,:) :: z_full, z_half real, intent(in), dimension(:,:) :: u_star, b_star real, intent(inout), dimension(:,:,:) :: k_m, k_t real, intent(out), dimension(:,:) :: h integer, intent(in), optional, dimension(:,:) :: kbot real, dimension(size(t,1),size(t,2),size(t,3)) :: svcp,z_full_ag, & k_m_save, k_t_save real, dimension(size(t,1),size(t,2),size(t,3)+1):: z_half_ag real, dimension(size(t,1),size(t,2)) :: z_surf integer :: i,j,k,nlev,nlat,nlon if(.not.module_is_initialized) call diffusivity_init nlev = size(t,3) k_m_save = k_m k_t_save = k_t !compute height of surface if (present(kbot)) then nlat = size(t,2) nlon = size(t,1) do j=1,nlat do i=1,nlon z_surf(i,j) = z_half(i,j,kbot(i,j)+1) enddo enddo else z_surf(:,:) = z_half(:,:,nlev+1) end if !compute density profile, and heights relative to surface do k = 1, nlev z_full_ag(:,:,k) = z_full(:,:,k) - z_surf(:,:) z_half_ag(:,:,k) = z_half(:,:,k) - z_surf(:,:) if(do_simple) then svcp(:,:,k) = t(:,:,k) + gcp*(z_full_ag(:,:,k)) else svcp(:,:,k) = t(:,:,k)*(1. + d608*q(:,:,k)) + gcp*(z_full_ag(:,:,k)) endif end do z_half_ag(:,:,nlev+1) = z_half(:,:,nlev+1) - z_surf(:,:) if(fixed_depth) then h = depth_0 else call pbl_depth(svcp,u,v,z_full_ag,u_star,b_star,h,kbot=kbot) end if if(pbl_mcm) then call diffusivity_pbl_mcm (u,v, t, p_full, p_half, & z_full_ag, z_half_ag, h, k_m, k_t) else call diffusivity_pbl (svcp, u, v, z_half_ag, h, u_star, b_star,& k_m, k_t, kbot=kbot) end if if(free_atm_diff) & call diffusivity_free (svcp, u, v, z_full_ag, z_half_ag, h, k_m, k_t) k_m = k_m + k_m_save k_t = k_t + k_t_save !NOTE THAT THIS LINE MUST FOLLOW DIFFUSIVITY_FREE SO THAT ENTRAINMENT !K's DO NOT GET OVERWRITTEN IN DIFFUSIVITY_FREE SUBROUTINE if(entr_ratio .gt. 0. .and. .not. fixed_depth .and. do_entrain) & call diffusivity_entr(svcp,z_full_ag,h,u_star,b_star,k_m,k_t) !set background diffusivities if(background_m.gt.0.0) k_m = max(k_m,background_m) if(background_t.gt.0.0) k_t = max(k_t,background_t) return end subroutine diffusivity !======================================================================= subroutine pbl_depth(t, u, v, z, u_star, b_star, h, kbot) real, intent(in) , dimension(:,:,:) :: t, u, v, z real, intent(in) , dimension(:,:) :: u_star,b_star real, intent(out), dimension(:,:) :: h integer,intent(in) , optional, dimension(:,:) :: kbot real, dimension(size(t,1),size(t,2),size(t,3)) :: rich real, dimension(size(t,1),size(t,2)) :: ws,k_t_ref,& h_inner,tbot real :: rich1, rich2,& h1,h2,svp,t1,t2 integer, dimension(size(t,1),size(t,2)) :: ibot integer :: i,j,k,nlon,& nlat, nlev nlev = size(t,3) nlat = size(t,2) nlon = size(t,1) !assign ibot, compute tbot (virtual temperature at lowest level) if (present(kbot)) then ibot(:,:) = kbot do j = 1,nlat do i = 1,nlon tbot(i,j) = t(i,j,ibot(i,j)) enddo enddo else ibot(:,:) = nlev tbot(:,:) = t(:,:,nlev) end if !compute richardson number for use in pbl depth of neutral/stable side do k = 1,nlev rich(:,:,k) = z(:,:,k)*grav*(t(:,:,k)-tbot(:,:))/tbot(:,:)& /(u(:,:,k)*u(:,:,k) + v(:,:,k)*v(:,:,k) + small ) end do !compute ws to be used in evaluating parcel buoyancy !ws = u_star / phi(h_inner,u_star,b_star) . To find phi !a call to mo_diff is made. h_inner(:,:)=frac_inner*znom call mo_diff(h_inner, u_star, b_star, ws, k_t_ref) ws = max(small,ws/vonkarm/h_inner) do j = 1, nlat do i = 1, nlon !do neutral or stable case if (b_star(i,j).le.0. .or. do_simple) then h1 = z(i,j,ibot(i,j)) h(i,j) = h1 rich1 = rich(i,j,ibot(i,j)) do k = ibot(i,j)-1, 1, -1 rich2 = rich(i,j,k) h2 = z(i,j,k) if(rich2.gt.rich_crit_pbl) then h(i,j) = h2 + (h1 - h2)*(rich2 - rich_crit_pbl)& /(rich2 - rich1 ) go to 10 endif rich1 = rich2 h1 = h2 enddo !do unstable case else svp = tbot(i,j)*(1.+ & (parcel_buoy*u_star(i,j)*b_star(i,j)/grav/ws(i,j)) ) h1 = z(i,j,ibot(i,j)) h(i,j) = h1 t1 = tbot(i,j) do k = ibot(i,j)-1 , 1, -1 h2 = z(i,j,k) t2 = t(i,j,k) if (t2.gt.svp) then h(i,j) = h2 + (h1 - h2)*(t2 - svp)/(t2 - t1 ) go to 10 end if h1 = h2 t1 = t2 enddo end if 10 continue enddo enddo return end subroutine pbl_depth !======================================================================= subroutine diffusivity_pbl(t, u, v, z_half, h, u_star, b_star, & k_m, k_t, kbot) real, intent(in) , dimension(:,:,:) :: t, u, v, z_half real, intent(in) , dimension(:,:) :: h, u_star, b_star real, intent(inout) , dimension(:,:,:) :: k_m, k_t integer, intent(in) , optional, dimension(:,:) :: kbot real, dimension(size(t,1),size(t,2)) :: h_inner, k_m_ref,& k_t_ref, factor real, dimension(size(t,1),size(t,2),size(t,3)+1) :: zm real :: h_inner_max integer :: k, kk, nlev nlev = size(t,3) !assign z_half to zm, and set to zero any values of zm < 0. !the setting to zero is necessary so that when using eta model !below ground half levels will have zero k_m and k_t zm = z_half if (present(kbot)) then where(zm < 0.) zm = 0. end where end if h_inner = frac_inner*h h_inner_max = maxval(h_inner) kk = nlev do k = 2, nlev if( minval(zm(:,:,k)) < h_inner_max) then kk = k exit end if end do k_m = 0.0 k_t = 0.0 call mo_diff(h_inner , u_star, b_star, k_m_ref , k_t_ref) call mo_diff(zm(:,:,kk:nlev), u_star, b_star, k_m(:,:,kk:nlev), k_t(:,:,kk:nlev)) do k = 2, nlev where(zm(:,:,k) >= h_inner .and. zm(:,:,k) < h) factor = (zm(:,:,k)/h_inner)* & (1.0 - (zm(:,:,k) - h_inner)/(h - h_inner))**2 k_m(:,:,k) = k_m_ref*factor k_t(:,:,k) = k_t_ref*factor end where end do return end subroutine diffusivity_pbl !======================================================================= subroutine diffusivity_pbl_mcm(u, v, t, p_full, p_half, z_full, z_half, & h, k_m, k_t) real, intent(in) , dimension(:,:,:) :: u, v, t, z_full, z_half real, intent(in) , dimension(:,:,:) :: p_full, p_half real, intent(in) , dimension(:,:) :: h real, intent(inout) , dimension(:,:,:) :: k_m, k_t integer :: k, nlev real, dimension(size(z_full,1),size(z_full,2)) :: delta_u, delta_v, delta_z real :: htcrit_ss real :: h_ss real, dimension(size(z_full,1),size(z_full,2)) :: sig_half, z_half_ss, elmix_ss ! htcrit_ss = height at which mixing length is a maximum (75m) ! h_ss = height at which mixing length vanishes (4900m) ! elmix_ss = mixing length ! Define some constants: ! salaps = standard atmospheric lapse rate (K/m) ! tsfc = idealized global mean surface temperature (15C) real :: tsfc = 288.16 real :: salaps = -6.5e-3 nlev = size(z_full,3) k_m = 0. h_ss = depth_0 htcrit_ss = frac_inner*h_ss do k = 2, nlev ! TK mods 8/13/01: (code derived from SS) ! Compute the height of each half level assuming a constant ! standard lapse rate using the above procedure. ! WARNING: These should be used with caution. They will ! have large errors above the tropopause. ! In order to determine the height, the layer mean temperature ! from the surface to that level is required. A surface ! temperature of 15 deg Celsius and a standard lapse rate of ! -6.5 deg/km will be used to estimate an average temperature ! profile. sig_half = p_half(:,:,k)/p_half(:,:,nlev+1) z_half_ss = -rdgas * .5*(tsfc+tsfc*(sig_half**(-rdgas*salaps/grav))) * alog(sig_half)/grav !compute mixing length as in SS (no geographical variation) elmix_ss = 0. where (z_half_ss < htcrit_ss .and. z_half_ss > 0.) elmix_ss = vonkarm*z_half_ss endwhere where (z_half_ss >= htcrit_ss .and. z_half_ss < h_ss) elmix_ss = vonkarm*htcrit_ss*(h_ss-z_half_ss)/(h_ss-htcrit_ss) endwhere delta_z = rdgas*0.5*(t(:,:,k)+t(:,:,k-1))*(p_full(:,:,k)-p_full(:,:,k-1))/& (grav*p_half(:,:,k)) delta_u = u(:,:,k-1) - u(:,:,k) delta_v = v(:,:,k-1) - v(:,:,k) k_m(:,:,k) = elmix_ss * elmix_ss *& sqrt(delta_u*delta_u + delta_v*delta_v)/delta_z end do k_t = k_m return end subroutine diffusivity_pbl_mcm !======================================================================= subroutine diffusivity_free(t, u, v, z, zz, h, k_m, k_t) real, intent(in) , dimension(:,:,:) :: t, u, v, z, zz real, intent(in) , dimension(:,:) :: h real, intent(inout) , dimension(:,:,:) :: k_m, k_t real, dimension(size(t,1),size(t,2)) :: dz, b, speed2, rich, fri, & alpz, fri2 integer :: k do k = 2, size(t,3) !---------------------------------------------------------------------- ! define the richardson number. set it to zero if it is negative. save ! a copy of it for later use (rich2). !---------------------------------------------------------------------- dz = z(:,:,k-1) - z(:,:,k) b = grav*(t(:,:,k-1)-t(:,:,k))/t(:,:,k) speed2 = (u(:,:,k-1) - u(:,:,k))**2 + (v(:,:,k-1) - v(:,:,k))**2 rich= b*dz/(speed2+small) rich = max(rich, 0.0) if (free_atm_skyhi_diff) then !--------------------------------------------------------------------- ! limit the standard richardson number to between 0 and the critical ! value (rich2). compute the richardson number factor needed in the ! eddy mixing coefficient using this standard richardson number. !--------------------------------------------------------------------- where (rich(:,:) >= rich_crit_diff) fri2(:,:) = 0.0 elsewhere fri2(:,:) = (1.0 - rich/rich_crit_diff)**2 endwhere endif !--------------------------------------------------------------------- ! if ampns is activated, compute the delta z factor. define rich ! including this factor. !--------------------------------------------------------------------- if (ampns) then alpz(:,:) = MIN ( (1. + 1.e-04*(dz(:,:)**1.5)), ampns_max) rich(:,:) = rich(:,:) / alpz(:,:) endif !--------------------------------------------------------------------- ! compute the richardson number factor to be used in the eddy ! mixing coefficient. if ampns is on, this value includes it; other- ! wise it does not. !--------------------------------------------------------------------- fri(:,:) = (1.0 - rich/rich_crit_diff)**2 !--------------------------------------------------------------------- ! compute the eddy mixing coefficients in the free atmosphere ( zz ! > h). in the non-ampns case, values are obtained only when the ! standard richardson number is sub-critical; in the ampns case values ! are obtained only when the richardson number computed with the ! ampns factor is sub critical. when the ampns factor is activated, ! it is also included in the mixing coefficient. the value of mixing ! for temperature, etc. is reduced dependent on the ri stability ! factor calculated without the ampns factor. !--------------------------------------------------------------------- if (free_atm_skyhi_diff) then !--------------------------------------------------------------------- ! this is the skyhi-like formulation -- possible ampns factor, ratio ! of k_m to k_t defined based on computed stability factor. !--------------------------------------------------------------------- if (ampns) then where (rich < rich_crit_diff .and. zz(:,:,k) > h) k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)* & ( 1. + 1.e-04*(dz(:,:)**1.5))/dz k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:)) end where else where (rich < rich_crit_diff .and. zz(:,:,k) > h) k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:)) end where endif else !--------------------------------------------------------------------- ! this is the non-skyhi-like formulation -- no ampns factor, ratio ! of k_m to k_t defined by rich_prandtl. !--------------------------------------------------------------------- where (rich < rich_crit_diff .and. zz(:,:,k) > h) k_t(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz k_m(:,:,k) = k_t(:,:,k)*rich_prandtl end where end if end do end subroutine diffusivity_free !======================================================================= subroutine molecular_diff ( temp, press, k_m, k_t) real, intent(in), dimension (:,:,:) :: temp, press real, intent(inout), dimension (:,:,:) :: k_m, k_t real, dimension (size(temp,1), size(temp,2)) :: temp_half, & rho_half, rbop2d integer :: k !--------------------------------------------------------------------- do k=2,size(temp,3) temp_half(:,:) = 0.5*(temp(:,:,k) + temp(:,:,k-1)) rho_half(:,:) = press(:,:,k)/(rdgas*temp_half(:,:) ) rbop2d(:,:) = beta*temp_half(:,:)*sqrt(temp_half(:,:))/ & (rho_half(:,:)*(temp_half(:,:)+rbop1)) k_m(:,:,k) = rbop2d(:,:) k_t(:,:,k) = rbop2d(:,:)*rbop2 end do k_m(:,:,1) = 0.0 k_t(:,:,1) = 0.0 end subroutine molecular_diff !======================================================================= subroutine diffusivity_entr(t, z, h, u_star, b_star, k_m, k_t) real, intent(in) , dimension(:,:,:) :: t, z real, intent(in) , dimension(:,:) :: h, u_star, b_star real, intent(inout) , dimension(:,:,:) :: k_m, k_t integer :: k, nlev nlev=size(t,3) do k = 2,nlev where (b_star .gt. 0. .and. z(:,:,k-1) .gt. h .and. & z(:,:,k) .le. h) k_t(:,:,k) = (z(:,:,k-1)-z(:,:,k))*entr_ratio*t(:,:,k)* & u_star*b_star/grav/max(small,t(:,:,k-1)-t(:,:,k)) k_m(:,:,k) = k_t(:,:,k) end where enddo end subroutine diffusivity_entr !======================================================================= end module diffusivity_mod