diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 1f778472..955f9b1f 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1026,8 +1026,15 @@ #endif + + + + + + + diff --git a/src/core_atmosphere/diagnostics/Makefile b/src/core_atmosphere/diagnostics/Makefile index 8925e1f2..29d84923 100644 --- a/src/core_atmosphere/diagnostics/Makefile +++ b/src/core_atmosphere/diagnostics/Makefile @@ -11,12 +11,13 @@ DIAGNOSTIC_MODULES = \ mpas_pv_diagnostics.o \ mpas_soundings.o \ mpas_vertical_interpolation.o \ + mpas_stability_indices.o mpas_isobaric_diagnostics.o: mpas_atm_diagnostics_utils.o mpas_vertical_interpolation.o mpas_cloud_diagnostics.o: mpas_atm_diagnostics_utils.o -mpas_convective_diagnostics.o: mpas_atm_diagnostics_utils.o +mpas_convective_diagnostics.o: mpas_atm_diagnostics_utils.o mpas_stability_indices.o mpas_pv_diagnostics.o: mpas_atm_diagnostics_utils.o @@ -24,6 +25,8 @@ mpas_soundings.o: mpas_vertical_interpolation.o: +mpas_stability_indices.o: + ################### Generally no need to modify below here ################### diff --git a/src/core_atmosphere/diagnostics/Registry_convective.xml b/src/core_atmosphere/diagnostics/Registry_convective.xml index 59ad1c2f..039c6abc 100644 --- a/src/core_atmosphere/diagnostics/Registry_convective.xml +++ b/src/core_atmosphere/diagnostics/Registry_convective.xml @@ -4,11 +4,34 @@ + + + + + + + + + + + description="Most Unstable Convective available potential energy"/> + description="Most Unstable Convective inhibition"/> + + + + diff --git a/src/core_atmosphere/diagnostics/Registry_isobaric.xml b/src/core_atmosphere/diagnostics/Registry_isobaric.xml index 0012f9c0..ff98a340 100644 --- a/src/core_atmosphere/diagnostics/Registry_isobaric.xml +++ b/src/core_atmosphere/diagnostics/Registry_isobaric.xml @@ -52,6 +52,9 @@ + + diff --git a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F index 163ee377..c23f2673 100644 --- a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F @@ -10,6 +10,7 @@ module mpas_convective_diagnostics use mpas_derived_types, only : MPAS_pool_type, MPAS_clock_type, MPAS_LOG_ERR, MPAS_LOG_CRIT use mpas_kind_types, only : RKIND, R8KIND use mpas_log, only : mpas_log_write + use stability_indices type (MPAS_pool_type), pointer :: mesh type (MPAS_pool_type), pointer :: state @@ -236,8 +237,12 @@ subroutine convective_diagnostics_compute() integer, pointer :: nCells, nCellsSolve, nVertLevels, nVertLevelsP1 ! Fields that are computed in this routine + real (kind=RKIND), dimension(:), pointer :: sfccape + real (kind=RKIND), dimension(:), pointer :: sfccin real (kind=RKIND), dimension(:), pointer :: cape real (kind=RKIND), dimension(:), pointer :: cin + real (kind=RKIND), dimension(:), pointer :: mlcape + real (kind=RKIND), dimension(:), pointer :: mlcin real (kind=RKIND), dimension(:), pointer :: lcl real (kind=RKIND), dimension(:), pointer :: lfc real (kind=RKIND), dimension(:), pointer :: srh_0_1km @@ -251,6 +256,10 @@ subroutine convective_diagnostics_compute() real (kind=RKIND), dimension(:), pointer :: temperature_surface real (kind=RKIND), dimension(:), pointer :: dewpoint_surface + real (kind=RKIND), dimension(:), pointer :: rli + real (kind=RKIND), dimension(:), pointer :: rki + real (kind=RKIND), dimension(:), pointer :: rwi + ! Other fields used in the computation of convective diagnostics ! defined above real (kind=RKIND), dimension(:,:), pointer :: height @@ -268,6 +277,9 @@ subroutine convective_diagnostics_compute() real (kind=RKIND) :: u_storm, v_storm, u_srh_bot, v_srh_bot, u_srh_top, v_srh_top real (kind=RKIND) :: u_mean, v_mean, u_shear, v_shear, shear_magnitude real (kind=RKIND) :: b_term, cape_out, cin_out + + real (kind=RKIND) :: rli_out, rki_out, rwi_out + real (kind=RKIND), dimension(:), allocatable :: dudz, dvdz, zp real (kind=RKIND), dimension(:), allocatable :: zrel, srh real (kind=RKIND), dimension(:), allocatable :: p_in, t_in, td_in @@ -276,16 +288,34 @@ subroutine convective_diagnostics_compute() real (kind=RKIND) :: evp - logical :: need_cape, need_cin, need_lcl, need_lfc, need_srh_01km, need_srh_03km, need_uzonal_sfc, need_uzonal_1km, & + logical :: need_sfccape, need_sfccin + logical :: need_cape, need_cin, need_mlcape, need_mlcin, need_lcl, need_lfc, need_srh_01km, need_srh_03km, need_uzonal_sfc, need_uzonal_1km, & need_uzonal_6km, need_umerid_sfc, need_umerid_1km, need_umerid_6km, need_tsfc, need_tdsfc logical :: need_any_diags + logical :: need_rli, need_rki, need_rwi need_any_diags = .false. + need_rli = MPAS_field_will_be_written('rli') + need_any_diags = need_any_diags .or. need_rli + need_rki = MPAS_field_will_be_written('rki') + need_any_diags = need_any_diags .or. need_rki + need_rwi = MPAS_field_will_be_written('rwi') + need_any_diags = need_any_diags .or. need_rwi + + need_sfccape = MPAS_field_will_be_written('sfccape') + need_any_diags = need_any_diags .or. need_sfccape + need_sfccin = MPAS_field_will_be_written('sfccin') + need_any_diags = need_any_diags .or. need_sfccin + need_cape = MPAS_field_will_be_written('cape') need_any_diags = need_any_diags .or. need_cape need_cin = MPAS_field_will_be_written('cin') need_any_diags = need_any_diags .or. need_cin + need_mlcape = MPAS_field_will_be_written('mlcape') + need_any_diags = need_any_diags .or. need_mlcape + need_mlcin = MPAS_field_will_be_written('mlcin') + need_any_diags = need_any_diags .or. need_mlcin need_lcl = MPAS_field_will_be_written('lcl') need_any_diags = need_any_diags .or. need_lcl need_lfc = MPAS_field_will_be_written('lfc') @@ -318,8 +348,13 @@ subroutine convective_diagnostics_compute() call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(mesh, 'nVertLevelsP1', nVertLevelsP1) + call mpas_pool_get_array(diag, 'sfccape', sfccape) + call mpas_pool_get_array(diag, 'sfccin', sfccin) + call mpas_pool_get_array(diag, 'cape', cape) call mpas_pool_get_array(diag, 'cin', cin) + call mpas_pool_get_array(diag, 'mlcape', mlcape) + call mpas_pool_get_array(diag, 'mlcin', mlcin) call mpas_pool_get_array(diag, 'lcl', lcl) call mpas_pool_get_array(diag, 'lfc', lfc) call mpas_pool_get_array(diag, 'srh_0_1km', srh_0_1km) @@ -333,6 +368,10 @@ subroutine convective_diagnostics_compute() call mpas_pool_get_array(diag, 'temperature_surface', temperature_surface) call mpas_pool_get_array(diag, 'dewpoint_surface', dewpoint_surface) + call mpas_pool_get_array(diag, 'rli', rli) + call mpas_pool_get_array(diag, 'rki', rki) + call mpas_pool_get_array(diag, 'rwi', rwi) + call mpas_pool_get_array(mesh, 'zgrid', height) call mpas_pool_get_array(diag, 'uReconstructMeridional', umeridional) call mpas_pool_get_array(diag, 'uReconstructZonal', uzonal) @@ -451,6 +490,21 @@ subroutine convective_diagnostics_compute() end do ! calculate cape and cin + + if (need_sfccape .or. need_sfccin) then + do iCell=1, nCellsSolve + p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND + t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND + td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND + + call getcape( 1, nVertLevels, p_in, t_in, td_in, cape_out, cin_out ) + + sfccape(iCell) = cape_out + sfccin(iCell) = cin_out + + end do + end if + if (need_cape .or. need_cin) then do iCell=1, nCellsSolve p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND @@ -463,7 +517,7 @@ subroutine convective_diagnostics_compute() ! /(17.625-log(relhum(k,iCell))-((17.625*t_in(k))/(243.04+t_in(k)))) ! end do - call getcape( nVertLevels, p_in, t_in, td_in, cape_out, cin_out ) + call getcape( 2, nVertLevels, p_in, t_in, td_in, cape_out, cin_out ) cape(iCell) = cape_out cin(iCell) = cin_out @@ -471,6 +525,40 @@ subroutine convective_diagnostics_compute() end do end if + if (need_mlcape .or. need_mlcin) then + do iCell=1, nCellsSolve + p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND + t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND + td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND + + call getcape( 3, nVertLevels, p_in, t_in, td_in, cape_out, cin_out ) + + mlcape(iCell) = cape_out + mlcin(iCell) = cin_out + + end do + end if + +! Stability indices + + if (need_rli .or. need_rki .or. need_rwi ) then + do iCell=1, nCellsSolve + p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) / 100.0_RKIND + t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) - 273.15_RKIND + td_in(1:nVertLevels) = dewpoint(1:nVertLevels,iCell) - 273.15_RKIND + +! print*,'Pressure (1:nVertLevels) ',p_in(1),p_in(nVertLevels) +! Enter variables from surface to top + + call INDX1 (nVertLevels, p_in, t_in, td_in, rli_out, rki_out, rwi_out ) + + rli(iCell) = rli_out + rki(iCell) = rki_out + rwi(iCell) = rwi_out + + end do + end if + deallocate(temperature) deallocate(dewpoint) @@ -592,10 +680,10 @@ end function integral_zpoint !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !----------------------------------------------------------------------- - subroutine getcape( nk , p_in , t_in , td_in, cape , cin ) + subroutine getcape( source, nk , p_in , t_in , td_in, cape , cin ) implicit none - integer, intent(in) :: nk + integer, intent(in) :: source, nk real (kind=RKIND), dimension(nk), intent(in) :: p_in,t_in,td_in real (kind=RKIND), intent(out) :: cape,cin @@ -639,7 +727,7 @@ subroutine getcape( nk , p_in , t_in , td_in, cape , cin ) ! results,larger number makes code ! go faster) - integer, parameter :: source = 2 ! Source parcel: +! integer, parameter :: source = 2 ! Source parcel: ! 1 = surface ! 2 = most unstable (max theta-e) ! 3 = mixed-layer (specify ml_depth) diff --git a/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F index 2701c25b..1d5ce0d3 100644 --- a/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F @@ -44,7 +44,7 @@ module mpas_isobaric_diagnostics need_kzm_isobaric, need_ni_isobaric, need_nr_isobaric, & need_qr_isobaric, need_qs_isobaric, need_qv_isobaric, & need_qc_isobaric, need_qg_isobaric, need_qi_isobaric, & - need_zgeo_isobaric + need_zgeo_isobaric, need_theta_isobaric contains @@ -152,6 +152,7 @@ subroutine isobaric_diagnostics_compute() need_zgeo_isobaric = .false. need_cldfrac_isobaric = .false. need_temperature_isobaric = .false. + need_theta_isobaric = .false. need_dewpoint_isobaric = .false. need_relhum_isobaric = .false. need_w_isobaric = .false. @@ -198,11 +199,12 @@ subroutine isobaric_diagnostics_compute() need_zgeo_isobaric = MPAS_field_will_be_written('zgeo_isobaric') need_any_diags = need_any_diags .or. need_zgeo_isobaric - need_cldfrac_isobaric = MPAS_field_will_be_written('cldfrac_isobaric') need_any_diags = need_any_diags .or. need_cldfrac_isobaric need_temperature_isobaric = MPAS_field_will_be_written('temperature_isobaric') need_any_diags = need_any_diags .or. need_temperature_isobaric + need_theta_isobaric = MPAS_field_will_be_written('theta_isobaric') + need_any_diags = need_any_diags .or. need_theta_isobaric need_dewpoint_isobaric = MPAS_field_will_be_written('dewpoint_isobaric') need_any_diags = need_any_diags .or. need_dewpoint_isobaric need_relhum_isobaric = MPAS_field_will_be_written('relhum_isobaric') @@ -325,6 +327,7 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph real (kind=RKIND), dimension(:,:), pointer :: cldfrac_isobaric real (kind=RKIND), dimension(:,:), pointer :: temperature_isobaric + real (kind=RKIND), dimension(:,:), pointer :: theta_isobaric real (kind=RKIND), dimension(:,:), pointer :: dewpoint_isobaric real (kind=RKIND), dimension(:,:), pointer :: relhum_isobaric real (kind=RKIND), dimension(:,:), pointer :: w_isobaric @@ -450,6 +453,7 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph call mpas_pool_get_array(diag, 'cldfrac_isobaric', cldfrac_isobaric) call mpas_pool_get_array(diag, 'temperature_isobaric', temperature_isobaric) + call mpas_pool_get_array(diag, 'theta_isobaric', theta_isobaric) call mpas_pool_get_array(diag, 'dewpoint_isobaric', dewpoint_isobaric) call mpas_pool_get_array(diag, 'relhum_isobaric', relhum_isobaric) call mpas_pool_get_array(diag, 'w_isobaric', w_isobaric) @@ -674,6 +678,20 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag, diag_physics, tend_ph ! call mpas_log_write('--- end interpolate temperature:') end if + if ( need_theta_isobaric ) then + + ! THETA = T * (1000./P1)**Rd_Cp + + do k = 1, nIsoLevelsT + do icell = 1, nCells + theta_isobaric(k,icell) = & + temperature_isobaric(k,icell)*(100000.0/t_iso_levels(k))**0.2857142 + enddo + enddo + + ! call mpas_log_write('--- end interpolate theta:') + end if + if (need_mp_thompson .eqv. .true. .or. need_bl_mynn .eqv. .true.) then if (NEED_NI_ISOBARIC) then diff --git a/src/core_atmosphere/diagnostics/mpas_stability_indices.F b/src/core_atmosphere/diagnostics/mpas_stability_indices.F new file mode 100644 index 00000000..603ab83d --- /dev/null +++ b/src/core_atmosphere/diagnostics/mpas_stability_indices.F @@ -0,0 +1,363 @@ +MODULE stability_indices + + use mpas_kind_types, only : RKIND + +! STABILITY ANALYSIS USING PARCEL METHOD + + IMPLICIT NONE + +! integer, parameter :: RKIND=8 + + PRIVATE + PUBLIC :: INDX1 + +CONTAINS + + SUBROUTINE INDX1 (JNO, P, TS, TSD, RLI, RKI, RWI) +!----------------------------------------------------------------------- +! +! This routine computes stability indices (Stone 1983, 1984) +! The following fields are computed by this routine: +! +! Lifted index +! K index +! Showalter index +! +! Stone, H.M., 1983: Stability Analysis Program, NOAA Eastern Region +! Computer Programs and Problems NWS ERCP-No. 9, National Weather +! ·Service, Garden City, NY. +! +! Stone, H.M., 1984: Stability Analysis Plot Program, NOAA Eastern Region +! Computer Programs and Problems NWS ERCP-No. 16, National Weather +! ·Service, Garden City, NY. +! +!----------------------------------------------------------------------- + +! COMPUTES LIFTED INDEX, K INDEX, AND SHOWALTER INDEX + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: JNO + INTEGER :: I,J,K + + REAL (kind=RKIND), DIMENSION(JNO), INTENT(IN) :: P,TS,TSD + REAL (kind=RKIND), DIMENSION(JNO) :: TW + REAL (kind=RKIND) :: PL(3),TL(3),TDL(3) + REAL (kind=RKIND), INTENT(OUT) :: RLI, RKI, RWI + + REAL (kind=RKIND) :: TC,TH,WTH,WTC,WSUM,THSUM,P1,P2,PFINISH + REAL (kind=RKIND) :: PLOG1,PLOG2,PLOG3,PLOG4 + REAL (kind=RKIND) :: DP1,E2,FACTOR,FACTORD,FACTORT,PPARCEL,TDE1,TDE2,TDPARCEL + REAL (kind=RKIND) :: TE1,TE2,TH1,TH2,THAVG,THW,TP,TP500,TPARCEL,W,W1,W2,WAVEG + REAL (kind=RKIND) :: WAVG,WFW,X + + PL(1) = 850.0_RKIND + PL(2) = 700.0_RKIND + PL(3) = 500.0_RKIND + +! COMPUTE WET BULB POTENTIAL TEMP AT ALL SGFNT LEVELS ABV SFC + + DO K = 1, JNO + TC=TCONOF(TS(K),TSD(K)) ! CONDENSATION TEMPERATURE + TH=THETA(TS(K)+273.16,1000.0_RKIND,P(K))-273.16_RKIND ! POT TEMP DEG C + WTH=WOBF(TH) + WTC=WOBF(TC) + TW(K)=TH-WTH+WTC ! WET BULB POT TEMP - DEG C + ENDDO + +! WRITE (*,65)" SIGNIFICANT LEVELS" +!65 FORMAT (20A) +! WRITE(*,70)" P T TD TW" +!70 FORMAT(20A) +! DO K = JNO,0,-1 +! WRITE(*,75) P(K),TS(K),TSD(K),TW(K) +! ENDDO +!75 FORMAT(2X,F6.1,4X,F6.1,4X,F6.1,4x,F6.1) + +! STOP + + DP1=50. ! AVERAGES OVER FIRST -DPI- MBS. + WSUM=0. + THSUM=0. + J=1 + P1=P(1) + TE1=TS(1) + TDE1=TSD(1) + TH1=THETA(TE1+273.16,1000.0_RKIND,P(1)) ! POT TEMP + W1=WMROF(P(1),TDE1) + PFINISH=P1-DP1 + 3 P2=PFINISH + IF(P(J+1)-P2) 1,2,2 ! <0, 0, >0 + 2 P2=P(J+1) + J=J+1 + 1 PLOG1=ALOG(P(J)/P(J+1)) + FACTORT=(TS(J)-TS(J+1))/PLOG1 + FACTORD=(TSD(J)-TSD(J+1))/PLOG1 + PLOG2=ALOG(P2/P(J+1)) + TE2=TS(J+1)+FACTORT*PLOG2 ! ENVIRONMENT TEMP AT P2 + TDE2=TSD(J+1)+FACTORD*PLOG2 ! ENVIRONMENT DEWPT AT P2 + TH2=THETA(TE2+273.16,1000.0_RKIND,P2) ! POT TEMP AT TE2,P2 + W2=WMROF(P2,TDE2) ! MIXING RATIO AT P2 + PLOG3=ALOG(P1/P2) + TH=.5*(TH1+TH2)*PLOG3 ! AVG POT TEMP IN LYR P1-P2 + W=.5*(W1+W2)*PLOG3 ! AVG MIX RATIO IN LYR P1-P2 + THSUM=THSUM+TH + WSUM=WSUM+W + P1=P2 + TH1=TH2 + W1=W2 + IF (P2.GT.PFINISH) GO TO 3 + +! COMPUTE AVG VALUES FOR FIRST -DPI- MBS. + + PLOG4=ALOG(P(1)/PFINISH) + THAVG=THSUM/PLOG4 + WAVG=WSUM/PLOG4 + PPARCEL=P(1)-.5*DP1 + TPARCEL=THETA(THAVG,PPARCEL,1000.0_RKIND)-273.16 ! DEG C + X=.0200*(TPARCEL-12.5+7500./PPARCEL) ! NON-IDEAL GAS CORRECTION + WFW=1.+.0000045*PPARCEL+.00140*X*X ! NON-IDEAL GAS CORRECTION + E2=.001*WAVG*PPARCEL/((WAVG*.001+.62197)*WFW) ! VAPOR PRES (MB) + TDPARCEL=DPTOF(E2) + TC=TCONOF(TPARCEL,TDPARCEL) + TH=THAVG-273.16 ! POT TEMP DEG C + WTH=WOBF(TH) + WTC=WOBF(TC) + THW=TH-WTH+WTC ! EQUIV WET BULB POT TEMP (DEG C) + TP500=SATLFT(THW,500.0_RKIND) + +! GET TEMP AND DEWPT AT 850,700,500 MBS + + DO 5 J=1,3 + DO 4 I=1,JNO + IF (PL(J)-P(I)) 4,6,7 +4 CONTINUE + PRINT*,"RAOB TERMINATES TOO SOON, P(JNO) = ",P(JNO) + RETURN +6 TL(J)=TS(I) + TDL(J)=TSD(I) + GO TO 5 +7 FACTOR=ALOG(PL(J)/P(I))/ALOG(P(I-1)/P(I)) + TL(J)=TS(I)+FACTOR*(TS(I-1)-TS(I)) + TDL(J)=TSD(I)+FACTOR*(TSD(I-1)-TSD(I)) +5 CONTINUE + +! COMPUTE LIFTED INDEX + + RLI=TL(3)-TP500 ! LIFTED INDEX + +! WRITE(*,*) +! WRITE(*,80) TL(3),TP500 + +! COMPUTE K INDEX + + RKI=(TL(1)-TL(3))+TDL(1)-(TL(2)-TDL(2)) ! K INDEX + +! COMPUTE SHOWALTER INDEX + + TC=TCONOF(TL(1),TDL(1)) + TH=THETA(TL(1)+273.16,1000.0_RKIND,850.0_RKIND)-273.16_RKIND ! DEG C + WTH=WOBF(TH) + WTC=WOBF(TC) + THW=TH-WTH+WTC ! EQUIV WET BULB POT TEMP + TP=SATLFT(THW,500.0_RKIND) + RWI=TL(3)-TP ! SHOWALTER INDEX + +! WRITE (*,*) +! WRITE (*,80) TL(3), TP +!80 FORMAT (10X,"TL - ",F6.2,4X,"TP - ",F6.2) +! WRITE (*,*)"Stability Indices:" +! WRITE (*,81) RLI,RKI,RWI +!81 FORMAT (10X,"LI - ",F4.0,4X,"KI - ",F4.0,4X,"SWI - ",F4.0) + + RETURN + END + +!-------------- + + REAL(kind=RKIND) FUNCTION WOBF(T) + + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: T + REAL (kind=RKIND) :: X, POL + +! COMPUTE BY DOUBLE ASYMPTOTIC APPROXIMATION +! CONSIDER SEPARATELY IF .GT. OR .LE. 20 DEG. +! CENT. FOR ALL TEMPS...THETU=THETA-UOBF(THETA)+UOBF(TEMPCON) +! CENT. FOR ALL TEMPS...THETM=THETA-U08F(THETA)+UOBF(TEMP) + X=T-20.0 + IF(X) 10,10,20 +10 CONTINUE +! CURVE FIG FOR COOL TEMPERATURE RANGE + POL=1.000+X*(-8.8416605E-3+X*(1.4714143E-4+X*(-9.6719890E-7 & + +X*(-3.2607217E-8+X*(-3.8598073E-10))))) + POL=POL*POL + WOBF=15.130/(POL*POL) + RETURN +20 CONTINUE +! CURVE FIT FOR WARMER TEMPERATURES + POL=1.000+X*(3.6182989E-3+X*(-1.3603273E-5+X*(4.9618922E-7 & + +X*(-6.1059365E-9+X*(3.9401551E-11+X*(-1.2588129E-13 & + +X*(1.6688280E-16) )))))) + POL=POL*POL + WOBF=29.930/(POL*POL)+0.9600*X-14.800 + RETURN + END + + REAL(kind=RKIND) FUNCTION SATLFT (THM,P) + +! COMPUTES TEMPERATURE (DEG C) WHERE THETA MOIST (DEG C) CROSSES P (MB) + + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: THM,P + +! CONSIDER THE EXPONENTIAL FOR POTENTIAL TEMPERATURE AS ROCP + REAL, PARAMETER:: ROCP=0.28571428 + REAL (kind=RKIND) :: PWRP,TONE,TTWO,EONE,ETWO,RATE,EOR + + IF(ABS(P-1000.0)-0.0010) 100,100,200 +100 SATLFT=THM + RETURN +200 PWRP=(P/1000.0)**ROCP +! COMPUTE TEMPERATURE OF DRY ADIABATIC LIFT FOR FIRST GUESS + TONE=(THM+273.16)*PWRP-273.16 +! CONSIDER PSEUDO-ADIABAT, EW1, THROUGH TONE AT P. +! COMPUTE EONE=EW1-THM + EONE=WOBF(TONE)-WOBF(THM) + RATE=1.0 + GO TO 330 +300 CONTINUE +! CONTRIBUTION TO ITERATION IS CHANGE IN T +! CORRESPONDING TO CHANGE IN E + RATE=(TTWO-TONE)/(ETWO-EONE) + TONE=TTWO + EONE=ETWO +330 CONTINUE +! COMPUTE ESTIMATED SATLIFT, TTWO + TTWO=TONE-EONE*RATE +! CONSIDER PSEUDO-ADIABAT, EW2, THROUGH TTWO AT P. +! COMPUTE ETW0=EW2-THM + ETWO=(TTWO+273.16)/PWRP-273.16 + ETWO=ETWO+WOBF(TTWO)-WOBF(ETWO )-THM +! CORRECTION TO TTWO IS EOR + EOR=ETWO*RATE + IF(ABS(EOR)-0.1000) 400,400,300 +400 SATLFT=TTWO-EOR + RETURN + END + + REAL(kind=RKIND) FUNCTION TCONOF(TEMP,DEWPT) +! COMPUTES CONDENSATION TEMPERATURE (DEGREES CENT) BY LIFTING + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: TEMP, DEWPT + REAL (kind=RKIND) :: S, T, DLT + + S=TEMP-DEWPT +! CONSIDER TEMP AND DEWPT TO BE LIKE UNITS (C OR K) + T=TEMP + IF(100.-TEMP) 4,5,5 +4 T=TEMP-273.16 +! COMPUTE CURVE FIT IN MOST EFFICIENT MANNER +5 DLT=S*(1.2185+0.001278*T+S*(-0.002190+11.73E-6*S-5.20E-6*T)) + TCONOF=T-DLT + RETURN + END + + REAL(kind=RKIND) FUNCTION WMROF(P,TD) +! COMPUTE MIXING RATIO (G/KG)...DEWPOINT (DEGREES C OR K)... PRESSURE (MB) + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: P, TD + REAL (kind=RKIND) :: T, X, WFW, FWESW + + T=TD + IF (100.-T) 3,4,4 +3 T=T-273.16 +! CURVE FIT CORRECTION FOR NON-IDEAL GAS +4 X=0.0200*(T-12.5+7500.0/P) + WFW=1.+0.0000045*P+0.00140*X*X +! COMPUTE ACCORDING TO STANDARD FORMULA + FWESW=WFW*VAPFW(T) + WMROF=621.97*(FWESW/(P-FWESW)) + RETURN + END + + REAL(kind=RKIND) FUNCTION THETA(T,P2,P1) + + REAL (kind=RKIND), INTENT(IN) :: T + REAL (kind=RKIND), INTENT(IN) :: P2 + REAL (kind=RKIND), INTENT(IN) :: P1 + + REAL (kind=RKIND) :: Rd, Cp, Rd_Cp + + Rd = 287.058_RKIND ! Specific gas constant for dry air (J/kg/K) + Cp = 1004.0_RKIND ! Specific heat at constant pressure (J/kg/K) + Rd_Cp = Rd/Cp + +! DRY ADIABATIC (T,P1) TO (THETA,P2) + + THETA = T * (P2/P1)**Rd_Cp ! DRY ADIABATIC (T,P1) TO (THETA,P2) + + RETURN + END + + REAL(kind=RKIND) FUNCTION DPTOF(EW) +! COMPUTE DEWPOINT, DPT, IN DEGREES C GIVEN WATER VAPOR PRESSURE (MB) + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: EW + REAL (kind=RKIND) :: TOL,X,BOT,DELTM, EDP,DTDE, DELT, DM + +! CREATE TOLERANCE TO DEGREE DESIRED + TOL=0.00010 + IF (EW-0.21382876E-09) 20,20,30 +20 DPTOF=-10000. + RETURN +30 IF (1013.0-EW) 20, 100, 100 +! CREATE GUESS BY INVERTING TETEN-S FORMULA +100 X=ALOG(EW/6.1078) + BOT=17.269388-X + DPTOF=(237.3*X)/BOT + BOT=BOT*EW + DELTM=0. +200 EDP=VAPFW(DPTOF) +! CORRECT GUESS BY DERIVATIVE OF TEMPERATURE WITH RESPECT TO VAPOR PRES. +! CALCULATED FROM INVERSE OF TETEN-S FORMULA + DTDE=(DPTOF+237.3)/BOT + DELT=DTDE*(EW-EDP) + DPTOF=DPTOF+DELT +! CHECK THAT ITERATION IS NOT IN AN ENDLESS CYCLE, A RARE SITUATION +! IF NEEDED. CHANGE -TOL- AND EXIT + DM=DELT-DELTM + IF(ABS(DM).GE.1.E-7) GO TO 10 ! IF DM VERY SMALL, ITERATION IS ENDLESS + TOL=ABS(DELT) + PRINT*,"TOLERANCE (TOL) IN DPTOF CHANGED TO ",TOL," (NORMAL TOL = .00010)" +10 DELTM=-DELT +! CHECK TO SEE IF ANSWER CLOSE ENOUGH, IF NOT ITERATE OVER CORRECTION + IF (ABS(DELT)-TOL) 300,300,200 +! CHANGE SO DEWPOINT IS ALWAYS LESS THAN THE TEMP. +! COMPATIBILITY WITH TOL IS FORCED +300 DPTOF=DPTOF-TOL + RETURN + END + + REAL(kind=RKIND) FUNCTION VAPFW(T) +! COMPUTE SATURATION VAPOR PRESSURE OVER WATER, VAPFW, IN MBS. + IMPLICIT NONE + REAL (kind=RKIND), INTENT(IN) :: T + REAL (kind=RKIND) :: X, POL +! CONSIDER T(TEMPERATURE) IN DEGREES C OR DEGREES K. + X=T + IF (100.0-X) 3,4,4 +3 X=X-273.16 +! CURVE FIT FOR RANGE -50 < T < 100 DEGREES C. +4 POL = 0.99999683E-00 + X *(-0.90826951E-02 + & + X *(0.70736169E-04 + X *(-0.61117958E-06 + & + X *(0.43884187E-08 + X *(-0.29883885E-10 + & + X *(0.21874425E-12 + X *(-0.17892321E-14 + & + X *(0.11112018E-16 + X *(-0.30994571E-19) )))))))) + POL=POL*POL + POL=POL*POL + VAPFW=6.107800/ (POL*POL) + RETURN + END + +END MODULE stability_indices