Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Soil temperatures (TP's) in units of Kelvin throughout #329

Closed
wants to merge 3 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Externals.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -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_2
protocol = git
sparse = ../../../../config/GEOSgcm_GridComp_ldas.sparse

Expand Down
2 changes: 1 addition & 1 deletion components.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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_2
develop: develop
13 changes: 7 additions & 6 deletions src/Applications/LDAS_App/tile_bin2nc4.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down Expand Up @@ -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'
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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, &
Expand All @@ -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 !!!

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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,:)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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))

Expand Down Expand Up @@ -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

Expand All @@ -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.

Expand Down Expand Up @@ -3474,7 +3478,10 @@ 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.) ! = 0.2 + MAPL_TICE ! [KELVIN]
!real, parameter :: tp1_threshold = -HUGE(1.) + MAPL_TICE ! = 0.2 + MAPL_TICE ! [KELVIN]

integer :: n, n_e, kk, ii

Expand Down Expand Up @@ -3615,7 +3622,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)), &
Expand Down Expand Up @@ -4282,7 +4291,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 ) )
Expand Down Expand Up @@ -4418,18 +4429,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 )

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down