From 1f7004b95588bd9cc0191b962790f7ec38cbad5a Mon Sep 17 00:00:00 2001 From: saraqzhang Date: Wed, 8 Mar 2023 17:11:47 -0500 Subject: [PATCH 01/10] add ensemble sdv to diag ( selected variables ) --- src/Applications/LDAS_App/GEOSldas_HIST.rc | 10 ++ .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 96 ++++++++++++++++++- .../GEOS_LandAssimGridComp.F90 | 86 ++++++++++++++++- 3 files changed, 183 insertions(+), 9 deletions(-) diff --git a/src/Applications/LDAS_App/GEOSldas_HIST.rc b/src/Applications/LDAS_App/GEOSldas_HIST.rc index 40299e69..b8e5ec40 100644 --- a/src/Applications/LDAS_App/GEOSldas_HIST.rc +++ b/src/Applications/LDAS_App/GEOSldas_HIST.rc @@ -419,11 +419,21 @@ COLLECTIONS: 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , + 'WCSFsdv' , 'ENSAVG' , 'SFMC_FCSTsdv' , + 'WCRZsdv' , 'ENSAVG' , 'RZMC_FCSTsdv' , + 'WCPRsdv' , 'ENSAVG' , 'PRMC_FCSTsdv' , + 'TPSURFsdv' , 'ENSAVG' , 'TSURF_FCSTsdv' , + 'TSOIL1TILEsdv' , 'ENSAVG' , 'TSOIL1_FCSTsdv' , 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , + 'WCSF_ANAsdv' , 'LANDASSIM' , 'SFMC_ANAsdv' , + 'WCRZ_ANAsdv' , 'LANDASSIM' , 'RZMC_ANAsdv' , + 'WCPR_ANAsdv' , 'LANDASSIM' , 'PRMC_ANAsdv' , + 'TPSURF_ANAsdv' , 'LANDASSIM' , 'TSURF_ANAsdv' , + 'TSOIL1_ANAsdv' , 'LANDASSIM' , 'TSOIL1_ANAsdv' :: # For lndfcstana, *.frequency and *.ref_time must be consistent with the LDAS.rc resource diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index b32e8044..e4ada4b7 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -656,6 +656,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'ave_catchment_temp_incl_snw ens.sdv',& + UNITS = 'K' ,& + SHORT_NAME = 'TPSURFsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'temperature_top_snow_layer',& UNITS = 'K' ,& @@ -800,6 +809,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_surface_layer ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCSFsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'water_root_zone' ,& UNITS = 'm3 m-3' ,& @@ -809,6 +827,16 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_root_zone ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCRZsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC, & LONG_NAME = 'water_ave_prof' ,& UNITS = 'm3 m-3' ,& @@ -818,6 +846,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_ave_prof ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCPRsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'soil_temperatures_layer_1' ,& UNITS = 'K' ,& @@ -827,6 +864,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'soil_temperatures_layer_1 ens.sdv' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TSOIL1TILEsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'soil_temperatures_layer_2' ,& UNITS = 'K' ,& @@ -2287,7 +2333,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: HLATN,HLATN_enavg real, dimension(:),pointer :: QINFIL,QINFIL_enavg real, dimension(:),pointer :: GHFLX,GHFLX_enavg - real, dimension(:),pointer :: TPSURF,TPSURF_enavg + real, dimension(:),pointer :: TPSURF,TPSURF_enavg,TPSURF_ensdv real, dimension(:),pointer :: TPSNOW,TPSNOW_enavg real, dimension(:),pointer :: TPUNST,TPUNST_enavg real, dimension(:),pointer :: TPSAT,TPSAT_enavg @@ -2304,10 +2350,10 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: WET1,WET1_enavg real, dimension(:),pointer :: WET2,WET2_enavg real, dimension(:),pointer :: WET3,WET3_enavg - real, dimension(:),pointer :: WCSF,WCSF_enavg - real, dimension(:),pointer :: WCRZ,WCRZ_enavg - real, dimension(:),pointer :: WCPR,WCPR_enavg - real, dimension(:),pointer :: TP1,TP1_enavg + real, dimension(:),pointer :: WCSF,WCSF_enavg,WCSF_ensdv + real, dimension(:),pointer :: WCRZ,WCRZ_enavg,WCRZ_ensdv + real, dimension(:),pointer :: WCPR,WCPR_enavg,WCPR_ensdv + real, dimension(:),pointer :: TP1,TP1_enavg,TP1_ensdv real, dimension(:),pointer :: TP2,TP2_enavg real, dimension(:),pointer :: TP3,TP3_enavg real, dimension(:),pointer :: TP4,TP4_enavg @@ -2859,6 +2905,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, TPSURF_enavg, 'TPSURF' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, TPSURF_ensdv, 'TPSURFsdv' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, TPSNOW_enavg, 'TPSNOW' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TPUNST_enavg, 'TPUNST' ,rc=status) @@ -2891,12 +2939,20 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, WCSF_enavg, 'WCSF' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, WCSF_ensdv, 'WCSFsdv' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, WCRZ_enavg, 'WCRZ' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, WCRZ_ensdv, 'WCRZsdv' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, WCPR_enavg, 'WCPR' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, WCPR_ensdv, 'WCPRsdv' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, TP1_enavg, 'TSOIL1TILE' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, TP1_ensdv, 'TSOIL1TILEsdv' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, TP2_enavg, 'TSOIL2TILE' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TP3_enavg, 'TSOIL3TILE' ,rc=status) @@ -3107,6 +3163,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(QINFIL_enavg)) QINFIL_enavg = 0.0 if(associated(GHFLX_enavg)) GHFLX_enavg = 0.0 if(associated(TPSURF_enavg)) TPSURF_enavg = 0.0 + if(associated(TPSURF_ensdv)) TPSURF_ensdv = 0.0 if(associated(TPSNOW_enavg)) TPSNOW_enavg = 0.0 if(associated(TPUNST_enavg)) TPUNST_enavg = 0.0 if(associated(TPSAT_enavg)) TPSAT_enavg = 0.0 @@ -3123,9 +3180,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = 0.0 if(associated(WET3_enavg)) WET3_enavg = 0.0 if(associated(WCSF_enavg)) WCSF_enavg = 0.0 + if(associated(WCSF_ensdv)) WCSF_ensdv = 0.0 if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 + if(associated(WCRZ_ensdv)) WCRZ_ensdv = 0.0 if(associated(WCPR_enavg)) WCPR_enavg = 0.0 + if(associated(WCPR_ensdv)) WCPR_ensdv = 0.0 if(associated(TP1_enavg)) TP1_enavg = 0.0 + if(associated(TP1_ensdv)) TP1_ensdv = 0.0 if(associated(TP2_enavg)) TP2_enavg = 0.0 if(associated(TP3_enavg)) TP3_enavg = 0.0 if(associated(TP4_enavg)) TP4_enavg = 0.0 @@ -3306,6 +3367,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) GHFLX_enavg = GHFLX_enavg + GHFLX if(associated(TPSURF_enavg) .and. associated(TPSURF)) & TPSURF_enavg = TPSURF_enavg + TPSURF + if(associated(TPSURF_ensdv) .and. associated(TPSURF)) & + TPSURF_ensdv = TPSURF_ensdv + TPSURF*TPSURF if(associated(TPSNOW_enavg) .and. associated(TPSNOW)) & TPSNOW_enavg = TPSNOW_enavg + TPSNOW if(associated(TPUNST_enavg) .and. associated(TPUNST)) & @@ -3338,12 +3401,20 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) WET3_enavg = WET3_enavg + WET3 if(associated(WCSF_enavg) .and. associated(WCSF)) & WCSF_enavg = WCSF_enavg + WCSF + if(associated(WCSF_ensdv) .and. associated(WCSF)) & + WCSF_ensdv = WCSF_ensdv + WCSF*WCSF if(associated(WCRZ_enavg) .and. associated(WCRZ)) & WCRZ_enavg = WCRZ_enavg + WCRZ + if(associated(WCRZ_ensdv) .and. associated(WCRZ)) & + WCRZ_ensdv = WCRZ_ensdv + WCRZ*WCRZ if(associated(WCPR_enavg) .and. associated(WCPR)) & WCPR_enavg = WCPR_enavg + WCPR + if(associated(WCPR_ensdv) .and. associated(WCPR)) & + WCPR_ensdv = WCPR_ensdv + WCPR*WCPR if(associated(TP1_enavg) .and. associated(TP1)) & TP1_enavg = TP1_enavg + TP1 + if(associated(TP1_ensdv) .and. associated(TP1)) & + TP1_ensdv = TP1_ensdv + TP1*TP1 if(associated(TP2_enavg) .and. associated(TP2)) & TP2_enavg = TP2_enavg + TP2 if(associated(TP3_enavg) .and. associated(TP3)) & @@ -3591,6 +3662,9 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) !if(associated(RZEQ_enavg)) RZEQ_enavg = RZEQ_enavg/NUM_ENSEMBLE if(associated(GHFLX_enavg)) GHFLX_enavg = GHFLX_enavg/NUM_ENSEMBLE if(associated(TPSURF_enavg)) TPSURF_enavg = TPSURF_enavg/NUM_ENSEMBLE + if(associated(TPSURF_ensdv)) TPSURF_ensdv = TPSURF_ensdv/NUM_ENSEMBLE + if(associated(TPSURF_ensdv) .and. associated(TPSURF_enavg)) & + TPSURF_ensdv = sqrt(TPSURF_ensdv - TPSURF_enavg*TPSURF_enavg) if(associated(TPSNOW_enavg)) TPSNOW_enavg = TPSNOW_enavg/NUM_ENSEMBLE if(associated(TPUNST_enavg)) TPUNST_enavg = TPUNST_enavg/NUM_ENSEMBLE if(associated(TPSAT_enavg)) TPSAT_enavg = TPSAT_enavg/NUM_ENSEMBLE @@ -3607,9 +3681,21 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = WET2_enavg/NUM_ENSEMBLE if(associated(WET3_enavg)) WET3_enavg = WET3_enavg/NUM_ENSEMBLE if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE + if(associated(WCSF_ensdv)) WCSF_ensdv = WCSF_ensdv/NUM_ENSEMBLE + if(associated(WCSF_ensdv) .and. associated(WCSF_enavg)) & + WCSF_ensdv = sqrt(WCSF_ensdv -WCSF_enavg*WCSF_enavg) if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE + if(associated(WCRZ_ensdv)) WCRZ_ensdv = WCRZ_ensdv/NUM_ENSEMBLE + if(associated(WCRZ_ensdv) .and. associated(WCRZ_enavg)) & + WCRZ_ensdv = sqrt(WCRZ_ensdv -WCRZ_enavg*WCRZ_enavg) if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE + if(associated(WCPR_ensdv)) WCPR_ensdv = WCPR_ensdv/NUM_ENSEMBLE + if(associated(WCPR_ensdv) .and. associated(WCPR_enavg)) & + WCPR_ensdv = sqrt(WCPR_ensdv -WCPR_enavg*WCPR_enavg) if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 + if(associated(TP1_ensdv)) TP1_ensdv = TP1_ensdv/NUM_ENSEMBLE + if(associated(TP1_ensdv) .and. associated(TP1_enavg)) & + TP1_ensdv = sqrt(TP1_ensdv - TP1_enavg*TP1_enavg) if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index d72d3521..512fb49c 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -777,7 +777,16 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) VERIFY_(STATUS) - + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_surface_analysis ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCSF_ANAsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& LONG_NAME = 'soil_moisture_rootzone_analysis' ,& UNITS = 'm3 m-3' ,& @@ -787,6 +796,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_rootzone_analysis ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCRZ_ANAsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& LONG_NAME = 'soil_moisture_profile_analysis' ,& UNITS = 'm3 m-3' ,& @@ -795,6 +813,15 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_profile_analysis ens.sdv' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCPR_ANAsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& LONG_NAME = 'ave_catchment_temp_incl_snw_analysis' ,& @@ -805,6 +832,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'ave_catchment_temp_incl_snw_analysisi ens.sdv' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TPSURF_ANAsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& LONG_NAME = 'soil_temperatures_layer_1_analysis' ,& UNITS = 'K' ,& @@ -814,6 +850,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_temperatures_layer_1_analysis ens.sdv' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TSOIL1_ANAsdv' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + ! Exports for microwave radiative transfer model (mwRTM) call MAPL_AddExportSpec(GC ,& @@ -1456,7 +1501,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) type(cat_diagS_type), dimension(:), allocatable :: cat_diagS type(cat_diagS_type), dimension(:), allocatable :: cat_diagS_ensavg - + type(cat_diagS_type), dimension(:), allocatable :: cat_diagS_enssdv type(obs_type), dimension(:), pointer :: Observations_l => null() logical :: fresh_incr @@ -1488,6 +1533,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) real, dimension(:),pointer :: PRMC_ana=>null() ! profile soil moisture real, dimension(:),pointer :: TPSURF_ana=>null() ! tpsurf real, dimension(:),pointer :: TSOIL1_ana=>null() ! tsoil1 + real, dimension(:),pointer :: SFMC_anasdv=>null() ! surface soil moisture + real, dimension(:),pointer :: RZMC_anasdv=>null() ! rootzone soil moisture + real, dimension(:),pointer :: PRMC_anasdv=>null() ! profile soil moisture + real, dimension(:),pointer :: TPSURF_anasdv=>null() ! tpsurf + real, dimension(:),pointer :: TSOIL1_anasdv=>null() ! tsoil1 !! export for microwave radiative transfer model (mwRTM) @@ -1655,6 +1705,16 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(status) call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, TPSURF_anasdv, 'TPSURF_ANAsdv' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, TSOIL1_anasdv, 'TSOIL1_ANAsdv' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, SFMC_anasdv, 'WCSF_ANAsdv' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, RZMC_anasdv, 'WCRZ_ANAsdv' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, PRMC_anasdv, 'WCPR_ANAsdv' ,rc=status) + VERIFY_(status) ! exports for microwave radiative transfer model (mwRTM) @@ -1912,9 +1972,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) allocate(cat_progn_tmp( N_catl)) allocate(cat_diagS( N_catl)) allocate(cat_diagS_ensavg(N_catl)) + allocate(cat_diagS_enssdv(N_catl)) do ii=1,N_catl cat_diagS_ensavg(ii) = 0.0 ! initialize ens average + cat_diagS_enssdv(ii) = 0.0 end do do n_e=1,NUM_ENSEMBLE @@ -1930,13 +1992,23 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii) + cat_diagS(ii) end do - + cat_diagS_enssdv(:)%sfmc = cat_diagS_enssdv(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc + cat_diagS_enssdv(:)%rzmc = cat_diagS_enssdv(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc + cat_diagS_enssdv(:)%prmc = cat_diagS_enssdv(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc + cat_diagS_enssdv(:)%tsurf = cat_diagS_enssdv(:)%tsurf + cat_diagS(:)%tsurf*cat_diagS(:)%tsurf + cat_diagS_enssdv(:)%tp(1) = cat_diagS_enssdv(:)%tp(1) + cat_diagS(:)%tp(1) * cat_diagS(:)%tp(1) end do do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize + cat_diagS_enssdv(ii) = cat_diagS_enssdv(ii)/real(NUM_ENSEMBLE) end do + cat_diagS_enssdv(:)%sfmc = sqrt(cat_diagS_enssdv(:)%sfmc - cat_diagS_ensavg(:)%sfmc * cat_diagS_ensavg(:)%sfmc ) + cat_diagS_enssdv(:)%rzmc = sqrt(cat_diagS_enssdv(:)%rzmc - cat_diagS_ensavg(:)%rzmc * cat_diagS_ensavg(:)%rzmc ) + cat_diagS_enssdv(:)%prmc = sqrt(cat_diagS_enssdv(:)%prmc - cat_diagS_ensavg(:)%prmc * cat_diagS_ensavg(:)%prmc ) + cat_diagS_enssdv(:)%tsurf = sqrt(cat_diagS_enssdv(:)%tsurf - cat_diagS_ensavg(:)%tsurf * cat_diagS_ensavg(:)%tsurf ) + cat_diagS_enssdv(:)%tp(1) = sqrt(cat_diagS_enssdv(:)%tp(1) - cat_diagS_ensavg(:)%tp(1) * cat_diagS_ensavg(:)%tp(1)) ! set export variables if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc @@ -1944,12 +2016,18 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K + if(associated(SFMC_anasdv)) SFMC_anasdv(:) = cat_diagS_enssdv(:)%sfmc + if(associated(RZMC_anasdv)) RZMC_anasdv(:) = cat_diagS_enssdv(:)%rzmc + if(associated(PRMC_anasdv)) PRMC_anasdv(:) = cat_diagS_enssdv(:)%prmc + if(associated(TPSURF_anasdv)) TPSURF_anasdv(:) = cat_diagS_enssdv(:)%tsurf + if(associated(TSOIL1_anasdv)) TSOIL1_anasdv(:) = cat_diagS_enssdv(:)%tp(1) if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY deallocate(cat_progn_tmp) deallocate(cat_diagS) - deallocate(cat_diagS_ensavg) + deallocate(cat_diagS_ensavg) + deallocate(cat_diagS_enssdv) ! write analysis fields into SMAP L4_SM aup file ! whenever it was time for assimilation (regardless From 3c360e8bf746de0b1555100dcd9bdcdfe6882c3a Mon Sep 17 00:00:00 2001 From: saraqzhang Date: Wed, 15 Mar 2023 14:26:10 -0400 Subject: [PATCH 02/10] (1) make "std" in variable names for ens. standard deviation, (2) make variance factor (N-1) --- src/Applications/LDAS_App/GEOSldas_HIST.rc | 20 ++-- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 109 ++++++++++-------- .../GEOS_LandAssimGridComp.F90 | 90 ++++++++------- 3 files changed, 116 insertions(+), 103 deletions(-) diff --git a/src/Applications/LDAS_App/GEOSldas_HIST.rc b/src/Applications/LDAS_App/GEOSldas_HIST.rc index b8e5ec40..e0235c44 100644 --- a/src/Applications/LDAS_App/GEOSldas_HIST.rc +++ b/src/Applications/LDAS_App/GEOSldas_HIST.rc @@ -419,21 +419,21 @@ COLLECTIONS: 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , - 'WCSFsdv' , 'ENSAVG' , 'SFMC_FCSTsdv' , - 'WCRZsdv' , 'ENSAVG' , 'RZMC_FCSTsdv' , - 'WCPRsdv' , 'ENSAVG' , 'PRMC_FCSTsdv' , - 'TPSURFsdv' , 'ENSAVG' , 'TSURF_FCSTsdv' , - 'TSOIL1TILEsdv' , 'ENSAVG' , 'TSOIL1_FCSTsdv' , + 'WCSFstd' , 'ENSAVG' , 'SFMC_FCSTstd' , + 'WCRZstd' , 'ENSAVG' , 'RZMC_FCSTstd' , + 'WCPRstd' , 'ENSAVG' , 'PRMC_FCSTstd' , + 'TPSURFstd' , 'ENSAVG' , 'TSURF_FCSTstd' , + 'TSOIL1TILEstd' , 'ENSAVG' , 'TSOIL1_FCSTstd' , 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , - 'WCSF_ANAsdv' , 'LANDASSIM' , 'SFMC_ANAsdv' , - 'WCRZ_ANAsdv' , 'LANDASSIM' , 'RZMC_ANAsdv' , - 'WCPR_ANAsdv' , 'LANDASSIM' , 'PRMC_ANAsdv' , - 'TPSURF_ANAsdv' , 'LANDASSIM' , 'TSURF_ANAsdv' , - 'TSOIL1_ANAsdv' , 'LANDASSIM' , 'TSOIL1_ANAsdv' + 'WCSF_ANAstd' , 'LANDASSIM' , 'SFMC_ANAstd' , + 'WCRZ_ANAstd' , 'LANDASSIM' , 'RZMC_ANAstd' , + 'WCPR_ANAstd' , 'LANDASSIM' , 'PRMC_ANAstd' , + 'TPSURF_ANAstd' , 'LANDASSIM' , 'TSURF_ANAstd' , + 'TSOIL1_ANAstd' , 'LANDASSIM' , 'TSOIL1_ANAstd' :: # For lndfcstana, *.frequency and *.ref_time must be consistent with the LDAS.rc resource diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index e4ada4b7..6b72502b 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -657,9 +657,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'ave_catchment_temp_incl_snw ens.sdv',& + LONG_NAME = 'ave_catchment_temp_incl_snw ensstd',& UNITS = 'K' ,& - SHORT_NAME = 'TPSURFsdv' ,& + SHORT_NAME = 'TPSURFstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -810,9 +810,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_surface_layer ens.sdv' ,& + LONG_NAME = 'water_surface_layer ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCSFsdv' ,& + SHORT_NAME = 'WCSFstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -828,9 +828,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_root_zone ens.sdv' ,& + LONG_NAME = 'water_root_zone ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCRZsdv' ,& + SHORT_NAME = 'WCRZstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -839,7 +839,7 @@ subroutine SetServices(gc, rc) call MAPL_AddExportSpec(GC, & LONG_NAME = 'water_ave_prof' ,& - UNITS = 'm3 m-3' ,& + UNITS = 'm3 m-3' ,& SHORT_NAME = 'WCPR' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& @@ -847,9 +847,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_ave_prof ens.sdv' ,& + LONG_NAME = 'water_ave_prof ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCPRsdv' ,& + SHORT_NAME = 'WCPRstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -865,9 +865,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'soil_temperatures_layer_1 ens.sdv' ,& + LONG_NAME = 'soil_temperatures_layer_1 ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TSOIL1TILEsdv' ,& + SHORT_NAME = 'TSOIL1TILEstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -2333,7 +2333,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: HLATN,HLATN_enavg real, dimension(:),pointer :: QINFIL,QINFIL_enavg real, dimension(:),pointer :: GHFLX,GHFLX_enavg - real, dimension(:),pointer :: TPSURF,TPSURF_enavg,TPSURF_ensdv + real, dimension(:),pointer :: TPSURF,TPSURF_enavg,TPSURF_ensstd real, dimension(:),pointer :: TPSNOW,TPSNOW_enavg real, dimension(:),pointer :: TPUNST,TPUNST_enavg real, dimension(:),pointer :: TPSAT,TPSAT_enavg @@ -2350,10 +2350,10 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: WET1,WET1_enavg real, dimension(:),pointer :: WET2,WET2_enavg real, dimension(:),pointer :: WET3,WET3_enavg - real, dimension(:),pointer :: WCSF,WCSF_enavg,WCSF_ensdv - real, dimension(:),pointer :: WCRZ,WCRZ_enavg,WCRZ_ensdv - real, dimension(:),pointer :: WCPR,WCPR_enavg,WCPR_ensdv - real, dimension(:),pointer :: TP1,TP1_enavg,TP1_ensdv + real, dimension(:),pointer :: WCSF,WCSF_enavg,WCSF_ensstd + real, dimension(:),pointer :: WCRZ,WCRZ_enavg,WCRZ_ensstd + real, dimension(:),pointer :: WCPR,WCPR_enavg,WCPR_ensstd + real, dimension(:),pointer :: TP1,TP1_enavg,TP1_ensstd real, dimension(:),pointer :: TP2,TP2_enavg real, dimension(:),pointer :: TP3,TP3_enavg real, dimension(:),pointer :: TP4,TP4_enavg @@ -2447,6 +2447,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:), pointer :: CNBURN, CNBURN_enavg real, dimension(:), pointer :: CNFSEL, CNFSEL_enavg + real :: nn1 + ! Get my name and setup traceback handle call ESMF_GridCompget(gc, name=comp_name, rc=status) VERIFY_(status) @@ -2905,7 +2907,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, TPSURF_enavg, 'TPSURF' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TPSURF_ensdv, 'TPSURFsdv' ,rc=status) + call MAPL_GetPointer(export, TPSURF_ensstd, 'TPSURFstd' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TPSNOW_enavg, 'TPSNOW' ,rc=status) VERIFY_(status) @@ -2939,19 +2941,19 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, WCSF_enavg, 'WCSF' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCSF_ensdv, 'WCSFsdv' ,rc=status) + call MAPL_GetPointer(export, WCSF_ensstd, 'WCSFstd' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, WCRZ_enavg, 'WCRZ' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCRZ_ensdv, 'WCRZsdv' ,rc=status) + call MAPL_GetPointer(export, WCRZ_ensstd, 'WCRZstd' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, WCPR_enavg, 'WCPR' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCPR_ensdv, 'WCPRsdv' ,rc=status) + call MAPL_GetPointer(export, WCPR_ensstd, 'WCPRstd' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TP1_enavg, 'TSOIL1TILE' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TP1_ensdv, 'TSOIL1TILEsdv' ,rc=status) + call MAPL_GetPointer(export, TP1_ensstd, 'TSOIL1TILEstd' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TP2_enavg, 'TSOIL2TILE' ,rc=status) VERIFY_(status) @@ -3163,7 +3165,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(QINFIL_enavg)) QINFIL_enavg = 0.0 if(associated(GHFLX_enavg)) GHFLX_enavg = 0.0 if(associated(TPSURF_enavg)) TPSURF_enavg = 0.0 - if(associated(TPSURF_ensdv)) TPSURF_ensdv = 0.0 + if(associated(TPSURF_ensstd)) TPSURF_ensstd = 0.0 if(associated(TPSNOW_enavg)) TPSNOW_enavg = 0.0 if(associated(TPUNST_enavg)) TPUNST_enavg = 0.0 if(associated(TPSAT_enavg)) TPSAT_enavg = 0.0 @@ -3180,13 +3182,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = 0.0 if(associated(WET3_enavg)) WET3_enavg = 0.0 if(associated(WCSF_enavg)) WCSF_enavg = 0.0 - if(associated(WCSF_ensdv)) WCSF_ensdv = 0.0 + if(associated(WCSF_ensstd)) WCSF_ensstd = 0.0 if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 - if(associated(WCRZ_ensdv)) WCRZ_ensdv = 0.0 + if(associated(WCRZ_ensstd)) WCRZ_ensstd = 0.0 if(associated(WCPR_enavg)) WCPR_enavg = 0.0 - if(associated(WCPR_ensdv)) WCPR_ensdv = 0.0 + if(associated(WCPR_ensstd)) WCPR_ensstd = 0.0 if(associated(TP1_enavg)) TP1_enavg = 0.0 - if(associated(TP1_ensdv)) TP1_ensdv = 0.0 + if(associated(TP1_ensstd)) TP1_ensstd = 0.0 if(associated(TP2_enavg)) TP2_enavg = 0.0 if(associated(TP3_enavg)) TP3_enavg = 0.0 if(associated(TP4_enavg)) TP4_enavg = 0.0 @@ -3367,8 +3369,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) GHFLX_enavg = GHFLX_enavg + GHFLX if(associated(TPSURF_enavg) .and. associated(TPSURF)) & TPSURF_enavg = TPSURF_enavg + TPSURF - if(associated(TPSURF_ensdv) .and. associated(TPSURF)) & - TPSURF_ensdv = TPSURF_ensdv + TPSURF*TPSURF + if(associated(TPSURF_ensstd) .and. associated(TPSURF)) & + TPSURF_ensstd = TPSURF_ensstd + TPSURF*TPSURF if(associated(TPSNOW_enavg) .and. associated(TPSNOW)) & TPSNOW_enavg = TPSNOW_enavg + TPSNOW if(associated(TPUNST_enavg) .and. associated(TPUNST)) & @@ -3401,20 +3403,20 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) WET3_enavg = WET3_enavg + WET3 if(associated(WCSF_enavg) .and. associated(WCSF)) & WCSF_enavg = WCSF_enavg + WCSF - if(associated(WCSF_ensdv) .and. associated(WCSF)) & - WCSF_ensdv = WCSF_ensdv + WCSF*WCSF + if(associated(WCSF_ensstd) .and. associated(WCSF)) & + WCSF_ensstd = WCSF_ensstd + WCSF*WCSF if(associated(WCRZ_enavg) .and. associated(WCRZ)) & WCRZ_enavg = WCRZ_enavg + WCRZ - if(associated(WCRZ_ensdv) .and. associated(WCRZ)) & - WCRZ_ensdv = WCRZ_ensdv + WCRZ*WCRZ + if(associated(WCRZ_ensstd) .and. associated(WCRZ)) & + WCRZ_ensstd = WCRZ_ensstd + WCRZ*WCRZ if(associated(WCPR_enavg) .and. associated(WCPR)) & WCPR_enavg = WCPR_enavg + WCPR - if(associated(WCPR_ensdv) .and. associated(WCPR)) & - WCPR_ensdv = WCPR_ensdv + WCPR*WCPR + if(associated(WCPR_ensstd) .and. associated(WCPR)) & + WCPR_ensstd = WCPR_ensstd + WCPR*WCPR if(associated(TP1_enavg) .and. associated(TP1)) & TP1_enavg = TP1_enavg + TP1 - if(associated(TP1_ensdv) .and. associated(TP1)) & - TP1_ensdv = TP1_ensdv + TP1*TP1 + if(associated(TP1_ensstd) .and. associated(TP1)) & + TP1_ensstd = TP1_ensstd + TP1*TP1 if(associated(TP2_enavg) .and. associated(TP2)) & TP2_enavg = TP2_enavg + TP2 if(associated(TP3_enavg) .and. associated(TP3)) & @@ -3615,6 +3617,11 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(collect_land_counter == NUM_ENSEMBLE) then + if (NUM_ENSEMBLE >1 ) then + nn1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) + else + nn1 = real(NUM_ENSEMBLE) + endif collect_land_counter = 0 if(associated(TC_enavg)) TC_enavg = TC_enavg/NUM_ENSEMBLE if(associated(QC_enavg)) QC_enavg = QC_enavg/NUM_ENSEMBLE @@ -3662,9 +3669,9 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) !if(associated(RZEQ_enavg)) RZEQ_enavg = RZEQ_enavg/NUM_ENSEMBLE if(associated(GHFLX_enavg)) GHFLX_enavg = GHFLX_enavg/NUM_ENSEMBLE if(associated(TPSURF_enavg)) TPSURF_enavg = TPSURF_enavg/NUM_ENSEMBLE - if(associated(TPSURF_ensdv)) TPSURF_ensdv = TPSURF_ensdv/NUM_ENSEMBLE - if(associated(TPSURF_ensdv) .and. associated(TPSURF_enavg)) & - TPSURF_ensdv = sqrt(TPSURF_ensdv - TPSURF_enavg*TPSURF_enavg) + if(associated(TPSURF_ensstd)) TPSURF_ensstd = TPSURF_ensstd/(NUM_ENSEMBLE-1) + if(associated(TPSURF_ensstd) .and. associated(TPSURF_enavg)) & + TPSURF_ensstd = sqrt(TPSURF_ensstd - nn1*(TPSURF_enavg)**2) if(associated(TPSNOW_enavg)) TPSNOW_enavg = TPSNOW_enavg/NUM_ENSEMBLE if(associated(TPUNST_enavg)) TPUNST_enavg = TPUNST_enavg/NUM_ENSEMBLE if(associated(TPSAT_enavg)) TPSAT_enavg = TPSAT_enavg/NUM_ENSEMBLE @@ -3681,21 +3688,21 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = WET2_enavg/NUM_ENSEMBLE if(associated(WET3_enavg)) WET3_enavg = WET3_enavg/NUM_ENSEMBLE if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE - if(associated(WCSF_ensdv)) WCSF_ensdv = WCSF_ensdv/NUM_ENSEMBLE - if(associated(WCSF_ensdv) .and. associated(WCSF_enavg)) & - WCSF_ensdv = sqrt(WCSF_ensdv -WCSF_enavg*WCSF_enavg) + if(associated(WCSF_ensstd)) WCSF_ensstd = WCSF_ensstd/(NUM_ENSEMBLE-1) + if(associated(WCSF_ensstd) .and. associated(WCSF_enavg)) & + WCSF_ensstd = sqrt(WCSF_ensstd - nn1*(WCSF_enavg)**2) if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE - if(associated(WCRZ_ensdv)) WCRZ_ensdv = WCRZ_ensdv/NUM_ENSEMBLE - if(associated(WCRZ_ensdv) .and. associated(WCRZ_enavg)) & - WCRZ_ensdv = sqrt(WCRZ_ensdv -WCRZ_enavg*WCRZ_enavg) + if(associated(WCRZ_ensstd)) WCRZ_ensstd = WCRZ_ensstd/(NUM_ENSEMBLE-1) + if(associated(WCRZ_ensstd) .and. associated(WCRZ_enavg)) & + WCRZ_ensstd = sqrt(WCRZ_ensstd - nn1*(WCRZ_enavg)**2) if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE - if(associated(WCPR_ensdv)) WCPR_ensdv = WCPR_ensdv/NUM_ENSEMBLE - if(associated(WCPR_ensdv) .and. associated(WCPR_enavg)) & - WCPR_ensdv = sqrt(WCPR_ensdv -WCPR_enavg*WCPR_enavg) + if(associated(WCPR_ensstd)) WCPR_ensstd = WCPR_ensstd/(NUM_ENSEMBLE-1) + if(associated(WCPR_ensstd) .and. associated(WCPR_enavg)) & + WCPR_ensstd = sqrt(WCPR_ensstd - nn1*(WCPR_enavg)**2) if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 - if(associated(TP1_ensdv)) TP1_ensdv = TP1_ensdv/NUM_ENSEMBLE - if(associated(TP1_ensdv) .and. associated(TP1_enavg)) & - TP1_ensdv = sqrt(TP1_ensdv - TP1_enavg*TP1_enavg) + if(associated(TP1_ensstd)) TP1_ensstd = TP1_ensstd/(NUM_ENSEMBLE-1) + if(associated(TP1_ensstd) .and. associated(TP1_enavg)) & + TP1_ensstd = sqrt(TP1_ensstd - nn1*(TP1_enavg)**2) if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 512fb49c..0f143c6a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -779,9 +779,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_surface_analysis ens.sdv' ,& + LONG_NAME = 'soil_moisture_surface_analysis ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCSF_ANAsdv' ,& + SHORT_NAME = 'WCSF_ANAstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -797,9 +797,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_rootzone_analysis ens.sdv' ,& + LONG_NAME = 'soil_moisture_rootzone_analysis ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCRZ_ANAsdv' ,& + SHORT_NAME = 'WCRZ_ANAstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -815,9 +815,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_profile_analysis ens.sdv' ,& + LONG_NAME = 'soil_moisture_profile_analysis ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCPR_ANAsdv' ,& + SHORT_NAME = 'WCPR_ANAstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -833,9 +833,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'ave_catchment_temp_incl_snw_analysisi ens.sdv' ,& + LONG_NAME = 'ave_catchment_temp_incl_snw_analysisi ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TPSURF_ANAsdv' ,& + SHORT_NAME = 'TPSURF_ANAstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -851,9 +851,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_temperatures_layer_1_analysis ens.sdv' ,& + LONG_NAME = 'soil_temperatures_layer_1_analysis ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TSOIL1_ANAsdv' ,& + SHORT_NAME = 'TSOIL1_ANAstd' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -1478,7 +1478,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) type(date_time_type) :: date_time_new character(len=14) :: datestamp - integer :: N_catl, N_catg,N_obsl_max, n_e, ii + integer :: N_catl, N_catg,N_obsl_max, n_e, ii character(len=300) :: out_path character(len=ESMF_MAXSTR) :: exp_id @@ -1501,7 +1501,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) type(cat_diagS_type), dimension(:), allocatable :: cat_diagS type(cat_diagS_type), dimension(:), allocatable :: cat_diagS_ensavg - type(cat_diagS_type), dimension(:), allocatable :: cat_diagS_enssdv + type(cat_diagS_type), dimension(:), allocatable :: cat_diagS_ensstd type(obs_type), dimension(:), pointer :: Observations_l => null() logical :: fresh_incr @@ -1533,11 +1533,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) real, dimension(:),pointer :: PRMC_ana=>null() ! profile soil moisture real, dimension(:),pointer :: TPSURF_ana=>null() ! tpsurf real, dimension(:),pointer :: TSOIL1_ana=>null() ! tsoil1 - real, dimension(:),pointer :: SFMC_anasdv=>null() ! surface soil moisture - real, dimension(:),pointer :: RZMC_anasdv=>null() ! rootzone soil moisture - real, dimension(:),pointer :: PRMC_anasdv=>null() ! profile soil moisture - real, dimension(:),pointer :: TPSURF_anasdv=>null() ! tpsurf - real, dimension(:),pointer :: TSOIL1_anasdv=>null() ! tsoil1 + real, dimension(:),pointer :: SFMC_anastd=>null() ! surface soil moisture + real, dimension(:),pointer :: RZMC_anastd=>null() ! rootzone soil moisture + real, dimension(:),pointer :: PRMC_anastd=>null() ! profile soil moisture + real, dimension(:),pointer :: TPSURF_anastd=>null() ! tpsurf + real, dimension(:),pointer :: TSOIL1_anastd=>null() ! tsoil1 !! export for microwave radiative transfer model (mwRTM) @@ -1551,6 +1551,8 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) character(len=ESMF_MAXSTR) :: ensid_string integer :: ens, nymd, nhms, ens_id_width integer :: LandassimDTstep + real :: nn1 + #ifdef DBG_LANDASSIM_INPUTS ! vars for debugging purposes type(ESMF_Grid) :: TILEGRID @@ -1705,15 +1707,15 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) VERIFY_(status) call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TPSURF_anasdv, 'TPSURF_ANAsdv' ,rc=status) + call MAPL_GetPointer(export, TPSURF_anastd, 'TPSURF_ANAstd' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TSOIL1_anasdv, 'TSOIL1_ANAsdv' ,rc=status) + call MAPL_GetPointer(export, TSOIL1_anastd, 'TSOIL1_ANAstd' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, SFMC_anasdv, 'WCSF_ANAsdv' ,rc=status) + call MAPL_GetPointer(export, SFMC_anastd, 'WCSF_ANAstd' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, RZMC_anasdv, 'WCRZ_ANAsdv' ,rc=status) + call MAPL_GetPointer(export, RZMC_anastd, 'WCRZ_ANAstd' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, PRMC_anasdv, 'WCPR_ANAsdv' ,rc=status) + call MAPL_GetPointer(export, PRMC_anastd, 'WCPR_ANAstd' ,rc=status) VERIFY_(status) ! exports for microwave radiative transfer model (mwRTM) @@ -1972,11 +1974,11 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) allocate(cat_progn_tmp( N_catl)) allocate(cat_diagS( N_catl)) allocate(cat_diagS_ensavg(N_catl)) - allocate(cat_diagS_enssdv(N_catl)) + allocate(cat_diagS_ensstd(N_catl)) do ii=1,N_catl cat_diagS_ensavg(ii) = 0.0 ! initialize ens average - cat_diagS_enssdv(ii) = 0.0 + cat_diagS_ensstd(ii) = 0.0 end do do n_e=1,NUM_ENSEMBLE @@ -1992,23 +1994,27 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii) + cat_diagS(ii) end do - cat_diagS_enssdv(:)%sfmc = cat_diagS_enssdv(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc - cat_diagS_enssdv(:)%rzmc = cat_diagS_enssdv(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc - cat_diagS_enssdv(:)%prmc = cat_diagS_enssdv(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc - cat_diagS_enssdv(:)%tsurf = cat_diagS_enssdv(:)%tsurf + cat_diagS(:)%tsurf*cat_diagS(:)%tsurf - cat_diagS_enssdv(:)%tp(1) = cat_diagS_enssdv(:)%tp(1) + cat_diagS(:)%tp(1) * cat_diagS(:)%tp(1) + cat_diagS_ensstd(:)%sfmc = cat_diagS_ensstd(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc + cat_diagS_ensstd(:)%rzmc = cat_diagS_ensstd(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc + cat_diagS_ensstd(:)%prmc = cat_diagS_ensstd(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc + cat_diagS_ensstd(:)%tsurf = cat_diagS_ensstd(:)%tsurf + cat_diagS(:)%tsurf*cat_diagS(:)%tsurf + cat_diagS_ensstd(:)%tp(1) = cat_diagS_ensstd(:)%tp(1) + cat_diagS(:)%tp(1) * cat_diagS(:)%tp(1) end do do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize - cat_diagS_enssdv(ii) = cat_diagS_enssdv(ii)/real(NUM_ENSEMBLE) + cat_diagS_ensstd(ii) = cat_diagS_ensstd(ii)/real(NUM_ENSEMBLE-1) end do - - cat_diagS_enssdv(:)%sfmc = sqrt(cat_diagS_enssdv(:)%sfmc - cat_diagS_ensavg(:)%sfmc * cat_diagS_ensavg(:)%sfmc ) - cat_diagS_enssdv(:)%rzmc = sqrt(cat_diagS_enssdv(:)%rzmc - cat_diagS_ensavg(:)%rzmc * cat_diagS_ensavg(:)%rzmc ) - cat_diagS_enssdv(:)%prmc = sqrt(cat_diagS_enssdv(:)%prmc - cat_diagS_ensavg(:)%prmc * cat_diagS_ensavg(:)%prmc ) - cat_diagS_enssdv(:)%tsurf = sqrt(cat_diagS_enssdv(:)%tsurf - cat_diagS_ensavg(:)%tsurf * cat_diagS_ensavg(:)%tsurf ) - cat_diagS_enssdv(:)%tp(1) = sqrt(cat_diagS_enssdv(:)%tp(1) - cat_diagS_ensavg(:)%tp(1) * cat_diagS_ensavg(:)%tp(1)) + if (NUM_ENSEMBLE > 1 ) then + nn1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) + else + nn1 = real(NUM_ENSEMBLE) + endif + cat_diagS_ensstd(:)%sfmc = sqrt(cat_diagS_ensstd(:)%sfmc - nn1*(cat_diagS_ensavg(:)%sfmc)**2 ) + cat_diagS_ensstd(:)%rzmc = sqrt(cat_diagS_ensstd(:)%rzmc - nn1*(cat_diagS_ensavg(:)%rzmc)**2 ) + cat_diagS_ensstd(:)%prmc = sqrt(cat_diagS_ensstd(:)%prmc - nn1*(cat_diagS_ensavg(:)%prmc)**2 ) + cat_diagS_ensstd(:)%tsurf = sqrt(cat_diagS_ensstd(:)%tsurf - nn1*(cat_diagS_ensavg(:)%tsurf)**2 ) + cat_diagS_ensstd(:)%tp(1) = sqrt(cat_diagS_ensstd(:)%tp(1) - nn1*(cat_diagS_ensavg(:)%tp(1))**2 ) ! set export variables if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc @@ -2016,18 +2022,18 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K - if(associated(SFMC_anasdv)) SFMC_anasdv(:) = cat_diagS_enssdv(:)%sfmc - if(associated(RZMC_anasdv)) RZMC_anasdv(:) = cat_diagS_enssdv(:)%rzmc - if(associated(PRMC_anasdv)) PRMC_anasdv(:) = cat_diagS_enssdv(:)%prmc - if(associated(TPSURF_anasdv)) TPSURF_anasdv(:) = cat_diagS_enssdv(:)%tsurf - if(associated(TSOIL1_anasdv)) TSOIL1_anasdv(:) = cat_diagS_enssdv(:)%tp(1) + if(associated(SFMC_anastd)) SFMC_anastd(:) = cat_diagS_ensstd(:)%sfmc + if(associated(RZMC_anastd)) RZMC_anastd(:) = cat_diagS_ensstd(:)%rzmc + if(associated(PRMC_anastd)) PRMC_anastd(:) = cat_diagS_ensstd(:)%prmc + if(associated(TPSURF_anastd)) TPSURF_anastd(:) = cat_diagS_ensstd(:)%tsurf + if(associated(TSOIL1_anastd)) TSOIL1_anastd(:) = cat_diagS_ensstd(:)%tp(1) if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY deallocate(cat_progn_tmp) deallocate(cat_diagS) deallocate(cat_diagS_ensavg) - deallocate(cat_diagS_enssdv) + deallocate(cat_diagS_ensstd) ! write analysis fields into SMAP L4_SM aup file ! whenever it was time for assimilation (regardless From 5d085077ee3717f8703038a4464a1a3fe27c5e64 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 16 Mar 2023 09:44:53 -0400 Subject: [PATCH 03/10] additional changes in variables names for ENSSTD output --- src/Applications/LDAS_App/GEOSldas_HIST.rc | 70 +++++----- src/Applications/LDAS_App/tile_bin2nc4.F90 | 33 +++-- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 125 +++++++++--------- .../GEOS_LandAssimGridComp.F90 | 109 +++++++-------- 4 files changed, 183 insertions(+), 154 deletions(-) diff --git a/src/Applications/LDAS_App/GEOSldas_HIST.rc b/src/Applications/LDAS_App/GEOSldas_HIST.rc index e0235c44..ba9fb5c4 100644 --- a/src/Applications/LDAS_App/GEOSldas_HIST.rc +++ b/src/Applications/LDAS_App/GEOSldas_HIST.rc @@ -414,26 +414,26 @@ COLLECTIONS: inst3_1d_lndfcstana_Nt.mode: 'instantaneous', inst3_1d_lndfcstana_Nt.frequency: 030000, inst3_1d_lndfcstana_Nt.ref_time: 000000, - inst3_1d_lndfcstana_Nt.fields: 'WCSF' , 'ENSAVG' , 'SFMC_FCST' , - 'WCRZ' , 'ENSAVG' , 'RZMC_FCST' , - 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , - 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , - 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , - 'WCSFstd' , 'ENSAVG' , 'SFMC_FCSTstd' , - 'WCRZstd' , 'ENSAVG' , 'RZMC_FCSTstd' , - 'WCPRstd' , 'ENSAVG' , 'PRMC_FCSTstd' , - 'TPSURFstd' , 'ENSAVG' , 'TSURF_FCSTstd' , - 'TSOIL1TILEstd' , 'ENSAVG' , 'TSOIL1_FCSTstd' , - 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , - 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , - 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , - 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , - 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , - 'WCSF_ANAstd' , 'LANDASSIM' , 'SFMC_ANAstd' , - 'WCRZ_ANAstd' , 'LANDASSIM' , 'RZMC_ANAstd' , - 'WCPR_ANAstd' , 'LANDASSIM' , 'PRMC_ANAstd' , - 'TPSURF_ANAstd' , 'LANDASSIM' , 'TSURF_ANAstd' , - 'TSOIL1_ANAstd' , 'LANDASSIM' , 'TSOIL1_ANAstd' + inst3_1d_lndfcstana_Nt.fields: 'WCSF' , 'ENSAVG' , 'SFMC_FCST' , + 'WCRZ' , 'ENSAVG' , 'RZMC_FCST' , + 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , + 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , + 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , + 'WCSF_ENSSTD' , 'ENSAVG' , 'SFMC_FCST_ENSSTD' , + 'WCRZ_ENSSTD' , 'ENSAVG' , 'RZMC_FCST_ENSSTD' , + 'WCPR_ENSSTD' , 'ENSAVG' , 'PRMC_FCST_ENSSTD' , + 'TPSURF_ENSSTD' , 'ENSAVG' , 'TSURF_FCST_ENSSTD' , + 'TSOIL1TILE_ENSSTD' , 'ENSAVG' , 'TSOIL1_FCST_ENSSTD' , + 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , + 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , + 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , + 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , + 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , + 'WCSF_ANA_ENSSTD' , 'LANDASSIM' , 'SFMC_ANA_ENSSTD' , + 'WCRZ_ANA_ENSSTD' , 'LANDASSIM' , 'RZMC_ANA_ENSSTD' , + 'WCPR_ANA_ENSSTD' , 'LANDASSIM' , 'PRMC_ANA_ENSSTD' , + 'TPSURF_ANA_ENSSTD' , 'LANDASSIM' , 'TSURF_ANA_ENSSTD' , + 'TSOIL1_ANA_ENSSTD' , 'LANDASSIM' , 'TSOIL1_ANA_ENSSTD' :: # For lndfcstana, *.frequency and *.ref_time must be consistent with the LDAS.rc resource @@ -451,16 +451,26 @@ COLLECTIONS: inst3_2d_lndfcstana_Nx.regrid_name: 'GRIDNAME', inst3_2d_lndfcstana_Nx.grid_label: PC720x361-DC, inst3_2d_lndfcstana_Nx.deflate: 2, - inst3_2d_lndfcstana_Nx.fields: 'WCSF' , 'ENSAVG' , 'SFMC_FCST' , - 'WCRZ' , 'ENSAVG' , 'RZMC_FCST' , - 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , - 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , - 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , - 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , - 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , - 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , - 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , - 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , + inst3_2d_lndfcstana_Nx.fields: 'WCSF' , 'ENSAVG' , 'SFMC_FCST' , + 'WCRZ' , 'ENSAVG' , 'RZMC_FCST' , + 'WCPR' , 'ENSAVG' , 'PRMC_FCST' , + 'TPSURF' , 'ENSAVG' , 'TSURF_FCST' , + 'TSOIL1TILE' , 'ENSAVG' , 'TSOIL1_FCST' , + 'WCSF_ENSSTD' , 'ENSAVG' , 'SFMC_FCST_ENSSTD' , + 'WCRZ_ENSSTD' , 'ENSAVG' , 'RZMC_FCST_ENSSTD' , + 'WCPR_ENSSTD' , 'ENSAVG' , 'PRMC_FCST_ENSSTD' , + 'TPSURF_ENSSTD' , 'ENSAVG' , 'TSURF_FCST_ENSSTD' , + 'TSOIL1TILE_ENSSTD' , 'ENSAVG' , 'TSOIL1_FCST_ENSSTD' , + 'WCSF_ANA' , 'LANDASSIM' , 'SFMC_ANA' , + 'WCRZ_ANA' , 'LANDASSIM' , 'RZMC_ANA' , + 'WCPR_ANA' , 'LANDASSIM' , 'PRMC_ANA' , + 'TPSURF_ANA' , 'LANDASSIM' , 'TSURF_ANA' , + 'TSOIL1_ANA' , 'LANDASSIM' , 'TSOIL1_ANA' , + 'WCSF_ANA_ENSSTD' , 'LANDASSIM' , 'SFMC_ANA_ENSSTD' , + 'WCRZ_ANA_ENSSTD' , 'LANDASSIM' , 'RZMC_ANA_ENSSTD' , + 'WCPR_ANA_ENSSTD' , 'LANDASSIM' , 'PRMC_ANA_ENSSTD' , + 'TPSURF_ANA_ENSSTD' , 'LANDASSIM' , 'TSURF_ANA_ENSSTD' , + 'TSOIL1_ANA_ENSSTD' , 'LANDASSIM' , 'TSOIL1_ANA_ENSSTD' :: # ========================== EOF ============================================================== diff --git a/src/Applications/LDAS_App/tile_bin2nc4.F90 b/src/Applications/LDAS_App/tile_bin2nc4.F90 index a9d6c6a1..e50a9eb5 100644 --- a/src/Applications/LDAS_App/tile_bin2nc4.F90 +++ b/src/Applications/LDAS_App/tile_bin2nc4.F90 @@ -399,16 +399,29 @@ FUNCTION getAttribute (SHORT_NAME, LNAME, UNT) result (str_atr) ! land assimilation forecast and analysis for Catchment model diagnostics - case ('SFMC_FCST'); LONG_NAME = 'soil_moisture_surface_forecast'; UNITS = 'm3 m-3' - case ('RZMC_FCST'); LONG_NAME = 'soil_moisture_rootzone_forecast'; UNITS = 'm3 m-3' - case ('PRMC_FCST'); LONG_NAME = 'soil_moisture_profile_forecast'; UNITS = 'm3 m-3' - case ('TSURF_FCST'); LONG_NAME = 'ave_catchment_temp_incl_snw_forecast'; UNITS = 'K' - case ('TSOIL1_FCST'); LONG_NAME = 'soil_temperatures_layer_1_forecast'; UNITS = 'K' - case ('SFMC_ANA'); LONG_NAME = 'soil_moisture_surface_analysis'; UNITS = 'm3 m-3' - case ('RZMC_ANA'); LONG_NAME = 'soil_moisture_rootzone_analysis'; UNITS = 'm3 m-3' - case ('PRMC_ANA'); LONG_NAME = 'soil_moisture_profile_analysis'; UNITS = 'm3 m-3' - case ('TSURF_ANA'); LONG_NAME = 'ave_catchment_temp_incl_snw_analysis'; UNITS = 'K' - case ('TSOIL1_ANA'); LONG_NAME = 'soil_temperatures_layer_1_analysis'; UNITS = 'K' + case ('SFMC_FCST'); LONG_NAME = 'soil_moisture_surface_forecast'; UNITS = 'm3 m-3' + case ('RZMC_FCST'); LONG_NAME = 'soil_moisture_rootzone_forecast'; UNITS = 'm3 m-3' + case ('PRMC_FCST'); LONG_NAME = 'soil_moisture_profile_forecast'; UNITS = 'm3 m-3' + case ('TSURF_FCST'); LONG_NAME = 'ave_catchment_temp_incl_snw_forecast'; UNITS = 'K' + case ('TSOIL1_FCST'); LONG_NAME = 'soil_temperatures_layer_1_forecast'; UNITS = 'K' + + case ('SFMC_FCST_ENSSTD'); LONG_NAME = 'soil_moisture_surface_forecast_ensstd'; UNITS = 'm3 m-3' + case ('RZMC_FCST_ENSSTD'); LONG_NAME = 'soil_moisture_rootzone_forecast_ensstd'; UNITS = 'm3 m-3' + case ('PRMC_FCST_ENSSTD'); LONG_NAME = 'soil_moisture_profile_forecast_ensstd'; UNITS = 'm3 m-3' + case ('TSURF_FCST_ENSSTD'); LONG_NAME = 'ave_catchment_temp_incl_snw_forecast_ensstd'; UNITS = 'K' + case ('TSOIL1_FCST_ENSSTD'); LONG_NAME = 'soil_temperatures_layer_1_forecast_ensstd'; UNITS = 'K' + + case ('SFMC_ANA'); LONG_NAME = 'soil_moisture_surface_analysis'; UNITS = 'm3 m-3' + case ('RZMC_ANA'); LONG_NAME = 'soil_moisture_rootzone_analysis'; UNITS = 'm3 m-3' + case ('PRMC_ANA'); LONG_NAME = 'soil_moisture_profile_analysis'; UNITS = 'm3 m-3' + case ('TSURF_ANA'); LONG_NAME = 'ave_catchment_temp_incl_snw_analysis'; UNITS = 'K' + case ('TSOIL1_ANA'); LONG_NAME = 'soil_temperatures_layer_1_analysis'; UNITS = 'K' + + case ('SFMC_ANA_ENSSTD'); LONG_NAME = 'soil_moisture_surface_analysis_ensstd'; UNITS = 'm3 m-3' + case ('RZMC_ANA_ENSSTD'); LONG_NAME = 'soil_moisture_rootzone_analysis_ensstd'; UNITS = 'm3 m-3' + case ('PRMC_ANA_ENSSTD'); LONG_NAME = 'soil_moisture_profile_analysis_ensstd'; UNITS = 'm3 m-3' + case ('TSURF_ANA_ENSSTD'); LONG_NAME = 'ave_catchment_temp_incl_snw_analysis_ensstd'; UNITS = 'K' + case ('TSOIL1_ANA_ENSSTD'); LONG_NAME = 'soil_temperatures_layer_1_analysis_ensstd'; UNITS = 'K' ! other land assimilation fields diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index 6b72502b..833247df 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -21,7 +21,8 @@ module GEOS_EnsGridCompMod public :: catch_progn public :: catch_param - ! !DESCRIPTION: This GridComp collect ensemble member and then averages the vaiables form catchment + ! !DESCRIPTION: This GridComp collects ensemble members and then averages the variables from Catchment. + ! For select variables, the ensemble standard deviation is also computed. !EOP integer :: NUM_ENSEMBLE @@ -126,7 +127,7 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) -!! export for ens average +!! exports for ens average (and, for a few variables, the ens std) call MAPL_AddExportSpec(GC ,& LONG_NAME = 'canopy_temperature' ,& @@ -439,14 +440,14 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'sublimation' ,& - UNITS = 'kg m-2 s-1' ,& - SHORT_NAME = 'SUBLIM' ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'sublimation' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'SUBLIM' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) - VERIFY_(STATUS) + VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & LONG_NAME = 'upward_sensible_heat_flux' ,& @@ -657,9 +658,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'ave_catchment_temp_incl_snw ensstd',& + LONG_NAME = 'ave_catchment_temp_incl_snw_ensstd',& UNITS = 'K' ,& - SHORT_NAME = 'TPSURFstd' ,& + SHORT_NAME = 'TPSURF_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -810,9 +811,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_surface_layer ensstd' ,& + LONG_NAME = 'water_surface_layer_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCSFstd' ,& + SHORT_NAME = 'WCSF_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -828,9 +829,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_root_zone ensstd' ,& + LONG_NAME = 'water_root_zone_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCRZstd' ,& + SHORT_NAME = 'WCRZ_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -847,9 +848,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'water_ave_prof ensstd' ,& + LONG_NAME = 'water_ave_prof_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCPRstd' ,& + SHORT_NAME = 'WCPR_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -865,9 +866,9 @@ subroutine SetServices(gc, rc) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & - LONG_NAME = 'soil_temperatures_layer_1 ensstd' ,& + LONG_NAME = 'soil_temperatures_layer_1_ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TSOIL1TILEstd' ,& + SHORT_NAME = 'TSOIL1TILE_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -2333,7 +2334,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: HLATN,HLATN_enavg real, dimension(:),pointer :: QINFIL,QINFIL_enavg real, dimension(:),pointer :: GHFLX,GHFLX_enavg - real, dimension(:),pointer :: TPSURF,TPSURF_enavg,TPSURF_ensstd + real, dimension(:),pointer :: TPSURF,TPSURF_enavg,TPSURF_enstd real, dimension(:),pointer :: TPSNOW,TPSNOW_enavg real, dimension(:),pointer :: TPUNST,TPUNST_enavg real, dimension(:),pointer :: TPSAT,TPSAT_enavg @@ -2350,10 +2351,10 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:),pointer :: WET1,WET1_enavg real, dimension(:),pointer :: WET2,WET2_enavg real, dimension(:),pointer :: WET3,WET3_enavg - real, dimension(:),pointer :: WCSF,WCSF_enavg,WCSF_ensstd - real, dimension(:),pointer :: WCRZ,WCRZ_enavg,WCRZ_ensstd - real, dimension(:),pointer :: WCPR,WCPR_enavg,WCPR_ensstd - real, dimension(:),pointer :: TP1,TP1_enavg,TP1_ensstd + real, dimension(:),pointer :: WCSF,WCSF_enavg,WCSF_enstd + real, dimension(:),pointer :: WCRZ,WCRZ_enavg,WCRZ_enstd + real, dimension(:),pointer :: WCPR,WCPR_enavg,WCPR_enstd + real, dimension(:),pointer :: TP1,TP1_enavg,TP1_enstd real, dimension(:),pointer :: TP2,TP2_enavg real, dimension(:),pointer :: TP3,TP3_enavg real, dimension(:),pointer :: TP4,TP4_enavg @@ -2447,7 +2448,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:), pointer :: CNBURN, CNBURN_enavg real, dimension(:), pointer :: CNFSEL, CNFSEL_enavg - real :: nn1 + real :: NdivNm1 ! Get my name and setup traceback handle call ESMF_GridCompget(gc, name=comp_name, rc=status) @@ -2907,7 +2908,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, TPSURF_enavg, 'TPSURF' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TPSURF_ensstd, 'TPSURFstd' ,rc=status) + call MAPL_GetPointer(export, TPSURF_enstd, 'TPSURF_ENSSTD' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TPSNOW_enavg, 'TPSNOW' ,rc=status) VERIFY_(status) @@ -2941,19 +2942,19 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) VERIFY_(status) call MAPL_GetPointer(export, WCSF_enavg, 'WCSF' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCSF_ensstd, 'WCSFstd' ,rc=status) + call MAPL_GetPointer(export, WCSF_enstd, 'WCSF_ENSSTD' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, WCRZ_enavg, 'WCRZ' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCRZ_ensstd, 'WCRZstd' ,rc=status) + call MAPL_GetPointer(export, WCRZ_enstd, 'WCRZ_ENSSTD' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, WCPR_enavg, 'WCPR' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, WCPR_ensstd, 'WCPRstd' ,rc=status) + call MAPL_GetPointer(export, WCPR_enstd, 'WCPR_ENSSTD' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TP1_enavg, 'TSOIL1TILE' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TP1_ensstd, 'TSOIL1TILEstd' ,rc=status) + call MAPL_GetPointer(export, TP1_enstd, 'TSOIL1TILE_ENSSTD' ,rc=status) VERIFY_(status) call MAPL_GetPointer(export, TP2_enavg, 'TSOIL2TILE' ,rc=status) VERIFY_(status) @@ -3165,7 +3166,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(QINFIL_enavg)) QINFIL_enavg = 0.0 if(associated(GHFLX_enavg)) GHFLX_enavg = 0.0 if(associated(TPSURF_enavg)) TPSURF_enavg = 0.0 - if(associated(TPSURF_ensstd)) TPSURF_ensstd = 0.0 + if(associated(TPSURF_enstd)) TPSURF_enstd = 0.0 if(associated(TPSNOW_enavg)) TPSNOW_enavg = 0.0 if(associated(TPUNST_enavg)) TPUNST_enavg = 0.0 if(associated(TPSAT_enavg)) TPSAT_enavg = 0.0 @@ -3182,13 +3183,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = 0.0 if(associated(WET3_enavg)) WET3_enavg = 0.0 if(associated(WCSF_enavg)) WCSF_enavg = 0.0 - if(associated(WCSF_ensstd)) WCSF_ensstd = 0.0 + if(associated(WCSF_enstd)) WCSF_enstd = 0.0 if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 - if(associated(WCRZ_ensstd)) WCRZ_ensstd = 0.0 + if(associated(WCRZ_enstd)) WCRZ_enstd = 0.0 if(associated(WCPR_enavg)) WCPR_enavg = 0.0 - if(associated(WCPR_ensstd)) WCPR_ensstd = 0.0 + if(associated(WCPR_enstd)) WCPR_enstd = 0.0 if(associated(TP1_enavg)) TP1_enavg = 0.0 - if(associated(TP1_ensstd)) TP1_ensstd = 0.0 + if(associated(TP1_enstd)) TP1_enstd = 0.0 if(associated(TP2_enavg)) TP2_enavg = 0.0 if(associated(TP3_enavg)) TP3_enavg = 0.0 if(associated(TP4_enavg)) TP4_enavg = 0.0 @@ -3369,8 +3370,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) GHFLX_enavg = GHFLX_enavg + GHFLX if(associated(TPSURF_enavg) .and. associated(TPSURF)) & TPSURF_enavg = TPSURF_enavg + TPSURF - if(associated(TPSURF_ensstd) .and. associated(TPSURF)) & - TPSURF_ensstd = TPSURF_ensstd + TPSURF*TPSURF + if(associated(TPSURF_enstd) .and. associated(TPSURF)) & + TPSURF_enstd = TPSURF_enstd + TPSURF*TPSURF if(associated(TPSNOW_enavg) .and. associated(TPSNOW)) & TPSNOW_enavg = TPSNOW_enavg + TPSNOW if(associated(TPUNST_enavg) .and. associated(TPUNST)) & @@ -3403,20 +3404,20 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) WET3_enavg = WET3_enavg + WET3 if(associated(WCSF_enavg) .and. associated(WCSF)) & WCSF_enavg = WCSF_enavg + WCSF - if(associated(WCSF_ensstd) .and. associated(WCSF)) & - WCSF_ensstd = WCSF_ensstd + WCSF*WCSF + if(associated(WCSF_enstd) .and. associated(WCSF)) & + WCSF_enstd = WCSF_enstd + WCSF*WCSF if(associated(WCRZ_enavg) .and. associated(WCRZ)) & WCRZ_enavg = WCRZ_enavg + WCRZ - if(associated(WCRZ_ensstd) .and. associated(WCRZ)) & - WCRZ_ensstd = WCRZ_ensstd + WCRZ*WCRZ + if(associated(WCRZ_enstd) .and. associated(WCRZ)) & + WCRZ_enstd = WCRZ_enstd + WCRZ*WCRZ if(associated(WCPR_enavg) .and. associated(WCPR)) & WCPR_enavg = WCPR_enavg + WCPR - if(associated(WCPR_ensstd) .and. associated(WCPR)) & - WCPR_ensstd = WCPR_ensstd + WCPR*WCPR + if(associated(WCPR_enstd) .and. associated(WCPR)) & + WCPR_enstd = WCPR_enstd + WCPR*WCPR if(associated(TP1_enavg) .and. associated(TP1)) & TP1_enavg = TP1_enavg + TP1 - if(associated(TP1_ensstd) .and. associated(TP1)) & - TP1_ensstd = TP1_ensstd + TP1*TP1 + if(associated(TP1_enstd) .and. associated(TP1)) & + TP1_enstd = TP1_enstd + TP1*TP1 if(associated(TP2_enavg) .and. associated(TP2)) & TP2_enavg = TP2_enavg + TP2 if(associated(TP3_enavg) .and. associated(TP3)) & @@ -3618,9 +3619,9 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(collect_land_counter == NUM_ENSEMBLE) then if (NUM_ENSEMBLE >1 ) then - nn1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) + NdivNm1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) else - nn1 = real(NUM_ENSEMBLE) + NdivNm1 = real(NUM_ENSEMBLE) endif collect_land_counter = 0 if(associated(TC_enavg)) TC_enavg = TC_enavg/NUM_ENSEMBLE @@ -3669,9 +3670,9 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) !if(associated(RZEQ_enavg)) RZEQ_enavg = RZEQ_enavg/NUM_ENSEMBLE if(associated(GHFLX_enavg)) GHFLX_enavg = GHFLX_enavg/NUM_ENSEMBLE if(associated(TPSURF_enavg)) TPSURF_enavg = TPSURF_enavg/NUM_ENSEMBLE - if(associated(TPSURF_ensstd)) TPSURF_ensstd = TPSURF_ensstd/(NUM_ENSEMBLE-1) - if(associated(TPSURF_ensstd) .and. associated(TPSURF_enavg)) & - TPSURF_ensstd = sqrt(TPSURF_ensstd - nn1*(TPSURF_enavg)**2) + if(associated(TPSURF_enstd)) TPSURF_enstd = TPSURF_enstd/(NUM_ENSEMBLE-1) + if(associated(TPSURF_enstd) .and. associated(TPSURF_enavg)) & + TPSURF_enstd = sqrt(TPSURF_enstd - NdivNm1*(TPSURF_enavg)**2) if(associated(TPSNOW_enavg)) TPSNOW_enavg = TPSNOW_enavg/NUM_ENSEMBLE if(associated(TPUNST_enavg)) TPUNST_enavg = TPUNST_enavg/NUM_ENSEMBLE if(associated(TPSAT_enavg)) TPSAT_enavg = TPSAT_enavg/NUM_ENSEMBLE @@ -3688,21 +3689,21 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = WET2_enavg/NUM_ENSEMBLE if(associated(WET3_enavg)) WET3_enavg = WET3_enavg/NUM_ENSEMBLE if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE - if(associated(WCSF_ensstd)) WCSF_ensstd = WCSF_ensstd/(NUM_ENSEMBLE-1) - if(associated(WCSF_ensstd) .and. associated(WCSF_enavg)) & - WCSF_ensstd = sqrt(WCSF_ensstd - nn1*(WCSF_enavg)**2) + if(associated(WCSF_enstd)) WCSF_enstd = WCSF_enstd/(NUM_ENSEMBLE-1) + if(associated(WCSF_enstd) .and. associated(WCSF_enavg)) & + WCSF_enstd = sqrt(WCSF_enstd - NdivNm1*(WCSF_enavg)**2) if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE - if(associated(WCRZ_ensstd)) WCRZ_ensstd = WCRZ_ensstd/(NUM_ENSEMBLE-1) - if(associated(WCRZ_ensstd) .and. associated(WCRZ_enavg)) & - WCRZ_ensstd = sqrt(WCRZ_ensstd - nn1*(WCRZ_enavg)**2) + if(associated(WCRZ_enstd)) WCRZ_enstd = WCRZ_enstd/(NUM_ENSEMBLE-1) + if(associated(WCRZ_enstd) .and. associated(WCRZ_enavg)) & + WCRZ_enstd = sqrt(WCRZ_enstd - NdivNm1*(WCRZ_enavg)**2) if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE - if(associated(WCPR_ensstd)) WCPR_ensstd = WCPR_ensstd/(NUM_ENSEMBLE-1) - if(associated(WCPR_ensstd) .and. associated(WCPR_enavg)) & - WCPR_ensstd = sqrt(WCPR_ensstd - nn1*(WCPR_enavg)**2) + if(associated(WCPR_enstd)) WCPR_enstd = WCPR_enstd/(NUM_ENSEMBLE-1) + if(associated(WCPR_enstd) .and. associated(WCPR_enavg)) & + WCPR_enstd = sqrt(WCPR_enstd - NdivNm1*(WCPR_enavg)**2) if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 - if(associated(TP1_ensstd)) TP1_ensstd = TP1_ensstd/(NUM_ENSEMBLE-1) - if(associated(TP1_ensstd) .and. associated(TP1_enavg)) & - TP1_ensstd = sqrt(TP1_ensstd - nn1*(TP1_enavg)**2) + if(associated(TP1_enstd)) TP1_enstd = TP1_enstd/(NUM_ENSEMBLE-1) + if(associated(TP1_enstd) .and. associated(TP1_enavg)) & + TP1_enstd = sqrt(TP1_enstd - NdivNm1*(TP1_enavg)**2) if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 0f143c6a..08f26842 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -779,9 +779,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_surface_analysis ensstd' ,& + LONG_NAME = 'soil_moisture_surface_analysis_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCSF_ANAstd' ,& + SHORT_NAME = 'WCSF_ANA_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -797,9 +797,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_rootzone_analysis ensstd' ,& + LONG_NAME = 'soil_moisture_rootzone_analysis_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCRZ_ANAstd' ,& + SHORT_NAME = 'WCRZ_ANA_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -815,9 +815,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_moisture_profile_analysis ensstd' ,& + LONG_NAME = 'soil_moisture_profile_analysis_ensstd' ,& UNITS = 'm3 m-3' ,& - SHORT_NAME = 'WCPR_ANAstd' ,& + SHORT_NAME = 'WCPR_ANA_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -833,9 +833,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'ave_catchment_temp_incl_snw_analysisi ensstd' ,& + LONG_NAME = 'ave_catchment_temp_incl_snw_analysis_ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TPSURF_ANAstd' ,& + SHORT_NAME = 'TPSURF_ANA_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -851,9 +851,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& - LONG_NAME = 'soil_temperatures_layer_1_analysis ensstd' ,& + LONG_NAME = 'soil_temperatures_layer_1_analysis_ensstd' ,& UNITS = 'K' ,& - SHORT_NAME = 'TSOIL1_ANAstd' ,& + SHORT_NAME = 'TSOIL1_ANA_ENSSTD' ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) @@ -1528,16 +1528,17 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) !! export for analysis model diagnostics - real, dimension(:),pointer :: SFMC_ana=>null() ! surface soil moisture - real, dimension(:),pointer :: RZMC_ana=>null() ! rootzone soil moisture - real, dimension(:),pointer :: PRMC_ana=>null() ! profile soil moisture - real, dimension(:),pointer :: TPSURF_ana=>null() ! tpsurf - real, dimension(:),pointer :: TSOIL1_ana=>null() ! tsoil1 - real, dimension(:),pointer :: SFMC_anastd=>null() ! surface soil moisture - real, dimension(:),pointer :: RZMC_anastd=>null() ! rootzone soil moisture - real, dimension(:),pointer :: PRMC_anastd=>null() ! profile soil moisture - real, dimension(:),pointer :: TPSURF_anastd=>null() ! tpsurf - real, dimension(:),pointer :: TSOIL1_anastd=>null() ! tsoil1 + real, dimension(:),pointer :: SFMC_ana=>null() ! surface soil moisture + real, dimension(:),pointer :: RZMC_ana=>null() ! rootzone soil moisture + real, dimension(:),pointer :: PRMC_ana=>null() ! profile soil moisture + real, dimension(:),pointer :: TPSURF_ana=>null() ! tpsurf + real, dimension(:),pointer :: TSOIL1_ana=>null() ! tsoil1 + + real, dimension(:),pointer :: SFMC_ana_ensstd=>null() ! surface soil moisture + real, dimension(:),pointer :: RZMC_ana_ensstd=>null() ! rootzone soil moisture + real, dimension(:),pointer :: PRMC_ana_ensstd=>null() ! profile soil moisture + real, dimension(:),pointer :: TPSURF_ana_ensstd=>null() ! tpsurf + real, dimension(:),pointer :: TSOIL1_ana_ensstd=>null() ! tsoil1 !! export for microwave radiative transfer model (mwRTM) @@ -1551,7 +1552,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) character(len=ESMF_MAXSTR) :: ensid_string integer :: ens, nymd, nhms, ens_id_width integer :: LandassimDTstep - real :: nn1 + real :: NdivNm1 #ifdef DBG_LANDASSIM_INPUTS ! vars for debugging purposes @@ -1697,25 +1698,25 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! exports for analysis model diagnostics - call MAPL_GetPointer(export, TPSURF_ana, 'TPSURF_ANA' ,rc=status) + call MAPL_GetPointer(export, TPSURF_ana, 'TPSURF_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TSOIL1_ana, 'TSOIL1_ANA' ,rc=status) + call MAPL_GetPointer(export, TSOIL1_ana, 'TSOIL1_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, SFMC_ana, 'WCSF_ANA' ,rc=status) + call MAPL_GetPointer(export, SFMC_ana, 'WCSF_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, RZMC_ana, 'WCRZ_ANA' ,rc=status) + call MAPL_GetPointer(export, RZMC_ana, 'WCRZ_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) + call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TPSURF_anastd, 'TPSURF_ANAstd' ,rc=status) + call MAPL_GetPointer(export, TPSURF_ana_ensstd, 'TPSURF_ANA_ENSSTD' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, TSOIL1_anastd, 'TSOIL1_ANAstd' ,rc=status) + call MAPL_GetPointer(export, TSOIL1_ana_ensstd, 'TSOIL1_ANA_ENSSTD' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, SFMC_anastd, 'WCSF_ANAstd' ,rc=status) + call MAPL_GetPointer(export, SFMC_ana_ensstd, 'WCSF_ANA_ENSSTD' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, RZMC_anastd, 'WCRZ_ANAstd' ,rc=status) + call MAPL_GetPointer(export, RZMC_ana_ensstd, 'WCRZ_ANA_ENSSTD' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, PRMC_anastd, 'WCPR_ANAstd' ,rc=status) + call MAPL_GetPointer(export, PRMC_ana_ensstd, 'WCPR_ANA_ENSSTD' ,rc=status) VERIFY_(status) ! exports for microwave radiative transfer model (mwRTM) @@ -1994,39 +1995,43 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii) + cat_diagS(ii) end do - cat_diagS_ensstd(:)%sfmc = cat_diagS_ensstd(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc - cat_diagS_ensstd(:)%rzmc = cat_diagS_ensstd(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc - cat_diagS_ensstd(:)%prmc = cat_diagS_ensstd(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc + cat_diagS_ensstd(:)%sfmc = cat_diagS_ensstd(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc + cat_diagS_ensstd(:)%rzmc = cat_diagS_ensstd(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc + cat_diagS_ensstd(:)%prmc = cat_diagS_ensstd(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc cat_diagS_ensstd(:)%tsurf = cat_diagS_ensstd(:)%tsurf + cat_diagS(:)%tsurf*cat_diagS(:)%tsurf - cat_diagS_ensstd(:)%tp(1) = cat_diagS_ensstd(:)%tp(1) + cat_diagS(:)%tp(1) * cat_diagS(:)%tp(1) + cat_diagS_ensstd(:)%tp(1) = cat_diagS_ensstd(:)%tp(1) + cat_diagS(:)%tp(1)*cat_diagS(:)%tp(1) end do do ii=1,N_catl cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize cat_diagS_ensstd(ii) = cat_diagS_ensstd(ii)/real(NUM_ENSEMBLE-1) end do + if (NUM_ENSEMBLE > 1 ) then - nn1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) + NdivNm1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) else - nn1 = real(NUM_ENSEMBLE) + NdivNm1 = real(NUM_ENSEMBLE) endif - cat_diagS_ensstd(:)%sfmc = sqrt(cat_diagS_ensstd(:)%sfmc - nn1*(cat_diagS_ensavg(:)%sfmc)**2 ) - cat_diagS_ensstd(:)%rzmc = sqrt(cat_diagS_ensstd(:)%rzmc - nn1*(cat_diagS_ensavg(:)%rzmc)**2 ) - cat_diagS_ensstd(:)%prmc = sqrt(cat_diagS_ensstd(:)%prmc - nn1*(cat_diagS_ensavg(:)%prmc)**2 ) - cat_diagS_ensstd(:)%tsurf = sqrt(cat_diagS_ensstd(:)%tsurf - nn1*(cat_diagS_ensavg(:)%tsurf)**2 ) - cat_diagS_ensstd(:)%tp(1) = sqrt(cat_diagS_ensstd(:)%tp(1) - nn1*(cat_diagS_ensavg(:)%tp(1))**2 ) + + cat_diagS_ensstd(:)%sfmc = sqrt(cat_diagS_ensstd(:)%sfmc - NdivNm1*(cat_diagS_ensavg(:)%sfmc)**2 ) + cat_diagS_ensstd(:)%rzmc = sqrt(cat_diagS_ensstd(:)%rzmc - NdivNm1*(cat_diagS_ensavg(:)%rzmc)**2 ) + cat_diagS_ensstd(:)%prmc = sqrt(cat_diagS_ensstd(:)%prmc - NdivNm1*(cat_diagS_ensavg(:)%prmc)**2 ) + cat_diagS_ensstd(:)%tsurf = sqrt(cat_diagS_ensstd(:)%tsurf - NdivNm1*(cat_diagS_ensavg(:)%tsurf)**2 ) + cat_diagS_ensstd(:)%tp(1) = sqrt(cat_diagS_ensstd(:)%tp(1) - NdivNm1*(cat_diagS_ensavg(:)%tp(1))**2 ) + ! set export variables - if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc - if(associated(RZMC_ana)) RZMC_ana(:) = cat_diagS_ensavg(:)%rzmc - if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc - if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf - if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K - if(associated(SFMC_anastd)) SFMC_anastd(:) = cat_diagS_ensstd(:)%sfmc - if(associated(RZMC_anastd)) RZMC_anastd(:) = cat_diagS_ensstd(:)%rzmc - if(associated(PRMC_anastd)) PRMC_anastd(:) = cat_diagS_ensstd(:)%prmc - if(associated(TPSURF_anastd)) TPSURF_anastd(:) = cat_diagS_ensstd(:)%tsurf - if(associated(TSOIL1_anastd)) TSOIL1_anastd(:) = cat_diagS_ensstd(:)%tp(1) + if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc + if(associated(RZMC_ana)) RZMC_ana(:) = cat_diagS_ensavg(:)%rzmc + if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc + if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf + if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K + + if(associated(SFMC_ana_ensstd)) SFMC_ana_ensstd(:) = cat_diagS_ensstd(:)%sfmc + if(associated(RZMC_ana_ensstd)) RZMC_ana_ensstd(:) = cat_diagS_ensstd(:)%rzmc + if(associated(PRMC_ana_ensstd)) PRMC_ana_ensstd(:) = cat_diagS_ensstd(:)%prmc + if(associated(TPSURF_ana_ensstd)) TPSURF_ana_ensstd(:) = cat_diagS_ensstd(:)%tsurf + if(associated(TSOIL1_ana_ensstd)) TSOIL1_ana_ensstd(:) = cat_diagS_ensstd(:)%tp(1) if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY From 4ef74731d9d92da7d6c1c90ea7c08eb452267d31 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 16 Mar 2023 14:42:14 -0400 Subject: [PATCH 04/10] added new operators for cat_diagS type; protected against division by zero in ensstd calculations --- .../GEOS_LandAssimGridComp.F90 | 52 ++-- .../GEOSldas_GridComp/Shared/catch_types.F90 | 286 +++++++++++++++++- 2 files changed, 302 insertions(+), 36 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 08f26842..648df3c3 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -49,7 +49,8 @@ module GEOS_LandAssimGridCompMod use catch_types, only: cat_progn_type use catch_types, only: cat_param_type use catch_types, only: cat_diagS_type - use catch_types, only: assignment(=), operator (+), operator (/) + use catch_types, only: cat_diagS_sqrt + use catch_types, only: assignment(=), operator (+), operator (-), operator (*), operator (/) use clsm_bias_routines, only: initialize_obs_bias use clsm_bias_routines, only: read_cat_bias_inputs @@ -1552,7 +1553,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) character(len=ESMF_MAXSTR) :: ensid_string integer :: ens, nymd, nhms, ens_id_width integer :: LandassimDTstep - real :: NdivNm1 + real :: Nm1, NdivNm1 #ifdef DBG_LANDASSIM_INPUTS ! vars for debugging purposes @@ -1978,10 +1979,12 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) allocate(cat_diagS_ensstd(N_catl)) do ii=1,N_catl - cat_diagS_ensavg(ii) = 0.0 ! initialize ens average - cat_diagS_ensstd(ii) = 0.0 + cat_diagS_ensavg(ii) = 0.0 ! initialize sum for ens average + cat_diagS_ensstd(ii) = 0.0 ! initialize sum of squares for ensemble standard deviation end do + ! compute sum (and sum of squares) of ensemble members + do n_e=1,NUM_ENSEMBLE ! make a copy of cat_progn to ensure 0-diff (recompute_diagS() potentially alters its input cat_progn) @@ -1993,32 +1996,33 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) call recompute_diagS( N_catl, cat_param, cat_progn_tmp, cat_diagS ) do ii=1,N_catl - cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii) + cat_diagS(ii) + cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii) + cat_diagS(ii) ! sum + cat_diagS_ensstd(ii) = cat_diagS_ensstd(ii) + ( cat_diagS(ii) * cat_diagS(ii) ) ! sum of squares end do - cat_diagS_ensstd(:)%sfmc = cat_diagS_ensstd(:)%sfmc + cat_diagS(:)%sfmc*cat_diagS(:)%sfmc - cat_diagS_ensstd(:)%rzmc = cat_diagS_ensstd(:)%rzmc + cat_diagS(:)%rzmc*cat_diagS(:)%rzmc - cat_diagS_ensstd(:)%prmc = cat_diagS_ensstd(:)%prmc + cat_diagS(:)%prmc*cat_diagS(:)%rzmc - cat_diagS_ensstd(:)%tsurf = cat_diagS_ensstd(:)%tsurf + cat_diagS(:)%tsurf*cat_diagS(:)%tsurf - cat_diagS_ensstd(:)%tp(1) = cat_diagS_ensstd(:)%tp(1) + cat_diagS(:)%tp(1)*cat_diagS(:)%tp(1) + end do + ! finalize ensemble average and standard deviation + + Nm1 = real(NUM_ENSEMBLE-1) + + if (NUM_ENSEMBLE > 1) NdivNm1 = real(NUM_ENSEMBLE)/Nm1 + do ii=1,N_catl - cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize - cat_diagS_ensstd(ii) = cat_diagS_ensstd(ii)/real(NUM_ENSEMBLE-1) + + cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg + + if (NUM_ENSEMBLE > 1 ) then + + cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii)) ) + + else ! NUM_ENSEMBLE = 1 + + cat_diagS_ensstd(ii) = MAPL_UNDEF + + end if end do - if (NUM_ENSEMBLE > 1 ) then - NdivNm1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) - else - NdivNm1 = real(NUM_ENSEMBLE) - endif - - cat_diagS_ensstd(:)%sfmc = sqrt(cat_diagS_ensstd(:)%sfmc - NdivNm1*(cat_diagS_ensavg(:)%sfmc)**2 ) - cat_diagS_ensstd(:)%rzmc = sqrt(cat_diagS_ensstd(:)%rzmc - NdivNm1*(cat_diagS_ensavg(:)%rzmc)**2 ) - cat_diagS_ensstd(:)%prmc = sqrt(cat_diagS_ensstd(:)%prmc - NdivNm1*(cat_diagS_ensavg(:)%prmc)**2 ) - cat_diagS_ensstd(:)%tsurf = sqrt(cat_diagS_ensstd(:)%tsurf - NdivNm1*(cat_diagS_ensavg(:)%tsurf)**2 ) - cat_diagS_ensstd(:)%tp(1) = sqrt(cat_diagS_ensstd(:)%tp(1) - NdivNm1*(cat_diagS_ensavg(:)%tp(1))**2 ) - ! set export variables if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc diff --git a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 index 0c97fc9b..9ebcf89c 100644 --- a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 @@ -36,7 +36,9 @@ module catch_types public :: cat_progn_type, cat_diagS_type, cat_diagF_type public :: cat_param_type, cat_force_type - public :: assignment (=), operator (/), operator (+) + public :: assignment (=), operator (*), operator (/), operator (+), operator (-) + + public :: cat_diagS_sqrt public :: catprogn2wesn, catprogn2htsn, catprogn2sndz, catprogn2ghtcnt @@ -373,6 +375,12 @@ module catch_types module procedure scalar2cat_force end interface + interface operator (*) + module procedure cat_diagS_mult_scalar + module procedure cat_diagF_mult_scalar + module procedure cat_diagS_mult_cat_diagS + end interface + interface operator (/) module procedure cat_progn_div_scalar module procedure cat_diagS_div_scalar @@ -387,6 +395,11 @@ module catch_types module procedure add_cat_force end interface + interface operator (-) + module procedure subtract_cat_diagS + module procedure subtract_cat_diagF + end interface + ! ---------------------------------------------------------------- contains @@ -466,34 +479,145 @@ end subroutine scalar2cat_diagF ! ----------------------------------------------------------- - function cat_diagS_div_scalar( cat_diagS, scalar ) + function cat_diagS_mult_scalar( cat_diagS, scalar ) implicit none - type(cat_diagS_type) :: cat_diagS_div_scalar + type(cat_diagS_type) :: cat_diagS_mult_scalar type(cat_diagS_type), intent(in) :: cat_diagS real, intent(in) :: scalar integer :: i ! local - cat_diagS_div_scalar%ar1 = cat_diagS%ar1 / scalar - cat_diagS_div_scalar%ar2 = cat_diagS%ar2 / scalar + cat_diagS_mult_scalar%ar1 = cat_diagS%ar1 * scalar + cat_diagS_mult_scalar%ar2 = cat_diagS%ar2 * scalar + + cat_diagS_mult_scalar%asnow = cat_diagS%asnow * scalar + + cat_diagS_mult_scalar%sfmc = cat_diagS%sfmc * scalar + cat_diagS_mult_scalar%rzmc = cat_diagS%rzmc * scalar + cat_diagS_mult_scalar%prmc = cat_diagS%prmc * scalar + + cat_diagS_mult_scalar%tsurf = cat_diagS%tsurf * scalar + + do i=1,N_gt + cat_diagS_mult_scalar%tp(i) = cat_diagS%tp(i) * scalar + end do + + do i=1,N_snow + cat_diagS_mult_scalar%tpsn(i) = cat_diagS%tpsn(i)* scalar + end do + + end function cat_diagS_mult_scalar + + ! ----------------------------------------------------------- + + function cat_diagS_mult_cat_diagS( cat_diagS_1, cat_diagS_2 ) + + implicit none + + type(cat_diagS_type) :: cat_diagS_mult_cat_diagS + type(cat_diagS_type), intent(in) :: cat_diagS_1, cat_diagS_2 + + integer :: i ! local + + cat_diagS_mult_cat_diagS%ar1 = cat_diagS_1%ar1 * cat_diagS_2%ar1 + cat_diagS_mult_cat_diagS%ar2 = cat_diagS_1%ar2 * cat_diagS_2%ar2 + + cat_diagS_mult_cat_diagS%asnow = cat_diagS_1%asnow * cat_diagS_2%asnow + + cat_diagS_mult_cat_diagS%sfmc = cat_diagS_1%sfmc * cat_diagS_2%sfmc + cat_diagS_mult_cat_diagS%rzmc = cat_diagS_1%rzmc * cat_diagS_2%rzmc + cat_diagS_mult_cat_diagS%prmc = cat_diagS_1%prmc * cat_diagS_2%prmc + + cat_diagS_mult_cat_diagS%tsurf = cat_diagS_1%tsurf * cat_diagS_2%tsurf + + do i=1,N_gt + cat_diagS_mult_cat_diagS%tp(i) = cat_diagS_1%tp(i) * cat_diagS_2%tp(i) + end do + + do i=1,N_snow + cat_diagS_mult_cat_diagS%tpsn(i) = cat_diagS_1%tpsn(i) * cat_diagS_2%tpsn(i) + end do + + end function cat_diagS_mult_cat_diagS - cat_diagS_div_scalar%asnow = cat_diagS%asnow / scalar + ! ----------------------------------------------------------- - cat_diagS_div_scalar%sfmc = cat_diagS%sfmc / scalar - cat_diagS_div_scalar%rzmc = cat_diagS%rzmc / scalar - cat_diagS_div_scalar%prmc = cat_diagS%prmc / scalar + function cat_diagF_mult_scalar( cat_diagF, scalar ) - cat_diagS_div_scalar%tsurf = cat_diagS%tsurf / scalar + implicit none + + type(cat_diagF_type) :: cat_diagF_mult_scalar + type(cat_diagF_type), intent(in) :: cat_diagF + + real, intent(in) :: scalar + + cat_diagF_mult_scalar%shflux = cat_diagF%shflux * scalar + cat_diagF_mult_scalar%lhflux = cat_diagF%lhflux * scalar + cat_diagF_mult_scalar%ghflux = cat_diagF%ghflux * scalar + + cat_diagF_mult_scalar%evap = cat_diagF%evap * scalar + cat_diagF_mult_scalar%eint = cat_diagF%eint * scalar + cat_diagF_mult_scalar%esoi = cat_diagF%esoi * scalar + cat_diagF_mult_scalar%eveg = cat_diagF%eveg * scalar + cat_diagF_mult_scalar%esno = cat_diagF%esno * scalar + + + cat_diagF_mult_scalar%runoff = cat_diagF%runoff * scalar + cat_diagF_mult_scalar%runsrf = cat_diagF%runsrf * scalar + cat_diagF_mult_scalar%bflow = cat_diagF%bflow * scalar + + cat_diagF_mult_scalar%snmelt = cat_diagF%snmelt * scalar + + cat_diagF_mult_scalar%lwup = cat_diagF%lwup * scalar + cat_diagF_mult_scalar%swup = cat_diagF%swup * scalar + + cat_diagF_mult_scalar%qinfil = cat_diagF%qinfil * scalar + + cat_diagF_mult_scalar%hsnacc = cat_diagF%hsnacc * scalar + cat_diagF_mult_scalar%evacc = cat_diagF%evacc * scalar + cat_diagF_mult_scalar%shacc = cat_diagF%shacc * scalar + cat_diagF_mult_scalar%lhacc = cat_diagF%lhacc * scalar + cat_diagF_mult_scalar%eacc_0 = cat_diagF%eacc_0 * scalar + + cat_diagF_mult_scalar%t2m = cat_diagF%t2m * scalar + cat_diagF_mult_scalar%q2m = cat_diagF%q2m * scalar + + + end function cat_diagF_mult_scalar + + ! ----------------------------------------------------------- + + function cat_diagS_div_scalar( cat_diagS, scalar ) + + implicit none + + type(cat_diagS_type) :: cat_diagS_div_scalar + type(cat_diagS_type), intent(in) :: cat_diagS + + real, intent(in) :: scalar + + integer :: i ! local + + cat_diagS_div_scalar%ar1 = cat_diagS%ar1 / scalar + cat_diagS_div_scalar%ar2 = cat_diagS%ar2 / scalar + + cat_diagS_div_scalar%asnow = cat_diagS%asnow / scalar + + cat_diagS_div_scalar%sfmc = cat_diagS%sfmc / scalar + cat_diagS_div_scalar%rzmc = cat_diagS%rzmc / scalar + cat_diagS_div_scalar%prmc = cat_diagS%prmc / scalar + + cat_diagS_div_scalar%tsurf = cat_diagS%tsurf / scalar do i=1,N_gt - cat_diagS_div_scalar%tp(i) = cat_diagS%tp(i) / scalar + cat_diagS_div_scalar%tp(i) = cat_diagS%tp(i) / scalar end do do i=1,N_snow - cat_diagS_div_scalar%tpsn(i)= cat_diagS%tpsn(i)/ scalar + cat_diagS_div_scalar%tpsn(i) = cat_diagS%tpsn(i)/ scalar end do end function cat_diagS_div_scalar @@ -617,6 +741,144 @@ function add_cat_diagF( cat_diagF_1, cat_diagF_2 ) end function add_cat_diagF + ! ----------------------------------------------------------- + + function subtract_cat_diagS( cat_diagS_1, cat_diagS_2 ) + + implicit none + + type(cat_diagS_type) :: subtract_cat_diagS + type(cat_diagS_type), intent(in) :: cat_diagS_1, cat_diagS_2 + + integer :: i ! local + + subtract_cat_diagS%ar1 = cat_diagS_1%ar1 - cat_diagS_2%ar1 + subtract_cat_diagS%ar2 = cat_diagS_1%ar2 - cat_diagS_2%ar2 + + subtract_cat_diagS%asnow = cat_diagS_1%asnow - cat_diagS_2%asnow + + subtract_cat_diagS%sfmc = cat_diagS_1%sfmc - cat_diagS_2%sfmc + subtract_cat_diagS%rzmc = cat_diagS_1%rzmc - cat_diagS_2%rzmc + subtract_cat_diagS%prmc = cat_diagS_1%prmc - cat_diagS_2%prmc + + subtract_cat_diagS%tsurf = cat_diagS_1%tsurf - cat_diagS_2%tsurf + + do i=1,N_gt + subtract_cat_diagS%tp(i) = cat_diagS_1%tp(i) - cat_diagS_2%tp(i) + end do + + do i=1,N_snow + subtract_cat_diagS%tpsn(i) = cat_diagS_1%tpsn(i) - cat_diagS_2%tpsn(i) + end do + + end function subtract_cat_diagS + + ! ----------------------------------------------------------- + + function subtract_cat_diagF( cat_diagF_1, cat_diagF_2 ) + + implicit none + + type(cat_diagF_type) :: subtract_cat_diagF + type(cat_diagF_type), intent(in) :: cat_diagF_1, cat_diagF_2 + + subtract_cat_diagF%shflux = cat_diagF_1%shflux - cat_diagF_2%shflux + subtract_cat_diagF%lhflux = cat_diagF_1%lhflux - cat_diagF_2%lhflux + subtract_cat_diagF%ghflux = cat_diagF_1%ghflux - cat_diagF_2%ghflux + + subtract_cat_diagF%evap = cat_diagF_1%evap - cat_diagF_2%evap + subtract_cat_diagF%eint = cat_diagF_1%eint - cat_diagF_2%eint + subtract_cat_diagF%esoi = cat_diagF_1%esoi - cat_diagF_2%esoi + subtract_cat_diagF%eveg = cat_diagF_1%eveg - cat_diagF_2%eveg + subtract_cat_diagF%esno = cat_diagF_1%esno - cat_diagF_2%esno + + subtract_cat_diagF%runoff = cat_diagF_1%runoff - cat_diagF_2%runoff + subtract_cat_diagF%runsrf = cat_diagF_1%runsrf - cat_diagF_2%runsrf + subtract_cat_diagF%bflow = cat_diagF_1%bflow - cat_diagF_2%bflow + + subtract_cat_diagF%snmelt = cat_diagF_1%snmelt - cat_diagF_2%snmelt + + subtract_cat_diagF%lwup = cat_diagF_1%lwup - cat_diagF_2%lwup + subtract_cat_diagF%swup = cat_diagF_1%swup - cat_diagF_2%swup + + subtract_cat_diagF%qinfil = cat_diagF_1%qinfil - cat_diagF_2%qinfil + + subtract_cat_diagF%hsnacc = cat_diagF_1%hsnacc - cat_diagF_2%hsnacc + subtract_cat_diagF%evacc = cat_diagF_1%evacc - cat_diagF_2%evacc + subtract_cat_diagF%shacc = cat_diagF_1%shacc - cat_diagF_2%shacc + subtract_cat_diagF%lhacc = cat_diagF_1%lhacc - cat_diagF_2%lhacc + subtract_cat_diagF%eacc_0 = cat_diagF_1%eacc_0 - cat_diagF_2%eacc_0 + + subtract_cat_diagF%t2m = cat_diagF_1%t2m - cat_diagF_2%t2m + subtract_cat_diagF%q2m = cat_diagF_1%q2m - cat_diagF_2%q2m + + + end function subtract_cat_diagF + + ! ----------------------------------------------------------- + + function cat_diagS_square( cat_diagS ) + + implicit none + + type(cat_diagS_type) :: cat_diagS_square + type(cat_diagS_type), intent(in) :: cat_diagS + + integer :: i ! local + + cat_diagS_square%ar1 = cat_diagS%ar1 * cat_diagS%ar1 + cat_diagS_square%ar2 = cat_diagS%ar2 * cat_diagS%ar2 + + cat_diagS_square%asnow = cat_diagS%asnow * cat_diagS%asnow + + cat_diagS_square%sfmc = cat_diagS%sfmc * cat_diagS%sfmc + cat_diagS_square%rzmc = cat_diagS%rzmc * cat_diagS%rzmc + cat_diagS_square%prmc = cat_diagS%prmc * cat_diagS%prmc + + cat_diagS_square%tsurf = cat_diagS%tsurf * cat_diagS%tsurf + + do i=1,N_gt + cat_diagS_square%tp(i) = cat_diagS%tp(i) * cat_diagS%tp(i) + end do + + do i=1,N_snow + cat_diagS_square%tpsn(i) = cat_diagS%tpsn(i) * cat_diagS%tpsn(i) + end do + + end function cat_diagS_square + + ! ----------------------------------------------------------- + + function cat_diagS_sqrt( cat_diagS ) + + implicit none + + type(cat_diagS_type) :: cat_diagS_sqrt + type(cat_diagS_type), intent(in) :: cat_diagS + + integer :: i ! local + + cat_diagS_sqrt%ar1 = sqrt( cat_diagS%ar1 ) + cat_diagS_sqrt%ar2 = sqrt( cat_diagS%ar2 ) + + cat_diagS_sqrt%asnow = sqrt( cat_diagS%asnow ) + + cat_diagS_sqrt%sfmc = sqrt( cat_diagS%sfmc ) + cat_diagS_sqrt%rzmc = sqrt( cat_diagS%rzmc ) + cat_diagS_sqrt%prmc = sqrt( cat_diagS%prmc ) + + cat_diagS_sqrt%tsurf = sqrt( cat_diagS%tsurf ) + + do i=1,N_gt + cat_diagS_sqrt%tp(i) = sqrt( cat_diagS%tp(i) ) + end do + + do i=1,N_snow + cat_diagS_sqrt%tpsn(i) = sqrt( cat_diagS%tpsn(i) ) + end do + + end function cat_diagS_sqrt + ! ******************************************************************* subroutine scalar2cat_progn( cat_progn, scalar ) From 3a16724d12b9282d67e1063a89f3123d5c2cfbcf Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 16 Mar 2023 17:27:03 -0400 Subject: [PATCH 05/10] additional protections against division by zero for calc of ens std --- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 53 +++++++++++-------- 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index 833247df..d4b48e26 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -2448,7 +2448,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:), pointer :: CNBURN, CNBURN_enavg real, dimension(:), pointer :: CNFSEL, CNFSEL_enavg - real :: NdivNm1 + real :: Nm1, NdivNm1 ! Get my name and setup traceback handle call ESMF_GridCompget(gc, name=comp_name, rc=status) @@ -3618,11 +3618,10 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(collect_land_counter == NUM_ENSEMBLE) then - if (NUM_ENSEMBLE >1 ) then - NdivNm1 = real(NUM_ENSEMBLE)/real(NUM_ENSEMBLE-1) - else - NdivNm1 = real(NUM_ENSEMBLE) - endif + + Nm1 = real(NUM_ENSEMBLE-1) + if (NUM_ENSEMBLE>1) NdivNm1 = real(NUM_ENSEMBLE)/Nm1 + collect_land_counter = 0 if(associated(TC_enavg)) TC_enavg = TC_enavg/NUM_ENSEMBLE if(associated(QC_enavg)) QC_enavg = QC_enavg/NUM_ENSEMBLE @@ -3670,9 +3669,11 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) !if(associated(RZEQ_enavg)) RZEQ_enavg = RZEQ_enavg/NUM_ENSEMBLE if(associated(GHFLX_enavg)) GHFLX_enavg = GHFLX_enavg/NUM_ENSEMBLE if(associated(TPSURF_enavg)) TPSURF_enavg = TPSURF_enavg/NUM_ENSEMBLE - if(associated(TPSURF_enstd)) TPSURF_enstd = TPSURF_enstd/(NUM_ENSEMBLE-1) - if(associated(TPSURF_enstd) .and. associated(TPSURF_enavg)) & - TPSURF_enstd = sqrt(TPSURF_enstd - NdivNm1*(TPSURF_enavg)**2) + if((NUM_ENSEMBLE>1) .and. associated(TPSURF_enstd) .and. associated(TPSURF_enavg)) then + TPSURF_enstd = sqrt( TPSURF_enstd/Nm1 - NdivNm1*(TPSURF_enavg**2) ) + else if (associated(TPSURF_enstd)) then + TPSURF_enstd = MAPL_UNDEF + end if if(associated(TPSNOW_enavg)) TPSNOW_enavg = TPSNOW_enavg/NUM_ENSEMBLE if(associated(TPUNST_enavg)) TPUNST_enavg = TPUNST_enavg/NUM_ENSEMBLE if(associated(TPSAT_enavg)) TPSAT_enavg = TPSAT_enavg/NUM_ENSEMBLE @@ -3689,21 +3690,29 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET2_enavg)) WET2_enavg = WET2_enavg/NUM_ENSEMBLE if(associated(WET3_enavg)) WET3_enavg = WET3_enavg/NUM_ENSEMBLE if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE - if(associated(WCSF_enstd)) WCSF_enstd = WCSF_enstd/(NUM_ENSEMBLE-1) - if(associated(WCSF_enstd) .and. associated(WCSF_enavg)) & - WCSF_enstd = sqrt(WCSF_enstd - NdivNm1*(WCSF_enavg)**2) + if((NUM_ENSEMBLE>1) .and. associated(WCSF_enstd) .and. associated(WCSF_enavg)) then + WCSF_enstd = sqrt( WCSF_enstd/Nm1 - NdivNm1*(WCSF_enavg**2) ) + else if (associated(WCSF_enstd)) then + WCSF_enstd = MAPL_UNDEF + end if if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE - if(associated(WCRZ_enstd)) WCRZ_enstd = WCRZ_enstd/(NUM_ENSEMBLE-1) - if(associated(WCRZ_enstd) .and. associated(WCRZ_enavg)) & - WCRZ_enstd = sqrt(WCRZ_enstd - NdivNm1*(WCRZ_enavg)**2) + if((NUM_ENSEMBLE>1) .and. associated(WCRZ_enstd) .and. associated(WCRZ_enavg)) then + WCRZ_enstd = sqrt( WCRZ_enstd/Nm1 - NdivNm1*(WCRZ_enavg**2) ) + else if (associated(WCRZ_enstd)) then + WCRZ_enstd = MAPL_UNDEF + end if if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE - if(associated(WCPR_enstd)) WCPR_enstd = WCPR_enstd/(NUM_ENSEMBLE-1) - if(associated(WCPR_enstd) .and. associated(WCPR_enavg)) & - WCPR_enstd = sqrt(WCPR_enstd - NdivNm1*(WCPR_enavg)**2) - if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 - if(associated(TP1_enstd)) TP1_enstd = TP1_enstd/(NUM_ENSEMBLE-1) - if(associated(TP1_enstd) .and. associated(TP1_enavg)) & - TP1_enstd = sqrt(TP1_enstd - NdivNm1*(TP1_enavg)**2) + if((NUM_ENSEMBLE>1) .and. associated(WCPR_enstd) .and. associated(WCPR_enavg)) then + WCPR_enstd = sqrt( WCPR_enstd/Nm1 - NdivNm1*(WCPR_enavg**2) ) + else if (associated(WCPR_enstd)) then + WCPR_enstd = MAPL_UNDEF + end if + if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units K + if((NUM_ENSEMBLE>1) .and. associated(TP1_enstd) .and. associated(TP1_enavg)) then + TP1_enstd = sqrt( TP1_enstd/Nm1 - NdivNm1*(TP1_enavg**2) ) + else if (associated(TP1_enstd)) then + TP1_enstd = MAPL_UNDEF + end if if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 From 435066f0f3a891b8288aa99bc0f225cae7cb17d6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 16 Mar 2023 18:17:31 -0400 Subject: [PATCH 06/10] fixed new multiplication operators for cat_diagS (order of args) --- .../GEOSldas_GridComp/Shared/catch_types.F90 | 146 +++++++----------- 1 file changed, 57 insertions(+), 89 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 index 9ebcf89c..df22a2b3 100644 --- a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 @@ -376,8 +376,8 @@ module catch_types end interface interface operator (*) - module procedure cat_diagS_mult_scalar - module procedure cat_diagF_mult_scalar + module procedure scalar_mult_cat_diagS + module procedure scalar_mult_cat_diagF module procedure cat_diagS_mult_cat_diagS end interface @@ -479,37 +479,37 @@ end subroutine scalar2cat_diagF ! ----------------------------------------------------------- - function cat_diagS_mult_scalar( cat_diagS, scalar ) + function scalar_mult_cat_diagS( scalar, cat_diagS ) implicit none - type(cat_diagS_type) :: cat_diagS_mult_scalar + type(cat_diagS_type) :: scalar_mult_cat_diagS type(cat_diagS_type), intent(in) :: cat_diagS real, intent(in) :: scalar integer :: i ! local - cat_diagS_mult_scalar%ar1 = cat_diagS%ar1 * scalar - cat_diagS_mult_scalar%ar2 = cat_diagS%ar2 * scalar - - cat_diagS_mult_scalar%asnow = cat_diagS%asnow * scalar - - cat_diagS_mult_scalar%sfmc = cat_diagS%sfmc * scalar - cat_diagS_mult_scalar%rzmc = cat_diagS%rzmc * scalar - cat_diagS_mult_scalar%prmc = cat_diagS%prmc * scalar - - cat_diagS_mult_scalar%tsurf = cat_diagS%tsurf * scalar - - do i=1,N_gt - cat_diagS_mult_scalar%tp(i) = cat_diagS%tp(i) * scalar - end do - - do i=1,N_snow - cat_diagS_mult_scalar%tpsn(i) = cat_diagS%tpsn(i)* scalar + scalar_mult_cat_diagS%ar1 = scalar * cat_diagS%ar1 + scalar_mult_cat_diagS%ar2 = scalar * cat_diagS%ar2 + + scalar_mult_cat_diagS%asnow = scalar * cat_diagS%asnow + + scalar_mult_cat_diagS%sfmc = scalar * cat_diagS%sfmc + scalar_mult_cat_diagS%rzmc = scalar * cat_diagS%rzmc + scalar_mult_cat_diagS%prmc = scalar * cat_diagS%prmc + + scalar_mult_cat_diagS%tsurf = scalar * cat_diagS%tsurf + + do i=1,N_gt + scalar_mult_cat_diagS%tp(i) = scalar * cat_diagS%tp(i) + end do + + do i=1,N_snow + scalar_mult_cat_diagS%tpsn(i) = scalar * cat_diagS%tpsn(i) end do - end function cat_diagS_mult_scalar + end function scalar_mult_cat_diagS ! ----------------------------------------------------------- @@ -545,48 +545,48 @@ end function cat_diagS_mult_cat_diagS ! ----------------------------------------------------------- - function cat_diagF_mult_scalar( cat_diagF, scalar ) + function scalar_mult_cat_diagF( scalar, cat_diagF ) implicit none - type(cat_diagF_type) :: cat_diagF_mult_scalar + type(cat_diagF_type) :: scalar_mult_cat_diagF type(cat_diagF_type), intent(in) :: cat_diagF real, intent(in) :: scalar - cat_diagF_mult_scalar%shflux = cat_diagF%shflux * scalar - cat_diagF_mult_scalar%lhflux = cat_diagF%lhflux * scalar - cat_diagF_mult_scalar%ghflux = cat_diagF%ghflux * scalar - - cat_diagF_mult_scalar%evap = cat_diagF%evap * scalar - cat_diagF_mult_scalar%eint = cat_diagF%eint * scalar - cat_diagF_mult_scalar%esoi = cat_diagF%esoi * scalar - cat_diagF_mult_scalar%eveg = cat_diagF%eveg * scalar - cat_diagF_mult_scalar%esno = cat_diagF%esno * scalar - - - cat_diagF_mult_scalar%runoff = cat_diagF%runoff * scalar - cat_diagF_mult_scalar%runsrf = cat_diagF%runsrf * scalar - cat_diagF_mult_scalar%bflow = cat_diagF%bflow * scalar - - cat_diagF_mult_scalar%snmelt = cat_diagF%snmelt * scalar - - cat_diagF_mult_scalar%lwup = cat_diagF%lwup * scalar - cat_diagF_mult_scalar%swup = cat_diagF%swup * scalar - - cat_diagF_mult_scalar%qinfil = cat_diagF%qinfil * scalar - - cat_diagF_mult_scalar%hsnacc = cat_diagF%hsnacc * scalar - cat_diagF_mult_scalar%evacc = cat_diagF%evacc * scalar - cat_diagF_mult_scalar%shacc = cat_diagF%shacc * scalar - cat_diagF_mult_scalar%lhacc = cat_diagF%lhacc * scalar - cat_diagF_mult_scalar%eacc_0 = cat_diagF%eacc_0 * scalar - - cat_diagF_mult_scalar%t2m = cat_diagF%t2m * scalar - cat_diagF_mult_scalar%q2m = cat_diagF%q2m * scalar - - - end function cat_diagF_mult_scalar + scalar_mult_cat_diagF%shflux = scalar * cat_diagF%shflux + scalar_mult_cat_diagF%lhflux = scalar * cat_diagF%lhflux + scalar_mult_cat_diagF%ghflux = scalar * cat_diagF%ghflux + + scalar_mult_cat_diagF%evap = scalar * cat_diagF%evap + scalar_mult_cat_diagF%eint = scalar * cat_diagF%eint + scalar_mult_cat_diagF%esoi = scalar * cat_diagF%esoi + scalar_mult_cat_diagF%eveg = scalar * cat_diagF%eveg + scalar_mult_cat_diagF%esno = scalar * cat_diagF%esno + + + scalar_mult_cat_diagF%runoff = scalar * cat_diagF%runoff + scalar_mult_cat_diagF%runsrf = scalar * cat_diagF%runsrf + scalar_mult_cat_diagF%bflow = scalar * cat_diagF%bflow + + scalar_mult_cat_diagF%snmelt = scalar * cat_diagF%snmelt + + scalar_mult_cat_diagF%lwup = scalar * cat_diagF%lwup + scalar_mult_cat_diagF%swup = scalar * cat_diagF%swup + + scalar_mult_cat_diagF%qinfil = scalar * cat_diagF%qinfil + + scalar_mult_cat_diagF%hsnacc = scalar * cat_diagF%hsnacc + scalar_mult_cat_diagF%evacc = scalar * cat_diagF%evacc + scalar_mult_cat_diagF%shacc = scalar * cat_diagF%shacc + scalar_mult_cat_diagF%lhacc = scalar * cat_diagF%lhacc + scalar_mult_cat_diagF%eacc_0 = scalar * cat_diagF%eacc_0 + + scalar_mult_cat_diagF%t2m = scalar * cat_diagF%t2m + scalar_mult_cat_diagF%q2m = scalar * cat_diagF%q2m + + + end function scalar_mult_cat_diagF ! ----------------------------------------------------------- @@ -817,38 +817,6 @@ end function subtract_cat_diagF ! ----------------------------------------------------------- - function cat_diagS_square( cat_diagS ) - - implicit none - - type(cat_diagS_type) :: cat_diagS_square - type(cat_diagS_type), intent(in) :: cat_diagS - - integer :: i ! local - - cat_diagS_square%ar1 = cat_diagS%ar1 * cat_diagS%ar1 - cat_diagS_square%ar2 = cat_diagS%ar2 * cat_diagS%ar2 - - cat_diagS_square%asnow = cat_diagS%asnow * cat_diagS%asnow - - cat_diagS_square%sfmc = cat_diagS%sfmc * cat_diagS%sfmc - cat_diagS_square%rzmc = cat_diagS%rzmc * cat_diagS%rzmc - cat_diagS_square%prmc = cat_diagS%prmc * cat_diagS%prmc - - cat_diagS_square%tsurf = cat_diagS%tsurf * cat_diagS%tsurf - - do i=1,N_gt - cat_diagS_square%tp(i) = cat_diagS%tp(i) * cat_diagS%tp(i) - end do - - do i=1,N_snow - cat_diagS_square%tpsn(i) = cat_diagS%tpsn(i) * cat_diagS%tpsn(i) - end do - - end function cat_diagS_square - - ! ----------------------------------------------------------- - function cat_diagS_sqrt( cat_diagS ) implicit none From 37dbf00b0edbf0ad4a284aa7f67bf2abc6e7ab52 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 16 Mar 2023 22:15:29 -0400 Subject: [PATCH 07/10] ensure non-negative ens std output --- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 10 +++---- .../GEOS_LandAssimGridComp.F90 | 26 +++++++++---------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index d4b48e26..40f8170d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -3670,7 +3670,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(GHFLX_enavg)) GHFLX_enavg = GHFLX_enavg/NUM_ENSEMBLE if(associated(TPSURF_enavg)) TPSURF_enavg = TPSURF_enavg/NUM_ENSEMBLE if((NUM_ENSEMBLE>1) .and. associated(TPSURF_enstd) .and. associated(TPSURF_enavg)) then - TPSURF_enstd = sqrt( TPSURF_enstd/Nm1 - NdivNm1*(TPSURF_enavg**2) ) + TPSURF_enstd = max( sqrt( TPSURF_enstd/Nm1 - NdivNm1*(TPSURF_enavg**2) ), 0. ) else if (associated(TPSURF_enstd)) then TPSURF_enstd = MAPL_UNDEF end if @@ -3691,25 +3691,25 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET3_enavg)) WET3_enavg = WET3_enavg/NUM_ENSEMBLE if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE if((NUM_ENSEMBLE>1) .and. associated(WCSF_enstd) .and. associated(WCSF_enavg)) then - WCSF_enstd = sqrt( WCSF_enstd/Nm1 - NdivNm1*(WCSF_enavg**2) ) + WCSF_enstd = max( sqrt( WCSF_enstd/Nm1 - NdivNm1*(WCSF_enavg**2) ), 0. ) else if (associated(WCSF_enstd)) then WCSF_enstd = MAPL_UNDEF end if if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE if((NUM_ENSEMBLE>1) .and. associated(WCRZ_enstd) .and. associated(WCRZ_enavg)) then - WCRZ_enstd = sqrt( WCRZ_enstd/Nm1 - NdivNm1*(WCRZ_enavg**2) ) + WCRZ_enstd = max( sqrt( WCRZ_enstd/Nm1 - NdivNm1*(WCRZ_enavg**2) ), 0. ) else if (associated(WCRZ_enstd)) then WCRZ_enstd = MAPL_UNDEF end if if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE if((NUM_ENSEMBLE>1) .and. associated(WCPR_enstd) .and. associated(WCPR_enavg)) then - WCPR_enstd = sqrt( WCPR_enstd/Nm1 - NdivNm1*(WCPR_enavg**2) ) + WCPR_enstd = max( sqrt( WCPR_enstd/Nm1 - NdivNm1*(WCPR_enavg**2) ), 0. ) else if (associated(WCPR_enstd)) then WCPR_enstd = MAPL_UNDEF end if if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units K if((NUM_ENSEMBLE>1) .and. associated(TP1_enstd) .and. associated(TP1_enavg)) then - TP1_enstd = sqrt( TP1_enstd/Nm1 - NdivNm1*(TP1_enavg**2) ) + TP1_enstd = max( sqrt( TP1_enstd/Nm1 - NdivNm1*(TP1_enavg**2) ), 0. ) else if (associated(TP1_enstd)) then TP1_enstd = MAPL_UNDEF end if diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 648df3c3..b93b1980 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -2025,19 +2025,19 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! set export variables - if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc - if(associated(RZMC_ana)) RZMC_ana(:) = cat_diagS_ensavg(:)%rzmc - if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc - if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf - if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K - - if(associated(SFMC_ana_ensstd)) SFMC_ana_ensstd(:) = cat_diagS_ensstd(:)%sfmc - if(associated(RZMC_ana_ensstd)) RZMC_ana_ensstd(:) = cat_diagS_ensstd(:)%rzmc - if(associated(PRMC_ana_ensstd)) PRMC_ana_ensstd(:) = cat_diagS_ensstd(:)%prmc - if(associated(TPSURF_ana_ensstd)) TPSURF_ana_ensstd(:) = cat_diagS_ensstd(:)%tsurf - if(associated(TSOIL1_ana_ensstd)) TSOIL1_ana_ensstd(:) = cat_diagS_ensstd(:)%tp(1) - - if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY + if(associated(SFMC_ana)) SFMC_ana(:) = cat_diagS_ensavg(:)%sfmc + if(associated(RZMC_ana)) RZMC_ana(:) = cat_diagS_ensavg(:)%rzmc + if(associated(PRMC_ana)) PRMC_ana(:) = cat_diagS_ensavg(:)%prmc + if(associated(TPSURF_ana)) TPSURF_ana(:) = cat_diagS_ensavg(:)%tsurf + if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) + MAPL_TICE ! convert to K + + if(associated(SFMC_ana_ensstd)) SFMC_ana_ensstd(:) = max( cat_diagS_ensstd(:)%sfmc , 0. ) + if(associated(RZMC_ana_ensstd)) RZMC_ana_ensstd(:) = max( cat_diagS_ensstd(:)%rzmc , 0. ) + if(associated(PRMC_ana_ensstd)) PRMC_ana_ensstd(:) = max( cat_diagS_ensstd(:)%prmc , 0. ) + if(associated(TPSURF_ana_ensstd)) TPSURF_ana_ensstd(:) = max( cat_diagS_ensstd(:)%tsurf , 0. ) + if(associated(TSOIL1_ana_ensstd)) TSOIL1_ana_ensstd(:) = max( cat_diagS_ensstd(:)%tp(1) , 0. ) + + if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY deallocate(cat_progn_tmp) deallocate(cat_diagS) From 42f5e71f2766370e7fddd1679662f77c0d074b24 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 20 Mar 2023 15:49:54 -0400 Subject: [PATCH 08/10] facilitate HISTORY output of ensstd even when ensavg is not written out; rearrange do loop in GEOS_LandAssimGridComp.F90 --- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 10 +++---- .../GEOS_LandAssimGridComp.F90 | 29 ++++++++++++------- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index 40f8170d..6a769f74 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -3165,7 +3165,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(HLATN_enavg)) HLATN_enavg = 0.0 if(associated(QINFIL_enavg)) QINFIL_enavg = 0.0 if(associated(GHFLX_enavg)) GHFLX_enavg = 0.0 - if(associated(TPSURF_enavg)) TPSURF_enavg = 0.0 + if(associated(TPSURF_enavg) .or. associated(TPSURF_enstd)) TPSURF_enavg = 0.0 if(associated(TPSURF_enstd)) TPSURF_enstd = 0.0 if(associated(TPSNOW_enavg)) TPSNOW_enavg = 0.0 if(associated(TPUNST_enavg)) TPUNST_enavg = 0.0 @@ -3182,13 +3182,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET1_enavg)) WET1_enavg = 0.0 if(associated(WET2_enavg)) WET2_enavg = 0.0 if(associated(WET3_enavg)) WET3_enavg = 0.0 - if(associated(WCSF_enavg)) WCSF_enavg = 0.0 + if(associated(WCSF_enavg) .or. associated(WCSF_enstd)) WCSF_enavg = 0.0 if(associated(WCSF_enstd)) WCSF_enstd = 0.0 - if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 + if(associated(WCRZ_enavg) .or. associated(WCRZ_enstd)) WCRZ_enavg = 0.0 if(associated(WCRZ_enstd)) WCRZ_enstd = 0.0 - if(associated(WCPR_enavg)) WCPR_enavg = 0.0 + if(associated(WCPR_enavg) .or. associated(WCPR_enstd)) WCPR_enavg = 0.0 if(associated(WCPR_enstd)) WCPR_enstd = 0.0 - if(associated(TP1_enavg)) TP1_enavg = 0.0 + if(associated(TP1_enavg) .or. associated(TP1_enstd)) TP1_enavg = 0.0 if(associated(TP1_enstd)) TP1_enstd = 0.0 if(associated(TP2_enavg)) TP2_enavg = 0.0 if(associated(TP3_enavg)) TP3_enavg = 0.0 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index b93b1980..fd84a8f8 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -2004,24 +2004,31 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! finalize ensemble average and standard deviation - Nm1 = real(NUM_ENSEMBLE-1) + if (NUM_ENSEMBLE > 1) then - if (NUM_ENSEMBLE > 1) NdivNm1 = real(NUM_ENSEMBLE)/Nm1 + Nm1 = real(NUM_ENSEMBLE-1) + + NdivNm1 = real(NUM_ENSEMBLE)/Nm1 - do ii=1,N_catl - - cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg + do ii=1,N_catl + + cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg - if (NUM_ENSEMBLE > 1 ) then - cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii)) ) + + end do + + else ! NUM_ENSEMBLE = 1 - else ! NUM_ENSEMBLE = 1 - + ! no need to normalize ens avg, set ens std to undef + + do ii=1,N_catl + cat_diagS_ensstd(ii) = MAPL_UNDEF + + end do - end if - end do + end if ! set export variables From f6a1b71be51578ab73e32d7ac54b637e60899338 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 22 Mar 2023 11:34:52 -0400 Subject: [PATCH 09/10] undoing previous commit which failed to facilitate ensstd HISTORY output when ensavg HISTORY output is not also requested --- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index 6a769f74..40f8170d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -3165,7 +3165,7 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(HLATN_enavg)) HLATN_enavg = 0.0 if(associated(QINFIL_enavg)) QINFIL_enavg = 0.0 if(associated(GHFLX_enavg)) GHFLX_enavg = 0.0 - if(associated(TPSURF_enavg) .or. associated(TPSURF_enstd)) TPSURF_enavg = 0.0 + if(associated(TPSURF_enavg)) TPSURF_enavg = 0.0 if(associated(TPSURF_enstd)) TPSURF_enstd = 0.0 if(associated(TPSNOW_enavg)) TPSNOW_enavg = 0.0 if(associated(TPUNST_enavg)) TPUNST_enavg = 0.0 @@ -3182,13 +3182,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WET1_enavg)) WET1_enavg = 0.0 if(associated(WET2_enavg)) WET2_enavg = 0.0 if(associated(WET3_enavg)) WET3_enavg = 0.0 - if(associated(WCSF_enavg) .or. associated(WCSF_enstd)) WCSF_enavg = 0.0 + if(associated(WCSF_enavg)) WCSF_enavg = 0.0 if(associated(WCSF_enstd)) WCSF_enstd = 0.0 - if(associated(WCRZ_enavg) .or. associated(WCRZ_enstd)) WCRZ_enavg = 0.0 + if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 if(associated(WCRZ_enstd)) WCRZ_enstd = 0.0 - if(associated(WCPR_enavg) .or. associated(WCPR_enstd)) WCPR_enavg = 0.0 + if(associated(WCPR_enavg)) WCPR_enavg = 0.0 if(associated(WCPR_enstd)) WCPR_enstd = 0.0 - if(associated(TP1_enavg) .or. associated(TP1_enstd)) TP1_enavg = 0.0 + if(associated(TP1_enavg)) TP1_enavg = 0.0 if(associated(TP1_enstd)) TP1_enstd = 0.0 if(associated(TP2_enavg)) TP2_enavg = 0.0 if(associated(TP3_enavg)) TP3_enavg = 0.0 From 89a8cc5c34b8a8e21b91a06b9bdafc7ac0bdd2fb Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 22 Mar 2023 11:58:17 -0400 Subject: [PATCH 10/10] added rudimentary documentation of ENSSTD output into GEOSldas_HIST.rc --- src/Applications/LDAS_App/GEOSldas_HIST.rc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Applications/LDAS_App/GEOSldas_HIST.rc b/src/Applications/LDAS_App/GEOSldas_HIST.rc index ba9fb5c4..adb52654 100644 --- a/src/Applications/LDAS_App/GEOSldas_HIST.rc +++ b/src/Applications/LDAS_App/GEOSldas_HIST.rc @@ -407,8 +407,10 @@ COLLECTIONS: # For lndfcstana, *.frequency and *.ref_time must be consistent with the LDAS.rc resource # parameters LANDASSIM_DT and LANDASSIM_T0. +# Output of the ensemble std-dev (ENSSTD) requires simultaneous output of the ensemble mean. If the +# ensemble mean is not written out for a given field, that field's ENSSTD output will be MAPL_UNDEF. - inst3_1d_lndfcstana_Nt.descr: 'Tile-space,3-Hourly,Instantaneous,Single-Level,Assimilation,Ensemble-Average Land Forecast and Analysis Diagnostics', + inst3_1d_lndfcstana_Nt.descr: 'Tile-space,3-Hourly,Instantaneous,Single-Level,Assimilation,Ensemble Land Forecast and Analysis Diagnostics', inst3_1d_lndfcstana_Nt.nbits: 12, inst3_1d_lndfcstana_Nt.template: '%y4%m2%d2_%h2%n2z.bin', inst3_1d_lndfcstana_Nt.mode: 'instantaneous', @@ -438,8 +440,10 @@ COLLECTIONS: # For lndfcstana, *.frequency and *.ref_time must be consistent with the LDAS.rc resource # parameters LANDASSIM_DT and LANDASSIM_T0. +# Output of the ensemble std-dev (ENSSTD) requires simultaneous output of the ensemble mean. If the +# ensemble mean is not written out for a given field, that field's ENSSTD output will be MAPL_UNDEF. - inst3_2d_lndfcstana_Nx.descr: '2d,3-Hourly,Instantaneous,Single-Level,Assimilation,Ensemble-Average Land Forecast and Analysis Diagnostics', + inst3_2d_lndfcstana_Nx.descr: '2d,3-Hourly,Instantaneous,Single-Level,Assimilation,Ensemble Land Forecast and Analysis Diagnostics', inst3_2d_lndfcstana_Nx.nbits: 12, inst3_2d_lndfcstana_Nx.template: '%y4%m2%d2_%h2%n2z.nc4', inst3_2d_lndfcstana_Nx.archive: '%c/Y%y4',