diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/CMakeLists.txt index cd547dbcc..fa38b61a1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/CMakeLists.txt @@ -6,6 +6,7 @@ set (srcs shoc.F90 edmf.F90 scm_surface.F90 + pblh_library.F90 ) esma_add_library (${this} diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 index 5f6fdc53c..55830f7b0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/GEOS_TurbulenceGridComp.F90 @@ -19,7 +19,8 @@ module GEOS_TurbulenceGridCompMod use shoc use edmf_mod, only: run_edmf,mfparams use scm_surface, only : surface_layer, surface - + use GEOSturbulence_PBLH_Library + #ifdef _CUDA use cudafor #endif @@ -1712,15 +1713,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'planetary_boundary_layer_height_horiz_tke', & - SHORT_NAME = 'ZPBLHTKE', & - UNITS = 'm', & - DIMS = MAPL_DimsHorzOnly, & - VLOCATION = MAPL_VLocationNone, & - RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & LONG_NAME = 'turbulent_kinetic_energy', & SHORT_NAME = 'TKE', & @@ -1775,6 +1767,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'wind_turning_angle_across_boundary_layer', & + SHORT_NAME = 'TURNANG', & + UNITS = '1', & + DIMS = MAPL_DimsHorzOnly, & + VLOCATION = MAPL_VLocationNone, & + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'surface_based_inversion_frequency', & SHORT_NAME = 'SBIFRQ', & @@ -2937,7 +2938,6 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(:,: ), pointer :: TCZPBL => null() real, dimension(:,: ), pointer :: ZPBL2 => null() real, dimension(:,: ), pointer :: ZPBL10P => null() - real, dimension(:,: ), pointer :: ZPBLHTKE => null() real, dimension(:,:,:), pointer :: TKE => null() real, dimension(:,: ), pointer :: ZPBLRI => null() real, dimension(:,: ), pointer :: ZPBLRI2 => null() @@ -2946,6 +2946,7 @@ subroutine REFRESH(IM,JM,LM,RC) real, dimension(:,: ), pointer :: ZPBLRFRCT => null() real, dimension(:,: ), pointer :: SBIFRQ => null() real, dimension(:,: ), pointer :: SBITOP => null() + real, dimension(:,: ), pointer :: TURNANG => null() real, dimension(:,: ), pointer :: KPBL => null() real, dimension(:,: ), pointer :: KPBL_SC => null() real, dimension(:,: ), pointer :: ZPBL_SC => null() @@ -3298,9 +3299,7 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, ZPBL10p, 'ZPBL10p', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, ZPBLHTKE, 'ZPBLHTKE', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, TKE, 'TKE', RC=STATUS) + call MAPL_GetPointer(EXPORT, TKE, 'TKE', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, ZPBLRI, 'ZPBLRI', RC=STATUS) VERIFY_(STATUS) @@ -3308,13 +3307,15 @@ subroutine REFRESH(IM,JM,LM,RC) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, ZPBLTHV, 'ZPBLTHV', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, ZPBLQV, 'ZPBLQV', RC=STATUS) + call MAPL_GetPointer(EXPORT, ZPBLQV, 'ZPBLQV', RC=STATUS) + VERIFY_(STATUS) + call MAPL_GetPointer(EXPORT, ZPBLRFRCT, 'ZPBLRFRCT', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, ZPBLRFRCT, 'ZPBLRFRCT', RC=STATUS) + call MAPL_GetPointer(EXPORT, SBIFRQ, 'SBIFRQ', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, SBIFRQ, 'SBIFRQ', RC=STATUS) + call MAPL_GetPointer(EXPORT, SBITOP, 'SBITOP', RC=STATUS) VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, SBITOP, 'SBITOP', RC=STATUS) + call MAPL_GetPointer(EXPORT, TURNANG, 'TURNANG', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, LWCRT, 'LWCRT', ALLOC=.TRUE., RC=STATUS) VERIFY_(STATUS) @@ -4501,290 +4502,39 @@ subroutine REFRESH(IM,JM,LM,RC) KPBLMIN = count(PREF < 50000.) - ZPBL = MAPL_UNDEF - if (associated(PPBL)) PPBL = MAPL_UNDEF + if (CALC_TCZPBL) call find_bulk_ri_pblh(im,jm,lm,u,v,z,thv,tczpbl,kpbltc,tcri_crit) - if (CALC_TCZPBL) then - TCZPBL = MAPL_UNDEF - thetavs = T(:,:,LM)*(1.0+MAPL_VIREPS*Q(:,:,LM)/(1.0-Q(:,:,LM)))*(TH(:,:,LM)/T(:,:,LM)) - tcrib(:,:,LM) = 0.0 - do I = 1, IM - do J = 1, JM - do L=LM-1,1,-1 - thetavh(I,J) = T(I,J,L)*(1.0+MAPL_VIREPS*Q(I,J,L)/(1.0-Q(I,J,L)))*(TH(I,J,L)/T(I,J,L)) - uv2h(I,J) = max(U(I,J,L)**2+V(I,J,L)**2,1.0E-8) - tcrib(I,J,L) = MAPL_GRAV*(thetavh(I,J)-thetavs(I,J))*Z(I,J,L)/(thetavs(I,J)*uv2h(I,J)) - if (tcrib(I,J,L) >= tcri_crit) then - TCZPBL(I,J) = Z(I,J,L+1)+(tcri_crit-tcrib(I,J,L+1))/(tcrib(I,J,L)-tcrib(I,J,L+1))*(Z(I,J,L)-Z(I,J,L+1)) - KPBLTC(I,J) = float(L) - exit - end if - end do - end do - end do - where (TCZPBL<0.) - TCZPBL = Z(:,:,LM) - KPBLTC = float(LM) - end where - end if ! CALC_TCZPBL - - if (CALC_ZPBL2) then - ZPBL2 = MAPL_UNDEF - - do I = 1, IM - do J = 1, JM - do L=LM,2,-1 - if ((KH(I,J,L-1) < 2.).and.(KH(I,J,L) >= 2.).and.(ZPBL2(I,J)==MAPL_UNDEF)) then - ZPBL2(I,J) = Z(I,J,L) - KPBL2(I,J) = float(L) - end if - end do - end do - end do + if (CALC_ZPBL2) call find_kh2_pblh(im,jm,lm,kpblmin,z,kh,zpbl2,kpbl2) - where ( ZPBL2 .eq. MAPL_UNDEF ) - ZPBL2 = Z(:,:,LM) - KPBL2 = float(LM) - end where - ZPBL2 = MIN(ZPBL2,Z(:,:,KPBLMIN)) - end if ! CALC_ZPBL2 + if (CALC_ZPBL10p) call find_kh10p_pblh(im,jm,lm,kpblmin,z,zl0,kh,zpbl10p,kpbl10p) - if (CALC_ZPBL10p) then - ZPBL10p = MAPL_UNDEF - - do I = 1, IM - do J = 1, JM - temparray(1:LM+1) = KH(I,J,0:LM) - do L = LM,2,-1 - locmax = maxloc(temparray,1) - minlval = max(0.001,0.0001*maxval(temparray)) - if(temparray(locmax-1)maxkh) maxkh = temparray(L) - if(temparray(L-1)= 0.1*maxkh) & - .and. (ZPBL10p(I,J) == MAPL_UNDEF ) ) then - ZPBL10p(I,J) = ZL0(I,J,L)+ & - ((ZL0(I,J,L-1)-ZL0(I,J,L))/(temparray(L)-temparray(L+1))) * (0.1*maxkh-temparray(L+1)) - KPBL10p(I,J) = float(L) - end if - end do - if ( ZPBL10p(I,J) .eq. MAPL_UNDEF .or. (maxkh.lt.1.)) then - ZPBL10p(I,J) = Z(I,J,LM) - KPBL10p(I,J) = float(LM) - endif - end do - end do + ! RI local diagnostic for pbl height threshold 0. + if (associated(ZPBLRI)) call find_ri_pblh(im,jm,lm,kpblmin,z,ri,zpblri,ri_crit) - ZPBL10p = MIN(ZPBL10p,Z(:,:,KPBLMIN)) - end if ! CALC_ZPBL10p + ! RI local diagnostic for pbl height threshold 0.2 + if (associated(ZPBLRI2)) call find_ri_pblh(im,jm,lm,kpblmin,z,ri,zpblri2,ri_crit2) - ! HTKE pbl height - if (associated(ZPBLHTKE)) then - ZPBLHTKE = MAPL_UNDEF - end if ! ZPBLHTKE + ! Thetav gradient based pbl height diagnostic + if (associated(ZPBLTHV)) call find_thv_pblh(im,jm,lm,kpblmin,z,thv,zpblthv) - ! RI local diagnostic for pbl height thresh 0. - if (associated(ZPBLRI)) then - ZPBLRI = MAPL_UNDEF - where (RI(:,:,LM-1)>ri_crit) ZPBLRI = Z(:,:,LM) + ! Refractivity-based PBLH + if (associated(ZPBLRFRCT)) call find_rfrct_pblh(im,jm,lm,z,plo,t,q,zpblrfrct) - do I = 1, IM - do J = 1, JM - do L=LM-1,1,-1 - if( (RI(I,J,L-1)>ri_crit) .and. (ZPBLRI(I,J) == MAPL_UNDEF) ) then - ZPBLRI(I,J) = Z(I,J,L+1)+(ri_crit-RI(I,J,L))/(RI(I,J,L-1)-RI(I,J,L))*(Z(I,J,L)-Z(I,J,L+1)) - end if - end do - end do - end do - - where ( ZPBLRI .eq. MAPL_UNDEF ) ZPBLRI = Z(:,:,LM) - ZPBLRI = MIN(ZPBLRI,Z(:,:,KPBLMIN)) - where ( ZPBLRI < 0.0 ) ZPBLRI = Z(:,:,LM) - end if ! ZPBLRI - - ! RI local diagnostic for pbl height thresh 0.2 - if (associated(ZPBLRI2)) then - ZPBLRI2 = MAPL_UNDEF - where (RI(:,:,LM-1) > ri_crit2) ZPBLRI2 = Z(:,:,LM) - - do I = 1, IM - do J = 1, JM - do L=LM-1,1,-1 - if( (RI(I,J,L-1)>ri_crit2) .and. (ZPBLRI2(I,J) == MAPL_UNDEF) ) then - ZPBLRI2(I,J) = Z(I,J,L+1)+(ri_crit2-RI(I,J,L))/(RI(I,J,L-1)-RI(I,J,L))*(Z(I,J,L)-Z(I,J,L+1)) - end if - end do - end do - end do - - where ( ZPBLRI2 .eq. MAPL_UNDEF ) ZPBLRI2 = Z(:,:,LM) - ZPBLRI2 = MIN(ZPBLRI2,Z(:,:,KPBLMIN)) - where ( ZPBLRI2 < 0.0 ) ZPBLRI2 = Z(:,:,LM) - end if ! ZPBLRI2 - - ! thetav gradient based pbl height diagnostic - if (associated(ZPBLTHV)) then - ZPBLTHV = MAPL_UNDEF - - do I = 1, IM - do J = 1, JM - - do L=LM,1,-1 - thetav(L) = TH(I,J,L)*(1.0+MAPL_VIREPS*Q(I,J,L)/(1.0-Q(I,J,L))) - end do - - maxdthvdz = 0 - - do L=LM-1,1,-1 - if(Z(I,J,L)<=Z(I,J,KPBLMIN)) then - dthvdz = (thetav(L+1)-thetav(L))/(Z(I,J,L+1)-Z(I,J,L)) - if(dthvdz>maxdthvdz) then - maxdthvdz = dthvdz - ZPBLTHV(I,J) = 0.5*(Z(I,J,L+1)+Z(I,J,L)) - end if - end if - end do - - end do - end do - end if ! ZPBLTHV - -!========================================================================= -! ZPBL defined by minimum in vertical gradient of refractivity. -! As shown in Ao, et al, 2012: "Planetary boundary layer heights from -! GPS radio occultation refractivity and humidity profiles", Climate and -! Dynamics. https://doi.org/10.1029/2012JD017598 -!========================================================================= - if (associated(ZPBLRFRCT)) then - - a1 = 0.776 ! K/Pa - a2 = 3.73e3 ! K2/Pa - - WVP = Q * PLO / (Q*(1.-0.622)+0.622) ! water vapor partial pressure - - ! Pressure gradient term - dum3d(:,:,2:LM-1) = (PLO(:,:,1:LM-2)-PLO(:,:,3:LM)) / (Z(:,:,1:LM-2)-Z(:,:,3:LM)) - dum3d(:,:,1) = (PLO(:,:,1)-PLO(:,:,2)) / (Z(:,:,1)-Z(:,:,2)) - dum3d(:,:,LM) = (PLO(:,:,LM-1)-PLO(:,:,LM)) / (Z(:,:,LM-1)-Z(:,:,LM)) - tmp3d = a1 * dum3d / T - - ! Add Temperature gradient term - dum3d(:,:,2:LM-1) = (T(:,:,1:LM-2)-T(:,:,3:LM)) / (Z(:,:,1:LM-2)-Z(:,:,3:LM)) - dum3d(:,:,1) = (T(:,:,1)-T(:,:,2)) / (Z(:,:,1)-Z(:,:,2)) - dum3d(:,:,LM) = (T(:,:,LM-1)-T(:,:,LM)) / (Z(:,:,LM-1)-Z(:,:,LM)) - tmp3d = tmp3d - (a1*plo/T**2 + 2.*a2*WVP/T**3)*dum3d - - ! Add vapor pressure gradient term - dum3d(:,:,2:LM-1) = (WVP(:,:,1:LM-2)-WVP(:,:,3:LM)) / (Z(:,:,1:LM-2)-Z(:,:,3:LM)) - dum3d(:,:,1) = (WVP(:,:,1)-WVP(:,:,2)) / (Z(:,:,1)-Z(:,:,2)) - dum3d(:,:,LM) = (WVP(:,:,LM-1)-WVP(:,:,LM)) / (Z(:,:,LM-1)-Z(:,:,LM)) - tmp3d = tmp3d + (a2/T**2)*dum3d - - ! ZPBL is height of minimum in refractivity (tmp3d) - do I = 1,IM - do J = 1,JM - K = MINLOC(tmp3d(I,J,:),DIM=1,BACK=.TRUE.) ! return last index, if multiple - ZPBLRFRCT(I,J) = Z(I,J,K) - end do - end do - - end if ! ZPBLRFRCT - - - ! PBL height diagnostic based on specific humidity gradient ! PBLH defined as level with minimum QV gradient - if (associated(ZPBLQV)) then - ZPBLQV = MAPL_UNDEF - - do I = 1, IM - do J = 1, JM - - maxdthvdz = 0 ! re-using variables from ZPBLTHV calc above - - do L=LM-1,1,-1 - if(Z(I,J,L)<=Z(I,J,KPBLMIN)) then - dthvdz = -1.*(Q(I,J,L+1)-Q(I,J,L))/(Z(I,J,L+1)-Z(I,J,L)) - if(dthvdz>maxdthvdz) then - maxdthvdz = dthvdz - ZPBLQV(I,J) = 0.5*(Z(I,J,L+1)+Z(I,J,L)) - end if - end if - end do - - end do - end do - end if ! ZPBLQV + if (associated(ZPBLQV)) call find_qv_pblh(im,jm,lm,kpblmin,z,q,zpblqv) + ! Surface-based inversion height, frequency + if (associated(SBITOP)) call return_surface_inversion_stats(im,jm,lm,z,t,sbitop,sbifrq) - if (associated(SBITOP) .or. associated(SBIFRQ) ) then + ! Trade inversion base height, temperature jump, frequency + if (associated(TRINVBS)) call return_trade_inversion_stats(im,jm,lm,plo,t,trinvbs,trinvdelt,trinvfrq) - SBIFRQ = 0. - SBITOP = MAPL_UNDEF + ! Wind turning angle from PBL top to surface + if (associated(TURNANG)) call find_turning_angle(im,jm,lm,u,v,z,tczpbl,kpbltc,turnang) - do I = 1, IM - do J = 1, JM - if (T(I,J,LM-1).gt.T(I,J,LM)) then - SBIFRQ(I,J) = 1. - do L = LM-1,1,-1 - if (T(I,J,L).gt.T(I,J,L+1)) then - SBITOP(I,J) = Z(I,J,L) - else - exit - end if - end do - end if - end do - end do - - end if ! SBITOP, SBIFRQ - - ! Trade inversion base height - if (associated(TRINVBS)) then - TRINVBS = MAPL_UNDEF - TRINVDELT = MAPL_UNDEF - TRINVFRQ = 0. - do I = 1,IM - do J = 1,JM - K = LM - - do while (PLO(I,J,K).gt.95000.) - K = K-1 - end do - do L = K,1,-1 ! K is first level above 950mb - if (PLO(I,J,L).lt.60000.) exit - - if (T(I,J,L-1).ge.T(I,J,L)) then ! if next level is warmer... - LTOP = L ! L is index of minimum T so far - do while (T(I,J,LTOP).ge.T(I,J,L)) ! find depth of warm layer - LTOP = LTOP-1 - end do - LTOP = LTOP+1 ! LTOP is index of highest level inside warm layer - - if ( MAXVAL(T(I,J,LTOP:L))-T(I,J,L).ge.0.5 .or. & - (MAXVAL(T(I,J,LTOP:L))-T(I,J,L).gt.0.01 .and. PLO(I,J,L)-PLO(I,J,LTOP)>2500.) ) then - - ! only save if DELTA-T exceeds any previous inversion - if ( TRINVFRQ(I,J).eq.0. .or. & - (TRINVFRQ(I,J).ne.0. .and. MAXVAL(T(I,J,LTOP:L))-T(I,J,L).gt.TRINVDELT(I,J)) ) then - TRINVBS(I,J) = PLO(I,J,L) - TRINVDELT(I,J) = MAXVAL(T(I,J,LTOP:L))-T(I,J,L) - TRINVFRQ(I,J) = 1. - end if - - end if - end if ! next level warmer - - end do ! L vert loop - - end do - end do - end if + + ZPBL = MAPL_UNDEF + if (associated(PPBL)) PPBL = MAPL_UNDEF SELECT CASE(PBLHT_OPTION) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/pblh_library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/pblh_library.F90 new file mode 100644 index 000000000..8ddcdc89e --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSturbulence_GridComp/pblh_library.F90 @@ -0,0 +1,446 @@ +module GEOSturbulence_PBLH_Library + + use ESMF + use MAPL + use MAPL_ConstantsMod, only: mapl_grav, mapl_cp, mapl_karman, mapl_rdry, & + mapl_pi, mapl_p00, mapl_pi + use MAPL_SatVaporMod, only: MAPL_EQsat + + implicit none + private + + character(len=ESMF_MAXSTR) :: IAm="GEOSturbulence_PBLH_Library" + + public :: find_bulk_ri_pblh + public :: find_kh2_pblh + public :: find_kh10p_pblh + public :: find_ri_pblh + public :: find_thv_pblh + public :: find_rfrct_pblh + public :: find_qv_pblh + public :: return_surface_inversion_stats + public :: return_trade_inversion_stats + public :: find_turning_angle + +contains + + !---------------------------------------------------------------------------- + ! Bulk-Richardson-Number PBLH (aka Transcom), described by Seidel et al 2012. + ! Scans from surface to find level at which bulk Ri first exceeds a critical + ! value specified by bulkri_crit parameter. + !---------------------------------------------------------------------------- + subroutine find_bulk_ri_pblh(im,jm,lm,u,v,z,thv,zpbl_bulkri,kpbl_bulkri,bulkri_crit) + + integer, intent(in) :: im, jm, lm + real , intent(in) :: bulkri_crit + real , dimension(im,jm,lm), intent(in) :: u, v, z, thv + real , dimension(im,jm) , intent(out) :: zpbl_bulkri + real , dimension(im,jm) , intent(out) :: kpbl_bulkri + + ! Local vars + integer :: i, j, k + real :: uv2h + real, dimension(lm) :: tcrib + + zpbl_bulkri = MAPL_UNDEF + tcrib(lm) = 0.0 + do i = 1, im + do j = 1, jm + do k=lm-1,1,-1 + uv2h = max(u(i,j,k)**2+v(i,j,k)**2,1.0E-8) + tcrib(k) = MAPL_GRAV * (thv(i,j,k)-thv(i,j,LM)) * z(i,j,k) / (thv(i,j,LM)*uv2h) + if (tcrib(k) >= bulkri_crit) then + zpbl_bulkri(i,j) = z(i,j,k+1)+(bulkri_crit-tcrib(k+1))/(tcrib(k)-tcrib(k+1))*(z(i,j,k)-z(i,j,k+1)) + kpbl_bulkri(i,j) = float(k) + exit + end if + end do + end do + end do + where (zpbl_bulkri<0.) + zpbl_bulkri = z(:,:,lm) + kpbl_bulkri = float(lm) + end where + + end subroutine find_bulk_ri_pblh + + + !---------------------------------------------------------------------------- + ! Diffusivity-based PBLH with threshold of 2 m2/s2 + ! Scans from surface to find level at which KH drops below threshold. + !---------------------------------------------------------------------------- + subroutine find_kh2_pblh(im,jm,lm,kpblmin,z,kh,zpbl_kh2,kpbl_kh2) + + integer, intent(in) :: im, jm, lm, kpblmin + real , dimension(im,jm,lm) , intent(in) :: z ! surface-relative heights + real , dimension(im,jm,0:lm), intent(in) :: kh + real , dimension(im,jm) , intent(out) :: zpbl_kh2 + real , dimension(im,jm) , intent(out) :: kpbl_kh2 + + ! Local vars + integer :: i, j, k + + zpbl_kh2 = MAPL_UNDEF + + do i = 1,im + do j = 1,jm + do k=lm,2,-1 + if ((kh(i,j,k-1) < 2.).and.(kh(i,j,k) >= 2.).and.(zpbl_kh2(i,j)==MAPL_UNDEF)) then + zpbl_kh2(i,j) = z(i,j,k) + kpbl_kh2(i,j) = float(k) + end if + end do + end do + end do + + where ( zpbl_kh2 .eq. MAPL_UNDEF ) + zpbl_kh2 = z(:,:,lm) + kpbl_kh2 = float(lm) + end where + zpbl_kh2 = min(zpbl_kh2,z(:,:,kpblmin)) + + end subroutine find_kh2_pblh + + + !---------------------------------------------------------------------------- + ! Diffusivity-based PBLH with threshold of 10 percent of maximum + ! Scans from surface to find level at which KH drops below threshold. + !---------------------------------------------------------------------------- + subroutine find_kh10p_pblh(im,jm,lm,kpblmin,z,ze,kh,zpbl_kh10p,kpbl_kh10p) + + integer, intent(in) :: im, jm, lm, kpblmin + real , dimension(im,jm,lm) , intent(in) :: z ! surface-relative heights + real , dimension(im,jm,0:lm), intent(in) :: ze ! surface-relative edge heights + real , dimension(im,jm,0:lm), intent(in) :: kh + real , dimension(im,jm) , intent(out) :: zpbl_kh10p + real , dimension(im,jm) , intent(out) :: kpbl_kh10p + + ! Local vars + real, dimension(lm+1) :: temparray + integer :: i, j, k, locmax + real :: minlval, maxkh + + zpbl_kh10p = MAPL_UNDEF + + do i = 1,im + do j = 1,jm + + temparray(1:lm+1) = kh(I,J,0:lm) + do k = lm,2,-1 + locmax = maxloc(temparray,1) + minlval = max(0.001,0.0001*maxval(temparray)) + if(temparray(locmax-1)maxkh) maxkh = temparray(k) + if(temparray(k-1)= 0.1*maxkh) & + .and. (zpbl_kh10p(i,j) == MAPL_UNDEF ) ) then + zpbl_kh10p(i,j) = ze(i,j,k)+ & + ((ze(i,j,k-1)-ze(i,j,k))/(temparray(k)-temparray(k+1))) * (0.1*maxkh-temparray(k+1)) + kpbl_kh10p(i,j) = float(k) + end if + end do + if ( zpbl_kh10p(i,j) .eq. MAPL_UNDEF .or. (maxkh.lt.1.)) then + zpbl_kh10p(i,j) = z(i,j,lm) + kpbl_kh10p(i,j) = float(lm) + endif + end do + end do + + zpbl_kh10p = min(zpbl_kh10p,ze(:,:,kpblmin)) + + end subroutine find_kh10p_pblh + + + !---------------------------------------------------------------------------- + ! Local Richardson number PBLH with a threshold of ri_crit. + ! Scans from surface to find level at which RI drops below threshold. + !---------------------------------------------------------------------------- + subroutine find_ri_pblh(im,jm,lm,kpblmin,z,ri,zpbl_ri,ri_crit) + + integer, intent(in) :: im, jm, lm, kpblmin + real , intent(in) :: ri_crit + real , dimension(im,jm,lm) , intent(in) :: z ! surface-relative heights + real , dimension(im,jm,0:lm), intent(in) :: ri ! local Richardson number + real , dimension(im,jm) , intent(out) :: zpbl_ri + + ! local variables + integer :: i, j, k + + zpbl_ri = MAPL_UNDEF + where (ri(:,:,lm-1)>ri_crit) zpbl_ri = z(:,:,lm) + + do i = 1, im + do j = 1, jm + do k=lm-1,1,-1 + if( (ri(i,j,k-1)>ri_crit) .and. (zpbl_ri(i,j) == MAPL_UNDEF) ) then + zpbl_ri(i,j) = z(i,j,k+1)+(ri_crit-ri(i,j,k))/(ri(i,j,k-1)-ri(i,j,k))*(z(i,j,k)-z(i,j,k+1)) + end if + end do + end do + end do + + where ( zpbl_ri .eq. MAPL_UNDEF ) zpbl_ri = z(:,:,lm) + zpbl_ri = min(zpbl_ri,z(:,:,kpblmin)) + where ( zpbl_ri < 0.0 ) zpbl_ri = z(:,:,lm) + + end subroutine find_ri_pblh + + + !----------------------------------------------------------------------------- + ! Potential temperature gradient PBLH. + ! Scans from surface to find height of maximum virtual potential temperature + ! gradient. + !----------------------------------------------------------------------------- + subroutine find_thv_pblh(im,jm,lm,kpblmin,z,thv,zpbl_thv) + + integer, intent(in) :: im, jm, lm, kpblmin + real , dimension(im,jm,lm), intent(in) :: z ! surface-relative heights + real , dimension(im,jm,lm), intent(in) :: thv ! virtual potential temperature + real , dimension(im,jm) , intent(out) :: zpbl_thv + + ! local variables + integer :: i, j, k + real :: dthvdz, maxdthvdz + + zpbl_thv = MAPL_UNDEF + + do i = 1, im + do j = 1, jm + + maxdthvdz = 0 + + do k=lm-1,1,-1 + if(z(i,j,k)<=z(i,j,kpblmin)) then + dthvdz = (thv(i,j,k+1)-thv(i,j,k))/(z(i,j,k+1)-z(i,j,k)) + if(dthvdz>maxdthvdz) then + maxdthvdz = dthvdz + zpbl_thv(i,j) = 0.5*(z(i,j,k+1)+z(i,j,k)) + end if + end if + end do + + end do + end do + + end subroutine find_thv_pblh + + + !----------------------------------------------------------------------------- + ! Specific humidity gradient PBLH. + ! Scans from surface to find height of maximum gradient in specific humidity. + !----------------------------------------------------------------------------- + subroutine find_qv_pblh(im,jm,lm,kpblmin,z,qv,zpbl_qv) + + integer, intent(in) :: im, jm, lm, kpblmin + real , dimension(im,jm,lm), intent(in) :: z ! surface-relative heights + real , dimension(im,jm,lm), intent(in) :: qv ! specific humidity + real , dimension(im,jm) , intent(out) :: zpbl_qv + + ! local variables + integer :: i, j, k + real :: dqvdz, maxdqvdz + + zpbl_qv = MAPL_UNDEF + + do i = 1, im + do j = 1, jm + + maxdqvdz = 0. + do k=lm-1,1,-1 + if(z(i,j,k)<=z(i,j,kpblmin)) then + dqvdz = -1.*(qv(i,j,k+1)-qv(i,j,k))/(z(i,j,k+1)-z(i,j,k)) + if(dqvdz>maxdqvdz) then + maxdqvdz = dqvdz + zpbl_qv(i,j) = 0.5*(z(i,j,k+1)+z(i,j,k)) + end if + end if + end do + + end do + end do + + end subroutine find_qv_pblh + + + !========================================================================= + ! ZPBL defined by minimum in vertical gradient of refractivity. + ! As shown in Ao, et al, 2012: "Planetary boundary layer heights from + ! GPS radio occultation refractivity and humidity profiles", Climate and + ! Dynamics. https://doi.org/10.1029/2012JD017598 + !========================================================================= + subroutine find_rfrct_pblh(im,jm,lm,z,p,t,qv,zpbl_rfrct) + + integer, intent(in) :: im, jm, lm + real , dimension(im,jm,lm), intent(in) :: z, & ! surface-relative heights + p, & + t, & + qv + real , dimension(im,jm) , intent(out) :: zpbl_rfrct + + ! local variables + integer :: i, j, k + real :: a1, a2 + real, dimension(im,jm,lm) :: wvp, dum3d, tmp3d + + a1 = 0.776 ! K/Pa + a2 = 3.73e3 ! K2/Pa + + wvp = qv * p / (qv*(1.-0.622)+0.622) ! water vapor partial pressure + + ! pressure gradient term + dum3d(:,:,2:lm-1) = (p(:,:,1:lm-2)-p(:,:,3:lm)) / (z(:,:,1:lm-2)-z(:,:,3:lm)) + dum3d(:,:,1) = (p(:,:,1)-p(:,:,2)) / (z(:,:,1)-z(:,:,2)) + dum3d(:,:,lm) = (p(:,:,lm-1)-p(:,:,lm)) / (z(:,:,lm-1)-z(:,:,lm)) + tmp3d = a1 * dum3d / t + + ! add temperature gradient term + dum3d(:,:,2:lm-1) = (t(:,:,1:lm-2)-t(:,:,3:lm)) / (z(:,:,1:lm-2)-z(:,:,3:lm)) + dum3d(:,:,1) = (t(:,:,1)-t(:,:,2)) / (z(:,:,1)-z(:,:,2)) + dum3d(:,:,lm) = (t(:,:,lm-1)-t(:,:,lm)) / (z(:,:,lm-1)-z(:,:,lm)) + tmp3d = tmp3d - (a1*p/t**2 + 2.*a2*wvp/t**3)*dum3d + + ! add vapor pressure gradient term + dum3d(:,:,2:lm-1) = (wvp(:,:,1:lm-2)-wvp(:,:,3:lm)) / (z(:,:,1:lm-2)-z(:,:,3:lm)) + dum3d(:,:,1) = (wvp(:,:,1)-wvp(:,:,2)) / (z(:,:,1)-z(:,:,2)) + dum3d(:,:,lm) = (wvp(:,:,lm-1)-wvp(:,:,lm)) / (z(:,:,lm-1)-z(:,:,lm)) + tmp3d = tmp3d + (a2/t**2)*dum3d + + ! zpbl is height of minimum in refractivity (tmp3d) + do i = 1,im + do j = 1,jm + k = minloc(tmp3d(i,j,:),dim=1,back=.true.) ! return last index, if multiple + zpbl_rfrct(i,j) = z(i,j,k) + end do + end do + + end subroutine find_rfrct_pblh + + + !----------------------------------------------------------------------------- + ! Surface-based inversion (top) height and frequency + ! Scans up from surface while absolute temperature increases with height. + !----------------------------------------------------------------------------- + subroutine return_surface_inversion_stats(im,jm,lm,z,t,sbitop,sbifrq) + + integer, intent(in) :: im, jm, lm + real , dimension(im,jm,lm), intent(in) :: z ! surface-relative heights + real , dimension(im,jm,lm), intent(in) :: t ! temperature + real , dimension(im,jm) , intent(out) :: sbitop ! sbi top height + real , dimension(:,:) , pointer :: sbifrq ! sbi frequency + + ! local variables + integer :: i, j, k + + sbitop = MAPL_UNDEF + + do i = 1, im + do j = 1, jm + if (t(i,j,lm-1).gt.t(i,j,lm)) then + do k = lm-1,1,-1 + if (t(i,j,k).gt.t(i,j,k+1)) then + sbitop(i,j) = z(i,j,k) + else + exit + end if + end do + end if + end do + end do + if (associated(sbifrq)) then + sbifrq = 0. + where(sbitop.ne.MAPL_UNDEF) sbifrq=1. + end if + + end subroutine return_surface_inversion_stats + + + !----------------------------------------------------------------------------- + ! Trade inversion height, magnitude and frequency + ! Searches upward from 950mb to 600mb for either an absolute temperature + ! increase of at least 0.5 K, or a segment of at least 25mb with no + ! temperature decrease. + !----------------------------------------------------------------------------- + subroutine return_trade_inversion_stats(im,jm,lm,p,t,trinvbs,trinvdelt,trinvfrq) + + integer, intent(in) :: im, jm, lm + real , dimension(im,jm,lm), intent(in) :: p ! pressure + real , dimension(im,jm,lm), intent(in) :: t ! temperature + real , dimension(im,jm) , intent(out) :: trinvbs ! inversion base height + real , dimension(:,:) , pointer :: trinvdelt ! inversion temperature jump + real , dimension(:,:) , pointer :: trinvfrq ! inversion frequency + + ! local variables + integer :: i, j, k, ktop + real, dimension(im,jm) :: tmpdelt, tmpfrq + + trinvbs = MAPL_UNDEF + tmpdelt = MAPL_UNDEF + tmpfrq = 0. + do i = 1,im + do j = 1,jm + k = lm + + do while (p(i,j,k).gt.95000.) + k = k-1 + end do + do k = k,1,-1 ! k is first level above 950mb + if (p(i,j,k).lt.60000.) exit + if (t(i,j,k-1).ge.t(i,j,k)) then ! if next level is warmer... + ktop = k ! k is index of minimum t so far + do while (t(i,j,ktop).ge.t(i,j,k)) ! find depth of warm layer + ktop = ktop-1 + end do + ktop = ktop+1 ! ktop is index of highest level inside warm layer + + if ( maxval(t(i,j,ktop:k))-t(i,j,k).ge.0.5 .or. & + (maxval(t(i,j,ktop:k))-t(i,j,k).gt.0.01 .and. p(i,j,k)-p(i,j,ktop)>2500.) ) then + + ! only save if delta-t exceeds any previous inversion + if ( tmpfrq(i,j).eq.0. .or. & + (tmpfrq(i,j).ne.0. .and. maxval(t(i,j,ktop:k))-t(i,j,k).gt.tmpdelt(i,j)) ) then + trinvbs(i,j) = p(i,j,k) + tmpdelt(i,j) = maxval(t(i,j,ktop:k))-t(i,j,k) + tmpfrq(i,j) = 1. + end if + + end if + end if ! next level warmer + + end do ! k vert loop + end do ! j + end do ! i + if (associated(trinvfrq)) trinvfrq = tmpfrq + if (associated(trinvdelt)) trinvdelt = tmpdelt + + end subroutine return_trade_inversion_stats + + subroutine find_turning_angle(im,jm,lm,u,v,z,zpbl,kpbl,turnangle) + + integer, intent(in) :: im, jm, lm + real , dimension(im,jm,lm), intent(in) :: u, v, z + real , dimension(im,jm) , intent(in) :: zpbl, kpbl + real , dimension(im,jm) , intent(out) :: turnangle + + ! Local vars + integer :: i, j, k + real, dimension(im,jm) :: utop, vtop + + do i = 1,im + do j = 1,jm + k = nint(kpbl(i,j)) + utop(i,j) = u(i,j,k) + (zpbl(i,j)-z(i,j,k)) * (u(i,j,k)-u(i,j,k+1))/(z(i,j,k)-z(i,j,k+1)) + vtop(i,j) = v(i,j,k) + (zpbl(i,j)-z(i,j,k)) * (v(i,j,k)-v(i,j,k+1))/(z(i,j,k)-z(i,j,k+1)) + end do + end do + turnangle = (180./mapl_pi)*atan2( utop*v(:,:,lm)-vtop*u(:,:,lm), utop*u(:,:,lm) + vtop*v(:,:,lm) ) + + end subroutine find_turning_angle + +end module GEOSturbulence_PBLH_Library + +