From 9c6db55e6b6bd87128c379c8306eb68d9f261f68 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Sun, 25 Oct 2020 21:08:15 -0400 Subject: [PATCH 1/2] Changing TP's from C to K --- Externals.cfg | 2 +- components.yaml | 2 +- src/Applications/LDAS_App/tile_bin2nc4.F90 | 13 ++--- .../GEOSens_GridComp/GEOS_EnsGridComp.F90 | 13 ++--- .../GEOS_LandAssimGridComp.F90 | 11 ++-- .../clsm_ensdrv_drv_routines.F90 | 3 ++ .../clsm_ensupd_enkf_update.F90 | 4 +- .../clsm_ensupd_upd_routines.F90 | 53 +++++++++++++------ .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 3 +- 9 files changed, 70 insertions(+), 34 deletions(-) diff --git a/Externals.cfg b/Externals.cfg index e0632a1d..89c88100 100644 --- a/Externals.cfg +++ b/Externals.cfg @@ -32,7 +32,7 @@ protocol = git required = True repo_url = git@github.com:GEOS-ESM/GEOSgcm_GridComp.git local_path = ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp -branch = develop +branch = feature/borescan_TP_C_to_K protocol = git sparse = ../../../../config/GEOSgcm_GridComp_ldas.sparse diff --git a/components.yaml b/components.yaml index a81b09a2..1ff7dc33 100644 --- a/components.yaml +++ b/components.yaml @@ -32,5 +32,5 @@ GEOSgcm_GridComp: local: ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp remote: git@github.com:GEOS-ESM/GEOSgcm_GridComp.git sparse: ./config/GEOSgcm_GridComp_ldas.sparse - branch: develop + branch: feature/borescan_TP_C_to_K develop: develop diff --git a/src/Applications/LDAS_App/tile_bin2nc4.F90 b/src/Applications/LDAS_App/tile_bin2nc4.F90 index 9143ea38..b9ce71cf 100644 --- a/src/Applications/LDAS_App/tile_bin2nc4.F90 +++ b/src/Applications/LDAS_App/tile_bin2nc4.F90 @@ -246,6 +246,7 @@ FUNCTION getAttribute (SHORT_NAME, LNAME, UNT) result (str_atr) case ('TB_LAND_1410MHZ_40DEG_VPOL'); LONG_NAME = 'brightness_temperature_land_1410MHz_40deg_Vpol'; UNITS = 'K' ! Done for SM_L4 + ! Bilja temp marker so I know I have changed something here C to K units case ('Tair'); LONG_NAME = 'air_temperature_at_RefH'; UNITS = 'K' case ('TA'); LONG_NAME = 'air_temperature_at_RefH'; UNITS = 'K' @@ -286,12 +287,12 @@ FUNCTION getAttribute (SHORT_NAME, LNAME, UNT) result (str_atr) case ('TPSURF'); LONG_NAME = 'ave_catchment_temp_incl_snw'; UNITS = 'K' case ('GRN'); LONG_NAME = 'greeness_fraction'; UNITS = '1' case ('LAI'); LONG_NAME = 'leaf_area_index'; UNITS = '1' - case ('TP1'); LONG_NAME = 'soil_temperatures_layer_1'; UNITS = 'C' - case ('TP2'); LONG_NAME = 'soil_temperatures_layer_2'; UNITS = 'C' - case ('TP3'); LONG_NAME = 'soil_temperatures_layer_3'; UNITS = 'C' - case ('TP4'); LONG_NAME = 'soil_temperatures_layer_4'; UNITS = 'C' - case ('TP5'); LONG_NAME = 'soil_temperatures_layer_5'; UNITS = 'C' - case ('TP6'); LONG_NAME = 'soil_temperatures_layer_6'; UNITS = 'C' + case ('TP1'); LONG_NAME = 'soil_temperatures_layer_1'; UNITS = 'K' + case ('TP2'); LONG_NAME = 'soil_temperatures_layer_2'; UNITS = 'K' + case ('TP3'); LONG_NAME = 'soil_temperatures_layer_3'; UNITS = 'K' + case ('TP4'); LONG_NAME = 'soil_temperatures_layer_4'; UNITS = 'K' + case ('TP5'); LONG_NAME = 'soil_temperatures_layer_5'; UNITS = 'K' + case ('TP6'); LONG_NAME = 'soil_temperatures_layer_6'; UNITS = 'K' case ('PRECTOTLAND');LONG_NAME = 'Total_precipitation_land'; UNITS = 'kg m-2 s-1' case ('PRECSNOLAND');LONG_NAME = 'snowfall_land'; UNITS = 'kg m-2 s-1' case ('SNOWMASS') ;LONG_NAME = 'snow_mass'; UNITS = 'kg m-2' diff --git a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 index 46fd1115..3bdca6ca 100644 --- a/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSens_GridComp/GEOS_EnsGridComp.F90 @@ -3257,12 +3257,13 @@ subroutine Collect_land_ens(gc, import, export, clock, rc) if(associated(WCSF_enavg)) WCSF_enavg = WCSF_enavg/NUM_ENSEMBLE if(associated(WCRZ_enavg)) WCRZ_enavg = WCRZ_enavg/NUM_ENSEMBLE if(associated(WCPR_enavg)) WCPR_enavg = WCPR_enavg/NUM_ENSEMBLE - if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K - if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K - if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K - if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K - if(associated(TP5_enavg)) TP5_enavg = TP5_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K - if(associated(TP6_enavg)) TP6_enavg = TP6_enavg/NUM_ENSEMBLE + MAPL_TICE ! convert to K +! *** Bilja removed C to K conversion + if(associated(TP1_enavg)) TP1_enavg = TP1_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K + if(associated(TP2_enavg)) TP2_enavg = TP2_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K + if(associated(TP3_enavg)) TP3_enavg = TP3_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K + if(associated(TP4_enavg)) TP4_enavg = TP4_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K + if(associated(TP5_enavg)) TP5_enavg = TP5_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K + if(associated(TP6_enavg)) TP6_enavg = TP6_enavg/NUM_ENSEMBLE ! + MAPL_TICE ! convert to K if(associated(EMIS_enavg)) EMIS_enavg = EMIS_enavg/NUM_ENSEMBLE if(associated(ALBVR_enavg)) ALBVR_enavg = ALBVR_enavg/NUM_ENSEMBLE if(associated(ALBVF_enavg)) ALBVF_enavg = ALBVF_enavg/NUM_ENSEMBLE diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index fb1cb5a2..c3d82528 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1902,7 +1902,10 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) 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 + ! *** Bilja In recompute_diagS (few lines above) cat_diagS%tp is returned + ! here in K. This makes cat_diagS_ensavg(:)%tp be in K as well. So, no + ! need to transfer to K + if(associated(TSOIL1_ana)) TSOIL1_ana(:) = cat_diagS_ensavg(:)%tp(1) ! + MAPL_TICE ! convert to K deallocate(cat_progn_tmp) deallocate(cat_diagS) @@ -2199,7 +2202,9 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) _VERIFY(status) ! convert Catchment model variables into inputs suitable for the mwRTM - ! NOTE: input TP1 must be in degree Celsius! + ! Bilja 'TP1' now should be imported in K, meaning I it will be sent to + ! catch2mwRTM_vars as such. The catch2mwRTM_vars is adjusted to deal with + ! tp in Kelving allocate(sfmc_mwRTM(N_catl), tsoil_mwRTM (N_catl)) call catch2mwRTM_vars( & N_catl, & @@ -2208,7 +2213,7 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) mwRTM_param%poros, & WCSF, & TPSURF, & - TP1, & ! units deg C !!! + TP1, & ! units deg K !!! sfmc_mwRTM, & tsoil_mwRTM ) ! units Kelvin !!! diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensdrv_drv_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensdrv_drv_routines.F90 index b1ec9d14..22b81c25 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensdrv_drv_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensdrv_drv_routines.F90 @@ -239,6 +239,9 @@ subroutine recompute_diagS( N_catd, cat_param, cat_progn, cat_diagS ) call catch_calc_tp( N_catd, cat_param%poros, & catprogn2ghtcnt(N_catd,cat_progn), tp ) + ! Bilja *** if catch_calc_tp returns Kelvin, then the cat_diagS will + ! become Kelvin as well. I made catch_calc_tp to get me 'tp' in Kelvin, + ! so cat_diagS(:)%tp(ii) = tp(ii,:) will leave this subroutine in K do ii=1,N_gt cat_diagS(:)%tp(ii) = tp(ii,:) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index f317e3c7..f459fa9c 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -2853,11 +2853,13 @@ subroutine write_smapL4SMaup( option, date_time, work_path, exp_id, N_ens, & tsurf(:,n_e) ) ! NOTE: "tp" is returned in CELSIUS [for consistency w/ catchment.F90] + ! *** Bilja I've changed catch_calc_tp to return tp in Kelvin call catch_calc_tp( N_catl, cat_param%poros, & catprogn2ghtcnt(N_catl,cat_progn(:,n_e)), tp ) - tp1(:,n_e) = tp(1,:) + MAPL_TICE + ! *** Bilja tp returned in Kelvin, removing addition of TICE + tp1(:,n_e) = tp(1,:) ! + MAPL_TICE end do diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 3b0631a1..2e3a4b2a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -1339,6 +1339,7 @@ subroutine get_obs_pred( & if (get_tp_l) then ! NOTE: "tp" is returned in CELSIUS [for consistency w/ catchment.F90] + ! *** Bilja changed "tp" in catch_calc_tp to be retured in KELVIN ! updated to new interface - reichle, 3 Apr 2012 @@ -1358,6 +1359,8 @@ subroutine get_obs_pred( & ar1_l, ar2_l, ar4_l, tsurf_excl_snow ) ! catch_calc_FT() expects "tp" in CELSIUS + ! *** Bilja changed catch_calc_FT to expect tp in Kelvin on input + ! here tp_l is in Kelvin call catch_calc_FT( N_catl, asnow, tp_l(1,:), tsurf_excl_snow, FT_l(:,n_e)) @@ -2272,7 +2275,7 @@ subroutine qc_model_based_for_sat_tsurf( N_cat, precip, SWE, tp1, & real, dimension(N_cat), intent(in) :: precip ! Rainf+Snowf [kg/m2/s] real, dimension(N_cat), intent(in) :: SWE ! total SWE [kg/m2] - real, dimension(N_cat), intent(in) :: tp1 ! soil temperature [C] + real, dimension(N_cat), intent(in) :: tp1 ! soil temperature [K] real, dimension(N_cat), intent(inout) :: tsurf @@ -2290,20 +2293,21 @@ subroutine qc_model_based_for_sat_tsurf( N_cat, precip, SWE, tp1, & logical, dimension(N_cat) :: frozen - real, dimension(N_cat) :: tp1_in_Kelvin + !real, dimension(N_cat) :: tp1_in_Kelvin ! --------------------------------------- - + ! **** Bilja elimiated tp1_in_Kelvin variable; tp1 now arrives in [K] frozen = .false. if (avoid_frozen) then - tp1_in_Kelvin = tp1 + MAPL_TICE + !tp1_in_Kelvin = tp1 + MAPL_TICE do i=1,N_cat if ( (tsurf( i) < temp_threshold) .or. & - (tp1_in_Kelvin( i) < temp_threshold) & + (tp1 ( i) < temp_threshold) & + !(tp1_in_Kelvin( i) < temp_threshold) & ) & frozen(i) = .true. @@ -3474,7 +3478,9 @@ subroutine cat_enkf_increments( & real, parameter :: SWE_threshold = +HUGE(1.) ! = 1.e-4 ! [kg/m2] - real, parameter :: tp1_threshold = -HUGE(1.) ! = 0.2 ! [CELSIUS] + ! Bilja *** tp is in Kelvin, which makes tp1_ensavg to be in Kelvin, so + ! this threshold has to be in K too (adding MAPL_TICE) + real, parameter :: tp1_threshold = -HUGE(1.) + MAPL_TICE ! = 0.2 + MAPL_TICE ! [KELVIN] integer :: n, n_e, kk, ii @@ -3615,7 +3621,9 @@ subroutine cat_enkf_increments( & ! soil temperature ! - ! NOTE: "tp" is returned in CELSIUS + ! *** Bilja changed catch_calc_tp to return Kelvin + ! NOTE: "tp" is returned in KELVIN + ! NOTE: fice is fractional ice content call catch_calc_tp( N_catd, cat_param%poros, & catprogn2ghtcnt(N_catd,cat_progn(:,n_e)), & @@ -4282,7 +4290,9 @@ subroutine cat_enkf_increments( & do n_e=1,N_ens - ! tp must be in CELSIUS; FT_Teff is returned in Kelvin + ! *** Bilja changed catch_calc_FT to expect tp in Kelvin on input. + ! Here tp IS in Kelvin + ! tp must be in KELVIN; FT_Teff is returned in Kelvin call catch_calc_FT( N_catd, asnow(:,n_e), tp(1,:,n_e), & tsurf_excl_snow(:,n_e), FT_state(:), FT_Teff(:,n_e ) ) @@ -4418,18 +4428,31 @@ subroutine cat_enkf_increments( & fice_minus = fice(1,kk,n_e) ! model forecast - tp1_minus = tp(1,kk,n_e) ! model forecast [CELSIUS] - + ! *** Bilja tp is in Kelvin here. I need to + ! change the following lines to account for the 273 + ! difference between C and K + tp1_minus = tp(1,kk,n_e) ! model forecast [Kelvin] + ! [K] <---- [K] fice_plus = fice_minus ! ice fraction does not change - - tp1_plus = tp1_minus + deltaT ! tentative tp1 analysis [CELSIUS] - + ! [-] [-] + tp1_plus = tp1_minus + deltaT ! tentative tp1 analysis [Kelvin] + ! [K] <---- [K] [K or C doesnt matter] ! avoid phase change of soil temp - if ((tp1_minus*tp1_plus) < 0.) tp1_plus = 0. - + ! *** Bilja: the following line says: if model forcast (tp) + ! and tentative tp1 analysis are of opposite sign, then make + ! tentative analysis to be equal to 0 deg Celsius before sending it to catch_calc_ght + ! if ((tp1_minus*tp1_plus) < 0.) tp1_plus = 0. + ! To rewrite this condition using Kelvin I do the following: + if ((tp1_minus .gt. MAPL_TICE .and. tp1_plus .lt. MAPL_TICE) .or. & + (tp1_minus .lt. MAPL_TICE .and. tp1_plus .gt. MAPL_TICE)) & + tp1_plus = MAPL_TICE + ! compute ght1_plus from tp1_plus and fice_plus + ! *** Bilja I changed the catch_calc_ght to expect Kelvin + ! (due to its other calls). So, I can send tp1_plus in + ! Kelvins to catch_calc_ght. tp1_plus IS in Kelvin call catch_calc_ght( cat_param(kk)%dzgt(1), & cat_param(kk)%poros, tp1_plus, fice_plus, ght1_plus ) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index a93fd2c2..09406faf 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -222,9 +222,10 @@ subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & ! (change prompted by revision of Catchment model parameter CSOIL_2) ! - reichle, 23 Dec 2015 +! **** Bilja removed C to K conversion ! NOTE: "tp" is in deg Celsius - tsoil_mwRTM = tp1_catch + MAPL_TICE + tsoil_mwRTM = tp1_catch ! + MAPL_TICE end subroutine catch2mwRTM_vars From 935d7e4bedecc9413ddd64841b55c952b9b617f9 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 30 Oct 2020 11:58:40 -0400 Subject: [PATCH 2/2] Possible fix to tp1_threshold --- Externals.cfg | 2 +- components.yaml | 2 +- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Externals.cfg b/Externals.cfg index 89c88100..86673218 100644 --- a/Externals.cfg +++ b/Externals.cfg @@ -32,7 +32,7 @@ protocol = git required = True repo_url = git@github.com:GEOS-ESM/GEOSgcm_GridComp.git local_path = ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp -branch = feature/borescan_TP_C_to_K +branch = feature/borescan_TP_C_to_K_2 protocol = git sparse = ../../../../config/GEOSgcm_GridComp_ldas.sparse diff --git a/components.yaml b/components.yaml index 1ff7dc33..b037c648 100644 --- a/components.yaml +++ b/components.yaml @@ -32,5 +32,5 @@ GEOSgcm_GridComp: local: ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp remote: git@github.com:GEOS-ESM/GEOSgcm_GridComp.git sparse: ./config/GEOSgcm_GridComp_ldas.sparse - branch: feature/borescan_TP_C_to_K + branch: feature/borescan_TP_C_to_K_2 develop: develop diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 2e3a4b2a..dd449b8d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3480,7 +3480,8 @@ subroutine cat_enkf_increments( & ! Bilja *** tp is in Kelvin, which makes tp1_ensavg to be in Kelvin, so ! this threshold has to be in K too (adding MAPL_TICE) - real, parameter :: tp1_threshold = -HUGE(1.) + MAPL_TICE ! = 0.2 + MAPL_TICE ! [KELVIN] + real, parameter :: tp1_threshold = -HUGE(1.) ! = 0.2 + MAPL_TICE ! [KELVIN] + !real, parameter :: tp1_threshold = -HUGE(1.) + MAPL_TICE ! = 0.2 + MAPL_TICE ! [KELVIN] integer :: n, n_e, kk, ii