diff --git a/Makefile b/Makefile index 1184026d..ddf302fe 100755 --- a/Makefile +++ b/Makefile @@ -79,6 +79,7 @@ OBJS_MKSRFDATA = \ Aggregation_Topography.o \ Aggregation_TopographyFactors.o \ Aggregation_Urban.o \ + Aggregation_SoilTexture.o \ MOD_MeshFilter.o \ MOD_RegionClip.o \ MKSRFDATA.o @@ -143,6 +144,7 @@ OBJS_BASIC = \ MOD_DBedrockReadin.o \ MOD_SoilColorRefl.o \ MOD_SoilParametersReadin.o \ + MOD_SoilTextureReadin.o \ MOD_HtopReadin.o \ MOD_UrbanReadin.o \ MOD_BGC_CNSummary.o \ diff --git a/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 b/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 index 94bbd490..43b73491 100644 --- a/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 +++ b/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 @@ -20,13 +20,20 @@ MODULE MOD_Catch_SubsurfaceFlow IMPLICIT NONE real(r8), parameter :: e_ice = 6.0 ! soil ice impedance factor - real(r8), parameter :: raniso = 1. ! anisotropy ratio, unitless + + ! anisotropy ratio of lateral/vertical hydraulic conductivity (unitless) + ! for USDA soil texture class: + ! 0: undefined + ! 1: clay; 2: silty clay; 3: sandy clay; 4: clay loam; 5: silty clay loam; 6: sandy clay loam; & + ! 7: loam; 8: silty loam; 9: sandy loam; 10: silt; 11: loamy sand; 12: sand + real(r8), parameter :: raniso(0:12) = (/ 1., & + 48., 40., 28., 24., 20., 14., 12., 10., 4., 2., 3., 2. /) ! -- neighbour variables -- type(pointer_real8_1d), allocatable :: agwt_nb (:) ! ground water area (for patchtype <= 2) of neighbours [m^2] type(pointer_real8_1d), allocatable :: theta_a_nb (:) ! saturated volume content [-] type(pointer_real8_1d), allocatable :: zwt_nb (:) ! water table depth [m] - type(pointer_real8_1d), allocatable :: Ks_nb (:) ! saturated hydraulic conductivity [m/s] + type(pointer_real8_1d), allocatable :: Kl_nb (:) ! lateral hydraulic conductivity [m/s] type(pointer_real8_1d), allocatable :: wdsrf_nb (:) ! depth of surface water [m] type(pointer_logic_1d), allocatable :: islake_nb (:) ! whether a neighbour is water body @@ -59,7 +66,7 @@ SUBROUTINE basin_neighbour_init () CALL allocate_neighbour_data (agwt_nb ) CALL allocate_neighbour_data (theta_a_nb) CALL allocate_neighbour_data (zwt_nb ) - CALL allocate_neighbour_data (Ks_nb ) + CALL allocate_neighbour_data (Kl_nb ) CALL allocate_neighbour_data (wdsrf_nb ) CALL allocate_neighbour_data (islake_nb ) @@ -127,7 +134,7 @@ SUBROUTINE subsurface_flow (deltime) real(r8), allocatable :: theta_a_h (:) real(r8), allocatable :: zwt_h (:) - real(r8), allocatable :: Ks_h (:) ! [m/s] + real(r8), allocatable :: Kl_h (:) ! [m/s] real(r8), allocatable :: xsubs_h (:) ! [m/s] real(r8), allocatable :: xsubs_fc (:) ! [m/s] @@ -135,17 +142,17 @@ SUBROUTINE subsurface_flow (deltime) real(r8) :: theta_s_h, air_h, icefrac, imped, delp real(r8) :: sumwt, sumarea, zwt_mean real(r8) :: zsubs_h_up, zsubs_h_dn - real(r8) :: slope, bdamp, Ks_fc, Ks_in + real(r8) :: slope, bdamp, Kl_fc, Kl_in real(r8) :: ca, cb real(r8) :: alp real(r8), allocatable :: theta_a_bsn (:) real(r8), allocatable :: zwt_bsn (:) - real(r8), allocatable :: Ks_bsn (:) ! [m/s] + real(r8), allocatable :: Kl_bsn (:) ! [m/s] integer :: jnb - real(r8) :: zsubs_up, zwt_up, Ks_up, theta_a_up, area_up - real(r8) :: zsubs_dn, zwt_dn, Ks_dn, theta_a_dn, area_dn + real(r8) :: zsubs_up, zwt_up, Kl_up, theta_a_up, area_up + real(r8) :: zsubs_dn, zwt_dn, Kl_dn, theta_a_dn, area_dn real(r8) :: lenbdr, xsubs_nb logical :: iam_lake, nb_is_lake, has_river @@ -185,7 +192,7 @@ SUBROUTINE subsurface_flow (deltime) IF (numbasin > 0) THEN allocate (theta_a_bsn (numbasin)); theta_a_bsn = 0. allocate (zwt_bsn (numbasin)); zwt_bsn = 0. - allocate (Ks_bsn (numbasin)); Ks_bsn = 0. + allocate (Kl_bsn (numbasin)); Kl_bsn = 0. ENDIF DO ibasin = 1, numbasin @@ -199,7 +206,7 @@ SUBROUTINE subsurface_flow (deltime) allocate (theta_a_h (nhru)); theta_a_h = 0. allocate (zwt_h (nhru)); zwt_h = 0. - allocate (Ks_h (nhru)); Ks_h = 0. + allocate (Kl_h (nhru)); Kl_h = 0. DO i = 1, nhru @@ -253,23 +260,23 @@ SUBROUTINE subsurface_flow (deltime) ENDIF ENDIF - Ks_h(i) = 0. + Kl_h(i) = 0. sumwt = 0. DO ipatch = ps, pe IF (patchtype(ipatch) <= 2) THEN DO ilev = 1, nl_soil icefrac = min(1., wice_soisno(ilev,ipatch)/denice/dz_soi(ilev)/porsl(ilev,ipatch)) imped = 10.**(-e_ice*icefrac) - Ks_h(i) = Ks_h(i) + hru_patch%subfrc(ipatch) & + Kl_h(i) = Kl_h(i) + hru_patch%subfrc(ipatch) * raniso(soiltext(ipatch)) & * hksati(ilev,ipatch)/1.0e3 * imped * dz_soi(ilev)/zi_soi(nl_soil) ENDDO sumwt = sumwt + hru_patch%subfrc(ipatch) ENDIF ENDDO - IF (sumwt > 0) Ks_h(i) = Ks_h(i) / sumwt + IF (sumwt > 0) Kl_h(i) = Kl_h(i) / sumwt ELSE ! Frozen soil. - Ks_h(i) = 0. + Kl_h(i) = 0. ENDIF ENDDO @@ -285,11 +292,11 @@ SUBROUTINE subsurface_flow (deltime) j = hrus%inext(i) IF (j <= 0) CYCLE ! downstream is out of catchment - IF (Ks_h(i) == 0.) CYCLE ! this HRU is frozen + IF (Kl_h(i) == 0.) CYCLE ! this HRU is frozen j_is_river = (hrus%indx(j) == 0) - IF ((.not. j_is_river) .and. (Ks_h(j) == 0.)) CYCLE ! non-river downstream HRU is frozen + IF ((.not. j_is_river) .and. (Kl_h(j) == 0.)) CYCLE ! non-river downstream HRU is frozen zsubs_h_up = hrus%elva(i) - zwt_h(i) @@ -317,27 +324,27 @@ SUBROUTINE subsurface_flow (deltime) IF ((zsubs_h_up > zsubs_h_dn) .or. j_is_river) THEN IF (zwt_h(i) > 1.5) THEN ! from Fan et al., JGR 112(D10125) - Ks_fc = raniso * Ks_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp) + Kl_fc = Kl_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp) ELSE - Ks_fc = raniso * Ks_h(i) * ((1.5-zwt_h(i)) + bdamp) + Kl_fc = Kl_h(i) * ((1.5-zwt_h(i)) + bdamp) ENDIF ELSE IF (zwt_h(j) > 1.5) THEN - Ks_fc = raniso * Ks_h(j) * bdamp * exp(-(zwt_h(j)-1.5)/bdamp) + Kl_fc = Kl_h(j) * bdamp * exp(-(zwt_h(j)-1.5)/bdamp) ELSE - Ks_fc = raniso * Ks_h(j) * ((1.5-zwt_h(j)) + bdamp) + Kl_fc = Kl_h(j) * ((1.5-zwt_h(j)) + bdamp) ENDIF ENDIF - ca = hrus%flen(i) * Ks_fc / theta_a_h(i) / delp / hrus%agwt(i) * deltime + ca = hrus%flen(i) * Kl_fc / theta_a_h(i) / delp / hrus%agwt(i) * deltime IF (.not. j_is_river) THEN - cb = hrus%flen(i) * Ks_fc / theta_a_h(j) / delp / hrus%agwt(j) * deltime + cb = hrus%flen(i) * Kl_fc / theta_a_h(j) / delp / hrus%agwt(j) * deltime ELSE - cb = hrus%flen(i) * Ks_fc / delp / hrus%area(j) * deltime + cb = hrus%flen(i) * Kl_fc / delp / hrus%area(j) * deltime ENDIF - xsubs_fc(i) = (zsubs_h_up - zsubs_h_dn) * hrus%flen(i) * Ks_fc / (1+ca+cb) / delp + xsubs_fc(i) = (zsubs_h_up - zsubs_h_dn) * hrus%flen(i) * Kl_fc / (1+ca+cb) / delp xsubs_h(i) = xsubs_h(i) + xsubs_fc(i) / hrus%agwt(i) @@ -390,9 +397,9 @@ SUBROUTINE subsurface_flow (deltime) IF (zwt_h(i) > 1.5) THEN ! from Fan et al., JGR 112(D10125) - Ks_in = raniso * Ks_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp) + Kl_in = Kl_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp) ELSE - Ks_in = raniso * Ks_h(i) * ((1.5-zwt_h(i)) + bdamp) + Kl_in = Kl_h(i) * ((1.5-zwt_h(i)) + bdamp) ENDIF ps = hru_patch%substt(hrus%ihru(i)) @@ -403,7 +410,7 @@ SUBROUTINE subsurface_flow (deltime) DO ipatch = ps, pe IF (patchtype(ipatch) <= 2) THEN - xsubs_pch(ipatch) = - Ks_in * (zwt(ipatch) - zwt_mean) *6.0*pi/hrus%agwt(i) + xsubs_pch(ipatch) = - Kl_in * (zwt(ipatch) - zwt_mean) *6.0*pi/hrus%agwt(i) ! Update total subsurface lateral flow (2): Between patches xwsub(ipatch) = xwsub(ipatch) + xsubs_pch(ipatch) * 1.e3 ! m/s to mm/s ENDIF @@ -416,12 +423,12 @@ SUBROUTINE subsurface_flow (deltime) IF (sumarea > 0) THEN theta_a_bsn (ibasin) = sum(theta_a_h * hrus%agwt) / sumarea zwt_bsn (ibasin) = sum(zwt_h * hrus%agwt) / sumarea - Ks_bsn (ibasin) = sum(Ks_h * hrus%agwt) / sumarea + Kl_bsn (ibasin) = sum(Kl_h * hrus%agwt) / sumarea ENDIF deallocate (theta_a_h) deallocate (zwt_h ) - deallocate (Ks_h ) + deallocate (Kl_h ) deallocate (xsubs_h ) deallocate (xsubs_fc ) @@ -429,7 +436,7 @@ SUBROUTINE subsurface_flow (deltime) CALL retrieve_neighbour_data (theta_a_bsn, theta_a_nb) CALL retrieve_neighbour_data (zwt_bsn , zwt_nb ) - CALL retrieve_neighbour_data (Ks_bsn , Ks_nb ) + CALL retrieve_neighbour_data (Kl_bsn , Kl_nb ) CALL retrieve_neighbour_data (wdsrf_bsn , wdsrf_nb ) DO ibasin = 1, numbasin @@ -449,7 +456,7 @@ SUBROUTINE subsurface_flow (deltime) ENDIF IF (.not. iam_lake) THEN - Ks_up = Ks_bsn (ibasin) + Kl_up = Kl_bsn (ibasin) zwt_up = zwt_bsn (ibasin) theta_a_up = theta_a_bsn(ibasin) zsubs_up = elementneighbour(ibasin)%myelva - zwt_up @@ -461,7 +468,7 @@ SUBROUTINE subsurface_flow (deltime) ENDIF IF (.not. nb_is_lake) THEN - Ks_dn = Ks_nb(ibasin)%val(jnb) + Kl_dn = Kl_nb(ibasin)%val(jnb) zwt_dn = zwt_nb(ibasin)%val(jnb) theta_a_dn = theta_a_nb(ibasin)%val(jnb) zsubs_dn = elementneighbour(ibasin)%elva(jnb) - zwt_dn @@ -474,8 +481,8 @@ SUBROUTINE subsurface_flow (deltime) IF ((.not. iam_lake) .and. (area_up <= 0)) CYCLE IF ((.not. nb_is_lake) .and. (area_dn <= 0)) CYCLE - IF ((.not. iam_lake) .and. (Ks_up == 0.) ) CYCLE - IF ((.not. nb_is_lake) .and. (Ks_dn == 0.) ) CYCLE + IF ((.not. iam_lake) .and. (Kl_up == 0. )) CYCLE + IF ((.not. nb_is_lake) .and. (Kl_dn == 0. )) CYCLE ! water body is dry. IF (iam_lake .and. (zsubs_up > zsubs_dn) .and. (wdsrf_bsn(ibasin) == 0.)) THEN @@ -507,22 +514,22 @@ SUBROUTINE subsurface_flow (deltime) IF (nb_is_lake .or. ((.not. iam_lake) .and. (zsubs_up > zsubs_dn))) THEN IF (zwt_up > 1.5) THEN ! from Fan et al., JGR 112(D10125) - Ks_fc = raniso * Ks_up * bdamp * exp(-(zwt_up-1.5)/bdamp) + Kl_fc = Kl_up * bdamp * exp(-(zwt_up-1.5)/bdamp) ELSE - Ks_fc = raniso * Ks_up * ((1.5-zwt_up) + bdamp) + Kl_fc = Kl_up * ((1.5-zwt_up) + bdamp) ENDIF ELSE IF (zwt_dn > 1.5) THEN - Ks_fc = raniso * Ks_dn * bdamp * exp(-(zwt_dn-1.5)/bdamp) + Kl_fc = Kl_dn * bdamp * exp(-(zwt_dn-1.5)/bdamp) ELSE - Ks_fc = raniso * Ks_dn * ((1.5-zwt_dn) + bdamp) + Kl_fc = Kl_dn * ((1.5-zwt_dn) + bdamp) ENDIF ENDIF - ca = lenbdr * Ks_fc / theta_a_up / delp / area_up * deltime - cb = lenbdr * Ks_fc / theta_a_dn / delp / area_dn * deltime + ca = lenbdr * Kl_fc / theta_a_up / delp / area_up * deltime + cb = lenbdr * Kl_fc / theta_a_dn / delp / area_dn * deltime - xsubs_nb = (zsubs_up - zsubs_dn) * lenbdr * Ks_fc / (1+ca+cb) / delp + xsubs_nb = (zsubs_up - zsubs_dn) * lenbdr * Kl_fc / (1+ca+cb) / delp IF (.not. iam_lake) THEN xsubs_nb = xsubs_nb / sum(hrus%agwt) @@ -547,7 +554,7 @@ SUBROUTINE subsurface_flow (deltime) IF (allocated(theta_a_bsn)) deallocate(theta_a_bsn) IF (allocated(zwt_bsn )) deallocate(zwt_bsn ) - IF (allocated(Ks_bsn )) deallocate(Ks_bsn ) + IF (allocated(Kl_bsn )) deallocate(Kl_bsn ) ENDIF @@ -712,7 +719,7 @@ SUBROUTINE basin_neighbour_final () IF (allocated(theta_a_nb)) deallocate(theta_a_nb) IF (allocated(zwt_nb )) deallocate(zwt_nb ) - IF (allocated(Ks_nb )) deallocate(Ks_nb ) + IF (allocated(Kl_nb )) deallocate(Kl_nb ) IF (allocated(wdsrf_nb )) deallocate(wdsrf_nb ) IF (allocated(agwt_nb )) deallocate(agwt_nb ) IF (allocated(islake_nb )) deallocate(islake_nb ) diff --git a/main/MOD_Vars_TimeInvariants.F90 b/main/MOD_Vars_TimeInvariants.F90 index 8494298a..556e2ba5 100644 --- a/main/MOD_Vars_TimeInvariants.F90 +++ b/main/MOD_Vars_TimeInvariants.F90 @@ -216,6 +216,8 @@ MODULE MOD_Vars_TimeInvariants real(r8), allocatable :: fc_vgm (:,:) !a scaling factor by using air entry value in the Mualem model [-] #endif + integer, allocatable :: soiltext(:) ! USDA soil texture class + real(r8), allocatable :: fsatmax (:) ! maximum saturated area fraction [-] real(r8), allocatable :: fsatdcf (:) ! decay factor in calucation of saturated area fraction [1/m] @@ -336,6 +338,7 @@ SUBROUTINE allocate_TimeInvariants () allocate (sc_vgm (nl_soil,numpatch)) allocate (fc_vgm (nl_soil,numpatch)) #endif + allocate (soiltext(numpatch)) allocate (fsatmax (numpatch)) allocate (fsatdcf (numpatch)) @@ -455,6 +458,8 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart) CALL ncio_read_vector (file_restart, 'fc_vgm ' , nl_soil, landpatch, fc_vgm ) ! a scaling factor by using air entry value in the Mualem model [-] #endif + CALL ncio_read_vector (file_restart, 'soiltext', landpatch, soiltext, defval = 0 ) + CALL ncio_read_vector (file_restart, 'fsatmax', landpatch, fsatmax, defval = 0.38 ) CALL ncio_read_vector (file_restart, 'fsatdcf', landpatch, fsatdcf, defval = 0.125) @@ -635,6 +640,8 @@ SUBROUTINE WRITE_TimeInvariants (lc_year, casename, dir_restart) CALL ncio_write_vector (file_restart, 'fc_vgm ' , 'soil', nl_soil, 'patch', landpatch, fc_vgm , compress) ! a scaling factor by using air entry value in the Mualem model [-] #endif + CALL ncio_write_vector (file_restart, 'soiltext', 'patch', landpatch, soiltext) + CALL ncio_write_vector (file_restart, 'fsatmax', 'patch', landpatch, fsatmax) CALL ncio_write_vector (file_restart, 'fsatdcf', 'patch', landpatch, fsatdcf) @@ -781,6 +788,7 @@ SUBROUTINE deallocate_TimeInvariants () deallocate (sc_vgm ) deallocate (fc_vgm ) #endif + deallocate (soiltext ) deallocate (fsatmax ) deallocate (fsatdcf ) @@ -891,13 +899,15 @@ SUBROUTINE check_TimeInvariants () CALL check_vector_data ('BA_alpha [-] ', BA_alpha ) ! alpha in Balland and Arp(2005) thermal conductivity scheme CALL check_vector_data ('BA_beta [-] ', BA_beta ) ! beta in Balland and Arp(2005) thermal conductivity scheme + CALL check_vector_data ('soiltexture [-] ', soiltext , -1) ! + CALL check_vector_data ('htop [m] ', htop ) CALL check_vector_data ('hbot [m] ', hbot ) IF(DEF_USE_BEDROCK)THEN CALL check_vector_data ('dbedrock [m] ', dbedrock ) ! ENDIF - + CALL check_vector_data ('topoelv [m] ', topoelv ) ! CALL check_vector_data ('topostd [m] ', topostd ) ! CALL check_vector_data ('BVIC [-] ', BVIC ) ! diff --git a/mkinidata/MOD_Initialize.F90 b/mkinidata/MOD_Initialize.F90 index edc45a59..1492ec45 100644 --- a/mkinidata/MOD_Initialize.F90 +++ b/mkinidata/MOD_Initialize.F90 @@ -80,6 +80,7 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & USE MOD_LakeDepthReadin USE MOD_PercentagesPFTReadin USE MOD_SoilParametersReadin + USE MOD_SoilTextureReadin IMPLICIT NONE @@ -356,9 +357,6 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & ENDIF ENDIF - - - #ifdef vanGenuchten_Mualem_SOIL_MODEL IF (p_is_worker) THEN IF (numpatch > 0) THEN @@ -376,6 +374,8 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & ENDIF #endif + CALL soiltext_readin (dir_landdata, lc_year) + ! ............................................................... ! 1.4 Plant time-invariant variables ! ............................................................... diff --git a/mkinidata/MOD_SoilTextureReadin.F90 b/mkinidata/MOD_SoilTextureReadin.F90 new file mode 100644 index 00000000..db2462e6 --- /dev/null +++ b/mkinidata/MOD_SoilTextureReadin.F90 @@ -0,0 +1,55 @@ +#include + +MODULE MOD_SoilTextureReadin + + USE MOD_Precision + IMPLICIT NONE + SAVE + + ! PUBLIC MEMBER FUNCTIONS: + PUBLIC :: soiltext_readin + +CONTAINS + + SUBROUTINE soiltext_readin (dir_landdata, lc_year) + + USE MOD_Precision + USE MOD_SPMD_Task + USE MOD_Namelist + USE MOD_UserDefFun + USE MOD_LandPatch + USE MOD_NetCDFVector + USE MOD_Vars_TimeInvariants, only : soiltext +#ifdef SinglePoint + USE MOD_SingleSrfdata +#endif + + IMPLICIT NONE + + character(len=256), intent(in) :: dir_landdata + integer, intent(in) :: lc_year ! which year of land cover data used + + ! Local Variables + character(len=256) :: lndname + character(len=4) :: cyear + integer :: ipatch + +#ifdef SinglePoint + soiltext(:) = SITE_soil_texture +#else + write(cyear,'(i4.4)') lc_year + lndname = trim(dir_landdata)// '/soil/' // trim(cyear) // '/soiltexture_patches.nc' + CALL ncio_read_vector (lndname, 'soiltext_patches', landpatch, soiltext) +#endif + + IF (p_is_worker) THEN + IF (numpatch > 0) THEN + WHERE (soiltext < 0 ) soiltext = 0 + WHERE (soiltext > 12) soiltext = 0 + ENDIF + ENDIF + + END SUBROUTINE soiltext_readin + +END MODULE MOD_SoilTextureReadin + diff --git a/mksrfdata/Aggregation_SoilParameters.F90 b/mksrfdata/Aggregation_SoilParameters.F90 index cf0ac67f..a576ab79 100644 --- a/mksrfdata/Aggregation_SoilParameters.F90 +++ b/mksrfdata/Aggregation_SoilParameters.F90 @@ -280,6 +280,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = vf_quartz_mineral_s_grid, data_r8_2d_out1 = vf_quartz_mineral_s_one) + CALL fillnan (vf_quartz_mineral_s_one) vf_quartz_mineral_s_patches (ipatch) = sum (vf_quartz_mineral_s_one * (area_one/sum(area_one))) ELSE vf_quartz_mineral_s_patches (ipatch) = -1.0e36_r8 @@ -357,6 +358,10 @@ SUBROUTINE Aggregation_SoilParameters ( & data_r8_2d_in2 = vf_sand_s_grid, data_r8_2d_out2 = vf_sand_s_one, & data_r8_2d_in3 = vf_om_s_grid, data_r8_2d_out3 = vf_om_s_one) + CALL fillnan (vf_gravels_s_one) + CALL fillnan (vf_sand_s_one ) + CALL fillnan (vf_om_s_one ) + vf_gravels_s_patches (ipatch) = sum (vf_gravels_s_one * (area_one/sum(area_one))) vf_sand_s_patches (ipatch) = sum (vf_sand_s_one * (area_one/sum(area_one))) vf_om_s_patches (ipatch) = sum (vf_om_s_one * (area_one/sum(area_one))) @@ -529,6 +534,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = wf_gravels_s_grid, data_r8_2d_out1 = wf_gravels_s_one) + CALL fillnan (wf_gravels_s_one) wf_gravels_s_patches (ipatch) = sum (wf_gravels_s_one * (area_one/sum(area_one))) ELSE wf_gravels_s_patches (ipatch) = -1.0e36_r8 @@ -591,6 +597,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = wf_sand_s_grid, data_r8_2d_out1 = wf_sand_s_one) + CALL fillnan (wf_sand_s_one) wf_sand_s_patches (ipatch) = sum (wf_sand_s_one * (area_one/sum(area_one))) ELSE wf_sand_s_patches (ipatch) = -1.0e36_r8 @@ -690,6 +697,12 @@ SUBROUTINE Aggregation_SoilParameters ( & data_r8_2d_in4 = theta_s_grid, data_r8_2d_out4 = theta_s_one, & data_r8_2d_in5 = k_s_grid, data_r8_2d_out5 = k_s_one, & data_r8_2d_in6 = L_vgm_grid, data_r8_2d_out6 = L_vgm_one ) + CALL fillnan (theta_r_one ) + CALL fillnan (alpha_vgm_one) + CALL fillnan (n_vgm_one ) + CALL fillnan (theta_s_one ) + CALL fillnan (k_s_one ) + CALL fillnan (L_vgm_one ) theta_r_patches (ipatch) = sum (theta_r_one * (area_one/sum(area_one))) alpha_vgm_patches (ipatch) = median (alpha_vgm_one, size(alpha_vgm_one), spval) n_vgm_patches (ipatch) = median (n_vgm_one, size(n_vgm_one), spval) @@ -951,6 +964,10 @@ SUBROUTINE Aggregation_SoilParameters ( & data_r8_2d_in2 = k_s_grid, data_r8_2d_out2 = k_s_one, & data_r8_2d_in3 = psi_s_grid, data_r8_2d_out3 = psi_s_one, & data_r8_2d_in4 = lambda_grid, data_r8_2d_out4 = lambda_one) + CALL fillnan (theta_s_one) + CALL fillnan (k_s_one ) + CALL fillnan (psi_s_one ) + CALL fillnan (lambda_one ) theta_s_patches (ipatch) = sum (theta_s_one * (area_one/sum(area_one))) k_s_patches (ipatch) = product(k_s_one**(area_one/sum(area_one))) psi_s_patches (ipatch) = median (psi_s_one, size(psi_s_one), spval) @@ -1133,6 +1150,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = csol_grid, data_r8_2d_out1 = csol_one) + CALL fillnan (csol_one) csol_patches (ipatch) = sum(csol_one*(area_one/sum(area_one))) ELSE csol_patches (ipatch) = -1.0e36_r8 @@ -1193,6 +1211,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = tksatu_grid, data_r8_2d_out1 = tksatu_one) + CALL fillnan (tksatu_one) tksatu_patches (ipatch) = product(tksatu_one**(area_one/sum(area_one))) ELSE tksatu_patches (ipatch) = -1.0e36_r8 @@ -1253,6 +1272,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = tksatf_grid, data_r8_2d_out1 = tksatf_one) + CALL fillnan (tksatf_one) tksatf_patches (ipatch) = product(tksatf_one**(area_one/sum(area_one))) ELSE tksatf_patches (ipatch) = -1.0e36_r8 @@ -1313,6 +1333,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = tkdry_grid, data_r8_2d_out1 = tkdry_one) + CALL fillnan (tkdry_one) tkdry_patches (ipatch) = product(tkdry_one**(area_one/sum(area_one))) ELSE tkdry_patches (ipatch) = -1.0e36_r8 @@ -1373,6 +1394,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = k_solids_grid, data_r8_2d_out1 = k_solids_one) + CALL fillnan (k_solids_one) k_solids_patches (ipatch) = product(k_solids_one**(area_one/sum(area_one))) ELSE k_solids_patches (ipatch) = -1.0e36_r8 @@ -1434,6 +1456,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = OM_density_s_grid, data_r8_2d_out1 = OM_density_s_one) + CALL fillnan (OM_density_s_one) OM_density_s_patches (ipatch) = sum (OM_density_s_one * (area_one/sum(area_one))) ELSE OM_density_s_patches (ipatch) = -1.0e36_r8 @@ -1495,6 +1518,7 @@ SUBROUTINE Aggregation_SoilParameters ( & IF (L /= 0) THEN CALL aggregation_request_data (landpatch, ipatch, gland, zip = USE_zip_for_aggregation, area = area_one, & data_r8_2d_in1 = BD_all_s_grid, data_r8_2d_out1 = BD_all_s_one) + CALL fillnan (BD_all_s_one) BD_all_s_patches (ipatch) = sum (BD_all_s_one * (area_one/sum(area_one))) ELSE BD_all_s_patches (ipatch) = -1.0e36_r8 diff --git a/mksrfdata/Aggregation_SoilTexture.F90 b/mksrfdata/Aggregation_SoilTexture.F90 new file mode 100644 index 00000000..6a9eb730 --- /dev/null +++ b/mksrfdata/Aggregation_SoilTexture.F90 @@ -0,0 +1,134 @@ +#include + +SUBROUTINE Aggregation_SoilTexture ( & + gland, dir_rawdata, dir_model_landdata, lc_year) + +! ---------------------------------------------------------------------- +! DESCRIPTION: +! Aggregate soil texture class within a patch. +! +! Use the USDA soil texture triangle (using the amount of sand, clay, and +! silt contents) to identify the soil texture in fine grid resolution and +! then finding the major soil type in a patch by counting number of fine +! grids with each type of soil and adopting the major one. +! +! Created by Shupeng Zhang, 01/2025 +! ---------------------------------------------------------------------- + + USE MOD_Precision + USE MOD_Namelist + USE MOD_SPMD_Task + USE MOD_Grid + USE MOD_LandPatch + USE MOD_NetCDFVector + USE MOD_NetCDFBlock +#ifdef RangeCheck + USE MOD_RangeCheck +#endif + + USE MOD_AggregationRequestData + + USE MOD_Utils, only: num_max_frequency +#ifdef SinglePoint + USE MOD_SingleSrfdata +#endif + +#ifdef SrfdataDiag + USE MOD_SrfdataDiag +#endif + + IMPLICIT NONE + ! arguments: + + integer, intent(in) :: lc_year + type(grid_type), intent(in) :: gland + character(len=*), intent(in) :: dir_rawdata + character(len=*), intent(in) :: dir_model_landdata + + ! local variables: + ! --------------------------------------------------------------- + character(len=256) :: landdir, lndname, cyear + integer :: ipatch + + type(block_data_int32_2d) :: soiltext + integer, allocatable :: soiltext_patches(:), soiltext_one(:) +#ifdef SrfdataDiag + integer :: typpatch(N_land_classification+1), ityp +#endif + + write(cyear,'(i4.4)') lc_year + landdir = trim(dir_model_landdata) // '/soil/' // trim(cyear) + +#ifdef USEMPI + CALL mpi_barrier (p_comm_glb, p_err) +#endif + IF (p_is_master) THEN + write(*,'(/, A)') 'Aggregate soil texture ...' + CALL system('mkdir -p ' // trim(adjustl(landdir))) + ENDIF +#ifdef USEMPI + CALL mpi_barrier (p_comm_glb, p_err) +#endif + +#ifdef SinglePoint + IF (USE_SITE_soilparameters) THEN + RETURN + ENDIF +#endif + + lndname = trim(dir_rawdata)//'/soil/soil_type.nc' + + IF (p_is_io) THEN + CALL allocate_block_data (gland, soiltext) + CALL ncio_read_block (lndname, 'soil_type_l6', gland, soiltext) + +#ifdef USEMPI + CALL aggregation_data_daemon (gland, data_i4_2d_in1 = soiltext) +#endif + ENDIF + + IF (p_is_worker) THEN + + IF (numpatch > 0) allocate (soiltext_patches (numpatch)) + + DO ipatch = 1, numpatch + CALL aggregation_request_data (landpatch, ipatch, gland, & + zip = USE_zip_for_aggregation, & + data_i4_2d_in1 = soiltext, data_i4_2d_out1 = soiltext_one) + soiltext_patches(ipatch) = num_max_frequency (soiltext_one) + ENDDO + +#ifdef USEMPI + CALL aggregation_worker_done () +#endif + ENDIF + +#ifdef USEMPI + CALL mpi_barrier (p_comm_glb, p_err) +#endif + +#ifdef RangeCheck + CALL check_vector_data ('soiltext_patches ', soiltext_patches) +#endif + +#ifndef SinglePoint + lndname = trim(landdir)//'/soiltexture_patches.nc' + CALL ncio_create_file_vector (lndname, landpatch) + CALL ncio_define_dimension_vector (lndname, landpatch, 'patch') + CALL ncio_write_vector (lndname, 'soiltext_patches', 'patch', landpatch, soiltext_patches, DEF_Srfdata_CompressLevel) + +#ifdef SrfdataDiag + typpatch = (/(ityp, ityp = 0, N_land_classification)/) + lndname = trim(dir_model_landdata)//'/diag/soiltexture_'//trim(cyear)//'.nc' + CALL srfdata_map_and_write (real(soiltext_patches,r8), landpatch%settyp, typpatch, & + m_patch2diag, -1., lndname, 'soiltexture', compress = 1, write_mode = 'one') +#endif +#else + SITE_soil_texture = soiltext_patches(1) +#endif + + IF (p_is_worker) THEN + IF (numpatch > 0) deallocate ( soiltext_patches ) + ENDIF + +END SUBROUTINE Aggregation_SoilTexture diff --git a/mksrfdata/MKSRFDATA.F90 b/mksrfdata/MKSRFDATA.F90 index 94814495..536d6587 100644 --- a/mksrfdata/MKSRFDATA.F90 +++ b/mksrfdata/MKSRFDATA.F90 @@ -285,8 +285,8 @@ PROGRAM MKSRFDATA CALL mesh_build () CALL landelm_build -#ifdef GRIDBASED - IF (.not. read_mesh_from_file) THEN +#ifndef CATCHMENT + IF (DEF_LANDONLY) THEN !TODO: distinguish USGS and IGBP land cover #ifndef LULC_USGS CALL mesh_filter (gpatch, trim(DEF_dir_rawdata)//'/landtype_update.nc', 'landtype') @@ -390,6 +390,8 @@ PROGRAM MKSRFDATA grid_urban_5km, grid_urban_500m) #endif + CALL Aggregation_SoilTexture (gsoil, dir_rawdata, dir_landdata, lc_year) + ! ................................................................ ! 4. Free memories. ! ................................................................ diff --git a/mksrfdata/MOD_AggregationRequestData.F90 b/mksrfdata/MOD_AggregationRequestData.F90 index b8e01ca7..030de100 100644 --- a/mksrfdata/MOD_AggregationRequestData.F90 +++ b/mksrfdata/MOD_AggregationRequestData.F90 @@ -661,4 +661,35 @@ END SUBROUTINE aggregation_worker_done #endif + + SUBROUTINE fillnan (vec) + + USE MOD_Precision + USE MOD_UserDefFun, only : isnan_ud + IMPLICIT NONE + + real(r8), intent(inout) :: vec(:) + + ! local variables + integer :: i, n + real(r8) :: s + + n = 0 + s = 0. + DO i = lbound(vec,1), ubound(vec,1) + IF (.not. isnan_ud(vec(i))) THEN + n = n + 1 + s = s + vec(i) + ENDIF + ENDDO + + IF ((n > 0) .and. (n < size(vec))) THEN + s = s/n + DO i = lbound(vec,1), ubound(vec,1) + IF (isnan_ud(vec(i))) vec(i) = s + ENDDO + ENDIF + + END SUBROUTINE fillnan + END MODULE MOD_AggregationRequestData diff --git a/mksrfdata/MOD_LandHRU.F90 b/mksrfdata/MOD_LandHRU.F90 index 9f749b2c..91fe4d5f 100644 --- a/mksrfdata/MOD_LandHRU.F90 +++ b/mksrfdata/MOD_LandHRU.F90 @@ -110,6 +110,8 @@ SUBROUTINE landhru_build () CALL mpi_send (landelm%eindex, numelm, MPI_INTEGER8, p_address_master, mpi_tag_data, p_comm_glb, p_err) CALL mpi_recv (numhru, 1, MPI_INTEGER4, p_address_master, mpi_tag_size, p_comm_glb, p_stat, p_err) CALL mpi_recv (lakeid, numelm, MPI_INTEGER4, p_address_master, mpi_tag_data, p_comm_glb, p_stat, p_err) + ELSE + numhru = 0 ENDIF ENDIF #else @@ -131,11 +133,13 @@ SUBROUTINE landhru_build () IF (p_is_worker) THEN - allocate (landhru%eindex (numhru)) - allocate (landhru%settyp (numhru)) - allocate (landhru%ipxstt (numhru)) - allocate (landhru%ipxend (numhru)) - allocate (landhru%ielm (numhru)) + IF (numhru > 0) THEN + allocate (landhru%eindex (numhru)) + allocate (landhru%settyp (numhru)) + allocate (landhru%ipxstt (numhru)) + allocate (landhru%ipxend (numhru)) + allocate (landhru%ielm (numhru)) + ENDIF numhru = 0 diff --git a/mksrfdata/MOD_LandPatch.F90 b/mksrfdata/MOD_LandPatch.F90 index 5f7619bb..97269686 100644 --- a/mksrfdata/MOD_LandPatch.F90 +++ b/mksrfdata/MOD_LandPatch.F90 @@ -286,16 +286,6 @@ SUBROUTINE landpatch_build (lc_year) CALL landpatch%set_vecgs - IF (DEF_LANDONLY) THEN - IF ((p_is_worker) .and. (numpatch > 0)) THEN - allocate(msk(numpatch)) - msk = (landpatch%settyp /= 0) - ENDIF - - CALL landpatch%pset_pack (msk, numpatch) - - IF (allocated(msk)) deallocate(msk) - ENDIF #if (!defined(URBAN_MODEL) && !defined(CROP)) #ifdef USEMPI diff --git a/mksrfdata/MOD_SingleSrfdata.F90 b/mksrfdata/MOD_SingleSrfdata.F90 index 89068641..cedd1cb6 100644 --- a/mksrfdata/MOD_SingleSrfdata.F90 +++ b/mksrfdata/MOD_SingleSrfdata.F90 @@ -75,6 +75,8 @@ MODULE MOD_SingleSrfdata #endif real(r8), allocatable :: SITE_soil_BA_alpha (:) real(r8), allocatable :: SITE_soil_BA_beta (:) + + integer, allocatable :: SITE_soil_texture (:) real(r8) :: SITE_dbedrock = 0. @@ -279,6 +281,7 @@ SUBROUTINE read_surface_data_single (fsrfdata, mksrfdata) #endif CALL ncio_read_serial (fsrfdata, 'soil_BA_alpha ', SITE_soil_BA_alpha ) CALL ncio_read_serial (fsrfdata, 'soil_BA_beta ', SITE_soil_BA_beta ) + CALL ncio_read_serial (fsrfdata, 'soil_texture ', SITE_soil_texture ) ENDIF IF (DEF_USE_BEDROCK) THEN @@ -460,6 +463,7 @@ SUBROUTINE read_urban_surface_data_single (fsrfdata, mksrfdata, mkrun) #endif CALL ncio_read_serial (fsrfdata, 'soil_BA_alpha ', SITE_soil_BA_alpha ) CALL ncio_read_serial (fsrfdata, 'soil_BA_beta ', SITE_soil_BA_beta ) + CALL ncio_read_serial (fsrfdata, 'soil_texture ', SITE_soil_texture ) ENDIF IF (DEF_USE_BEDROCK) THEN @@ -639,6 +643,9 @@ SUBROUTINE write_surface_data_single (numpatch, numpft) CALL ncio_write_serial (fsrfdata, 'soil_BA_beta ', SITE_soil_BA_beta , 'soil') CALL ncio_put_attr (fsrfdata, 'soil_BA_alpha', 'source', source) CALL ncio_put_attr (fsrfdata, 'soil_BA_beta ', 'source', source) + + CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture, 'soil') + CALL ncio_put_attr (fsrfdata, 'soil_texture ', 'source', source) IF(DEF_USE_BEDROCK)THEN CALL ncio_write_serial (fsrfdata, 'depth_to_bedrock', SITE_dbedrock) @@ -835,6 +842,9 @@ SUBROUTINE write_urban_surface_data_single (numurban) CALL ncio_put_attr (fsrfdata, 'soil_BA_alpha', 'source', source) CALL ncio_put_attr (fsrfdata, 'soil_BA_beta ', 'source', source) + CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture, 'soil') + CALL ncio_put_attr (fsrfdata, 'soil_texture ', 'source', source) + IF(DEF_USE_BEDROCK)THEN CALL ncio_write_serial (fsrfdata, 'depth_to_bedrock', SITE_dbedrock) CALL ncio_put_attr (fsrfdata, 'depth_to_bedrock', 'source', datasource(USE_SITE_dbedrock)) @@ -936,6 +946,8 @@ SUBROUTINE single_srfdata_final () #endif IF (allocated(SITE_soil_BA_alpha )) deallocate(SITE_soil_BA_alpha ) IF (allocated(SITE_soil_BA_beta )) deallocate(SITE_soil_BA_beta ) + + IF (allocated(SITE_soil_texture )) deallocate(SITE_soil_texture ) IF (allocated(SITE_sf_lut )) deallocate(SITE_sf_lut ) IF (allocated(SITE_slp_type )) deallocate(SITE_slp_type )