diff --git a/src/Applications/LDAS_App/GEOSldas_HIST.rc b/src/Applications/LDAS_App/GEOSldas_HIST.rc index 40299e69..adb52654 100644 --- a/src/Applications/LDAS_App/GEOSldas_HIST.rc +++ b/src/Applications/LDAS_App/GEOSldas_HIST.rc @@ -407,29 +407,43 @@ 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', 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' , - '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_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 # 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', @@ -441,16 +455,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 b32e8044..40f8170d 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' ,& @@ -656,6 +657,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'ave_catchment_temp_incl_snw_ensstd',& + UNITS = 'K' ,& + SHORT_NAME = 'TPSURF_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'temperature_top_snow_layer',& UNITS = 'K' ,& @@ -800,6 +810,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_surface_layer_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCSF_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'water_root_zone' ,& UNITS = 'm3 m-3' ,& @@ -809,15 +828,34 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_root_zone_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCRZ_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + + 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 ,& RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'water_ave_prof_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCPR_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'soil_temperatures_layer_1' ,& UNITS = 'K' ,& @@ -827,6 +865,15 @@ subroutine SetServices(gc, rc) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'soil_temperatures_layer_1_ensstd' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TSOIL1TILE_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + call MAPL_AddExportSpec(GC, & LONG_NAME = 'soil_temperatures_layer_2' ,& UNITS = 'K' ,& @@ -2287,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 + 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 @@ -2304,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 - 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_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 @@ -2401,6 +2448,8 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) real, dimension(:), pointer :: CNBURN, CNBURN_enavg real, dimension(:), pointer :: CNFSEL, CNFSEL_enavg + real :: Nm1, NdivNm1 + ! Get my name and setup traceback handle call ESMF_GridCompget(gc, name=comp_name, rc=status) VERIFY_(status) @@ -2859,6 +2908,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_enstd, 'TPSURF_ENSSTD' ,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 +2942,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_enstd, 'WCSF_ENSSTD' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, WCRZ_enavg, 'WCRZ' ,rc=status) VERIFY_(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_enstd, 'WCPR_ENSSTD' ,rc=status) + VERIFY_(status) call MAPL_GetPointer(export, TP1_enavg, 'TSOIL1TILE' ,rc=status) VERIFY_(status) + call MAPL_GetPointer(export, TP1_enstd, 'TSOIL1TILE_ENSSTD' ,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 +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_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 @@ -3123,9 +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_enstd)) WCSF_enstd = 0.0 if(associated(WCRZ_enavg)) WCRZ_enavg = 0.0 + if(associated(WCRZ_enstd)) WCRZ_enstd = 0.0 if(associated(WCPR_enavg)) WCPR_enavg = 0.0 + if(associated(WCPR_enstd)) WCPR_enstd = 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 if(associated(TP4_enavg)) TP4_enavg = 0.0 @@ -3306,6 +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_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)) & @@ -3338,12 +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_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_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_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_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)) & @@ -3544,6 +3618,10 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(collect_land_counter == NUM_ENSEMBLE) then + + 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 @@ -3591,6 +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((NUM_ENSEMBLE>1) .and. associated(TPSURF_enstd) .and. associated(TPSURF_enavg)) then + TPSURF_enstd = max( sqrt( TPSURF_enstd/Nm1 - NdivNm1*(TPSURF_enavg**2) ), 0. ) + 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 @@ -3607,9 +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((NUM_ENSEMBLE>1) .and. associated(WCSF_enstd) .and. associated(WCSF_enavg)) then + 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 = 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(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! units now K, rreichle & borescan, 6 Nov 2020 + if((NUM_ENSEMBLE>1) .and. associated(WCPR_enstd) .and. associated(WCPR_enavg)) then + 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 = max( sqrt( TP1_enstd/Nm1 - NdivNm1*(TP1_enavg**2) ), 0. ) + 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 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 8fc69ef4..e7f4ac86 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 @@ -777,7 +778,16 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) VERIFY_(STATUS) - + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_surface_analysis_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCSF_ANA_ENSSTD' ,& + 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 +797,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_rootzone_analysis_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCRZ_ANA_ENSSTD' ,& + 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 +814,15 @@ subroutine SetServices ( GC, RC ) VLOCATION = MAPL_VLocationNone ,& RC=STATUS ) VERIFY_(STATUS) + + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_moisture_profile_analysis_ensstd' ,& + UNITS = 'm3 m-3' ,& + SHORT_NAME = 'WCPR_ANA_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddExportSpec(GC ,& LONG_NAME = 'ave_catchment_temp_incl_snw_analysis' ,& @@ -805,6 +833,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'ave_catchment_temp_incl_snw_analysis_ensstd' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TPSURF_ANA_ENSSTD' ,& + 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 +851,15 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddExportSpec(GC ,& + LONG_NAME = 'soil_temperatures_layer_1_analysis_ensstd' ,& + UNITS = 'K' ,& + SHORT_NAME = 'TSOIL1_ANA_ENSSTD' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RC=STATUS ) + VERIFY_(STATUS) + ! Exports for microwave radiative transfer model (mwRTM) call MAPL_AddExportSpec(GC ,& @@ -1431,7 +1477,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 @@ -1454,7 +1500,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_ensstd type(obs_type), dimension(:), pointer :: Observations_l => null() logical :: fresh_incr @@ -1481,11 +1527,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_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) @@ -1499,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 :: Nm1, NdivNm1 + #ifdef DBG_LANDASSIM_INPUTS ! vars for debugging purposes type(ESMF_Grid) :: TILEGRID @@ -1643,15 +1697,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) VERIFY_(status) - call MAPL_GetPointer(export, TSOIL1_ana, 'TSOIL1_ANA' ,rc=status) + call MAPL_GetPointer(export, SFMC_ana, 'WCSF_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, SFMC_ana, 'WCSF_ANA' ,rc=status) + call MAPL_GetPointer(export, RZMC_ana, 'WCRZ_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, RZMC_ana, 'WCRZ_ANA' ,rc=status) + call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) VERIFY_(status) - call MAPL_GetPointer(export, PRMC_ana, 'WCPR_ANA' ,rc=status) + call MAPL_GetPointer(export, TPSURF_ana_ensstd, 'TPSURF_ANA_ENSSTD' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, TSOIL1_ana_ensstd, 'TSOIL1_ANA_ENSSTD' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, SFMC_ana_ensstd, 'WCSF_ANA_ENSSTD' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, RZMC_ana_ensstd, 'WCRZ_ANA_ENSSTD' ,rc=status) + VERIFY_(status) + call MAPL_GetPointer(export, PRMC_ana_ensstd, 'WCPR_ANA_ENSSTD' ,rc=status) VERIFY_(status) ! exports for microwave radiative transfer model (mwRTM) @@ -1910,11 +1974,15 @@ 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_ensstd(N_catl)) do ii=1,N_catl - cat_diagS_ensavg(ii) = 0.0 ! initialize ens average + 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) @@ -1926,28 +1994,60 @@ 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 - + end do - do ii=1,N_catl - cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize - end do + ! finalize ensemble average and standard deviation + + if (NUM_ENSEMBLE > 1) then - ! set export variables + Nm1 = real(NUM_ENSEMBLE-1) + + NdivNm1 = real(NUM_ENSEMBLE)/Nm1 - 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 + do ii=1,N_catl + + cat_diagS_ensavg(ii) = cat_diagS_ensavg(ii)/real(NUM_ENSEMBLE) ! normalize --> ens avg + + 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 + + ! no need to normalize ens avg, set ens std to undef - if(associated(MWRTM_VEGOPACITY)) MWRTM_VEGOPACITY(:) = mwRTM_param(:)%VEGOPACITY + do ii=1,N_catl + + cat_diagS_ensstd(ii) = MAPL_UNDEF + + end do + + end if + + ! 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(:) = 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) - deallocate(cat_diagS_ensavg) + deallocate(cat_diagS_ensavg) + deallocate(cat_diagS_ensstd) ! write analysis fields into SMAP L4_SM aup file ! whenever it was time for assimilation (regardless diff --git a/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 b/src/Components/GEOSldas_GridComp/Shared/catch_types.F90 index 0c97fc9b..df22a2b3 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 scalar_mult_cat_diagS + module procedure scalar_mult_cat_diagF + 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 scalar_mult_cat_diagS( scalar, cat_diagS ) implicit none - type(cat_diagS_type) :: cat_diagS_div_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_div_scalar%ar1 = cat_diagS%ar1 / scalar - cat_diagS_div_scalar%ar2 = cat_diagS%ar2 / 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 scalar_mult_cat_diagS + + ! ----------------------------------------------------------- + + 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 + + ! ----------------------------------------------------------- + + function scalar_mult_cat_diagF( scalar, cat_diagF ) + + implicit none + + type(cat_diagF_type) :: scalar_mult_cat_diagF + type(cat_diagF_type), intent(in) :: cat_diagF + + real, intent(in) :: 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 - 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_diagS_div_scalar( cat_diagS, scalar ) - cat_diagS_div_scalar%tsurf = cat_diagS%tsurf / 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,112 @@ 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_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 )