From 8c116cbee0a21f836d7ca3fea34a22c693d36af4 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Fri, 5 May 2023 10:13:15 -0400 Subject: [PATCH 1/6] update mwRTM_get_Tb --- .../GEOS_LandAssimGridComp.F90 | 1 + .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 109 ++++++++++++++++-- 2 files changed, 101 insertions(+), 9 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index e7f4ac86..843c081d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -2386,6 +2386,7 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) SWE, & dummy_real, & ! intent(in), "Tair", not used as long as "incl_atm_terms=.false." incl_atm_terms, & + 'mironov', & Tb_h_tmp, Tb_v_tmp ) ! intent(out) 'TB_LAND_1410MHZ_40DEG_HPOL', 'TB_LAND_1410MHZ_40DEG_VPOL' deallocate(dummy_real) else diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index 9c299dab..8b472f03 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -193,7 +193,8 @@ module mwRTM_routines ! **************************************************************** subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & - sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, tp1_in_Kelvin ) + sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, & + sfmc_cat2rtm_scale,sfmc_cat2rtm_offset, tp1_in_Kelvin ) ! convert soil moisture, surface temperature, and soil temperature from the Catchment ! model into soil moisture and soil temperature inputs for the microwave radiative @@ -216,6 +217,8 @@ subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & real, dimension(N_tile), intent(out) :: sfmc_mwRTM, tsoil_mwRTM logical, intent(in), optional :: tp1_in_Kelvin + real, dimension(N_tile), intent(in), optional :: sfmc_cat2rtm_scale + real, dimension(N_tile), intent(in), optional :: sfmc_cat2rtm_offset ! local variables @@ -230,6 +233,12 @@ subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & sfmc_mwRTM = sfmc_catch * poros_mwRTM / poros_catch + if (present(sfmc_cat2rtm_scale)) & + sfmc_mwRTM = sfmc_mwRTM * sfmc_cat2rtm_scale + sfmc_cat2rtm_offset + + where (sfmc_mwRTM < 1.e-2) sfmc_mwRTM = 1.e-2 + where (sfmc_mwRTM >= poros_mwRTM) sfmc_mwRTM = poros_mwRTM + ! diagnose soil temperature to be used with mwRTM ! (change prompted by revision of Catchment model parameter CSOIL_2) ! - reichle, 23 Dec 2015 @@ -385,10 +394,12 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & soiltemp_in_C = soiltemp(n) - MAPL_TICE - CALL DIELWANG (freq, soilmoist(n), soiltemp_in_C, & - mwp(n)%wang_wt, mwp(n)%wang_wp, mwp(n)%poros, & - mwp(n)%sand, mwp(n)%clay, c_er ) - + !CALL DIELWANG (freq, soilmoist(n), soiltemp_in_C, & + ! mwp(n)%wang_wt, mwp(n)%wang_wp, mwp(n)%poros, & + ! mwp(n)%sand, mwp(n)%clay, c_er ) + + CALL mironov(freq,soilmoist(n),mwp(n)%clay,c_er) + ! soil reflectivity for smooth surface based on dielect const. (Fresnel) tmpc1 = SQRT(c_er - sin_inc**2) @@ -431,8 +442,8 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & if (freq < 2.e9) then - Q = 0. ! Q is assumed zero at low frequency - + !Q = 0. ! Q is assumed zero at low frequency + Q = 0.1771 * h_mc else Q = 0.35 * (1.0 - exp(-0.6 * (mwp(n)%rgh_polmix**2) * (freq/1.e9) )) @@ -441,9 +452,11 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & ! rough surface reflectivity - rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrh) - rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrv) + !rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrh) + !rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrv) + rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**2) + rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**2) ! ------------------------------------------------------------- ! @@ -854,6 +867,84 @@ END SUBROUTINE DIEL_WAT ! ********************************************************************** + SUBROUTINE mironov(freq,mv,clayfrac,er_r) + + + IMPLICIT NONE + REAL, INTENT(IN) :: freq ! microwave frequency [Hz] + REAL, INTENT(IN) :: mv ! volumetric soil water content [m3/m3] + REAL, INTENT(IN) :: clayfrac ! clay fraction [0-1] + + COMPLEX, INTENT(OUT) :: er_r + + !REAL, PARAMETER :: PI = acos(-1.0) + REAL :: PI, f, C + REAL :: nd, kd, mvt, eps0b, eps0u, taub, tauu, sigb, sigu + REAL :: eps0, epsinf, epsb_real, epsu_real, epsb_imag, epsu_imag + REAL :: nb, kb, nu, ku + REAL :: nm, km + REAL :: er_r_real, er_r_imag + + PI = MAPL_PI + + f = freq !*1e9 ! Section IV + C = clayfrac*100 ! Section VI + + !! Mironov's regression expressions based on Curtis, Dobson, and Hallikainen datasets + + nd = 1.634 - 0.539e-2 * C + 0.2748e-4 * C**2 ! Eqn 17 + kd = 0.03952 - 0.04038e-2 * C ! Eqn 18 + mvt = 0.02863 + 0.30673e-2 * C ! Eqn 19 + eps0b = 79.8 - 85.4e-2 * C + 32.7e-4 * C**2 ! Eqn 20 + taub = 1.062e-11 + 3.450e-12 * 1e-2 * C ! Eqn 21 + sigb = 0.3112 + 0.467e-2 * C ! Eqn 22 + sigu = 0.3631 + 1.217e-2 * C ! Eqn 23 + eps0u = 100 ! Eqn 24 + tauu = 8.5e-12 ! Eqn 25 + + !! Debye relaxation equations for water as a function of frequency + + eps0 = 8.854e-12 ! Vacuum permittivity + epsinf = 4.9 ! Section IV + epsb_real = epsinf + ( (eps0b - epsinf)/(1 + (2*PI*f*taub)**2) ) ! Eqn 16 + epsb_imag = (eps0b - epsinf)/(1 + (2*PI*f*taub)**2) * (2*PI*f*taub) + sigb/(2*PI*eps0*f) ! Eqn 16 + epsu_real = epsinf + ( (eps0u - epsinf)/(1 + (2*PI*f*tauu)**2) ) ! Eqn 16 + epsu_imag = (eps0u - epsinf)/(1 + (2*PI*f*tauu)**2) * (2*PI*f*tauu) + sigu/(2*PI*eps0*f) ! Eqn 16 + + !! Refractive indices + + nb = 1/sqrt(2.0) * sqrt( sqrt(epsb_real**2 + epsb_imag**2) + epsb_real ) ! Eqn 14 + kb = 1/sqrt(2.0) * sqrt( sqrt(epsb_real**2 + epsb_imag**2) - epsb_real ) ! Eqn 14 + nu = 1/sqrt(2.0) * sqrt( sqrt(epsu_real**2 + epsu_imag**2) + epsu_real ) ! Eqn 14 + ku = 1/sqrt(2.0) * sqrt( sqrt(epsu_real**2 + epsu_imag**2) - epsu_real ) ! Eqn 14 + + !! n(*) are refractive indices, k(*) are normalized attenuation coefficients + !! m: moist soil + !! d: dry soil + !! b: bound soil water (BSW) + !! u: unbound (free) soil water (FSW) + + IF (mv <= mvt) THEN + + nm = nd + (nb - 1) * mv ! Eqn 12 + km = kd + kb * mv ! Eqn 13 + er_r_real = nm**2 - km**2 ! Eqn 11 + er_r_imag = 2 * nm * km ! Eqn 11 + + ELSE + + nm = nd + (nb - 1) * mvt + (nu - 1) * (mv - mvt) ! Eqn 12 + km = kd + kb * mvt + ku * (mv - mvt) ! Eqn 13 + er_r_real = nm**2 - km**2 ! Eqn 11 + er_r_imag = 2 * nm * km ! Eqn 11 + + ENDIF + + er_r = CMPLX(er_r_real,er_r_imag) + + END SUBROUTINE mironov + + ! ********************************************************************** end module mwRTM_routines From 24ae9a9501de83a36ac132bb7381aeea281442f3 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Fri, 5 May 2023 10:15:13 -0400 Subject: [PATCH 2/6] update mwRTM_get_Tb --- .../clsm_ensupd_upd_routines.F90 | 12 +++- .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 68 +++++++++++++------ 2 files changed, 57 insertions(+), 23 deletions(-) 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 2010d18b..7e3a3681 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 @@ -1375,14 +1375,22 @@ subroutine get_obs_pred( & call mwRTM_get_Tb( & N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .true., & - Tb_h_vec, Tb_v_vec ) + 'wang',Tb_h_vec, Tb_v_vec ) case (2) call mwRTM_get_Tb( & N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .false., & - Tb_h_vec, Tb_v_vec ) + 'wang',Tb_h_vec, Tb_v_vec ) + + case (3) + + call mwRTM_get_Tb( & + N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & + lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .false., & + 'mironov',Tb_h_vec, Tb_v_vec ) + case default diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index 8b472f03..ef8da3ed 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -193,8 +193,8 @@ module mwRTM_routines ! **************************************************************** subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & - sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, & - sfmc_cat2rtm_scale,sfmc_cat2rtm_offset, tp1_in_Kelvin ) + sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, & + tp1_in_Kelvin ) ! convert soil moisture, surface temperature, and soil temperature from the Catchment ! model into soil moisture and soil temperature inputs for the microwave radiative @@ -217,8 +217,6 @@ subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & real, dimension(N_tile), intent(out) :: sfmc_mwRTM, tsoil_mwRTM logical, intent(in), optional :: tp1_in_Kelvin - real, dimension(N_tile), intent(in), optional :: sfmc_cat2rtm_scale - real, dimension(N_tile), intent(in), optional :: sfmc_cat2rtm_offset ! local variables @@ -232,13 +230,7 @@ subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & ! proper functioning of mwRTM calibration to SMOS obs sfmc_mwRTM = sfmc_catch * poros_mwRTM / poros_catch - - if (present(sfmc_cat2rtm_scale)) & - sfmc_mwRTM = sfmc_mwRTM * sfmc_cat2rtm_scale + sfmc_cat2rtm_offset - where (sfmc_mwRTM < 1.e-2) sfmc_mwRTM = 1.e-2 - where (sfmc_mwRTM >= poros_mwRTM) sfmc_mwRTM = poros_mwRTM - ! diagnose soil temperature to be used with mwRTM ! (change prompted by revision of Catchment model parameter CSOIL_2) ! - reichle, 23 Dec 2015 @@ -261,7 +253,7 @@ end subroutine catch2mwRTM_vars subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & LAI, soilmoist, soiltemp, SWE, Tair, incl_atm_terms, & - Tb_h, Tb_v ) + soildiel_model, Tb_h, Tb_v ) !--------------------------------------------------- !RTM adapted from Steven Chan @@ -319,6 +311,7 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & real, dimension(N_tile), intent(in) :: Tair ! [K] logical, intent(in) :: incl_atm_terms + character(len=*), intent(in) :: soildiel_model real, dimension(N_tile), intent(out) :: Tb_h, Tb_v ! [K] @@ -393,12 +386,21 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & ! soil dielectric constant soiltemp_in_C = soiltemp(n) - MAPL_TICE - - !CALL DIELWANG (freq, soilmoist(n), soiltemp_in_C, & - ! mwp(n)%wang_wt, mwp(n)%wang_wp, mwp(n)%poros, & - ! mwp(n)%sand, mwp(n)%clay, c_er ) + select case (soildiel_model) + + case('wang') + + CALL DIELWANG (freq, soilmoist(n), soiltemp_in_C, & + mwp(n)%wang_wt, mwp(n)%wang_wp, mwp(n)%poros, & + mwp(n)%sand, mwp(n)%clay, c_er ) + + case('mironov') CALL mironov(freq,soilmoist(n),mwp(n)%clay,c_er) + + case default + CALL mironov(freq,soilmoist(n),mwp(n)%clay,c_er) + end select ! soil reflectivity for smooth surface based on dielect const. (Fresnel) @@ -441,9 +443,21 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & ! 2) polarization mixing, Q as defined in CMEM: if (freq < 2.e9) then - - !Q = 0. ! Q is assumed zero at low frequency + + select case (soildiel_model) + + case('wang') + + Q = 0. ! Q is assumed zero at low frequency + + case('mironov') Q = 0.1771 * h_mc + + case default + Q = 0.1771 * h_mc + + end select + else Q = 0.35 * (1.0 - exp(-0.6 * (mwp(n)%rgh_polmix**2) * (freq/1.e9) )) @@ -452,12 +466,24 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & ! rough surface reflectivity - !rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrh) - !rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrv) - + select case (soildiel_model) + + case('wang') + + rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrh) + rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**mwp(n)%rgh_nrv) + + case('mironov') + rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**2) rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**2) - + + case default + + rsh = ( (1-Q) * roh + Q * rov) * EXP(-h_mc*cos_inc**2) + rsv = ( (1-Q) * rov + Q * roh) * EXP(-h_mc*cos_inc**2) + end select + ! ------------------------------------------------------------- ! ! Tb at top of vegetation (excl atmos contribution) (tau-omega model) From bbe4e02a8e6a709371050615c027759bdf65c467 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 8 May 2023 18:06:22 -0400 Subject: [PATCH 3/6] use RTM_ID alone to identify mwRTM config; cleanup of Mironov subroutine --- .../LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml | 7 +- .../GEOS_LandAssimGridComp.F90 | 29 +- .../clsm_ensupd_upd_routines.F90 | 57 +--- .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 295 +++++++++++------- 4 files changed, 221 insertions(+), 167 deletions(-) diff --git a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml index f8d77eaf..259eeb52 100644 --- a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -103,9 +103,10 @@ fcsterr_inflation_fac = -9999. ! %RTM_ID = ID of radiative transfer model to use for Tb forward modeling ! (subroutine get_obs_pred()) ! 0 = none -! 1 = tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) -! 2 = same as 1 but without Pellarin atmospheric corrections -! 3 = ... +! 1 = L-band tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) (SMOS) +! 2 = same as 1 but without Pellarin atm corr (SMAP) +! 3 = same as 1 but with Mironov and SMAP L2_SM pol mixing (SMOS) +! 4 = same as 3 but without Pellarin atm corr (targeted for SMAP L4_SM Version 8) ! %bias_Npar = number of obs bias states tracked per day (integer) ! %bias_trel = e-folding time scale of obs bias memory [s] ! %bias_tcut = cutoff time for confident obs bias estimate [s] diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index 843c081d..e90e09d3 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -115,6 +115,8 @@ module GEOS_LandAssimGridCompMod logical :: mwRTM logical, allocatable :: tb_nodata(:) + + character(len=400) :: err_msg contains @@ -2248,7 +2250,7 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) ! hard-coded SMAP Tb parameters real, parameter :: freq = 1.41e9 ! microwave frequency [Hz] real, parameter :: inc_angle = 40. ! incidence angle [deg] - logical, parameter :: incl_atm_terms = .false. ! no atmospheric correction, ie, get Tb at top-of-vegetation + integer, parameter :: RTM_ID = 4 ! config of RTM - see obs_param (LDAS_DEFAULT_inputs_ensupd.nml) integer :: status character(len=ESMF_MAXSTR) :: Iam='CALC_LAND_TB' @@ -2375,23 +2377,34 @@ subroutine CALC_LAND_TB(gc, import, export, clock, rc) allocate(TB_h_tmp(N_catl), TB_v_tmp(N_catl)) - if (.not. incl_atm_terms) then + select case (RTM_ID) + + case(2,4) + allocate(dummy_real(N_catl)) ! allocate needed for GNU compiler + call mwRTM_get_Tb( & N_catl, freq, inc_angle, mwRTM_param, & - dummy_real, & ! intent(in), "elev", not used as long as "incl_atm_terms=.false." + dummy_real, & ! intent(in), "elev", not used as long as RTM_ID=4 (formerly "incl_atm_terms=.false.") LAI, & sfmc_mwRTM, & tsoil_mwRTM, & SWE, & dummy_real, & ! intent(in), "Tair", not used as long as "incl_atm_terms=.false." - incl_atm_terms, & - 'mironov', & + RTM_ID, & Tb_h_tmp, Tb_v_tmp ) ! intent(out) 'TB_LAND_1410MHZ_40DEG_HPOL', 'TB_LAND_1410MHZ_40DEG_VPOL' deallocate(dummy_real) - else - _ASSERT(.false., "top-of-atmosphere Tb calculation not yet implemented (incl_atm_terms=.true.)") - end if + + case(1,3) + + _ASSERT(.false., "top-of-atmosphere Tb calculation (requested per RTM_ID) not yet implemented") + + case default + + err_msg = 'unknown RTM_ID (during CALC_LAND_TB)' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end select if (collect_tb_counter == 0) then TB_H_enavg = 0. 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 7e3a3681..6b8726ad 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 @@ -1350,57 +1350,26 @@ subroutine get_obs_pred( & do j=1,N_TbuniqFreqAngRTMid - freq = Tb_freq_ang_RTMid(j,1) - inc_angle = Tb_freq_ang_RTMid(j,2) - RTM_id = Tb_freq_ang_RTMid(j,3) + freq = Tb_freq_ang_RTMid(j,1) + inc_angle = Tb_freq_ang_RTMid(j,2) + RTM_id = nint(Tb_freq_ang_RTMid(j,3)) ! Select a specific configuration of the RTM via the field ! "RTM_ID" in the "obs_param" type. ! ! %RTM_ID = ID of radiative transfer model to use for Tb forward modeling ! (subroutine get_obs_pred()) - ! 0 = none - ! 1 = tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) - ! 2 = same as 1 but without Pellarin atmospheric corrections - ! 3 = ... + ! 0 = none + ! 1 = L-band tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) (SMOS) + ! 2 = same as 1 but without Pellarin atm corr (SMAP) + ! 3 = same as 1 but with Mironov and SMAP L2_SM pol mixing (SMOS) + ! 4 = same as 3 but without Pellarin atm corr (targeted for SMAP L4_SM Version 8) + + call mwRTM_get_Tb( & + N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & + lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, RTM_ID, & + Tb_h_vec, Tb_v_vec ) - select case (RTM_id) - - case (1) - - ! bug fix: previously, mwRTM_get_Tb() was called without specifying the - ! sub-array of "stemp_l" that corresponds to ensemble member n_e - ! - reichle, 11 Dec 2013 - - call mwRTM_get_Tb( & - N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & - lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .true., & - 'wang',Tb_h_vec, Tb_v_vec ) - - case (2) - - call mwRTM_get_Tb( & - N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & - lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .false., & - 'wang',Tb_h_vec, Tb_v_vec ) - - case (3) - - call mwRTM_get_Tb( & - N_catl, freq, inc_angle, mwRTM_param, tile_coord_l%elev, & - lai, smoist, stemp_l(:,n_e), SWE, met_force%Tair, .false., & - 'mironov',Tb_h_vec, Tb_v_vec ) - - - case default - - write (tmpstring10,*) RTM_ID - - err_msg = 'unknown or inconsistent RTM_ID=' // tmpstring10 - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end select - Tb_h_l(:,j,n_e) = Tb_h_vec Tb_v_l(:,j,n_e) = Tb_v_vec diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index ef8da3ed..20fe664b 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -9,9 +9,10 @@ module mwRTM_routines ! %RTM_ID = ID of radiative transfer model to use for Tb forward modeling ! (subroutine get_obs_pred()) ! 0 = none - ! 1 = tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) - ! 2 = same as 1 but without Pellarin atmospheric corrections - ! 3 = ... + ! 1 = L-band tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) (SMOS) + ! 2 = same as 1 but without Pellarin atm corr (SMAP) + ! 3 = same as 1 but with Mironov and SMAP L2_SM pol mixing (SMOS) + ! 4 = same as 3 but without Pellarin atm corr (targeted for SMAP L4_SM Version 8) ! ! reichle, 16 May 2011 ! reichle, 31 Mar 2015 - added RTM_ID @@ -193,8 +194,7 @@ module mwRTM_routines ! **************************************************************** subroutine catch2mwRTM_vars( N_tile, vegcls_catch, poros_catch, poros_mwRTM, & - sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, & - tp1_in_Kelvin ) + sfmc_catch, tsurf_catch, tp1_catch, sfmc_mwRTM, tsoil_mwRTM, tp1_in_Kelvin ) ! convert soil moisture, surface temperature, and soil temperature from the Catchment ! model into soil moisture and soil temperature inputs for the microwave radiative @@ -252,8 +252,7 @@ end subroutine catch2mwRTM_vars ! **************************************************************** subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & - LAI, soilmoist, soiltemp, SWE, Tair, incl_atm_terms, & - soildiel_model, Tb_h, Tb_v ) + LAI, soilmoist, soiltemp, SWE, Tair, RTM_ID, Tb_h, Tb_v ) !--------------------------------------------------- !RTM adapted from Steven Chan @@ -310,8 +309,7 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & real, dimension(N_tile), intent(in) :: SWE ! [kg/m2] "mm" real, dimension(N_tile), intent(in) :: Tair ! [K] - logical, intent(in) :: incl_atm_terms - character(len=*), intent(in) :: soildiel_model + integer, intent(in) :: RTM_ID real, dimension(N_tile), intent(out) :: Tb_h, Tb_v ! [K] @@ -351,14 +349,27 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & !if (logit) write(logunit,*) 'entering mwRTM_get_Tb...' ! check first element of elevation against no-data-value - ! (elevation is needed only when incl_atm_terms=.true.) + ! (elevation is only needed for RTMs with Pellarin atm corr) - if (incl_atm_terms) then + select case (RTM_ID) + + case(1,3) + if ( abs(elev(1)-nodata_generic) Date: Tue, 9 May 2023 11:41:56 -0400 Subject: [PATCH 4/6] fixed typo in variable declaration (mwRTM_routines.F90) --- .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index 20fe664b..ac50968b 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -952,8 +952,8 @@ SUBROUTINE mironov(freq,mv,clayfrac,er_r) REAL :: tmptaub, tmptaub2plus1, tmpdiffepsb REAL :: tmptauu, tmptauu2plus1, tmpdiffepsu - REAL :: tmpreal, tmpepsbreal2, tmpepbsimag2 - REAL :: tmpepsureal2, tmpepbuimag2 + REAL :: tmpreal, tmpepsbreal2, tmpepsbimag2 + REAL :: tmpepsureal2, tmpepsuimag2 ! -------------------------------------------------------------------------------------- From 3f3bb9912fc2b234bbc1bf7aa6cf87910e266487 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 9 May 2023 12:03:56 -0400 Subject: [PATCH 5/6] adding module use statement for LDAS_exceptionsMod (GEOS_LandAssimGridComp.F90) --- .../GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index e90e09d3..f515ad8d 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -17,6 +17,8 @@ module GEOS_LandAssimGridCompMod use ESMF_CFIOMOD, only: ESMF_CFIOstrTemplate use GEOS_ExportCatchIncrGridCompMod, only: ExportCatchIncrSetServices=>SetServices use MAPL_ConstantsMod, only: MAPL_TICE + + use LDAS_exceptionsMod, only: ldas_abort, LDAS_GENERIC_ERROR use LDAS_TileCoordType, only: tile_coord_type use LDAS_TileCoordType, only: grid_def_type @@ -68,7 +70,6 @@ module GEOS_LandAssimGridCompMod use, intrinsic :: ieee_arithmetic - implicit none include 'mpif.h' From 1333041c0ca24770e38d5461e22ca469b142cc9c Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 11 May 2023 09:27:01 -0400 Subject: [PATCH 6/6] additional comments and clean-up for new mwRTM (GEOS_LandAssimGridComp.F90; mwRTM_routines.F90) --- .../GEOS_LandAssimGridComp.F90 | 4 +- .../GEOSlandassim_GridComp/mwRTM_routines.F90 | 62 +++++++++---------- 2 files changed, 32 insertions(+), 34 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index f515ad8d..f13c1bfc 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -2237,7 +2237,9 @@ end subroutine UPDATE_ASSIM ! ****************************************************************************** ! subroutine to calculate Tb for HISTORY output - + ! + ! IMPORTANT: hardwired mwRTM configuration for SMAP L-band Tb w/o Pellarin atm correction (RTM_ID=4) + subroutine CALC_LAND_TB(gc, import, export, clock, rc) type(ESMF_GridComp), intent(inout) :: gc ! Gridded component diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 index ac50968b..8a2fc5b6 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/mwRTM_routines.F90 @@ -408,7 +408,7 @@ subroutine mwRTM_get_Tb( N_tile, freq, inc_angle, mwp, elev, & case(3,4) - CALL mironov(freq,soilmoist(n),mwp(n)%clay,c_er) + CALL MIRONOV( freq, soilmoist(n), mwp(n)%clay, c_er ) case default @@ -924,11 +924,11 @@ END SUBROUTINE DIEL_WAT ! ********************************************************************** - SUBROUTINE mironov(freq,mv,clayfrac,er_r) + SUBROUTINE MIRONOV( freq, mv, clayfrac, er_r ) ! Soil dielectric mixing model by Mironov et al IEEE TGRS 2009, doi:10.1109/TGRS.2008.2011631 ! - ! 8 May 2023: Implementation taken from SMAP L2_SM_P retrieval model; clean-up by reichle + ! 8 May 2023: Implementation taken from SMAP L2_SM_P retrieval model by qliu; clean-up by reichle IMPLICIT NONE @@ -952,13 +952,12 @@ SUBROUTINE mironov(freq,mv,clayfrac,er_r) REAL :: tmptaub, tmptaub2plus1, tmpdiffepsb REAL :: tmptauu, tmptauu2plus1, tmpdiffepsu - REAL :: tmpreal, tmpepsbreal2, tmpepsbimag2 - REAL :: tmpepsureal2, tmpepsuimag2 + REAL :: tmpreal, tmprealb, tmprealu ! -------------------------------------------------------------------------------------- - f = freq ! Section IV - C = clayfrac*100 ! Section VI + f = freq ! Section IV + C = clayfrac*100 ! Section VI !! Mironov's regression expressions based on Curtis, Dobson, and Hallikainen datasets !! @@ -975,17 +974,17 @@ SUBROUTINE mironov(freq,mv,clayfrac,er_r) !! b: bound soil water (BSW) !! u: unbound (free) soil water (FSW) - nd = 1.634 - 0.539e-2 * C + 0.2748e-4 * C**2 ! Eqn 17 - kd = 0.03952 - 0.04038e-2 * C ! Eqn 18 - mvt = 0.02863 + 0.30673e-2 * C ! Eqn 19 - eps0b = 79.8 - 85.4e-2 * C + 32.7e-4 * C**2 ! Eqn 20 - taub = 1.062e-11 + 3.450e-12 * 1e-2 * C ! Eqn 21 - sigb = 0.3112 + 0.467e-2 * C ! Eqn 22 - sigu = 0.3631 + 1.217e-2 * C ! Eqn 23 - eps0u = 100. ! Eqn 24 - tauu = 8.5e-12 ! Eqn 25 + nd = 1.634 - 0.539e-2 * C + 0.2748e-4 * C**2 ! Eqn 17 + kd = 0.03952 - 0.04038e-2 * C ! Eqn 18 + mvt = 0.02863 + 0.30673e-2 * C ! Eqn 19 + eps0b = 79.8 - 85.4e-2 * C + 32.7e-4 * C**2 ! Eqn 20 + taub = 1.062e-11 + 3.450e-12 * 1e-2 * C ! Eqn 21 + sigb = 0.3112 + 0.467e-2 * C ! Eqn 22 + sigu = 0.3631 + 1.217e-2 * C ! Eqn 23 + eps0u = 100. ! Eqn 24 + tauu = 8.5e-12 ! Eqn 25 - !! Debye relaxation equations for water as a function of frequency ! Eqn 16 + !! Debye relaxation equations for water as a function of frequency ! Eqn 16 tmp2PIf = 2.*MAPL_PI*f @@ -1007,35 +1006,32 @@ SUBROUTINE mironov(freq,mv,clayfrac,er_r) !! Refractive indices and normalized attenuation coefficients - tmpreal = 1/sqrt(2.0) + tmpreal = 1/sqrt(2.0) - tmpepsbreal2 = epsb_real**2 - tmpepsbimag2 = epsb_imag**2 + tmprealb = sqrt( epsb_real**2 + epsb_imag**2 ) + tmprealu = sqrt( epsu_real**2 + epsu_imag**2 ) - tmpepsureal2 = epsu_real**2 - tmpepsuimag2 = epsu_imag**2 - - nb = tmpreal * sqrt( sqrt(tmpepsbreal2 + tmpepsbimag2) + epsb_real ) ! Eqn 14 - kb = tmpreal * sqrt( sqrt(tmpepsbreal2 + tmpepsbimag2) - epsb_real ) ! Eqn 15 - nu = tmpreal * sqrt( sqrt(tmpepsureal2 + tmpepsuimag2) + epsu_real ) ! Eqn 14 - ku = tmpreal * sqrt( sqrt(tmpepsureal2 + tmpepsuimag2) - epsu_real ) ! Eqn 15 + nb = tmpreal * sqrt( tmprealb + epsb_real ) ! Eqn 14 + kb = tmpreal * sqrt( tmprealb - epsb_real ) ! Eqn 15 + nu = tmpreal * sqrt( tmprealu + epsu_real ) ! Eqn 14 + ku = tmpreal * sqrt( tmprealu - epsu_real ) ! Eqn 15 IF (mv <= mvt) THEN - nm = nd + (nb - 1.) * mv ! Eqn 12 - km = kd + kb * mv ! Eqn 13 + nm = nd + (nb - 1.) * mv ! Eqn 12 + km = kd + kb * mv ! Eqn 13 ELSE - nm = nd + (nb - 1.) * mvt + (nu - 1.) * (mv - mvt) ! Eqn 12 - km = kd + kb * mvt + ku * (mv - mvt) ! Eqn 13 + nm = nd + (nb - 1.) * mvt + (nu - 1.) * (mv - mvt) ! Eqn 12 + km = kd + kb * mvt + ku * (mv - mvt) ! Eqn 13 ENDIF ! complex dielectric constant of moist soil - er_r_real = nm**2 - km**2 ! Eqn 11 - er_r_imag = 2. * nm * km ! Eqn 11 + er_r_real = nm**2 - km**2 ! Eqn 11 + er_r_imag = 2. * nm * km ! Eqn 11 er_r = CMPLX(er_r_real,er_r_imag)