From 374a9fc3a3f945fdaacd06802aaadde2e3e4499c Mon Sep 17 00:00:00 2001 From: Julio P R Fernandez Date: Thu, 26 Mar 2026 09:12:03 -0300 Subject: [PATCH 1/3] =?UTF-8?q?Atualizando=20o=20topico=20History=20com=20?= =?UTF-8?q?informa=C3=A7=C3=B5es=20de=20meu=20fork?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index bae3a13e..2ca1c02c 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,9 @@ The MONAN Model is managed by a scientific committee appointed by INPE's directo History ==== + +Versão do MONAN com novas variáveis diagnósticas no meu fork (Pablo) + - Version 1.4.2-rc (Release Candidate) - Speed up by about 7%. Cleanup to become the initial version of the C3P Community Cloud-Convection Parameterization. Joint development between INPE and NOAA/GSL. Removed files not needed anymore. Adding effects of PCW (Neelin et al. 2009) and vertical shear of horizontal wind on the entrainment rate. This should improve model simulations of MCSs. Additional trigger function based on Xie et al 2019. Additional comments and references. - Version 1.4.1-rc (Release Candidate) - Terrain height (ter, calculated previously in the pre processing with init_atmosphere_model) included in the Registry.xml's input section so that can be read from 'init' file and post processed. kubota relhum evalute modification on mpas_isobaric_diagnostics.F. - Version 1.4.0-rc (Release Candidate) - Changing the number of isobaric levels from 22 to 18. From 47294426fb203f5cfaf8d86eec7f45619779448e Mon Sep 17 00:00:00 2001 From: Julio P R Fernandez Date: Thu, 26 Mar 2026 15:04:39 -0300 Subject: [PATCH 2/3] =?UTF-8?q?Atualizando=20Registry=20e=20codigos=20das?= =?UTF-8?q?=20novas=20vari=C3=A1veis=20diagnosticas?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/core_atmosphere/Registry.xml | 5 + src/core_atmosphere/diagnostics/Makefile | 5 +- .../diagnostics/Registry_convective.xml | 21 +- .../diagnostics/Registry_isobaric.xml | 3 + .../diagnostics/mpas_convective_diagnostics.F | 72 +++- .../diagnostics/mpas_isobaric_diagnostics.F | 22 +- .../diagnostics/mpas_stability_indices.F | 363 ++++++++++++++++++ 7 files changed, 481 insertions(+), 10 deletions(-) create mode 100644 src/core_atmosphere/diagnostics/mpas_stability_indices.F diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 1f778472..10e7f7e8 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1026,8 +1026,13 @@ #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..40691133 100644 --- a/src/core_atmosphere/diagnostics/Registry_convective.xml +++ b/src/core_atmosphere/diagnostics/Registry_convective.xml @@ -4,11 +4,28 @@ + + + + + + + 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..ea8f6791 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 @@ -238,6 +239,8 @@ subroutine convective_diagnostics_compute() ! Fields that are computed in this routine 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 +254,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 +275,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 +286,28 @@ 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_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_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') @@ -320,6 +342,8 @@ subroutine convective_diagnostics_compute() 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 +357,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) @@ -463,7 +491,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 +499,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 +654,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 +701,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 From e93cfeb339c6896c42d9d2c38e61d8ac97d7e195 Mon Sep 17 00:00:00 2001 From: Julio P R Fernandez Date: Thu, 26 Mar 2026 22:44:08 -0300 Subject: [PATCH 3/3] =?UTF-8?q?Adiciona=20vari=C3=A1vel=20Surface=20CAPE?= =?UTF-8?q?=20&=20CIN?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/core_atmosphere/Registry.xml | 2 ++ .../diagnostics/Registry_convective.xml | 6 +++++ .../diagnostics/mpas_convective_diagnostics.F | 26 +++++++++++++++++++ 3 files changed, 34 insertions(+) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 10e7f7e8..955f9b1f 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1029,6 +1029,8 @@ + + diff --git a/src/core_atmosphere/diagnostics/Registry_convective.xml b/src/core_atmosphere/diagnostics/Registry_convective.xml index 40691133..039c6abc 100644 --- a/src/core_atmosphere/diagnostics/Registry_convective.xml +++ b/src/core_atmosphere/diagnostics/Registry_convective.xml @@ -13,6 +13,12 @@ + + + + diff --git a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F index ea8f6791..c23f2673 100644 --- a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F @@ -237,6 +237,8 @@ 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 @@ -286,6 +288,7 @@ subroutine convective_diagnostics_compute() real (kind=RKIND) :: evp + 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 @@ -300,6 +303,11 @@ subroutine convective_diagnostics_compute() 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') @@ -340,6 +348,9 @@ 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) @@ -479,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