From e67ae45b21d552b9e7b4c467bc4dee34dffdcbfe Mon Sep 17 00:00:00 2001 From: thomasuhoh Date: Fri, 9 Jun 2023 12:48:30 +0200 Subject: [PATCH 1/3] Adding additional pressure level output --- src/core_atmosphere/Registry.xml | 43 ++ .../diagnostics/Registry_isobaric.xml | 245 ++++++--- .../diagnostics/mpas_isobaric_diagnostics.F | 479 +++++++++++++++--- 3 files changed, 634 insertions(+), 133 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 087e8fb620..45e218108b 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1011,38 +1011,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ diff --git a/src/core_atmosphere/diagnostics/Registry_isobaric.xml b/src/core_atmosphere/diagnostics/Registry_isobaric.xml index 853be6cde3..97b2be01bc 100644 --- a/src/core_atmosphere/diagnostics/Registry_isobaric.xml +++ b/src/core_atmosphere/diagnostics/Registry_isobaric.xml @@ -16,78 +16,102 @@ - + + - - - - + + + + - + + - + + - - + + + + - - - - - - - - - - - - - + - - + + @@ -100,114 +124,219 @@ + + + + + + + + - + + + + + + - - - - + + + + + - + + - + + - - + + + + + + - - - - + + + + - + + - + + - - + + + + - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F index e52c71b125..ef887a2be5 100644 --- a/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_isobaric_diagnostics.F @@ -35,7 +35,14 @@ module mpas_isobaric_diagnostics need_w_50, need_w_100, need_w_200, need_w_250, need_w_500, need_w_700, need_w_850, need_w_925, & need_vorticity_50, need_vorticity_100, need_vorticity_200, need_vorticity_250, need_vorticity_500, need_vorticity_700, need_vorticity_850, need_vorticity_925, & need_t_isobaric, need_z_isobaric, need_meanT_500_300 - logical :: need_temp, need_relhum, need_dewpoint, need_w, need_uzonal, need_umeridional, need_vorticity, need_height + LOGICAL :: need_qv_200, need_qv_250, need_qv_500, need_qv_700, need_qv_850, need_qv_925, need_qv, need_temp_300, need_umeridional_300, & + need_uzonal_300 , need_qv_300, need_temp_400, need_umeridional_400,need_uzonal_400 , need_qv_400, & + need_temp_600, need_umeridional_600,need_uzonal_600 , need_qv_600, need_height_300, need_height_400, need_height_600, & + need_temp_1000, need_umeridional_1000,need_uzonal_1000 , need_qv_1000, need_height_1000, & + need_temp_20, need_umeridional_20,need_uzonal_20 , need_qv_20, need_height_20, & + need_temp_800, need_umeridional_800,need_uzonal_800 , need_qv_800, need_height_800, & + need_temp_950, need_umeridional_950,need_uzonal_950 , need_qv_950, need_height_950 + logical :: need_temp, need_relhum, need_dewpoint, need_w, need_uzonal, need_umeridional, need_vorticity, need_height, need_qv contains @@ -100,6 +107,7 @@ subroutine isobaric_diagnostics_compute() need_umeridional = .false. need_vorticity = .false. need_height = .false. + need_qv = .false. need_mslp = MPAS_field_will_be_written('mslp') need_any_diags = need_any_diags .or. need_mslp @@ -302,6 +310,135 @@ subroutine isobaric_diagnostics_compute() need_meanT_500_300 = MPAS_field_will_be_written('meanT_500_300') need_any_diags = need_any_diags .or. need_meanT_500_300 +! Additional pressure levels + + need_qv_200 = MPAS_field_will_be_written('qv_200hPa') + need_qv = need_qv .or. need_qv_200 + need_any_diags = need_any_diags .or. need_qv_200 + need_qv_250 = MPAS_field_will_be_written('qv_250hPa') + need_qv = need_qv .or. need_qv_250 + need_any_diags = need_any_diags .or. need_qv_250 + need_qv_300 = MPAS_field_will_be_written('qv_300hPa') + need_qv = need_qv .or. need_qv_300 + need_any_diags = need_any_diags .or. need_qv_300 + need_qv_500 = MPAS_field_will_be_written('qv_500hPa') + need_qv = need_qv .or. need_qv_500 + need_any_diags = need_any_diags .or. need_qv_500 + need_qv_700 = MPAS_field_will_be_written('qv_700hPa') + need_qv = need_qv .or. need_qv_700 + need_any_diags = need_any_diags .or. need_qv_700 + need_qv_850 = MPAS_field_will_be_written('qv_850hPa') + need_qv = need_qv .or. need_qv_850 + need_any_diags = need_any_diags .or. need_qv_850 + need_qv_925 = MPAS_field_will_be_written('qv_925hPa') + need_qv = need_qv .or. need_qv_925 + need_any_diags = need_any_diags .or. need_qv_925 + + need_temp_20 = MPAS_field_will_be_written('temperature_20hPa') + need_temp = need_temp .or. need_temp_20 + need_any_diags = need_any_diags .or. need_temp_20 + need_umeridional_20 = MPAS_field_will_be_written('umeridional_20hPa') + need_umeridional = need_umeridional .or. need_umeridional_20 + need_any_diags = need_any_diags .or. need_umeridional_20 + need_uzonal_20 = MPAS_field_will_be_written('uzonal_20hPa') + need_uzonal = need_uzonal .or. need_uzonal_20 + need_any_diags = need_any_diags .or. need_uzonal_20 + need_qv_20 = MPAS_field_will_be_written('qv_20hPa') + need_qv = need_qv .or. need_qv_20 + need_any_diags = need_any_diags .or. need_qv_20 + + need_temp_300 = MPAS_field_will_be_written('temperature_300hPa') + need_temp = need_temp .or. need_temp_300 + need_any_diags = need_any_diags .or. need_temp_300 + need_umeridional_300 = MPAS_field_will_be_written('umeridional_300hPa') + need_umeridional = need_umeridional .or. need_umeridional_300 + need_any_diags = need_any_diags .or. need_umeridional_300 + need_uzonal_300 = MPAS_field_will_be_written('uzonal_300hPa') + need_uzonal = need_uzonal .or. need_uzonal_300 + need_any_diags = need_any_diags .or. need_uzonal_300 + + need_temp_400 = MPAS_field_will_be_written('temperature_400hPa') + need_temp = need_temp .or. need_temp_400 + need_any_diags = need_any_diags .or. need_temp_400 + need_umeridional_400 = MPAS_field_will_be_written('umeridional_400hPa') + need_umeridional = need_umeridional .or. need_umeridional_400 + need_any_diags = need_any_diags .or. need_umeridional_400 + need_uzonal_400 = MPAS_field_will_be_written('uzonal_400hPa') + need_uzonal = need_uzonal .or. need_uzonal_400 + need_any_diags = need_any_diags .or. need_uzonal_400 + need_qv_400 = MPAS_field_will_be_written('qv_400hPa') + need_qv = need_qv .or. need_qv_400 + need_any_diags = need_any_diags .or. need_qv_400 + + need_temp_600 = MPAS_field_will_be_written('temperature_600hPa') + need_temp = need_temp .or. need_temp_600 + need_any_diags = need_any_diags .or. need_temp_600 + need_umeridional_600 = MPAS_field_will_be_written('umeridional_600hPa') + need_umeridional = need_umeridional .or. need_umeridional_600 + need_any_diags = need_any_diags .or. need_umeridional_600 + need_uzonal_600 = MPAS_field_will_be_written('uzonal_600hPa') + need_uzonal = need_uzonal .or. need_uzonal_600 + need_any_diags = need_any_diags .or. need_uzonal_600 + need_qv_600 = MPAS_field_will_be_written('qv_600hPa') + need_qv = need_qv .or. need_qv_600 + need_any_diags = need_any_diags .or. need_qv_600 + + need_temp_800 = MPAS_field_will_be_written('temperature_800hPa') + need_temp = need_temp .or. need_temp_800 + need_any_diags = need_any_diags .or. need_temp_800 + need_umeridional_800 = MPAS_field_will_be_written('umeridional_800hPa') + need_umeridional = need_umeridional .or. need_umeridional_800 + need_any_diags = need_any_diags .or. need_umeridional_800 + need_uzonal_800 = MPAS_field_will_be_written('uzonal_800hPa') + need_uzonal = need_uzonal .or. need_uzonal_800 + need_any_diags = need_any_diags .or. need_uzonal_800 + need_qv_800 = MPAS_field_will_be_written('qv_800hPa') + need_qv = need_qv .or. need_qv_800 + need_any_diags = need_any_diags .or. need_qv_800 + + need_temp_950 = MPAS_field_will_be_written('temperature_950hPa') + need_temp = need_temp .or. need_temp_950 + need_any_diags = need_any_diags .or. need_temp_950 + need_umeridional_950 = MPAS_field_will_be_written('umeridional_950hPa') + need_umeridional = need_umeridional .or. need_umeridional_950 + need_any_diags = need_any_diags .or. need_umeridional_950 + need_uzonal_950 = MPAS_field_will_be_written('uzonal_950hPa') + need_uzonal = need_uzonal .or. need_uzonal_950 + need_any_diags = need_any_diags .or. need_uzonal_950 + need_qv_950 = MPAS_field_will_be_written('qv_950hPa') + need_qv = need_qv .or. need_qv_950 + need_any_diags = need_any_diags .or. need_qv_950 + + need_temp_1000 = MPAS_field_will_be_written('temperature_1000hPa') + need_temp = need_temp .or. need_temp_1000 + need_any_diags = need_any_diags .or. need_temp_1000 + need_umeridional_1000 = MPAS_field_will_be_written('umeridional_1000hPa') + need_umeridional = need_umeridional .or. need_umeridional_1000 + need_any_diags = need_any_diags .or. need_umeridional_1000 + need_uzonal_1000 = MPAS_field_will_be_written('uzonal_1000hPa') + need_uzonal = need_uzonal .or. need_uzonal_1000 + need_any_diags = need_any_diags .or. need_uzonal_1000 + need_qv_1000 = MPAS_field_will_be_written('qv_1000hPa') + need_qv = need_qv .or. need_qv_1000 + need_any_diags = need_any_diags .or. need_qv_1000 + + + need_height_100 = MPAS_field_will_be_written('height_100hPa') + need_height = need_height .or. need_height_100 + need_any_diags = need_any_diags .or. need_height_100 + need_height_300 = MPAS_field_will_be_written('height_300hPa') + need_height = need_height .or. need_height_300 + need_any_diags = need_any_diags .or. need_height_300 + need_height_400 = MPAS_field_will_be_written('height_400hPa') + need_height = need_height .or. need_height_400 + need_any_diags = need_any_diags .or. need_height_400 + need_height_600 = MPAS_field_will_be_written('height_600hPa') + need_height = need_height .or. need_height_600 + need_any_diags = need_any_diags .or. need_height_600 + need_height_1000 = MPAS_field_will_be_written('height_1000hPa') + need_height = need_height .or. need_height_1000 + need_any_diags = need_any_diags .or. need_height_1000 + if (need_any_diags) then call interp_diagnostics(mesh, state, 1, diag) end if @@ -344,15 +481,22 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) real (kind=RKIND), dimension(:,:), pointer :: t_isobaric real (kind=RKIND), dimension(:,:), pointer :: z_isobaric real (kind=RKIND), dimension(:), pointer :: meanT_500_300 - + + real (kind=RKIND), dimension(:), pointer :: temperature_20hPa real (kind=RKIND), dimension(:), pointer :: temperature_50hPa real (kind=RKIND), dimension(:), pointer :: temperature_100hPa real (kind=RKIND), dimension(:), pointer :: temperature_200hPa real (kind=RKIND), dimension(:), pointer :: temperature_250hPa + real (kind=RKIND), dimension(:), pointer :: temperature_300hPa + real (kind=RKIND), dimension(:), pointer :: temperature_400hPa real (kind=RKIND), dimension(:), pointer :: temperature_500hPa + real (kind=RKIND), dimension(:), pointer :: temperature_600hPa real (kind=RKIND), dimension(:), pointer :: temperature_700hPa + real (kind=RKIND), dimension(:), pointer :: temperature_800hPa real (kind=RKIND), dimension(:), pointer :: temperature_850hPa real (kind=RKIND), dimension(:), pointer :: temperature_925hPa + real (kind=RKIND), dimension(:), pointer :: temperature_950hPa + real (kind=RKIND), dimension(:), pointer :: temperature_1000hPa real (kind=RKIND), dimension(:), pointer :: relhum_50hPa real (kind=RKIND), dimension(:), pointer :: relhum_100hPa @@ -371,34 +515,55 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) real (kind=RKIND), dimension(:), pointer :: dewpoint_700hPa real (kind=RKIND), dimension(:), pointer :: dewpoint_850hPa real (kind=RKIND), dimension(:), pointer :: dewpoint_925hPa - + + real (kind=RKIND), dimension(:), pointer :: uzonal_20hPa real (kind=RKIND), dimension(:), pointer :: uzonal_50hPa real (kind=RKIND), dimension(:), pointer :: uzonal_100hPa real (kind=RKIND), dimension(:), pointer :: uzonal_200hPa real (kind=RKIND), dimension(:), pointer :: uzonal_250hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_300hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_400hPa real (kind=RKIND), dimension(:), pointer :: uzonal_500hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_600hPa real (kind=RKIND), dimension(:), pointer :: uzonal_700hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_800hPa real (kind=RKIND), dimension(:), pointer :: uzonal_850hPa real (kind=RKIND), dimension(:), pointer :: uzonal_925hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_950hPa + real (kind=RKIND), dimension(:), pointer :: uzonal_1000hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_20hPa real (kind=RKIND), dimension(:), pointer :: umeridional_50hPa real (kind=RKIND), dimension(:), pointer :: umeridional_100hPa real (kind=RKIND), dimension(:), pointer :: umeridional_200hPa real (kind=RKIND), dimension(:), pointer :: umeridional_250hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_300hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_400hPa real (kind=RKIND), dimension(:), pointer :: umeridional_500hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_600hPa real (kind=RKIND), dimension(:), pointer :: umeridional_700hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_800hPa real (kind=RKIND), dimension(:), pointer :: umeridional_850hPa real (kind=RKIND), dimension(:), pointer :: umeridional_925hPa - + real (kind=RKIND), dimension(:), pointer :: umeridional_950hPa + real (kind=RKIND), dimension(:), pointer :: umeridional_1000hPa + + real (kind=RKIND), dimension(:), pointer :: height_20hPa real (kind=RKIND), dimension(:), pointer :: height_50hPa real (kind=RKIND), dimension(:), pointer :: height_100hPa real (kind=RKIND), dimension(:), pointer :: height_200hPa real (kind=RKIND), dimension(:), pointer :: height_250hPa + real (kind=RKIND), dimension(:), pointer :: height_300hPa + real (kind=RKIND), dimension(:), pointer :: height_400hPa real (kind=RKIND), dimension(:), pointer :: height_500hPa + real (kind=RKIND), dimension(:), pointer :: height_600hPa real (kind=RKIND), dimension(:), pointer :: height_700hPa + real (kind=RKIND), dimension(:), pointer :: height_800hPa real (kind=RKIND), dimension(:), pointer :: height_850hPa real (kind=RKIND), dimension(:), pointer :: height_925hPa - + real (kind=RKIND), dimension(:), pointer :: height_950hPa + real (kind=RKIND), dimension(:), pointer :: height_1000hPa + real (kind=RKIND), dimension(:), pointer :: w_50hPa real (kind=RKIND), dimension(:), pointer :: w_100hPa real (kind=RKIND), dimension(:), pointer :: w_200hPa @@ -416,6 +581,22 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) real (kind=RKIND), dimension(:), pointer :: vorticity_700hPa real (kind=RKIND), dimension(:), pointer :: vorticity_850hPa real (kind=RKIND), dimension(:), pointer :: vorticity_925hPa + + real (kind=RKIND), dimension(:), pointer :: qv_20hPa + real (kind=RKIND), dimension(:), pointer :: qv_50hPa + real (kind=RKIND), dimension(:), pointer :: qv_100hPa + real (kind=RKIND), dimension(:), pointer :: qv_200hPa + real (kind=RKIND), dimension(:), pointer :: qv_250hPa + real (kind=RKIND), dimension(:), pointer :: qv_300hPa + real (kind=RKIND), dimension(:), pointer :: qv_400hPa + real (kind=RKIND), dimension(:), pointer :: qv_500hPa + real (kind=RKIND), dimension(:), pointer :: qv_600hPa + real (kind=RKIND), dimension(:), pointer :: qv_700hPa + real (kind=RKIND), dimension(:), pointer :: qv_800hPa + real (kind=RKIND), dimension(:), pointer :: qv_850hPa + real (kind=RKIND), dimension(:), pointer :: qv_925hPa + real (kind=RKIND), dimension(:), pointer :: qv_950hPa + real (kind=RKIND), dimension(:), pointer :: qv_1000hPa real (kind=RKIND) :: evp @@ -474,16 +655,23 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) call mpas_pool_get_array(diag, 't_isobaric', t_isobaric) call mpas_pool_get_array(diag, 'z_isobaric', z_isobaric) call mpas_pool_get_array(diag, 'meanT_500_300', meanT_500_300) - + + call mpas_pool_get_array(diag, 'temperature_20hPa', temperature_20hPa) call mpas_pool_get_array(diag, 'temperature_50hPa', temperature_50hPa) call mpas_pool_get_array(diag, 'temperature_100hPa', temperature_100hPa) call mpas_pool_get_array(diag, 'temperature_200hPa', temperature_200hPa) call mpas_pool_get_array(diag, 'temperature_250hPa', temperature_250hPa) + call mpas_pool_get_array(diag, 'temperature_300hPa', temperature_300hPa) + call mpas_pool_get_array(diag, 'temperature_400hPa', temperature_400hPa) call mpas_pool_get_array(diag, 'temperature_500hPa', temperature_500hPa) + call mpas_pool_get_array(diag, 'temperature_600hPa', temperature_600hPa) call mpas_pool_get_array(diag, 'temperature_700hPa', temperature_700hPa) + call mpas_pool_get_array(diag, 'temperature_800hPa', temperature_3800hPa) call mpas_pool_get_array(diag, 'temperature_850hPa', temperature_850hPa) call mpas_pool_get_array(diag, 'temperature_925hPa', temperature_925hPa) - + call mpas_pool_get_array(diag, 'temperature_950hPa', temperature_950hPa) + call mpas_pool_get_array(diag, 'temperature_1000hPa', temperature_1000hPa) + call mpas_pool_get_array(diag, 'relhum_50hPa', relhum_50hPa) call mpas_pool_get_array(diag, 'relhum_100hPa', relhum_100hPa) call mpas_pool_get_array(diag, 'relhum_200hPa', relhum_200hPa) @@ -501,33 +689,54 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) call mpas_pool_get_array(diag, 'dewpoint_700hPa', dewpoint_700hPa) call mpas_pool_get_array(diag, 'dewpoint_850hPa', dewpoint_850hPa) call mpas_pool_get_array(diag, 'dewpoint_925hPa', dewpoint_925hPa) - + + call mpas_pool_get_array(diag, 'uzonal_20hPa', uzonal_20hPa) call mpas_pool_get_array(diag, 'uzonal_50hPa', uzonal_50hPa) call mpas_pool_get_array(diag, 'uzonal_100hPa', uzonal_100hPa) call mpas_pool_get_array(diag, 'uzonal_200hPa', uzonal_200hPa) call mpas_pool_get_array(diag, 'uzonal_250hPa', uzonal_250hPa) + call mpas_pool_get_array(diag, 'uzonal_300hPa', uzonal_300hPa) + call mpas_pool_get_array(diag, 'uzonal_400hPa', uzonal_400hPa) call mpas_pool_get_array(diag, 'uzonal_500hPa', uzonal_500hPa) + call mpas_pool_get_array(diag, 'uzonal_600hPa', uzonal_600hPa) call mpas_pool_get_array(diag, 'uzonal_700hPa', uzonal_700hPa) + call mpas_pool_get_array(diag, 'uzonal_800hPa', uzonal_800hPa) call mpas_pool_get_array(diag, 'uzonal_850hPa', uzonal_850hPa) call mpas_pool_get_array(diag, 'uzonal_925hPa', uzonal_925hPa) - + call mpas_pool_get_array(diag, 'uzonal_950hPa', uzonal_950hPa) + call mpas_pool_get_array(diag, 'uzonal_1000hPa', uzonal_1000hPa) + + call mpas_pool_get_array(diag, 'umeridional_20hPa', umeridional_20hPa) call mpas_pool_get_array(diag, 'umeridional_50hPa', umeridional_50hPa) call mpas_pool_get_array(diag, 'umeridional_100hPa', umeridional_100hPa) call mpas_pool_get_array(diag, 'umeridional_200hPa', umeridional_200hPa) call mpas_pool_get_array(diag, 'umeridional_250hPa', umeridional_250hPa) + call mpas_pool_get_array(diag, 'umeridional_300hPa', umeridional_300hPa) + call mpas_pool_get_array(diag, 'umeridional_400hPa', umeridional_400hPa) call mpas_pool_get_array(diag, 'umeridional_500hPa', umeridional_500hPa) + call mpas_pool_get_array(diag, 'umeridional_600hPa', umeridional_600hPa) call mpas_pool_get_array(diag, 'umeridional_700hPa', umeridional_700hPa) + call mpas_pool_get_array(diag, 'umeridional_800hPa', umeridional_800hPa) call mpas_pool_get_array(diag, 'umeridional_850hPa', umeridional_850hPa) call mpas_pool_get_array(diag, 'umeridional_925hPa', umeridional_925hPa) - + call mpas_pool_get_array(diag, 'umeridional_950hPa', umeridional_950hPa) + call mpas_pool_get_array(diag, 'umeridional_1000hPa', umeridional_1000hPa) + + call mpas_pool_get_array(diag, 'height_20hPa', height_20hPa) call mpas_pool_get_array(diag, 'height_50hPa', height_50hPa) call mpas_pool_get_array(diag, 'height_100hPa', height_100hPa) call mpas_pool_get_array(diag, 'height_200hPa', height_200hPa) call mpas_pool_get_array(diag, 'height_250hPa', height_250hPa) + call mpas_pool_get_array(diag, 'height_300hPa', height_300hPa) + call mpas_pool_get_array(diag, 'height_400hPa', height_400hPa) call mpas_pool_get_array(diag, 'height_500hPa', height_500hPa) + call mpas_pool_get_array(diag, 'height_600hPa', height_600hPa) call mpas_pool_get_array(diag, 'height_700hPa', height_700hPa) + call mpas_pool_get_array(diag, 'height_800hPa', height_800hPa) call mpas_pool_get_array(diag, 'height_850hPa', height_850hPa) call mpas_pool_get_array(diag, 'height_925hPa', height_925hPa) + call mpas_pool_get_array(diag, 'height_950hPa', height_950hPa) + call mpas_pool_get_array(diag, 'height_1000hPa', height_1000hPa) call mpas_pool_get_array(diag, 'w_50hPa', w_50hPa) call mpas_pool_get_array(diag, 'w_100hPa', w_100hPa) @@ -658,18 +867,25 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) end if !interpolation to fixed pressure levels for fields located at cells centers and at mass points: - nIntP = 8 + nIntP = 15 if(.not.allocated(field_interp)) allocate(field_interp(nCells,nIntP) ) if(.not.allocated(press_interp)) allocate(press_interp(nCells,nIntP) ) do iCell = 1, nCells - press_interp(iCell,1) = 50.0_RKIND - press_interp(iCell,2) = 100.0_RKIND - press_interp(iCell,3) = 200.0_RKIND - press_interp(iCell,4) = 250.0_RKIND - press_interp(iCell,5) = 500.0_RKIND - press_interp(iCell,6) = 700.0_RKIND - press_interp(iCell,7) = 850.0_RKIND - press_interp(iCell,8) = 925.0_RKIND + press_interp(iCell,1) = 20.0_RKIND + press_interp(iCell,2) = 50.0_RKIND + press_interp(iCell,3) = 100.0_RKIND + press_interp(iCell,4) = 200.0_RKIND + press_interp(iCell,5) = 250.0_RKIND + press_interp(iCell,6) = 300.0_RKIND + press_interp(iCell,7) = 400.0_RKIND + press_interp(iCell,8) = 500.0_RKIND + press_interp(iCell,9) = 600.0_RKIND + press_interp(iCell,10) = 700.0_RKIND + press_interp(iCell,11) = 800.0_RKIND + press_interp(iCell,12) = 850.0_RKIND + press_interp(iCell,13) = 925.0_RKIND + press_interp(iCell,14) = 950.0_RKIND + press_interp(iCell,15) = 1000.0_RKIND enddo if(.not.allocated(press_in)) allocate(press_in(nCells,nVertLevels)) @@ -691,14 +907,21 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) - temperature_50hPa(1:nCells) = field_interp(1:nCells,1) - temperature_100hPa(1:nCells) = field_interp(1:nCells,2) - temperature_200hPa(1:nCells) = field_interp(1:nCells,3) - temperature_250hPa(1:nCells) = field_interp(1:nCells,4) - temperature_500hPa(1:nCells) = field_interp(1:nCells,5) - temperature_700hPa(1:nCells) = field_interp(1:nCells,6) - temperature_850hPa(1:nCells) = field_interp(1:nCells,7) - temperature_925hPa(1:nCells) = field_interp(1:nCells,8) + temperature_20hPa(1:nCells) = field_interp(1:nCells,1) + temperature_50hPa(1:nCells) = field_interp(1:nCells,2) + temperature_100hPa(1:nCells) = field_interp(1:nCells,3) + temperature_200hPa(1:nCells) = field_interp(1:nCells,4) + temperature_250hPa(1:nCells) = field_interp(1:nCells,5) + temperature_300hPa(1:nCells) = field_interp(1:nCells,6) + temperature_400hPa(1:nCells) = field_interp(1:nCells,7) + temperature_500hPa(1:nCells) = field_interp(1:nCells,8) + temperature_600hPa(1:nCells) = field_interp(1:nCells,9) + temperature_700hPa(1:nCells) = field_interp(1:nCells,10) + temperature_800hPa(1:nCells) = field_interp(1:nCells,11) + temperature_850hPa(1:nCells) = field_interp(1:nCells,12) + temperature_925hPa(1:nCells) = field_interp(1:nCells,13) + temperature_950hPa(1:nCells) = field_interp(1:nCells,14) + temperature_1000hPa(1:nCells) = field_interp(1:nCells,15) ! call mpas_log_write('--- end interpolate temperature:') end if @@ -712,14 +935,14 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) - relhum_50hPa(1:nCells) = field_interp(1:nCells,1) - relhum_100hPa(1:nCells) = field_interp(1:nCells,2) - relhum_200hPa(1:nCells) = field_interp(1:nCells,3) - relhum_250hPa(1:nCells) = field_interp(1:nCells,4) - relhum_500hPa(1:nCells) = field_interp(1:nCells,5) - relhum_700hPa(1:nCells) = field_interp(1:nCells,6) - relhum_850hPa(1:nCells) = field_interp(1:nCells,7) - relhum_925hPa(1:nCells) = field_interp(1:nCells,8) + relhum_50hPa(1:nCells) = field_interp(1:nCells,2) + relhum_100hPa(1:nCells) = field_interp(1:nCells,3) + relhum_200hPa(1:nCells) = field_interp(1:nCells,4) + relhum_250hPa(1:nCells) = field_interp(1:nCells,5) + relhum_500hPa(1:nCells) = field_interp(1:nCells,8) + relhum_700hPa(1:nCells) = field_interp(1:nCells,10) + relhum_850hPa(1:nCells) = field_interp(1:nCells,12) + relhum_925hPa(1:nCells) = field_interp(1:nCells,13) ! call mpas_log_write('--- end interpolate relative humidity:') end if @@ -732,16 +955,43 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) - dewpoint_50hPa(1:nCells) = field_interp(1:nCells,1) - dewpoint_100hPa(1:nCells) = field_interp(1:nCells,2) - dewpoint_200hPa(1:nCells) = field_interp(1:nCells,3) - dewpoint_250hPa(1:nCells) = field_interp(1:nCells,4) - dewpoint_500hPa(1:nCells) = field_interp(1:nCells,5) - dewpoint_700hPa(1:nCells) = field_interp(1:nCells,6) - dewpoint_850hPa(1:nCells) = field_interp(1:nCells,7) - dewpoint_925hPa(1:nCells) = field_interp(1:nCells,8) + dewpoint_50hPa(1:nCells) = field_interp(1:nCells,2) + dewpoint_100hPa(1:nCells) = field_interp(1:nCells,3) + dewpoint_200hPa(1:nCells) = field_interp(1:nCells,4) + dewpoint_250hPa(1:nCells) = field_interp(1:nCells,5) + dewpoint_500hPa(1:nCells) = field_interp(1:nCells,8) + dewpoint_700hPa(1:nCells) = field_interp(1:nCells,10) + dewpoint_850hPa(1:nCells) = field_interp(1:nCells,12) + dewpoint_925hPa(1:nCells) = field_interp(1:nCells,13) ! call mpas_log_write('--- end interpolate relative humidity:') end if + + if (NEED_QV) then + !... QVAPOR + do iCell = 1, nCells + do k = 1, nVertLevels + kk = nVertLevels+1-k + field_in(iCell,kk) = scalars(index_qv,k,iCell) + enddo + enddo + call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) + qv_20hPa(1:nCells) = field_interp(1:nCells,1) + qv_50hPa(1:nCells) = field_interp(1:nCells,2) + qv_100hPa(1:nCells) = field_interp(1:nCells,3) + qv_200hPa(1:nCells) = field_interp(1:nCells,4) + qv_250hPa(1:nCells) = field_interp(1:nCells,5) + qv_300hPa(1:nCells) = field_interp(1:nCells,6) + qv_400hPa(1:nCells) = field_interp(1:nCells,7) + qv_500hPa(1:nCells) = field_interp(1:nCells,8) + qv_600hPa(1:nCells) = field_interp(1:nCells,9) + qv_700hPa(1:nCells) = field_interp(1:nCells,10) + qv_800hPa(1:nCells) = field_interp(1:nCells,11) + qv_850hPa(1:nCells) = field_interp(1:nCells,12) + qv_925hPa(1:nCells) = field_interp(1:nCells,13) + qv_950hPa(1:nCells) = field_interp(1:nCells,14) + qv_1000hPa(1:nCells) = field_interp(1:nCells,15) +! call mpas_log_write('--- end interpolate qv:') + end if if (NEED_UZONAL) then !... u zonal wind: @@ -752,14 +1002,21 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) - uzonal_50hPa(1:nCells) = field_interp(1:nCells,1) - uzonal_100hPa(1:nCells) = field_interp(1:nCells,2) - uzonal_200hPa(1:nCells) = field_interp(1:nCells,3) - uzonal_250hPa(1:nCells) = field_interp(1:nCells,4) - uzonal_500hPa(1:nCells) = field_interp(1:nCells,5) - uzonal_700hPa(1:nCells) = field_interp(1:nCells,6) - uzonal_850hPa(1:nCells) = field_interp(1:nCells,7) - uzonal_925hPa(1:nCells) = field_interp(1:nCells,8) + uzonal_20hPa(1:nCells) = field_interp(1:nCells,1) + uzonal_50hPa(1:nCells) = field_interp(1:nCells,2) + uzonal_100hPa(1:nCells) = field_interp(1:nCells,3) + uzonal_200hPa(1:nCells) = field_interp(1:nCells,4) + uzonal_250hPa(1:nCells) = field_interp(1:nCells,5) + uzonal_300hPa(1:nCells) = field_interp(1:nCells,6) + uzonal_400hPa(1:nCells) = field_interp(1:nCells,7) + uzonal_500hPa(1:nCells) = field_interp(1:nCells,8) + uzonal_600hPa(1:nCells) = field_interp(1:nCells,9) + uzonal_700hPa(1:nCells) = field_interp(1:nCells,10) + uzonal_800hPa(1:nCells) = field_interp(1:nCells,11) + uzonal_850hPa(1:nCells) = field_interp(1:nCells,12) + uzonal_925hPa(1:nCells) = field_interp(1:nCells,13) + uzonal_950hPa(1:nCells) = field_interp(1:nCells,14) + uzonal_1000hPa(1:nCells) = field_interp(1:nCells,15) ! call mpas_log_write('--- end interpolate zonal wind:') end if @@ -772,14 +1029,21 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp) - umeridional_50hPa(1:nCells) = field_interp(1:nCells,1) - umeridional_100hPa(1:nCells) = field_interp(1:nCells,2) - umeridional_200hPa(1:nCells) = field_interp(1:nCells,3) - umeridional_250hPa(1:nCells) = field_interp(1:nCells,4) - umeridional_500hPa(1:nCells) = field_interp(1:nCells,5) - umeridional_700hPa(1:nCells) = field_interp(1:nCells,6) - umeridional_850hPa(1:nCells) = field_interp(1:nCells,7) - umeridional_925hPa(1:nCells) = field_interp(1:nCells,8) + umeridional_20hPa(1:nCells) = field_interp(1:nCells,1) + umeridional_50hPa(1:nCells) = field_interp(1:nCells,2) + umeridional_100hPa(1:nCells) = field_interp(1:nCells,3) + umeridional_200hPa(1:nCells) = field_interp(1:nCells,4) + umeridional_250hPa(1:nCells) = field_interp(1:nCells,5) + umeridional_300hPa(1:nCells) = field_interp(1:nCells,6) + umeridional_400hPa(1:nCells) = field_interp(1:nCells,7) + umeridional_500hPa(1:nCells) = field_interp(1:nCells,8) + umeridional_600hPa(1:nCells) = field_interp(1:nCells,9) + umeridional_700hPa(1:nCells) = field_interp(1:nCells,10) + umeridional_800hPa(1:nCells) = field_interp(1:nCells,11) + umeridional_850hPa(1:nCells) = field_interp(1:nCells,12) + umeridional_925hPa(1:nCells) = field_interp(1:nCells,13) + umeridional_950hPa(1:nCells) = field_interp(1:nCells,14) + umeridional_1000hPa(1:nCells) = field_interp(1:nCells,15) ! call mpas_log_write('--- end interpolate meridional wind:') end if @@ -806,14 +1070,21 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevelsP1,nIntP,press_in,field_in,press_interp,field_interp) - height_50hPa(1:nCells) = field_interp(1:nCells,1) - height_100hPa(1:nCells) = field_interp(1:nCells,2) - height_200hPa(1:nCells) = field_interp(1:nCells,3) - height_250hPa(1:nCells) = field_interp(1:nCells,4) - height_500hPa(1:nCells) = field_interp(1:nCells,5) - height_700hPa(1:nCells) = field_interp(1:nCells,6) - height_850hPa(1:nCells) = field_interp(1:nCells,7) - height_925hPa(1:nCells) = field_interp(1:nCells,8) + height_20hPa(1:nCells) = field_interp(1:nCells,1) + height_50hPa(1:nCells) = field_interp(1:nCells,2) + height_100hPa(1:nCells) = field_interp(1:nCells,3) + height_200hPa(1:nCells) = field_interp(1:nCells,4) + height_250hPa(1:nCells) = field_interp(1:nCells,5) + height_300hPa(1:nCells) = field_interp(1:nCells,6) + height_400hPa(1:nCells) = field_interp(1:nCells,7) + height_500hPa(1:nCells) = field_interp(1:nCells,8) + height_600hPa(1:nCells) = field_interp(1:nCells,9) + height_700hPa(1:nCells) = field_interp(1:nCells,10) + height_800hPa(1:nCells) = field_interp(1:nCells,11) + height_850hPa(1:nCells) = field_interp(1:nCells,12) + height_925hPa(1:nCells) = field_interp(1:nCells,13) + height_950hPa(1:nCells) = field_interp(1:nCells,14) + height_1000hPa(1:nCells) = field_interp(1:nCells,15) ! call mpas_log_write('--- end interpolate height:') !... vertical velocity @@ -824,14 +1095,14 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo call interp_tofixed_pressure(nCells,nVertLevelsP1,nIntP,press_in,field_in,press_interp,field_interp) - w_50hPa(1:nCells) = field_interp(1:nCells,1) - w_100hPa(1:nCells) = field_interp(1:nCells,2) - w_200hPa(1:nCells) = field_interp(1:nCells,3) - w_250hPa(1:nCells) = field_interp(1:nCells,4) - w_500hPa(1:nCells) = field_interp(1:nCells,5) - w_700hPa(1:nCells) = field_interp(1:nCells,6) - w_850hPa(1:nCells) = field_interp(1:nCells,7) - w_925hPa(1:nCells) = field_interp(1:nCells,8) + w_50hPa(1:nCells) = field_interp(1:nCells,2) + w_100hPa(1:nCells) = field_interp(1:nCells,3) + w_200hPa(1:nCells) = field_interp(1:nCells,4) + w_250hPa(1:nCells) = field_interp(1:nCells,5) + w_500hPa(1:nCells) = field_interp(1:nCells,8) + w_700hPa(1:nCells) = field_interp(1:nCells,10) + w_850hPa(1:nCells) = field_interp(1:nCells,12) + w_925hPa(1:nCells) = field_interp(1:nCells,13) if(allocated(field_in)) deallocate(field_in) if(allocated(press_in)) deallocate(press_in) @@ -893,6 +1164,64 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) if(allocated(pressureCp1) ) deallocate(pressureCp1 ) if(allocated(pressure_v) ) deallocate(pressure_v ) + +! if levels are below surface, set them to missing value (defined in the registry) + DO iCell = 1, nCells + IF(height_1000hPa(icell) > 360.) qv_1000hPa(icell) = -9999. + IF(height_950hPa(icell) > 700. ) qv_950hPa(icell) = -9999. + IF(height_925hPa(icell) > 1000.) qv_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) qv_850hPa(icell) = -9999. + IF(height_800hPa(icell) > 2150.) qv_800hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) qv_700hPa(icell) = -9999. + IF(height_600hPa(icell) > 4550.) qv_600hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) qv_500hPa(icell) = -9999. + IF(height_400hPa(icell) > 7900.) qv_400hPa(icell) = -9999. + + IF(height_1000hPa(icell) > 360.) temperature_1000hPa(icell) = -9999. + IF(height_950hPa(icell) > 700. ) temperature_950hPa(icell) = -9999. + IF(height_925hPa(icell) > 1000.) temperature_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) temperature_850hPa(icell) = -9999. + IF(height_800hPa(icell) > 2150.) temperature_800hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) temperature_700hPa(icell) = -9999. + IF(height_600hPa(icell) > 4550.) temperature_600hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) temperature_500hPa(icell) = -9999. + IF(height_400hPa(icell) > 7900.) temperature_400hPa(icell) = -9999. + + IF(height_1000hPa(icell) > 360.) umeridional_1000hPa(icell) = -9999. + IF(height_950hPa(icell) > 700. ) umeridional_950hPa(icell) = -9999. + IF(height_925hPa(icell) > 1000.) umeridional_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) umeridional_850hPa(icell) = -9999. + IF(height_800hPa(icell) > 2150.) umeridional_800hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) umeridional_700hPa(icell) = -9999. + IF(height_600hPa(icell) > 4550.) umeridional_600hPa(icell) =-9999. + IF(height_500hPa(icell) > 6000.) umeridional_500hPa(icell) = -9999. + IF(height_400hPa(icell) > 7900.) umeridional_400hPa(icell) = -9999. + + IF(height_1000hPa(icell) > 360.) uzonal_1000hPa(icell) = -9999. + IF(height_950hPa(icell) > 700. ) uzonal_950hPa(icell) = -9999. + IF(height_925hPa(icell) > 1000.) uzonal_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) uzonal_850hPa(icell) = -9999. + IF(height_800hPa(icell) > 2150.) uzonal_800hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) uzonal_700hPa(icell) = -9999. + IF(height_600hPa(icell) > 4550.) uzonal_600hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) uzonal_500hPa(icell) = -9999. + IF(height_400hPa(icell) > 7900.) uzonal_400hPa(icell) = -9999. + + IF(height_925hPa(icell) > 1000.) w_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) w_850hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) w_700hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) w_500hPa(icell) = -9999. + + IF(height_925hPa(icell) > 1000.) vorticity_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) vorticity_850hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) vorticity_700hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) vorticity_500hPa(icell) = -9999. + + IF(height_925hPa(icell) > 1000.) dewpoint_925hPa(icell) = -9999. + IF(height_850hPa(icell) > 1700.) dewpoint_850hPa(icell) = -9999. + IF(height_700hPa(icell) > 3300.) dewpoint_700hPa(icell) = -9999. + IF(height_500hPa(icell) > 6000.) dewpoint_500hPa(icell) = -9999. + ENDDO if (need_mslp) then !... compute SLP (requires temp, height, pressure, qvapor) From 6ad1bbc3fc177d4890cf82f334bbd886b5f85d75 Mon Sep 17 00:00:00 2001 From: thomasuhoh Date: Fri, 9 Jun 2023 14:59:59 +0200 Subject: [PATCH 2/3] Adding additional 2d diagnostic variables --- src/core_atmosphere/Registry.xml | 12 ++ .../diagnostics/Registry_convective.xml | 23 +++ .../diagnostics/mpas_convective_diagnostics.F | 157 +++++++++++++++++- 3 files changed, 190 insertions(+), 2 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 45e218108b..01be80e53b 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -490,6 +490,7 @@ + #ifdef DO_PHYSICS @@ -857,6 +858,7 @@ + #endif @@ -1122,6 +1124,7 @@ + #endif @@ -1142,6 +1145,15 @@ + + + + + + + + + #ifdef DO_PHYSICS diff --git a/src/core_atmosphere/diagnostics/Registry_convective.xml b/src/core_atmosphere/diagnostics/Registry_convective.xml index 59ad1c2f58..f8d98f7e12 100644 --- a/src/core_atmosphere/diagnostics/Registry_convective.xml +++ b/src/core_atmosphere/diagnostics/Registry_convective.xml @@ -55,5 +55,28 @@ + + + + + + + + + + + + + + + diff --git a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F index 163ee3774f..2ff7940378 100644 --- a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F @@ -14,6 +14,8 @@ module mpas_convective_diagnostics type (MPAS_pool_type), pointer :: mesh type (MPAS_pool_type), pointer :: state type (MPAS_pool_type), pointer :: diag + type (MPAS_pool_type), pointer :: diag_physics + type (MPAS_pool_type), pointer :: sfc_input type (MPAS_clock_type), pointer :: clock @@ -65,6 +67,8 @@ subroutine convective_diagnostics_setup(all_pools, simulation_clock) call mpas_pool_get_subpool(all_pools, 'mesh', mesh) call mpas_pool_get_subpool(all_pools, 'state', state) call mpas_pool_get_subpool(all_pools, 'diag', diag) + call mpas_pool_get_subpool(all_pools, 'diag_physics', diag_physics) + call mpas_pool_get_subpool(all_pools, 'sfc_input', sfc_input) clock => simulation_clock @@ -240,6 +244,17 @@ subroutine convective_diagnostics_compute() real (kind=RKIND), dimension(:), pointer :: cin real (kind=RKIND), dimension(:), pointer :: lcl real (kind=RKIND), dimension(:), pointer :: lfc + real (kind=RKIND), dimension(:), pointer :: fzlev + real (kind=RKIND), dimension(:), pointer :: speed_surface + real (kind=RKIND), dimension(:), pointer :: heatidx + real (kind=RKIND), dimension(:), pointer :: ivt + real (kind=RKIND), dimension(:), pointer :: ztd + real (kind=RKIND), dimension(:), pointer :: latCell + real (kind=RKIND), dimension(:), pointer :: terrain + real (kind=RKIND), dimension(:), pointer :: gust_10m + real (kind=RKIND), dimension(:), pointer :: ust + real (kind=RKIND), dimension(:), pointer :: v10 + real (kind=RKIND), dimension(:), pointer :: u10 real (kind=RKIND), dimension(:), pointer :: srh_0_1km real (kind=RKIND), dimension(:), pointer :: srh_0_3km real (kind=RKIND), dimension(:), pointer :: uzonal_surface @@ -271,13 +286,16 @@ subroutine convective_diagnostics_compute() 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 + real (kind=RKIND) :: part1, part2 + real (kind=RKIND) :: zf, zhd, zdiff, zwd_sfc, term1, term2 real (kind=RKIND), dimension(:,:), allocatable :: temperature, dewpoint - real (kind=RKIND) :: evp + real (kind=RKIND) :: evp, tf2, tf, rh2, pressure_in logical :: need_cape, need_cin, 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 + need_uzonal_6km, need_umerid_sfc, need_umerid_1km, need_umerid_6km, need_tsfc, need_tdsfc, need_fzlev , & + need_speed_surface, need_heatidx, need_temp_wet, need_ivt, need_ztd logical :: need_any_diags need_any_diags = .false. @@ -311,17 +329,46 @@ subroutine convective_diagnostics_compute() need_tdsfc = MPAS_field_will_be_written('dewpoint_surface') need_any_diags = need_any_diags .or. need_tdsfc + need_fzlev = MPAS_field_will_be_written('fzlev') + need_any_diags = need_any_diags .or. need_fzlev + need_speed_surface = MPAS_field_will_be_written('speed_surface') + need_any_diags = need_any_diags .or. need_speed_surface + need_heatidx = MPAS_field_will_be_written('heatidx') + need_any_diags = need_any_diags .or. need_heatidx + need_temp_wet = MPAS_field_will_be_written('temp_wet') + need_any_diags = need_any_diags .or. need_temp_wet + need_ivt = MPAS_field_will_be_written('ivt') + need_any_diags = need_any_diags .or. need_ivt + need_ztd = MPAS_field_will_be_written('ztd') + need_any_diags = need_any_diags .or. need_ztd + if (need_any_diags) then call mpas_pool_get_dimension(mesh, 'nCells', nCells) call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve) call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(mesh, 'nVertLevelsP1', nVertLevelsP1) + call mpas_pool_get_array(mesh, 'latCell', latCell) call mpas_pool_get_array(diag, 'cape', cape) call mpas_pool_get_array(diag, 'cin', cin) call mpas_pool_get_array(diag, 'lcl', lcl) call mpas_pool_get_array(diag, 'lfc', lfc) + call mpas_pool_get_array(diag, 'fzlev', fzlev) + call mpas_pool_get_array(diag, 'speed_surface', speed_surface) + call mpas_pool_get_array(diag, 'heatidx', heatidx) + call mpas_pool_get_array(diag_physics, 't2m', t2m) + call mpas_pool_get_array(diag_physics, 'q2', q2) + call mpas_pool_get_array(diag_physics, 'ust', ust) + call mpas_pool_get_array(diag_physics, 'v10', v10) + call mpas_pool_get_array(diag_physics, 'u10', u10) + call mpas_pool_get_array(diag, 'surface_pressure', surface_pressure) + call mpas_pool_get_array(diag, 'relhum_2m', relhum_2m) + call mpas_pool_get_array(diag, 'temp_wet', temp_wet) + call mpas_pool_get_array(diag, 'ivt', ivt) + call mpas_pool_get_array(diag, 'ztd', ztd) + call mpas_pool_get_array(sfc_input, 'ter', terrain) + call mpas_pool_get_array(diag, 'gust_10m', gust_10m) call mpas_pool_get_array(diag, 'srh_0_1km', srh_0_1km) call mpas_pool_get_array(diag, 'srh_0_3km', srh_0_3km) call mpas_pool_get_array(diag, 'uzonal_surface', uzonal_surface) @@ -384,6 +431,7 @@ subroutine convective_diagnostics_compute() umeridional_surface(iCell) = umeridional(1,iCell) temperature_surface(iCell) = temperature(1,iCell) dewpoint_surface(iCell) = dewpoint(1,iCell) + speed_surface(iCell) = sqrt(uzonal_surface(iCell)**2. + umeridional_surface(iCell)**2.) if (need_uzonal_1km) then uzonal_1km(iCell) = column_height_value(uzonal(1:nVertLevels,iCell), zp, 1000.0_RKIND, nVertLevels) end if @@ -471,6 +519,111 @@ subroutine convective_diagnostics_compute() end do end if + if (need_ivt .or. need_ztd) then + do iCell=1, nCellsSolve + p_in(1:nVertLevels) = (pressure_p(1:nVertLevels,iCell) + pressure_base(1:nVertLevels,iCell)) + t_in(1:nVertLevels) = temperature(1:nVertLevels,iCell) + qvapor_in(1:nVertLevels) = scalars(index_qv,1:nVertLevels,iCell) + !for IVT + part1=0._RKIND + part2=0._RKIND + do k = 1,nVertLevels-1 + if(p_in(k+1) >= 25000.0_RKIND)THEN + part1 = part1 + uzonal(k,iCell)*scalars(index_qv,k,iCell) * (p_in(k+1) - p_in(k)) + part2 = part2 + umeridional(k,iCell)*scalars(index_qv,k,iCell) * (p_in(k+1) - p_in(k)) + ivt(iCell) = sqrt((part1/9.806)**2+(part2/9.806)**2) + endif + enddo + + !for ZTD + part1=0._RKIND + part2=0._RKIND + zwd_sfc = 0._RKIND +! Saastamoinen + zf = 1. - 2.66e-3*cos(2.*latCell(iCell))-2.8e-7*terrain(iCell) + zhd = 2.2768e-5*surface_pressure(iCell)/zf + + do k = 1,nVertLevels-1 + zdiff = (height(k+1,iCell) - height(k,iCell))/0.622_RKIND + part1 = p_in(k)*qvapor_in(k)/t_in(k) + term1 = part1 * zdiff * 2.21e-7 + term2 = part1 * zdiff * 3.739e-3 / t_in(k) + zwd_sfc = zwd_sfc + term1+term2 + enddo + ztd(iCell) = (zwd_sfc + zhd)*100.0_RKIND + gust_10m(iCell) = sqrt(u10(iCell)*u10(iCell)+v10(iCell)*v10(iCell))+ (3.*2.4*ust(iCell)) + end do ! cells + end if + +! from WRF + IF (need_fzlev) then + fzlev = -999._RKIND + do iCell = 1,nCells + do k = 1,nVertLevels + IF (temperature(k,icell) .lt. 273.15_RKIND) THEN + ! Freezing level calculation + ! -------------------------- + IF (fzlev (icell) .eq. -999._RKIND) THEN ! At freezing level + IF (k .ne. 1) THEN ! If not at surface, interpolate. + fzlev (icell) = height (k-1,icell) + & + ((273.15_RKIND - temperature (k-1,icell)) & + /(temperature (k,icell) - temperature (k-1,icell))) & + *(height (k,icell) - height (k-1,icell)) + ELSE ! If at surface, use first level. + fzlev(icell) = height (k,icell) + ENDIF + ENDIF + ENDIF + enddo + enddo + ENDIF + + DO iCell = 1, nCells + IF(t2m(iCell) > 273.15_RKIND )THEN + relhum_2m(iCell) = 0.263*surface_pressure(iCell)*q2(iCell)*(exp(17.625*(t2m(iCell)-273.15_RKIND)/(t2m(iCell)-30.1)))**(-1.) + ELSE IF( t2m(iCell) .le. 273.15_RKIND) THEN +! https://journals.ametsoc.org/view/journals/apme/57/6/jamc-d-17-0334.1.xml + relhum_2m(iCell) = 0.263*surface_pressure(iCell)*q2(iCell)*(exp(22.587*(t2m(iCell)-273.15_RKIND)/(t2m(iCell)+0.71)))**(-1.) + ENDIF + ENDDO + where(relhum_2m > 100.) relhum_2m = 100._RKIND + + IF (need_heatidx) then + do iCell = 1,nCells + IF ( t2m(icell) > 294.26111) THEN + tF = ( (t2m(icell) - 273.15_RKIND) * (REAL (9)/REAL (5)) ) + REAL ( 32 ) + tF2 = tF ** 2 + rh2 = relhum_2m(iCell)/100. ** 2 + heatidx(icell) = -42.379 & + + ( 2.04901523 * tF ) & + + ( 10.14333127 * relhum_2m(iCell)/100. ) & + - ( 0.22475541 * tF * relhum_2m(iCell)/100. ) & + - ( 6.83783E-03 * tF2 ) & + - ( 5.481717E-02 * rh2 ) & + + ( 1.22874E-03 * tF2 * relhum_2m(iCell)/100. ) & + + ( 8.5282E-04 * tF * rh2 ) & + - ( 1.99E-06 * tF2 * rh2 ) + heatidx(icell) = ((heatidx(icell) - REAL (32)) * (REAL (5)/REAL (9))) + 273.15_RKIND + ELSE + heatidx(icell) = t2m(icell) + END IF + enddo + ENDIF + + IF(need_temp_wet)THEN +! Stull 2011: doi:10.1175/JAMC-D-11-0143.1 + DO iCell =1, nCells + IF(t2m(iCell) > 253.15 .AND. relhum_2m(iCell) .ge. 5. .AND. relhum_2m(iCell) .le. 99.)THEN + temp_wet(iCell) = (t2m(iCell)-273.15_RKIND)*atan(0.151977*sqrt(relhum_2m(iCell)+8.313659)) & + + atan(t2m(iCell)-273.15_RKIND+relhum_2m(iCell))& + - atan(relhum_2m(iCell)-1.676331)+0.00391838*((relhum_2m(iCell))**(1.5))*atan(0.023101*relhum_2m(iCell)) & + - 4.686035 +273.15_RKIND + ELSE + temp_wet(iCell) = -9999._RKIND + ENDIF + ENDDO + ENDIF + deallocate(temperature) deallocate(dewpoint) From 200de9e9afed0893b8c63223294f01b9db253ceb Mon Sep 17 00:00:00 2001 From: thomasuhoh Date: Wed, 14 Jun 2023 07:41:14 +0200 Subject: [PATCH 3/3] zero out diagnostic variables before calculation --- .../diagnostics/mpas_convective_diagnostics.F | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F index 2ff7940378..6d7b0bd4b1 100644 --- a/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_convective_diagnostics.F @@ -405,6 +405,18 @@ subroutine convective_diagnostics_compute() allocate(t_in(nVertLevels)) allocate(td_in(nVertLevels)) + ctt = 0._RKIND + ctpres = 0._RKIND + ivt = 0._RKIND + ztd = 0._RKIND + cape = 0._RKIND + cin = 0._RKIND + gust_10m = 0._RKIND + fzlev = 0._RKIND + relhum_2m = 0._RKIND + heatidx = 0._RKIND + temp_wet = 0._RKIND + do iCell = 1,nCells do k = 1,nVertLevels temperature(k,iCell) = (theta_m(k,iCell)/(1._RKIND+rvord*scalars(index_qv,k,iCell)))*exner(k,iCell)