diff --git a/Makefile b/Makefile old mode 100644 new mode 100755 index d2cbaeff..5aa3e86d --- a/Makefile +++ b/Makefile @@ -77,6 +77,7 @@ OBJS_MKSRFDATA = \ Aggregation_SoilParameters.o \ Aggregation_DBedrock.o \ Aggregation_Topography.o \ + Aggregation_TopographyFactors.o \ Aggregation_Urban.o \ MOD_MeshFilter.o \ MOD_RegionClip.o \ @@ -127,6 +128,7 @@ OBJS_BASIC = \ MOD_NdepData.o \ MOD_FireData.o \ MOD_OrbCoszen.o \ + MOD_OrbCosazi.o \ MOD_3DCanopyRadiation.o \ MOD_Aerosol.o \ MOD_SnowSnicar.o \ diff --git a/main/MOD_Forcing.F90 b/main/MOD_Forcing.F90 index 4b8ab955..e307efb3 100644 --- a/main/MOD_Forcing.F90 +++ b/main/MOD_Forcing.F90 @@ -51,7 +51,10 @@ MODULE MOD_Forcing type(pointer_real8_1d), allocatable :: forc_prc_grid (:) type(pointer_real8_1d), allocatable :: forc_prl_grid (:) type(pointer_real8_1d), allocatable :: forc_lwrad_grid(:) + type(pointer_real8_1d), allocatable :: forc_swrad_grid(:) type(pointer_real8_1d), allocatable :: forc_hgt_grid (:) + type(pointer_real8_1d), allocatable :: forc_us_grid (:) + type(pointer_real8_1d), allocatable :: forc_vs_grid (:) type(pointer_real8_1d), allocatable :: forc_t_part (:) type(pointer_real8_1d), allocatable :: forc_th_part (:) @@ -61,6 +64,9 @@ MODULE MOD_Forcing type(pointer_real8_1d), allocatable :: forc_prc_part (:) type(pointer_real8_1d), allocatable :: forc_prl_part (:) type(pointer_real8_1d), allocatable :: forc_frl_part (:) + type(pointer_real8_1d), allocatable :: forc_swrad_part (:) + type(pointer_real8_1d), allocatable :: forc_us_part (:) + type(pointer_real8_1d), allocatable :: forc_vs_part (:) logical, allocatable :: glacierss (:) @@ -128,6 +134,8 @@ SUBROUTINE forcing_init (dir_forcing, deltatime, ststamp, lc_year, etstamp) real(r8) :: missing_value integer :: ielm, istt, iend + integer :: iblkme, xblk, yblk, xloc, yloc + CALL init_user_specified_forcing ! CO2 data initialization @@ -237,7 +245,10 @@ SUBROUTINE forcing_init (dir_forcing, deltatime, ststamp, lc_year, etstamp) CALL mg2p_forc%allocate_part (forc_prc_grid ) CALL mg2p_forc%allocate_part (forc_prl_grid ) CALL mg2p_forc%allocate_part (forc_lwrad_grid ) + CALL mg2p_forc%allocate_part (forc_swrad_grid ) CALL mg2p_forc%allocate_part (forc_hgt_grid ) + CALL mg2p_forc%allocate_part (forc_us_grid ) + CALL mg2p_forc%allocate_part (forc_vs_grid ) CALL mg2p_forc%allocate_part (forc_t_part ) CALL mg2p_forc%allocate_part (forc_th_part ) @@ -247,6 +258,9 @@ SUBROUTINE forcing_init (dir_forcing, deltatime, ststamp, lc_year, etstamp) CALL mg2p_forc%allocate_part (forc_prc_part ) CALL mg2p_forc%allocate_part (forc_prl_part ) CALL mg2p_forc%allocate_part (forc_frl_part ) + CALL mg2p_forc%allocate_part (forc_swrad_part ) + CALL mg2p_forc%allocate_part (forc_us_part ) + CALL mg2p_forc%allocate_part (forc_vs_part ) CALL mg2p_forc%grid2part (topo_grid, forc_topo_grid ) CALL mg2p_forc%grid2part (maxelv_grid, forc_maxelv_grid) @@ -326,6 +340,7 @@ SUBROUTINE forcing_final () deallocate (forc_prc_grid ) deallocate (forc_prl_grid ) deallocate (forc_lwrad_grid ) + deallocate (forc_swrad_grid ) deallocate (forc_hgt_grid ) deallocate (forc_t_part ) @@ -336,6 +351,7 @@ SUBROUTINE forcing_final () deallocate (forc_prc_part ) deallocate (forc_prl_part ) deallocate (forc_frl_part ) + deallocate (forc_swrad_part ) ENDIF ENDIF @@ -355,11 +371,12 @@ END SUBROUTINE forcing_reset !-------------------------------- SUBROUTINE read_forcing (idate, dir_forcing) - + USE MOD_OrbCosazi USE MOD_Precision USE MOD_Namelist USE MOD_Const_Physical, only: rgas, grav USE MOD_Vars_TimeInvariants + USE MOD_Vars_TimeVariables, only: alb USE MOD_Vars_1DForcing USE MOD_Vars_2DForcing USE MOD_Block @@ -369,146 +386,144 @@ SUBROUTINE read_forcing (idate, dir_forcing) USE MOD_LandPatch USE MOD_RangeCheck USE MOD_UserSpecifiedForcing - USE MOD_ForcingDownscaling, only : rair, cpair, downscale_forcings_1c + USE MOD_ForcingDownscaling, only : rair, cpair, downscale_forcings + USE MOD_NetCDFVector IMPLICIT NONE + integer, intent(in) :: idate(3) character(len=*), intent(in) :: dir_forcing ! local variables: - integer :: ivar, istt, iend + integer :: ivar, istt, iend, id(3) integer :: iblkme, ib, jb, i, j, ilon, ilat, np, ipart, ne - real(r8) :: calday ! Julian cal day (1.xx to 365.xx) + real(r8) :: calday ! Julian cal day (1.xx to 365.xx) real(r8) :: sunang, cloud, difrat, vnrat real(r8) :: a, hsolar, ratio_rvrf type(block_data_real8_2d) :: forc_xy_solarin - + integer :: ii + character(10) :: cyear = "2005" + character(256):: lndname + type(timestamp) :: mtstamp - integer :: id(3) integer :: dtLB, dtUB - real(r8) :: cosz - integer :: year, month, mday + real(r8) :: cosz, coszen(numpatch), cosa, cosazi(numpatch), balb + INTEGER :: year, month, mday logical :: has_u,has_v - real solar, frl, prcp, tm, us, vs, pres, qm - real(r8) :: pco2m - - IF (p_is_io) THEN - - !------------------------------------------------------------ - ! READ in THE ATMOSPHERIC FORCING - - ! read lower and upper boundary forcing data - CALL metreadLBUB(idate, dir_forcing) - - ! set model time stamp - id(:) = idate(:) - !CALL adj2end(id) - mtstamp = id - - has_u = .true. - has_v = .true. - ! loop for variables - DO ivar = 1, NVAR - - IF (ivar == 5 .and. trim(vname(ivar)) == 'NULL') has_u = .false. - IF (ivar == 6 .and. trim(vname(ivar)) == 'NULL') has_v = .false. - IF (trim(vname(ivar)) == 'NULL') CYCLE ! no data, CYCLE - IF (trim(tintalgo(ivar)) == 'NULL') CYCLE - - ! to make sure the forcing data calculated is in the range of time - ! interval [LB, UB] - IF ( (mtstamp < tstamp_LB(ivar)) .or. (tstamp_UB(ivar) < mtstamp) ) THEN - write(6, *) "the data required is out of range! STOP!"; CALL CoLM_stop() - ENDIF + real(r8) :: pco2m + + IF (p_is_io) THEN + !------------------------------------------------------------ + ! READ in THE ATMOSPHERIC FORCING + ! read lower and upper boundary forcing data + CALL metreadLBUB(idate, dir_forcing) + ! set model time stamp + id(:) = idate(:) + !CALL adj2end(id) + mtstamp = id + has_u = .true. + has_v = .true. + ! loop for variables + DO ivar = 1, NVAR + IF (ivar == 5 .and. trim(vname(ivar)) == 'NULL') has_u = .false. + IF (ivar == 6 .and. trim(vname(ivar)) == 'NULL') has_v = .false. + IF (trim(vname(ivar)) == 'NULL') CYCLE ! no data, CYCLE + IF (trim(tintalgo(ivar)) == 'NULL') CYCLE + + ! to make sure the forcing data calculated is in the range of time + ! interval [LB, UB] + IF ( (mtstamp < tstamp_LB(ivar)) .or. (tstamp_UB(ivar) < mtstamp) ) THEN + write(6, *) "the data required is out of range! STOP!"; CALL CoLM_stop() + ENDIF - ! calcualte distance to lower/upper boundary - dtLB = mtstamp - tstamp_LB(ivar) - dtUB = tstamp_UB(ivar) - mtstamp + ! calcualte distance to lower/upper boundary + dtLB = mtstamp - tstamp_LB(ivar) + dtUB = tstamp_UB(ivar) - mtstamp - ! nearest method, for precipitation - IF (tintalgo(ivar) == 'nearest') THEN - IF (dtLB <= dtUB) THEN - CALL block_data_copy (forcn_LB(ivar), forcn(ivar)) - ELSE - CALL block_data_copy (forcn_UB(ivar), forcn(ivar)) - ENDIF - ENDIF + ! nearest method, for precipitation + IF (tintalgo(ivar) == 'nearest') THEN + IF (dtLB <= dtUB) THEN + CALL block_data_copy (forcn_LB(ivar), forcn(ivar)) + ELSE + CALL block_data_copy (forcn_UB(ivar), forcn(ivar)) + ENDIF + ENDIF - ! linear method, for T, Pres, Q, W, LW - IF (tintalgo(ivar) == 'linear') THEN - IF ( (dtLB+dtUB) > 0 ) THEN - CALL block_data_linear_interp ( & - forcn_LB(ivar), real(dtUB,r8)/real(dtLB+dtUB,r8), & - forcn_UB(ivar), real(dtLB,r8)/real(dtLB+dtUB,r8), & - forcn(ivar)) - ELSE - CALL block_data_copy (forcn_LB(ivar), forcn(ivar)) - ENDIF - ENDIF + ! linear method, for T, Pres, Q, W, LW + IF (tintalgo(ivar) == 'linear') THEN + IF ( (dtLB+dtUB) > 0 ) THEN + CALL block_data_linear_interp ( & + forcn_LB(ivar), real(dtUB,r8)/real(dtLB+dtUB,r8), & + forcn_UB(ivar), real(dtLB,r8)/real(dtLB+dtUB,r8), & + forcn(ivar)) + ELSE + CALL block_data_copy (forcn_LB(ivar), forcn(ivar)) + ENDIF + ENDIF - ! coszen method, for SW - IF (tintalgo(ivar) == 'coszen') THEN - DO iblkme = 1, gblock%nblkme - ib = gblock%xblkme(iblkme) - jb = gblock%yblkme(iblkme) + ! coszen method, for SW + IF (tintalgo(ivar) == 'coszen') THEN + DO iblkme = 1, gblock%nblkme + ib = gblock%xblkme(iblkme) + jb = gblock%yblkme(iblkme) - DO j = 1, gforc%ycnt(jb) - DO i = 1, gforc%xcnt(ib) + DO j = 1, gforc%ycnt(jb) + DO i = 1, gforc%xcnt(ib) - ilat = gforc%ydsp(jb) + j - ilon = gforc%xdsp(ib) + i - IF (ilon > gforc%nlon) ilon = ilon - gforc%nlon + ilat = gforc%ydsp(jb) + j + ilon = gforc%xdsp(ib) + i + IF (ilon > gforc%nlon) ilon = ilon - gforc%nlon - calday = calendarday(mtstamp) - cosz = orb_coszen(calday, gforc%rlon(ilon), gforc%rlat(ilat)) - cosz = max(0.001, cosz) - forcn(ivar)%blk(ib,jb)%val(i,j) = & - cosz / avgcos%blk(ib,jb)%val(i,j) * forcn_LB(ivar)%blk(ib,jb)%val(i,j) + calday = calendarday(mtstamp) + cosz = orb_coszen(calday, gforc%rlon(ilon), gforc%rlat(ilat)) + cosz = max(0.001, cosz) + forcn(ivar)%blk(ib,jb)%val(i,j) = & + cosz / avgcos%blk(ib,jb)%val(i,j) * forcn_LB(ivar)%blk(ib,jb)%val(i,j) - ENDDO - ENDDO ENDDO - ENDIF - + ENDDO ENDDO + ENDIF - ! preprocess for forcing data, only for QIAN data right now? - CALL metpreprocess (gforc, forcn) - - CALL allocate_block_data (gforc, forc_xy_solarin) - - CALL block_data_copy (forcn(1), forc_xy_t ) - CALL block_data_copy (forcn(2), forc_xy_q ) - CALL block_data_copy (forcn(3), forc_xy_psrf ) - CALL block_data_copy (forcn(3), forc_xy_pbot ) - CALL block_data_copy (forcn(4), forc_xy_prl, sca = 2/3._r8) - CALL block_data_copy (forcn(4), forc_xy_prc, sca = 1/3._r8) - CALL block_data_copy (forcn(7), forc_xy_solarin) - CALL block_data_copy (forcn(8), forc_xy_frl ) - IF (DEF_USE_CBL_HEIGHT) THEN - CALL block_data_copy (forcn(9), forc_xy_hpbl ) - ENDIF - - IF (has_u .and. has_v) THEN - CALL block_data_copy (forcn(5), forc_xy_us ) - CALL block_data_copy (forcn(6), forc_xy_vs ) - ELSEif (has_u) THEN - CALL block_data_copy (forcn(5), forc_xy_us , sca = 1/sqrt(2.0_r8)) - CALL block_data_copy (forcn(5), forc_xy_vs , sca = 1/sqrt(2.0_r8)) - ELSEif (has_v) THEN - CALL block_data_copy (forcn(6), forc_xy_us , sca = 1/sqrt(2.0_r8)) - CALL block_data_copy (forcn(6), forc_xy_vs , sca = 1/sqrt(2.0_r8)) - ELSE - IF (.not.trim(DEF_forcing%dataset) == 'CPL7') THEN - write(6, *) "At least one of the wind components must be provided! STOP!"; - CALL CoLM_stop() - ENDIF - ENDIF + ENDDO + + ! preprocess for forcing data, only for QIAN data right now? + CALL metpreprocess (gforc, forcn) + + CALL allocate_block_data (gforc, forc_xy_solarin) + + CALL block_data_copy (forcn(1), forc_xy_t ) + CALL block_data_copy (forcn(2), forc_xy_q ) + CALL block_data_copy (forcn(3), forc_xy_psrf ) + CALL block_data_copy (forcn(3), forc_xy_pbot ) + CALL block_data_copy (forcn(4), forc_xy_prl, sca = 2/3._r8) + CALL block_data_copy (forcn(4), forc_xy_prc, sca = 1/3._r8) + CALL block_data_copy (forcn(7), forc_xy_solarin) + CALL block_data_copy (forcn(8), forc_xy_frl ) + IF (DEF_USE_CBL_HEIGHT) THEN + CALL block_data_copy (forcn(9), forc_xy_hpbl ) + ENDIF + + IF (has_u .and. has_v) THEN + CALL block_data_copy (forcn(5), forc_xy_us ) + CALL block_data_copy (forcn(6), forc_xy_vs ) + ELSEIF (has_u) THEN + CALL block_data_copy (forcn(5), forc_xy_us , sca = 1/sqrt(2.0_r8)) + CALL block_data_copy (forcn(5), forc_xy_vs , sca = 1/sqrt(2.0_r8)) + ELSEIF (has_v) THEN + CALL block_data_copy (forcn(6), forc_xy_us , sca = 1/sqrt(2.0_r8)) + CALL block_data_copy (forcn(6), forc_xy_vs , sca = 1/sqrt(2.0_r8)) + ELSE + IF (.not.trim(DEF_forcing%dataset) == 'CPL7') THEN + write(6, *) "At least one of the wind components must be provided! STOP!"; + CALL CoLM_stop() + ENDIF + ENDIF - CALL flush_block_data (forc_xy_hgt_u, real(HEIGHT_V,r8)) - CALL flush_block_data (forc_xy_hgt_t, real(HEIGHT_T,r8)) - CALL flush_block_data (forc_xy_hgt_q, real(HEIGHT_Q,r8)) + CALL flush_block_data (forc_xy_hgt_u, real(HEIGHT_V,r8)) + CALL flush_block_data (forc_xy_hgt_t, real(HEIGHT_T,r8)) + CALL flush_block_data (forc_xy_hgt_q, real(HEIGHT_Q,r8)) IF (solarin_all_band) THEN @@ -545,49 +560,48 @@ SUBROUTINE read_forcing (idate, dir_forcing) ENDDO ELSE - !--------------------------------------------------------------- - ! as the downward solar is in full band, an empirical expression - ! will be used to divide fractions of band and incident - ! (visible, near-infrad, dirct, diffuse) - ! Julian calday (1.xx to 365.xx) - !--------------------------------------------------------------- - DO iblkme = 1, gblock%nblkme - ib = gblock%xblkme(iblkme) - jb = gblock%yblkme(iblkme) - - DO j = 1, gforc%ycnt(jb) - DO i = 1, gforc%xcnt(ib) - - ilat = gforc%ydsp(jb) + j - ilon = gforc%xdsp(ib) + i - IF (ilon > gforc%nlon) ilon = ilon - gforc%nlon - - a = forc_xy_solarin%blk(ib,jb)%val(i,j) - calday = calendarday(idate) - sunang = orb_coszen (calday, gforc%rlon(ilon), gforc%rlat(ilat)) - - cloud = (1160.*sunang-a)/(963.*sunang) - cloud = max(cloud,0.) - cloud = min(cloud,1.) - cloud = max(0.58,cloud) - - difrat = 0.0604/(sunang-0.0223)+0.0683 - IF(difrat.lt.0.) difrat = 0. - IF(difrat.gt.1.) difrat = 1. - - difrat = difrat+(1.0-difrat)*cloud - vnrat = (580.-cloud*464.)/((580.-cloud*499.)+(580.-cloud*464.)) - - forc_xy_sols %blk(ib,jb)%val(i,j) = a*(1.0-difrat)*vnrat - forc_xy_soll %blk(ib,jb)%val(i,j) = a*(1.0-difrat)*(1.0-vnrat) - forc_xy_solsd%blk(ib,jb)%val(i,j) = a*difrat*vnrat - forc_xy_solld%blk(ib,jb)%val(i,j) = a*difrat*(1.0-vnrat) - ENDDO - ENDDO - ENDDO - ENDIF - - ENDIF + !--------------------------------------------------------------- + ! as the downward solar is in full band, an empirical expression + ! will be used to divide fractions of band and incident + ! (visible, near-infrad, dirct, diffuse) + ! Julian calday (1.xx to 365.xx) + !--------------------------------------------------------------- + DO iblkme = 1, gblock%nblkme + ib = gblock%xblkme(iblkme) + jb = gblock%yblkme(iblkme) + + DO j = 1, gforc%ycnt(jb) + DO i = 1, gforc%xcnt(ib) + + ilat = gforc%ydsp(jb) + j + ilon = gforc%xdsp(ib) + i + IF (ilon > gforc%nlon) ilon = ilon - gforc%nlon + + a = forc_xy_solarin%blk(ib,jb)%val(i,j) + calday = calendarday(idate) + sunang = orb_coszen (calday, gforc%rlon(ilon), gforc%rlat(ilat)) + + cloud = (1160.*sunang-a)/(963.*sunang) + cloud = max(cloud,0.) + cloud = min(cloud,1.) + cloud = max(0.58,cloud) + + difrat = 0.0604/(sunang-0.0223)+0.0683 + IF(difrat.lt.0.) difrat = 0. + IF(difrat.gt.1.) difrat = 1. + + difrat = difrat+(1.0-difrat)*cloud + vnrat = (580.-cloud*464.)/((580.-cloud*499.)+(580.-cloud*464.)) + + forc_xy_sols %blk(ib,jb)%val(i,j) = a*(1.0-difrat)*vnrat + forc_xy_soll %blk(ib,jb)%val(i,j) = a*(1.0-difrat)*(1.0-vnrat) + forc_xy_solsd%blk(ib,jb)%val(i,j) = a*difrat*vnrat + forc_xy_solld%blk(ib,jb)%val(i,j) = a*difrat*(1.0-vnrat) + ENDDO + ENDDO + ENDDO + ENDIF + ENDIF ! [GET ATMOSPHERE CO2 CONCENTRATION DATA] year = idate(1) @@ -598,7 +612,6 @@ SUBROUTINE read_forcing (idate, dir_forcing) ENDIF - IF (.not. DEF_USE_Forcing_Downscaling) THEN ! Mapping the 2d atmospheric fields [lon_points]x[lat_points] @@ -653,41 +666,60 @@ SUBROUTINE read_forcing (idate, dir_forcing) ENDIF ELSE - - ! Mapping the 2d atmospheric fields [lon_points]x[lat_points] - ! -> the 1d vector of subgrid points [numelm] - CALL mg2p_forc%grid2pset (forc_xy_pco2m, forc_pco2m) - CALL mg2p_forc%grid2pset (forc_xy_po2m , forc_po2m ) - CALL mg2p_forc%grid2pset (forc_xy_us , forc_us ) - CALL mg2p_forc%grid2pset (forc_xy_vs , forc_vs ) - - CALL mg2p_forc%grid2pset (forc_xy_psrf , forc_psrf ) - - CALL mg2p_forc%grid2pset (forc_xy_sols , forc_sols ) - CALL mg2p_forc%grid2pset (forc_xy_soll , forc_soll ) - CALL mg2p_forc%grid2pset (forc_xy_solsd, forc_solsd) - CALL mg2p_forc%grid2pset (forc_xy_solld, forc_solld) - - CALL mg2p_forc%grid2pset (forc_xy_hgt_t, forc_hgt_t) - CALL mg2p_forc%grid2pset (forc_xy_hgt_u, forc_hgt_u) - CALL mg2p_forc%grid2pset (forc_xy_hgt_q, forc_hgt_q) + ! ------------------------------------------------------ + ! Forcing downscaling module + ! ------------------------------------------------------ + ! init forcing on patches + CALL mg2p_forc%grid2pset (forc_xy_pco2m, forc_pco2m) + CALL mg2p_forc%grid2pset (forc_xy_po2m , forc_po2m ) + CALL mg2p_forc%grid2pset (forc_xy_us , forc_us ) + CALL mg2p_forc%grid2pset (forc_xy_vs , forc_vs ) + CALL mg2p_forc%grid2pset (forc_xy_psrf , forc_psrf ) + CALL mg2p_forc%grid2pset (forc_xy_sols , forc_sols ) + CALL mg2p_forc%grid2pset (forc_xy_soll , forc_soll ) + CALL mg2p_forc%grid2pset (forc_xy_solsd, forc_solsd) + CALL mg2p_forc%grid2pset (forc_xy_solld, forc_solld) + CALL mg2p_forc%grid2pset (forc_xy_solarin, forc_swrad) + CALL mg2p_forc%grid2pset (forc_xy_hgt_t, forc_hgt_t) + CALL mg2p_forc%grid2pset (forc_xy_hgt_u, forc_hgt_u) + CALL mg2p_forc%grid2pset (forc_xy_hgt_q, forc_hgt_q) IF (DEF_USE_CBL_HEIGHT) THEN CALL mg2p_forc%grid2pset (forc_xy_hpbl, forc_hpbl) ENDIF - CALL mg2p_forc%grid2part (forc_xy_t , forc_t_grid ) - CALL mg2p_forc%grid2part (forc_xy_q , forc_q_grid ) - CALL mg2p_forc%grid2part (forc_xy_prc , forc_prc_grid ) - CALL mg2p_forc%grid2part (forc_xy_prl , forc_prl_grid ) - CALL mg2p_forc%grid2part (forc_xy_pbot , forc_pbot_grid ) - CALL mg2p_forc%grid2part (forc_xy_frl , forc_lwrad_grid) - CALL mg2p_forc%grid2part (forc_xy_hgt_t, forc_hgt_grid ) + ! Mapping the 2d atmospheric fields [lon_points]x[lat_points] + ! -> the 1d vector of subgrid points [numelm] + ! by selected mapping methods + CALL mg2p_forc%grid2part (forc_xy_t , forc_t_grid ) + CALL mg2p_forc%grid2part (forc_xy_q , forc_q_grid ) + CALL mg2p_forc%grid2part (forc_xy_prc , forc_prc_grid ) + CALL mg2p_forc%grid2part (forc_xy_prl , forc_prl_grid ) + CALL mg2p_forc%grid2part (forc_xy_pbot , forc_pbot_grid ) + CALL mg2p_forc%grid2part (forc_xy_frl , forc_lwrad_grid) + CALL mg2p_forc%grid2part (forc_xy_hgt_t, forc_hgt_grid ) + CALL mg2p_forc%grid2part (forc_xy_solarin, forc_swrad_grid) + CALL mg2p_forc%grid2part (forc_xy_us, forc_us_grid ) + CALL mg2p_forc%grid2part (forc_xy_vs, forc_vs_grid ) + + calday = calendarday(idate) + write(*,*) 'calday', calday IF (p_is_worker) THEN + DO np = 1, numpatch ! patches + + ! calculate albedo of each patches + IF (forc_sols(np)+forc_solsd(np)+forc_soll(np)+forc_solld(np) == 0) THEN + balb = 0 + ELSE + balb = (alb(1,1,np)*forc_sols(np) & + +alb(1,2,np)*forc_solsd(np) & + +alb(2,1,np)*forc_soll(np) & + +alb(2,2,np)*forc_solld(np)) & + /(forc_sols(np)+forc_solsd(np)+forc_soll(np)+forc_solld(np)) + ENDIF - DO np = 1, numpatch - DO ipart = 1, mg2p_forc%npart(np) + DO ipart = 1, mg2p_forc%npart(np) ! part loop of each patch IF (mg2p_forc%areapart(np)%val(ipart) == 0.) CYCLE @@ -695,10 +727,10 @@ SUBROUTINE read_forcing (idate, dir_forcing) ! the ground. Scientists have measured the most frigid temperature ever ! recorded on the continent's eastern highlands: about (180K) colder than ! dry ice. - IF(forc_t_grid(np)%val(ipart) < 180.) forc_t_grid(np)%val(ipart) = 180. + IF (forc_t_grid(np)%val(ipart) < 180.) forc_t_grid(np)%val(ipart) = 180. ! the highest air temp was found in Kuwait 326 K, Sulaibya 2012-07-31; ! Pakistan, Sindh 2010-05-26; Iraq, Nasiriyah 2011-08-03 - IF(forc_t_grid(np)%val(ipart) > 326.) forc_t_grid(np)%val(ipart) = 326. + IF (forc_t_grid(np)%val(ipart) > 326.) forc_t_grid(np)%val(ipart) = 326. forc_rho_grid(np)%val(ipart) = (forc_pbot_grid(np)%val(ipart) & - 0.378*forc_q_grid(np)%val(ipart)*forc_pbot_grid(np)%val(ipart) & @@ -707,29 +739,46 @@ SUBROUTINE read_forcing (idate, dir_forcing) forc_th_grid(np)%val(ipart) = forc_t_grid(np)%val(ipart) & * (1.e5/forc_pbot_grid(np)%val(ipart)) ** (rair/cpair) + ! caculate sun zenith angle and sun azimuth angle and turn to degree + coszen(np) = orb_coszen(calday, patchlonr(np), patchlatr(np)) + cosazi(np) = orb_cosazi(calday, patchlonr(np), patchlatr(np), coszen(np)) + + ! downscale forcing from grid to part + CALL downscale_forcings ( & + glacierss(np), & - CALL downscale_forcings_1c ( glacierss(np), & - ! forcing in gridcells + ! non-adjusted forcing forc_topo_grid(np)%val(ipart), forc_maxelv_grid(np)%val(ipart), & forc_t_grid(np)%val(ipart), forc_th_grid(np)%val(ipart), & forc_q_grid(np)%val(ipart), forc_pbot_grid(np)%val(ipart), & forc_rho_grid(np)%val(ipart), forc_prc_grid(np)%val(ipart), & forc_prl_grid(np)%val(ipart), forc_lwrad_grid(np)%val(ipart), & - forc_hgt_grid(np)%val(ipart), & - ! forcing in part of patches + forc_hgt_grid(np)%val(ipart), forc_swrad_grid(np)%val(ipart), & + forc_us_grid(np)%val(ipart), forc_vs_grid(np)%val(ipart), & + + ! topography-based factor on patch + slp_type_patches(:,np), asp_type_patches(:,np), area_type_patches(:,np), & + svf_patches(np), cur_patches(np), sf_lut_patches(:,:,np), & + + ! other factors + calday, coszen(np), cosazi(np), balb, & + + ! adjusted forcing forc_topo(np), forc_t_part(np)%val(ipart), & forc_th_part(np)%val(ipart), forc_q_part(np)%val(ipart), & forc_pbot_part(np)%val(ipart), forc_rhoair_part(np)%val(ipart), & forc_prc_part(np)%val(ipart), forc_prl_part(np)%val(ipart), & - forc_frl_part(np)%val(ipart)) - + forc_frl_part(np)%val(ipart), forc_swrad_part(np)%val(ipart), & + forc_us_part(np)%val(ipart), forc_vs_part(np)%val(ipart)) ENDDO ENDDO + ENDIF - ENDIF - - CALL mg2p_forc%normalize (forc_xy_frl, forc_frl_part) - + ! Conservation of short- and long- waves radiation in the grid of forcing + CALL mg2p_forc%normalize (forc_xy_solarin, forc_swrad_part) + CALL mg2p_forc%normalize (forc_xy_frl, forc_frl_part ) + + ! mapping parts to patches CALL mg2p_forc%part2pset (forc_t_part, forc_t ) CALL mg2p_forc%part2pset (forc_q_part, forc_q ) CALL mg2p_forc%part2pset (forc_pbot_part, forc_pbot ) @@ -737,8 +786,40 @@ SUBROUTINE read_forcing (idate, dir_forcing) CALL mg2p_forc%part2pset (forc_prc_part, forc_prc ) CALL mg2p_forc%part2pset (forc_prl_part, forc_prl ) CALL mg2p_forc%part2pset (forc_frl_part, forc_frl ) + CALL mg2p_forc%part2pset (forc_swrad_part, forc_swrad ) + CALL mg2p_forc%part2pset (forc_us_part, forc_us ) + CALL mg2p_forc%part2pset (forc_vs_part, forc_vs ) - ENDIF + ! divide fractions of downscaled shortwave radiation + IF (p_is_worker) THEN + DO j = 1, numpatch + a = forc_swrad(j) + IF (isnan(a)) a = 0 + calday = calendarday(idate) + sunang = orb_coszen (calday, patchlonr(j), patchlatr(j)) + IF (sunang.eq.0) THEN + cloud = 0. + ELSE + cloud = (1160.*sunang-a)/(963.*sunang) + ENDIF + cloud = max(cloud,0.0001) + cloud = min(cloud,1.) + cloud = max(0.58,cloud) + + difrat = 0.0604/(sunang-0.0223)+0.0683 + IF(difrat.lt.0.) difrat = 0. + IF(difrat.gt.1.) difrat = 1. + + difrat = difrat+(1.0-difrat)*cloud + vnrat = (580.-cloud*464.)/((580.-cloud*499.)+(580.-cloud*464.)) + + forc_sols(j) = a*(1.0-difrat)*vnrat + forc_soll(j) = a*(1.0-difrat)*(1.0-vnrat) + forc_solsd(j) = a*difrat*vnrat + forc_solld(j) = a*difrat*(1.0-vnrat) + ENDDO + ENDIF + ENDIF #ifdef RangeCheck #ifdef USEMPI diff --git a/main/MOD_ForcingDownscaling.F90 b/main/MOD_ForcingDownscaling.F90 index ce37a511..b25fc780 100644 --- a/main/MOD_ForcingDownscaling.F90 +++ b/main/MOD_ForcingDownscaling.F90 @@ -4,8 +4,7 @@ MODULE MOD_ForcingDownscaling !----------------------------------------------------------------------------- ! DESCRIPTION: -! Downscaling of gridcell-level meteorological forcings into column-level forcings -! (not included wind and solar radiation) +! Downscaling meteorological forcings ! ! INITIAL: ! The Community Land Model version 5.0 (CLM5.0) @@ -13,22 +12,23 @@ MODULE MOD_ForcingDownscaling ! REVISIONS: ! Updated by Yongjiu Dai, January 16, 2023 ! Revised by Lu Li, January 24, 2024 -! +! Revised by Sisi Chen, Lu Li, June, 2024 +!----------------------------------------------------------------------------- USE MOD_Precision USE MOD_Qsadv USE MOD_Namelist USE MOD_Const_Physical + USE MOD_Vars_Global + IMPLICIT NONE - real(r8), parameter :: SHR_CONST_MWDAIR = 28.966_r8 ! molecular weight dry air [kg/kmole] - real(r8), parameter :: SHR_CONST_MWWV = 18.016_r8 ! molecular weight water vapor - real(r8), parameter :: SHR_CONST_AVOGAD = 6.02214e26_r8 ! Avogadro's number [molecules/kmole] - real(r8), parameter :: SHR_CONST_BOLTZ = 1.38065e-23_r8 ! Boltzmann's constant [J/K/molecule] + real(r8), parameter :: SHR_CONST_MWDAIR = 28.966_r8 ! molecular weight dry air [kg/kmole] + real(r8), parameter :: SHR_CONST_MWWV = 18.016_r8 ! molecular weight water vapor + real(r8), parameter :: SHR_CONST_AVOGAD = 6.02214e26_r8 ! Avogadro's number [molecules/kmole] + real(r8), parameter :: SHR_CONST_BOLTZ = 1.38065e-23_r8 ! Boltzmann's constant [J/K/molecule] real(r8), parameter :: SHR_CONST_RGAS = SHR_CONST_AVOGAD*SHR_CONST_BOLTZ ! Universal gas constant [J/K/kmole] - real(r8), parameter :: SHR_CONST_RDAIR = SHR_CONST_RGAS/SHR_CONST_MWDAIR ! Dry air gas constant [J/K/kg] - - real(r8) :: rair = SHR_CONST_RDAIR ! Dry air gas constant [J/K/kg] + real(r8), parameter :: rair = SHR_CONST_RGAS/SHR_CONST_MWDAIR ! Dry air gas constant [J/K/kg] ! On the windward side of the range, annual mean lapse rates of 3.9-5.2 (deg km-1), ! substantially smaller than the often-assumed 6.5 (deg km-1). @@ -40,251 +40,20 @@ MODULE MOD_ForcingDownscaling ! ! Kunkel, K. E., 1989: Simple procedures for extrapolation of humidity variables in the mountainous western United States. ! J. Climate, 2, 656-669. lapse_rate = /Jan - Dec/4.4,5.9,7.1,7.8,8.1,8.2,8.1,8.1,7.7,6.8,5.5,4.7/ (deg km-1) - real(r8), parameter :: lapse_rate = 0.006_r8 ! surface temperature lapse rate (deg m-1) + real(r8), parameter :: lapse_rate = 0.006_r8 ! surface temperature lapse rate (deg m-1) SAVE ! PUBLIC MEMBER FUNCTIONS: - PUBLIC :: downscale_forcings ! Downscale atm forcing fields from gridcell to column - PUBLIC :: downscale_forcings_1c ! Downscale atm forcing fields from gridcell to column + PUBLIC :: downscale_forcings ! Downscale atmospheric forcing ! PRIVATE MEMBER FUNCTIONS: PRIVATE :: rhos ! calculate atmospheric density - PRIVATE :: downscale_longwave ! Downscale longwave radiation from gridcell to column - PRIVATE :: downscale_longwave_1c ! Downscale longwave radiation from gridcell to column - !----------------------------------------------------------------------------- CONTAINS -!----------------------------------------------------------------------------- - - SUBROUTINE downscale_forcings(& - num_gridcells,num_columns,begc,endc,glaciers,wt_column,& - - !slp_c, asp_c, cur_c, svf_c, sf_c,& - - forc_topo_g ,forc_t_g ,forc_th_g ,forc_q_g ,forc_pbot_g ,& - forc_rho_g ,forc_prc_g ,forc_prl_g ,forc_lwrad_g ,forc_hgt_grc,& - !forc_us_g ,forc_vs_g ,forc_swrad_g,& - - forc_topo_c ,forc_t_c ,forc_th_c ,forc_q_c ,forc_pbot_c ,& - forc_rho_c ,forc_prc_c ,forc_prl_c ,forc_lwrad_c) - !forc_swrad_c,forc_us_c ,forc_vs_c) - -!----------------------------------------------------------------------------- -! DESCRIPTION: -! Downscale atmospheric forcing fields from gridcell to column. -! -! Downscaling is done based on the difference between each land model column's elevation and -! the atmosphere's surface elevation (which is the elevation at which the atmospheric -! forcings are valid). -! -! Note that the downscaling procedure can result in changes in grid cell mean values -! compared to what was provided by the atmosphere. We conserve fluxes of mass and -! energy, but allow states such as temperature to differ. -!----------------------------------------------------------------------------- - - IMPLICIT NONE - - ! ARGUMENTS: - integer, intent(in) :: num_gridcells ! number of gridcell (element) - integer, intent(in) :: num_columns ! number of column (patch) - integer, intent(in) :: begc (1:num_gridcells) ! beginning column index - integer, intent(in) :: endc (1:num_gridcells) ! ending column index - logical, intent(in) :: glaciers (1:num_columns) ! true: glacier column (itypwat = 3) - real(r8), intent(in) :: wt_column(1:num_columns) ! weight of the column relative to the grid cell - - ! topography factor - !real(r8), intent(in) :: forc_slp_c(1:num_columns) ! slope factor - !real(r8), intent(in) :: forc_asp_c(1:num_columns) ! aspect factor - !real(r8), intent(in) :: forc_cur_c(1:num_columns) ! curvature factor - !real(r8), intent(in) :: forc_svf_c(1:num_columns) ! sky view factor - !real(r8), intent(in) :: forc_sf_c(1:num_columns) ! shadow factor - - ! Gridcell-level non-downscaled fields: - real(r8), intent(in) :: forc_topo_g (1:num_gridcells) ! atmospheric surface height [m] - real(r8), intent(in) :: forc_t_g (1:num_gridcells) ! atmospheric temperature [Kelvin] - real(r8), intent(in) :: forc_th_g (1:num_gridcells) ! atmospheric potential temperature [Kelvin] - real(r8), intent(in) :: forc_q_g (1:num_gridcells) ! atmospheric specific humidity [kg/kg] - real(r8), intent(in) :: forc_pbot_g (1:num_gridcells) ! atmospheric pressure [Pa] - real(r8), intent(in) :: forc_rho_g (1:num_gridcells) ! atmospheric density [kg/m**3] - real(r8), intent(in) :: forc_prc_g (1:num_gridcells) ! convective precipitation in grid [mm/s] - real(r8), intent(in) :: forc_prl_g (1:num_gridcells) ! large-scale precipitation in grid [mm/s] - real(r8), intent(in) :: forc_lwrad_g (1:num_gridcells) ! grid downward longwave [W/m**2] - real(r8), intent(in) :: forc_hgt_grc (1:num_gridcells) ! atmospheric reference height [m] - !real(r8), intent(in) :: forc_us_g (1:num_gridcells) ! atmospheric u wind [m/s] - !real(r8), intent(in) :: forc_vs_g (1:num_gridcells) ! atmospheric v wind [m/s] - !real(r8), intent(in) :: forc_swrad_g (1:num_gridcells) ! grid downward shortwave [W/m**2] - - ! Column-level downscaled fields: - real(r8), intent(in) :: forc_topo_c (1:num_columns) ! column surface height [m] - real(r8), intent(out) :: forc_t_c (1:num_columns) ! atmospheric temperature [Kelvin] - real(r8), intent(out) :: forc_th_c (1:num_columns) ! atmospheric potential temperature [Kelvin] - real(r8), intent(out) :: forc_q_c (1:num_columns) ! atmospheric specific humidity [kg/kg] - real(r8), intent(out) :: forc_pbot_c (1:num_columns) ! atmospheric pressure [Pa] - real(r8), intent(out) :: forc_rho_c (1:num_columns) ! atmospheric density [kg/m**3] - real(r8), intent(out) :: forc_prc_c (1:num_columns) ! column convective precipitation [mm/s] - real(r8), intent(out) :: forc_prl_c (1:num_columns) ! column large-scale precipitation [mm/s] - real(r8), intent(out) :: forc_lwrad_c(1:num_columns) ! column downward longwave [W/m**2] - !real(r8), intent(out) :: forc_swrad_c(1:num_columns) ! column downward shortwave [W/m**2] - !real(r8), intent(out) :: forc_us_c (1:num_columns) ! column u wind [m/s] - !real(r8), intent(out) :: forc_vs_c (1:num_columns) ! column v wind [m/s] - - ! Local variables for topo downscaling: - integer :: g,c ! indices - - real(r8) :: hsurf_g, hsurf_c - real(r8) :: Hbot, zbot - real(r8) :: tbot_g, pbot_g, thbot_g, qbot_g, qs_g, es_g, rhos_g - real(r8) :: tbot_c, pbot_c, thbot_c, qbot_c, qs_c, es_c, rhos_c - real(r8) :: rhos_c_estimate, rhos_g_estimate - real(r8) :: dum1, dum2 - - real(r8) :: max_elev_c ! the maximum column level elevation value within the grid - real(r8) :: delta_prc_c ! deviation of the column convective precipitation from the grid level precipitation - real(r8) :: delta_prl_c ! deviation of the column large-scale precipitation from the grid level precipitation - - !-------------------------------------------------------------------------- - - ! Initialize column forcing (needs to be done for ALL active columns) - DO g = 1, num_gridcells - - DO c = begc(g), endc(g) - forc_t_c (c) = forc_t_g (g) - forc_th_c (c) = forc_th_g (g) - forc_q_c (c) = forc_q_g (g) - forc_pbot_c (c) = forc_pbot_g (g) - forc_rho_c (c) = forc_rho_g (g) - forc_prc_c (c) = forc_prc_g (g) - forc_prl_c (c) = forc_prl_g (g) - forc_lwrad_c(c) = forc_lwrad_g(g) - END DO - END DO - - ! Downscale forc_t, forc_th, forc_q, forc_pbot, and forc_rho to columns. - DO g = 1, num_gridcells - hsurf_g = forc_topo_g(g) ! gridcell sfc elevation - tbot_g = forc_t_g(g) ! atm sfc temp - thbot_g = forc_th_g(g) ! atm sfc pot temp - qbot_g = forc_q_g(g) ! atm sfc spec humid - pbot_g = forc_pbot_g(g) ! atm sfc pressure - rhos_g = forc_rho_g(g) ! atm density - zbot = forc_hgt_grc(g) ! atm ref height - - max_elev_c = maxval(forc_topo_c(begc(g):endc(g))) ! maximum column level elevation value within the grid - - !real(r8) :: wd(1:num_columns) ! wind direction (rad) - !real(r8) :: slp_wd(1:num_columns) ! the slope in the direction of wind - !real(r8) :: norm_slp_wd ! normalize the slope in the direction of wind - !real(r8) :: norm_cur_c ! normalize curvature factor - - ! wind adjust factor - !wd = atan(forc_vs_g(g)/forc_us_g(g)) ! cal wind direction (rad) - !slp_wd = forc_slp_c(begc(g):endc(g))*cos(wd-forc_asp_c(begc(g):endc(g))) ! the slope in the direction of wind - !norm_slp_wd = slp_wd/(2*maxval(slp_wd)) ! normalize the slope in the direction of wind - !norm_cur_c = forc_cur_c(begc(g):endc(g))/(2*maxval(forc_cur_c(begc(g):endc(g)))) ! normalize curvature factor - - DO c = begc(g), endc(g) - ! This is a simple downscaling procedure - ! Note that forc_hgt, forc_u, forc_v and solar radiation are not downscaled. - - !asp_c = forc_asp_c(c) - !cur_c = forc_cur_c(c) - - hsurf_c = forc_topo_c(c) ! column sfc elevation - tbot_c = tbot_g-lapse_rate*(hsurf_c-hsurf_g) ! adjust temp for column - Hbot = rair*0.5_r8*(tbot_g+tbot_c)/grav ! scale ht at avg temp - pbot_c = pbot_g*exp(-(hsurf_c-hsurf_g)/Hbot) ! adjust press for column - - ! Derivation of potential temperature calculation: - ! - ! The textbook definition would be: - ! thbot_c = tbot_c * (p0/pbot_c)^(rair/cpair) - ! - ! Note that pressure is related to scale height as: - ! pbot_c = p0 * exp(-zbot/Hbot) - ! - ! Plugging this in to the textbook definition, then manipulating, we get: - ! thbot_c = tbot_c * (p0/(p0*exp(-zbot/Hbot)))^(rair/cpair) - ! = tbot_c * (1/exp(-zbot/Hbot))^(rair/cpair) - ! = tbot_c * (exp(zbot/Hbot))^(rair/cpair) - ! = tbot_c * exp((zbot/Hbot) * (rair/cpair)) - - ! But we want everything expressed in delta form, resulting in: - thbot_c = thbot_g + (tbot_c - tbot_g)*exp((zbot/Hbot)*(rair/cpair)) ! adjust pot temp for column - - CALL Qsadv(tbot_g,pbot_g,es_g,dum1,qs_g,dum2) ! es, qs for gridcell - CALL Qsadv(tbot_c,pbot_c,es_c,dum1,qs_c,dum2) ! es, qs for column - qbot_c = qbot_g*(qs_c/qs_g) ! adjust q for column - - rhos_c_estimate = rhos(qbot=qbot_c, pbot=pbot_c, tbot=tbot_c) - rhos_g_estimate = rhos(qbot=qbot_g, pbot=pbot_g, tbot=tbot_g) - rhos_c = rhos_g * (rhos_c_estimate / rhos_g_estimate) ! adjust density for column - - forc_t_c(c) = tbot_c - forc_th_c(c) = thbot_c - forc_q_c(c) = qbot_c - forc_pbot_c(c) = pbot_c - forc_rho_c(c) = rhos_c - - ! adjust precipitation - IF (trim(DEF_DS_precipitation_adjust_scheme) == 'I') THEN - ! Tesfa et al, 2020: Exploring Topography-Based Methods for Downscaling - ! Subgrid Precipitation for Use in Earth System Models. Equation (5) - ! https://doi.org/ 10.1029/2019JD031456 - - delta_prc_c = forc_prc_g(g) * (forc_topo_c(c) - forc_topo_g(g)) / max_elev_c - forc_prc_c(c) = forc_prc_g(g) + delta_prc_c ! convective precipitation [mm/s] - - delta_prl_c = forc_prl_g(g) * (forc_topo_c(c) - forc_topo_g(g)) / max_elev_c - forc_prl_c(c) = forc_prl_g(g) + delta_prl_c ! large scale precipitation [mm/s] - - ELSEIF (trim(DEF_DS_precipitation_adjust_scheme) == 'II') THEN - ! Liston, G. E. and Elder, K.: A meteorological distribution system - ! for high-resolution terrestrial modeling (MicroMet), J. Hydrometeorol., 7, 217-234, 2006. - ! Equation (33) and Table 1: chi range from January to December: - ! [0.35,0.35,0.35,0.30,0.25,0.20,0.20,0.20,0.20,0.25,0.30,0.35] (1/m) - - delta_prc_c = forc_prc_g(g) * 2.0*0.27e-3*(forc_topo_c(c) - forc_topo_g(g)) & - /(1.0 - 0.27e-3*(forc_topo_c(c) - forc_topo_g(g))) - forc_prc_c(c) = forc_prc_g(g) + delta_prc_c ! large scale precipitation [mm/s] - - delta_prl_c = forc_prl_g(g) * 2.0*0.27e-3*(forc_topo_c(c) - forc_topo_g(g)) & - /(1.0 - 0.27e-3*(forc_topo_c(c) - forc_topo_g(g))) - forc_prl_c(c) = forc_prl_g(g) + delta_prl_c ! large scale precipitation [mm/s] - - ELSEIF (trim(DEF_DS_precipitation_adjust_scheme) == 'III') THEN - ! Mei, Y., Maggioni, V., Houser, P., Xue, Y., & Rouf, T. (2020). A nonparametric statistical - ! technique for spatial downscaling of precipitation over High Mountain Asia. Water Resources Research, - ! 56, e2020WR027472. https://doi.org/ 10.1029/2020WR027472 - ! Change Random forest model to AutoML model. - !TODO: Lu Li; Need to done after all other forcings are downscaled - END IF - - IF (forc_prl_c(c) < 0) THEN - write(*,*) 'negative prl', forc_prl_g(g), max_elev_c, forc_topo_c(c), forc_topo_g(g) - END IF - - IF (forc_prc_c(c) < 0) THEN - write(*,*) 'negative prc', forc_prc_g(g), max_elev_c, forc_topo_c(c), forc_topo_g(g) - END IF - - forc_prc_c(c) = max(forc_prc_c(c), 0.) - forc_prl_c(c) = max(forc_prl_c(c), 0.) - - END DO - END DO - - CALL downscale_longwave(num_gridcells, num_columns, begc, endc, glaciers, wt_column, & - forc_topo_g, forc_t_g, forc_q_g, forc_pbot_g, forc_lwrad_g, & - forc_topo_c, forc_t_c, forc_q_c, forc_pbot_c, forc_lwrad_c) - - END SUBROUTINE downscale_forcings - - - !----------------------------------------------------------------------------- PURE FUNCTION rhos(qbot, pbot, tbot) @@ -297,7 +66,7 @@ PURE FUNCTION rhos(qbot, pbot, tbot) IMPLICIT NONE ! ARGUMENTS: - real(r8) :: rhos ! function result: atmospheric density (kg/m**3) + real(r8) :: rhos ! function result: atmospheric density (kg/m**3) real(r8), intent(in) :: qbot ! atmospheric specific humidity (kg/kg) real(r8), intent(in) :: pbot ! atmospheric pressure (Pa) real(r8), intent(in) :: tbot ! atmospheric temperature (K) @@ -306,7 +75,6 @@ PURE FUNCTION rhos(qbot, pbot, tbot) real(r8) :: egcm real(r8) :: wv_to_dair_weight_ratio ! ratio of molecular weight of water vapor to that of dry air [-] - !-------------------------------------------------------------------------- wv_to_dair_weight_ratio = SHR_CONST_MWWV/SHR_CONST_MWDAIR egcm = qbot*pbot / (wv_to_dair_weight_ratio + (1._r8 - wv_to_dair_weight_ratio)*qbot) @@ -314,194 +82,30 @@ PURE FUNCTION rhos(qbot, pbot, tbot) END FUNCTION rhos - - -!----------------------------------------------------------------------------- - - SUBROUTINE downscale_longwave(& - num_gridcells, num_columns, begc, endc, glaciers, wt_column, & - forc_topo_g, forc_t_g, forc_q_g, forc_pbot_g, forc_lwrad_g, & - forc_topo_c, forc_t_c, forc_q_c, forc_pbot_c, forc_lwrad_c) - -!----------------------------------------------------------------------------- -! DESCRIPTION: -! Downscale longwave radiation from gridcell to column -! Must be done AFTER temperature downscaling !----------------------------------------------------------------------------- - - IMPLICIT NONE - - ! ARGUMENTS: - integer, intent(in) :: num_gridcells ! number of gridcell - integer, intent(in) :: num_columns ! number of column - integer, intent(in) :: begc (1:num_gridcells) ! beginning column index - integer, intent(in) :: endc (1:num_gridcells) ! ending column index - logical, intent(in) :: glaciers (1:num_columns) ! true: glacier column - real(r8), intent(in) :: wt_column (1:num_columns) ! weight of the column relative to the grid cell - - real(r8), intent(in) :: forc_topo_g (1:num_gridcells) ! atmospheric surface height (m) - real(r8), intent(in) :: forc_t_g (1:num_gridcells) ! atmospheric temperature [Kelvin] - real(r8), intent(in) :: forc_q_g (1:num_gridcells) ! atmospheric specific humidity [kg/kg] - real(r8), intent(in) :: forc_pbot_g (1:num_gridcells) ! atmospheric pressure [Pa] - real(r8), intent(in) :: forc_lwrad_g (1:num_gridcells) ! downward longwave (W/m**2) - real(r8), intent(in) :: forc_topo_c (1:num_columns) ! column surface height (m) - real(r8), intent(in) :: forc_t_c (1:num_columns) ! atmospheric temperature [Kelvin] - real(r8), intent(in) :: forc_q_c (1:num_columns) ! atmospheric specific humidity [kg/kg] - real(r8), intent(in) :: forc_pbot_c (1:num_columns) ! atmospheric pressure [Pa] - real(r8), intent(out) :: forc_lwrad_c(1:num_columns) ! downward longwave (W/m**2) - - ! LOCAL VARIABLES: - integer :: c,g ! indices - real(r8) :: hsurf_c ! column-level elevation (m) - real(r8) :: hsurf_g ! gridcell-level elevation (m) - real(r8) :: sum_lwrad_g (1:num_gridcells) ! weighted sum of column-level lwrad - real(r8) :: sum_wts_g (1:num_gridcells) ! sum of weights that contribute to sum_lwrad_g - real(r8) :: lwrad_norm_g (1:num_gridcells) ! normalization factors - real(r8) :: newsum_lwrad_g (1:num_gridcells) ! weighted sum of column-level lwrad after normalization - - real(r8) :: pv_g ! the water vapor pressure at grid cell (hPa) - real(r8) :: pv_c ! the water vapor pressure at column (hPa) - real(r8) :: emissivity_clearsky_g ! clear-sky emissivity at grid cell - real(r8) :: emissivity_clearsky_c ! clear-sky emissivity at grid column - real(r8) :: emissivity_allsky_g ! all-sky emissivity at grid cell - real(r8) :: es_g, es_c, dum1, dum2, dum3 - - real(r8), parameter :: lapse_rate_longwave = 0.032_r8 ! longwave radiation lapse rate (W m-2 m-1) - real(r8), parameter :: longwave_downscaling_limit = 0.5_r8 ! relative limit for how much longwave downscaling can be done (unitless) - - !-------------------------------------------------------------------------- - - ! Initialize column forcing (needs to be done for ALL active columns) - DO g = 1, num_gridcells - DO c = begc(g), endc(g) - forc_lwrad_c(c) = forc_lwrad_g(g) - END DO - END DO - - ! Downscale the longwave radiation, conserving energy - - ! Initialize variables related to normalization - DO g = 1, num_gridcells - sum_lwrad_g(g) = 0._r8 - sum_wts_g(g) = 0._r8 - newsum_lwrad_g(g) = 0._r8 - END DO - - ! Do the downscaling - DO g = 1, num_gridcells - DO c = begc(g), endc(g) - - hsurf_g = forc_topo_g(g) - hsurf_c = forc_topo_c(c) - - IF (trim(DEF_DS_longwave_adjust_scheme) == 'I') THEN - ! Fiddes and Gruber, 2014, TopoSCALE v.1.0: downscaling gridded climate data in - ! complex terrain. Geosci. Model Dev., 7, 387-405. doi:10.5194/gmd-7-387-2014. - ! Equation (1) (2) (3); here, the empirical parameters x1 and x2 are different from - ! Konzelmann et al. (1994) where x1 = 0.443 and x2 = 8 (optimal for measurements on the Greenland ice sheet) - - CALL Qsadv(forc_t_g(g) ,forc_pbot_g(g) ,es_g,dum1,dum2,dum3) - CALL Qsadv(forc_t_c(c),forc_pbot_c(c),es_c,dum1,dum2,dum3) - pv_g = forc_q_g(g) *es_g/100._r8 ! (hPa) - pv_c = forc_q_c(c)*es_c/100._r8 ! (hPa) - - emissivity_clearsky_g = 0.23_r8 + 0.43_r8*(pv_g/forc_t_g(g))**(1._r8/5.7_r8) - emissivity_clearsky_c = 0.23_r8 + 0.43_r8*(pv_c/forc_t_c(c))**(1._r8/5.7_r8) - emissivity_allsky_g = forc_lwrad_g(g) / (5.67e-8_r8*forc_t_g(g)**4) - - forc_lwrad_c(c) = (emissivity_clearsky_c + (emissivity_allsky_g - emissivity_clearsky_g)) & - * 5.67e-8_r8*forc_t_c(c)**4 - ELSE - ! Longwave radiation is downscaled by assuming a linear decrease in downwelling longwave radiation - ! with increasing elevation (0.032 W m-2 m-1, limited to 0.5 - 1.5 times the gridcell mean value, - ! then normalized to conserve gridcell total energy) (Van Tricht et al., 2016, TC) Figure 6, - ! doi:10.5194/tc-10-2379-2016 - - IF (glaciers(c)) THEN - forc_lwrad_c(c) = forc_lwrad_g(g) - lapse_rate_longwave * (hsurf_c-hsurf_g) - - ! Here we assume that deltaLW = (dLW/dT)*(dT/dz)*deltaz - ! We get dLW/dT = 4*eps*sigma*T^3 = 4*LW/T from the Stefan-Boltzmann law, - ! evaluated at the mean temp. We assume the same temperature lapse rate as above. - - ELSE - forc_lwrad_c(c) = forc_lwrad_g(g) & - - 4.0_r8 * forc_lwrad_g(g)/(0.5_r8*(forc_t_c(c)+forc_t_g(g))) & - * lapse_rate * (hsurf_c - hsurf_g) - END IF - END IF - - ! But ensure that we don't depart too far from the atmospheric forcing value: - ! negative values of lwrad are certainly bad, but small positive values might - ! also be bad. We can especially run into trouble due to the normalization: a - ! small lwrad value in one column can lead to a big normalization factor, - ! leading to huge lwrad values in other columns. - - forc_lwrad_c(c) = min(forc_lwrad_c(c), & - forc_lwrad_g(g) * (1._r8 + longwave_downscaling_limit)) - forc_lwrad_c(c) = max(forc_lwrad_c(c), & - forc_lwrad_g(g) * (1._r8 - longwave_downscaling_limit)) - - ! Keep track of the gridcell-level weighted sum for later normalization. - ! This gridcell-level weighted sum just includes points for which we do the - ! downscaling (e.g., glc_mec points). Thus the contributing weights - ! generally do not add to 1. So to do the normalization properly, we also - ! need to keep track of the weights that have contributed to this sum. - - sum_lwrad_g(g) = sum_lwrad_g(g) + wt_column(c)*forc_lwrad_c(c) - sum_wts_g(g) = sum_wts_g(g) + wt_column(c) - END DO - - ! Normalize forc_lwrad_c(c) to conserve energy - IF (sum_wts_g(g) == 0._r8) THEN - lwrad_norm_g(g) = 1.0_r8 - ELSE IF (sum_lwrad_g(g) == 0._r8) THEN - lwrad_norm_g(g) = 1.0_r8 - ELSE ! The standard case - lwrad_norm_g(g) = forc_lwrad_g(g) / (sum_lwrad_g(g) / sum_wts_g(g)) - END IF - - DO c = begc(g), endc(g) - forc_lwrad_c(c) = forc_lwrad_c(c) * lwrad_norm_g(g) - newsum_lwrad_g(g) = newsum_lwrad_g(g) + wt_column(c)*forc_lwrad_c(c) - END DO - - END DO - - ! Make sure that, after normalization, the grid cell mean is conserved - DO g = 1, num_gridcells - IF (sum_wts_g(g) > 0._r8) THEN - IF (abs((newsum_lwrad_g(g) / sum_wts_g(g)) - forc_lwrad_g(g)) > 1.e-8_r8) THEN - write(6,*) 'g, newsum_lwrad_g, sum_wts_g, forc_lwrad_g: ', & - g, newsum_lwrad_g(g), sum_wts_g(g), forc_lwrad_g(g) - CALL abort - END IF - END IF - END DO - - END SUBROUTINE downscale_longwave - - - - !----------------------------------------------------------------------------- - SUBROUTINE downscale_forcings_1c (& - glaciers, & + SUBROUTINE downscale_forcings (& + glaciers, & - !slp_c, asp_c, cur_c, svf_c, sf_c,& + ! non-adjusted forcing + forc_topo_g ,forc_maxelv_g ,forc_t_g ,forc_th_g ,forc_q_g ,& + forc_pbot_g ,forc_rho_g ,forc_prc_g ,forc_prl_g ,forc_lwrad_g ,& + forc_hgt_g ,forc_swrad_g ,forc_us_g ,forc_vs_g , & - forc_topo_g ,forc_maxelv_g ,forc_t_g ,forc_th_g ,forc_q_g ,& - forc_pbot_g ,forc_rho_g ,forc_prc_g ,forc_prl_g ,forc_lwrad_g ,& - forc_hgt_grc,& - !forc_us_g ,forc_vs_g ,forc_swrad_g,& + ! topography-based factor on patch + slp_type_c, asp_type_c, area_type_c, svf_c, cur_c, sf_lut_c, & - forc_topo_c ,forc_t_c ,forc_th_c ,forc_q_c ,forc_pbot_c ,& - forc_rho_c ,forc_prc_c ,forc_prl_c ,forc_lwrad_c) - !forc_swrad_c,forc_us_c ,forc_vs_c) + ! other factors + julian_day, coszen, cosazi, alb, & + + ! adjusted forcing + forc_topo_c ,forc_t_c ,forc_th_c ,forc_q_c ,forc_pbot_c ,& + forc_rho_c ,forc_prc_c ,forc_prl_c ,forc_lwrad_c, forc_swrad_c, & + forc_us_c ,forc_vs_c) !----------------------------------------------------------------------------- ! DESCRIPTION: -! Downscale atmospheric forcing fields from gridcell to column. +! Downscale atmospheric forcing fields. ! ! Downscaling is done based on the difference between each land model column's elevation and ! the atmosphere's surface elevation (which is the elevation at which the atmospheric @@ -514,10 +118,25 @@ SUBROUTINE downscale_forcings_1c (& IMPLICIT NONE - ! ARGUMENTS: - logical, intent(in) :: glaciers ! true: glacier column (itypwat = 3) + integer, parameter :: S = 1370 ! solar constant (W/m**2) + real(r8), parameter :: thr = 85*PI/180 ! threshold of ?? - ! Gridcell-level non-downscaled fields: + ! ARGUMENTS: + logical, intent(in) :: glaciers ! true: glacier column (itypwat = 3) + real(r8), intent(in) :: julian_day ! day of year + real(r8), intent(in) :: coszen ! cosine of sun zenith angle at an hour + real(r8), intent(in) :: cosazi ! cosine of sun azimuth angle at an hour + real(r8), intent(in) :: alb ! blue sky albedo + + ! topography-based factor + real(r8), intent(in) :: svf_c ! sky view factor + real(r8), intent(in) :: cur_c ! curvature + real(r8), intent(in) :: sf_lut_c (1:num_azimuth,1:num_zenith) ! look up table of shadow mask of a patch + real(r8), intent(in) :: asp_type_c (1:num_type) ! topographic aspect of each type of one patch + real(r8), intent(in) :: slp_type_c (1:num_type) ! topographic slope of each character of one patch + real(r8), intent(in) :: area_type_c(1:num_type) ! area percentage of each character of one patch + + ! non-downscaled fields: real(r8), intent(in) :: forc_topo_g ! atmospheric surface height [m] real(r8), intent(in) :: forc_maxelv_g ! max atmospheric surface height [m] real(r8), intent(in) :: forc_t_g ! atmospheric temperature [Kelvin] @@ -528,9 +147,12 @@ SUBROUTINE downscale_forcings_1c (& real(r8), intent(in) :: forc_prc_g ! convective precipitation in grid [mm/s] real(r8), intent(in) :: forc_prl_g ! large-scale precipitation in grid [mm/s] real(r8), intent(in) :: forc_lwrad_g ! grid downward longwave [W/m**2] - real(r8), intent(in) :: forc_hgt_grc ! atmospheric reference height [m] + real(r8), intent(in) :: forc_swrad_g ! grid downward shortwave [W/m**2] + real(r8), intent(in) :: forc_hgt_g ! atmospheric reference height [m] + real(r8), intent(in) :: forc_us_g ! eastward wind [m/s] + real(r8), intent(in) :: forc_vs_g ! northward wind [m/s] - ! Column-level downscaled fields: + ! downscaled fields: real(r8), intent(in) :: forc_topo_c ! column surface height [m] real(r8), intent(out) :: forc_t_c ! atmospheric temperature [Kelvin] real(r8), intent(out) :: forc_th_c ! atmospheric potential temperature [Kelvin] @@ -540,40 +162,38 @@ SUBROUTINE downscale_forcings_1c (& real(r8), intent(out) :: forc_prc_c ! column convective precipitation [mm/s] real(r8), intent(out) :: forc_prl_c ! column large-scale precipitation [mm/s] real(r8), intent(out) :: forc_lwrad_c ! column downward longwave [W/m**2] - - ! Local variables for topo downscaling: - + real(r8), intent(out) :: forc_swrad_c ! column downward shortwave [W/m**2] + real(r8), intent(out) :: forc_us_c ! column eastward wind [m/s] + real(r8), intent(out) :: forc_vs_c ! column northward wind [m/s] + + ! Local variables for topo downscaling: real(r8) :: hsurf_g, hsurf_c real(r8) :: Hbot, zbot real(r8) :: tbot_g, pbot_g, thbot_g, qbot_g, qs_g, es_g, rhos_g real(r8) :: tbot_c, pbot_c, thbot_c, qbot_c, qs_c, es_c, rhos_c real(r8) :: rhos_c_estimate, rhos_g_estimate real(r8) :: dum1, dum2 - real(r8) :: max_elev_c ! the maximum column level elevation value within the grid real(r8) :: delta_prc_c ! deviation of the column convective precipitation from the grid level precipitation real(r8) :: delta_prl_c ! deviation of the column large-scale precipitation from the grid level precipitation +!----------------------------------------------------------------------------- - ! ! Downscale forc_t, forc_th, forc_q, forc_pbot, and forc_rho to columns. + ! -------------------------------------------------------------------------------------- + ! 1. adjust air temperature, potential temperature, specific humidity, pressure, density + ! -------------------------------------------------------------------------------------- hsurf_g = forc_topo_g ! gridcell sfc elevation tbot_g = forc_t_g ! atm sfc temp thbot_g = forc_th_g ! atm sfc pot temp qbot_g = forc_q_g ! atm sfc spec humid pbot_g = forc_pbot_g ! atm sfc pressure rhos_g = forc_rho_g ! atm density - zbot = forc_hgt_grc ! atm ref height - - ! This is a simple downscaling procedure - ! Note that forc_hgt, forc_u, forc_v and solar radiation are not downscaled. + zbot = forc_hgt_g ! atm ref height - !asp_c = forc_asp_c(c) - !cur_c = forc_cur_c(c) - - hsurf_c = forc_topo_c ! column sfc elevation - tbot_c = tbot_g-lapse_rate*(hsurf_c-hsurf_g) ! adjust temp for column + hsurf_c = forc_topo_c ! column sfc elevation + tbot_c = tbot_g-lapse_rate*(hsurf_c-hsurf_g) ! adjust [temp] for column Hbot = rair*0.5_r8*(tbot_g+tbot_c)/grav ! scale ht at avg temp - pbot_c = pbot_g*exp(-(hsurf_c-hsurf_g)/Hbot) ! adjust press for column + pbot_c = pbot_g*exp(-(hsurf_c-hsurf_g)/Hbot) ! adjust [press] for column ! Derivation of potential temperature calculation: ! @@ -590,23 +210,46 @@ SUBROUTINE downscale_forcings_1c (& ! = tbot_c * exp((zbot/Hbot) * (rair/cpair)) ! But we want everything expressed in delta form, resulting in: - thbot_c = thbot_g + (tbot_c - tbot_g)*exp((zbot/Hbot)*(rair/cpair)) ! adjust pot temp for column + thbot_c = thbot_g + (tbot_c - tbot_g)*exp((zbot/Hbot)*(rair/cpair)) ! adjust [pot temp] for column CALL Qsadv(tbot_g,pbot_g,es_g,dum1,qs_g,dum2) ! es, qs for gridcell CALL Qsadv(tbot_c,pbot_c,es_c,dum1,qs_c,dum2) ! es, qs for column - qbot_c = qbot_g*(qs_c/qs_g) ! adjust q for column + qbot_c = qbot_g*(qs_c/qs_g) ! adjust [q] for column rhos_c_estimate = rhos(qbot=qbot_c, pbot=pbot_c, tbot=tbot_c) rhos_g_estimate = rhos(qbot=qbot_g, pbot=pbot_g, tbot=tbot_g) - rhos_c = rhos_g * (rhos_c_estimate / rhos_g_estimate) ! adjust density for column + rhos_c = rhos_g * (rhos_c_estimate / rhos_g_estimate) ! adjust [density] for column + ! save forc_t_c = tbot_c forc_th_c = thbot_c forc_q_c = qbot_c forc_pbot_c = pbot_c forc_rho_c = rhos_c - ! adjust precipitation + ! -------------------------------------------------------------------------------------- + ! 2. adjust wind speed + ! -------------------------------------------------------------------------------------- + CALL downscale_wind(forc_us_g, forc_vs_g, & + forc_us_c, forc_vs_c, & + slp_type_c, asp_type_c, area_type_c, cur_c) + + ! -------------------------------------------------------------------------------------- + ! 3. adjust longwave radiation and shortwave radiation + ! -------------------------------------------------------------------------------------- + CALL downscale_longwave (glaciers, & + forc_topo_g, forc_t_g, forc_q_g, forc_pbot_g, forc_lwrad_g, & + forc_topo_c, forc_t_c, forc_q_c, forc_pbot_c, forc_lwrad_c) + + CALL downscale_shortwave(& + forc_topo_g, forc_pbot_g, forc_swrad_g, & + forc_topo_c, forc_pbot_c, forc_swrad_c, & + julian_day, coszen, cosazi, alb, & + slp_type_c, asp_type_c, svf_c, sf_lut_c, area_type_c) + + ! -------------------------------------------------------------------------------------- + ! 4. adjust precipitation + ! -------------------------------------------------------------------------------------- IF (trim(DEF_DS_precipitation_adjust_scheme) == 'I') THEN ! Tesfa et al, 2020: Exploring Topography-Based Methods for Downscaling ! Subgrid Precipitation for Use in Earth System Models. Equation (5) @@ -636,8 +279,8 @@ SUBROUTINE downscale_forcings_1c (& ! Mei, Y., Maggioni, V., Houser, P., Xue, Y., & Rouf, T. (2020). A nonparametric statistical ! technique for spatial downscaling of precipitation over High Mountain Asia. Water Resources Research, ! 56, e2020WR027472. https://doi.org/ 10.1029/2020WR027472 - ! Change Random forest model to AutoML model. - !TODO: Lu Li; Need to done after all other forcings are downscaled + + ! We implement this scheme in MOD_Forcing.F90 END IF IF (forc_prl_c < 0) THEN @@ -650,25 +293,81 @@ SUBROUTINE downscale_forcings_1c (& forc_prc_c = 0. END IF + END SUBROUTINE downscale_forcings + +!----------------------------------------------------------------------------- + + SUBROUTINE downscale_wind(forc_us_g, forc_vs_g, & + forc_us_c, forc_vs_c, & + slp_type_c, asp_type_c, area_type_c, cur_c) + +!----------------------------------------------------------------------------- +! DESCRIPTION: +! Downscale wind speed +! +! Liston, G. E. and Elder, K.: A meteorological distribution system +! for high-resolution terrestrial modeling (MicroMet), J. Hydrometeorol., 7, 217-234, 2006. +!----------------------------------------------------------------------------- + + IMPLICIT NONE + + ! ARGUMENTS: + real(r8), intent(in) :: forc_us_g ! eastward wind (m/s) + real(r8), intent(in) :: forc_vs_g ! northward wind (m/s) + real(r8), intent(out) :: forc_us_c ! adjusted eastward wind (m/s) + real(r8), intent(out) :: forc_vs_c ! adjusted northward wind (m/s) + + real(r8), intent(in) :: cur_c ! curvature + real(r8), intent(in) :: asp_type_c (1:num_type) ! topographic aspect of each character of one patch + real(r8), intent(in) :: slp_type_c (1:num_type) ! topographic slope of each character of one patch + real(r8), intent(in) :: area_type_c (1:num_type) ! area percentage of each character of one patch + + ! local variables + real(r8) :: wind_dir ! wind direction + real(r8) :: ws_g ! non-downscaled wind speed + real(r8) :: wind_dir_slp (1:num_type) ! the slope in the direction of the wind + real(r8) :: ws_c_type(1:num_type) ! downscaled wind speed of each type in each patch + real(r8) :: ws_c ! downscaled wind speed + integer :: g, c, i + +!----------------------------------------------------------------------------- + + ! calculate wind direction + IF (forc_us_g == 0.) THEN + wind_dir = 0 + ELSE + wind_dir = atan(forc_vs_g /forc_us_g) + ENDIF + + ! non-adjusted wind speed + ws_g = sqrt(forc_vs_g *forc_vs_g +forc_us_g *forc_us_g ) - CALL downscale_longwave_1c (glaciers, & - forc_topo_g, forc_t_g, forc_q_g, forc_pbot_g, forc_lwrad_g, & - forc_topo_c, forc_t_c, forc_q_c, forc_pbot_c, forc_lwrad_c) + ! compute the slope in the direction of the wind + DO i = 1, num_type + wind_dir_slp(i) = slp_type_c(i)*cos(wind_dir-asp_type_c(i)*PI/180) + ENDDO - END SUBROUTINE downscale_forcings_1c + ! compute wind speed ajustment + DO i = 1, num_type + ws_c_type(i) = ws_g *(1+(0.58*wind_dir_slp(i))+0.42*cur_c)*area_type_c(i) + ENDDO + ! adjusted wind speed + ws_c = sum(ws_c_type(:)) + forc_us_c = ws_c*cos(wind_dir) + forc_vs_c = ws_c*sin(wind_dir) + END SUBROUTINE downscale_wind !----------------------------------------------------------------------------- - SUBROUTINE downscale_longwave_1c (glaciers, & + SUBROUTINE downscale_longwave (glaciers, & forc_topo_g, forc_t_g, forc_q_g, forc_pbot_g, forc_lwrad_g, & forc_topo_c, forc_t_c, forc_q_c, forc_pbot_c, forc_lwrad_c) !----------------------------------------------------------------------------- ! DESCRIPTION: -! Downscale longwave radiation from gridcell to column -! Must be done AFTER temperature downscaling +! Downscale longwave radiation !----------------------------------------------------------------------------- IMPLICIT NONE @@ -681,11 +380,12 @@ SUBROUTINE downscale_longwave_1c (glaciers, & real(r8), intent(in) :: forc_q_g ! atmospheric specific humidity [kg/kg] real(r8), intent(in) :: forc_pbot_g ! atmospheric pressure [Pa] real(r8), intent(in) :: forc_lwrad_g ! downward longwave (W/m**2) + real(r8), intent(in) :: forc_topo_c ! column surface height (m) real(r8), intent(in) :: forc_t_c ! atmospheric temperature [Kelvin] real(r8), intent(in) :: forc_q_c ! atmospheric specific humidity [kg/kg] real(r8), intent(in) :: forc_pbot_c ! atmospheric pressure [Pa] - real(r8), intent(out) :: forc_lwrad_c ! downward longwave (W/m**2) + real(r8), intent(out):: forc_lwrad_c ! downward longwave (W/m**2) ! LOCAL VARIABLES: real(r8) :: hsurf_c ! column-level elevation (m) @@ -701,13 +401,10 @@ SUBROUTINE downscale_longwave_1c (glaciers, & real(r8), parameter :: lapse_rate_longwave = 0.032_r8 ! longwave radiation lapse rate (W m-2 m-1) real(r8), parameter :: longwave_downscaling_limit = 0.5_r8 ! relative limit for how much longwave downscaling can be done (unitless) - !-------------------------------------------------------------------------- +!-------------------------------------------------------------------------- - ! Initialize column forcing (needs to be done for ALL active columns) + ! Initialize (needs to be done for ALL active columns) forc_lwrad_c = forc_lwrad_g - - ! Do the downscaling - hsurf_g = forc_topo_g hsurf_c = forc_topo_c @@ -746,8 +443,8 @@ SUBROUTINE downscale_longwave_1c (glaciers, & forc_lwrad_c = forc_lwrad_g & - 4.0_r8 * forc_lwrad_g/(0.5_r8*(forc_t_c+forc_t_g)) & * lapse_rate * (hsurf_c - hsurf_g) - END IF - END IF + ENDIF + ENDIF ! But ensure that we don't depart too far from the atmospheric forcing value: ! negative values of lwrad are certainly bad, but small positive values might @@ -758,7 +455,187 @@ SUBROUTINE downscale_longwave_1c (glaciers, & forc_lwrad_c = min(forc_lwrad_c, forc_lwrad_g * (1._r8 + longwave_downscaling_limit)) forc_lwrad_c = max(forc_lwrad_c, forc_lwrad_g * (1._r8 - longwave_downscaling_limit)) + END SUBROUTINE downscale_longwave + +!----------------------------------------------------------------------------- + + SUBROUTINE downscale_shortwave( & + forc_topo_g, forc_pbot_g, forc_swrad_g, & + forc_topo_c, forc_pbot_c, forc_swrad_c, & + julian_day, coszen, cosazi, alb, & + slp_type_c, asp_type_c, svf_c, sf_lut_c, area_type_c) + +!----------------------------------------------------------------------------- +! DESCRIPTION: +! +! Rouf, T., Mei, Y., Maggioni, V., Houser, P., & Noonan, M. (2020). A Physically Based +! Atmospheric Variables Downscaling Technique. Journal of Hydrometeorology, +! 21(1), 93–108. https://doi.org/10.1175/JHM-D-19-0109.1 +! +! Sisi Chen, Lu Li, Yongjiu Dai, et al. Exploring Topography Downscaling Methods for +! Hyper-Resolution Land Surface Modeling. Authorea. April 25, 2024. +! DOI: 10.22541/au.171403656.68476353/v1 +! +! Must be done after downscaling of surface pressure +!----------------------------------------------------------------------------- + + IMPLICIT NONE + + integer, parameter :: S = 1370 ! solar constant (W/m**2) + real(r8), parameter :: thr = 85*PI/180 ! threshold of ?? + real(r8), parameter :: shortwave_downscaling_limit = 0.5_r8 ! relative limit for how much shortwave downscaling can be done (unitless) + + ! ARGUMENTS: + real(r8), intent(in) :: julian_day ! day of year + real(r8), intent(in) :: coszen ! zenith angle at an hour + real(r8), intent(in) :: cosazi ! azimuth angle at an hour + real(r8), intent(in) :: alb ! blue sky albedo + + real(r8), intent(in) :: forc_topo_g ! atmospheric surface height (m) + real(r8), intent(in) :: forc_pbot_g ! atmospheric pressure [Pa] + real(r8), intent(in) :: forc_swrad_g ! downward shortwave (W/m**2) + + real(r8), intent(in) :: forc_topo_c ! column surface height (m) + real(r8), intent(in) :: forc_pbot_c ! atmospheric pressure [Pa] + real(r8), intent(out):: forc_swrad_c ! downward shortwave (W/m**2) + + real(r8), intent(in) :: svf_c ! sky view factor + real(r8), intent(in) :: sf_lut_c (1:num_azimuth,1:num_zenith) ! look up table of shadow factor + real(r8), intent(in) :: asp_type_c (1:num_type) ! topographic aspect of each character of one patch (°) + real(r8), intent(in) :: slp_type_c (1:num_type) ! topographic slope of each character of one patch + real(r8), intent(in) :: area_type_c(1:num_type) ! area percentage of each character of one patch + + ! LOCAL VARIABLES: + real(r8) :: zen_rad, azi_rad, zen_deg, azi_deg ! rad and deg of sun zenith and azimuth angles + integer :: idx_azi, idx_zen ! index used to cal shadow factor from look up table + real(r8) :: sf_c ! shadow factor + real(r8) :: rt_R ! The ratio of the current distance between the sun and the earth ! to the average distance between the sun and the earth + real(r8) :: toa_swrad ! top of atmosphere shortwave radiation + real(r8) :: clr_idx ! atmospheric transparency + real(r8) :: diff_wgt ! diffuse weight + real(r8) :: k_c ! column broadband attenuation coefficient [Pa^-1] + real(r8) :: opt_factor ! optical length factor + real(r8) :: a_p + real(r8) :: svf, balb + + real(r8) :: diff_swrad_g, beam_swrad_g ! diffuse and beam radiation + real(r8) :: diff_swrad_c, beam_swrad_c, refl_swrad_c! downscaled diffuse and beam radiation + real(r8) :: beam_swrad_type (1:num_type) ! beam radiation of one characterized patch + real(r8) :: refl_swrad_type (1:num_type) ! reflect radiation of one characterized patch + real(r8) :: tcf_type (1:num_type) ! terrain configure factor + real(r8) :: cosill_type (1:num_type) ! illumination angle (cos) at defined types + + integer :: i + +!----------------------------------------------------------------------------- + + ! calculate shadow factor according to sun zenith and azimuth angle + zen_rad = acos(coszen) + azi_rad = acos(cosazi) + zen_deg = zen_rad*180/PI ! turn deg + azi_deg = azi_rad*180.0/PI ! turn deg + + idx_azi = INT(azi_deg*num_azimuth/360) + idx_zen = INT(zen_deg*num_zenith/90) + IF (idx_azi==0) idx_azi = 1 + IF (idx_zen==0) idx_zen = 1 + + sf_c = sf_lut_c(idx_azi, idx_zen) + IF (sf_c<0) sf_c = 0 + IF (sf_c>1) sf_c = 1 + + ! calculate top-of-atmosphere incident shortwave radiation + rt_R = 1-0.01672*cos(0.9856*(julian_day-4)) + toa_swrad = S*(rt_R**2)*coszen + + ! calculate clearness index + IF (toa_swrad.le.0) THEN + clr_idx = 1 + ELSE + clr_idx = forc_swrad_g/toa_swrad + ENDIF + IF (clr_idx>1) clr_idx = 1 + + ! calculate diffuse weight + ! Ruiz-Arias, J. A., Alsamamra, H., Tovar-Pescador, J., & Pozo-Vázquez, D. (2010). + ! Proposal of a regressive model for the hourly diffuse solar radiation under all sky + ! conditions. Energy Conversion and Management, 51(5), 881–893. + ! https://doi.org/10.1016/j.enconman.2009.11.024 + diff_wgt = 0.952-1.041*exp(-1*exp(2.3-4.702*clr_idx)) + IF (diff_wgt>1) diff_wgt = 1 + IF (diff_wgt<0) diff_wgt = 0 + + ! calculate diffuse and beam radiation + diff_swrad_g = forc_swrad_g*diff_wgt + beam_swrad_g = forc_swrad_g*(1-diff_wgt) + + ! calcualte broadband attenuation coefficient [Pa^-1] + IF (clr_idx.le.0) THEN + k_c = 0 + ELSE + k_c = log(clr_idx)/forc_pbot_c + ENDIF + + ! calculate factor to account for the difference of optical path length due to pressure difference + opt_factor = exp(k_c*(forc_pbot_g-forc_pbot_c)) + ! Control the boundary of optical path length + IF ((opt_factor>10000).or.(opt_factor<-10000)) opt_factor = 0 + + ! Adjust the zenith angle so that the range of zenith angles is less than 85° + IF (zen_rad>thr) zen_rad=thr + + ! loop for four defined types to downscale beam radiation + DO i = 1, num_type + ! calculate the cosine of solar illumination angle, cos(θ), + ! ranging between −1 and 1, indicates if the sun is below or + ! above the local horizon (note that values lower than 0 are set to 0 indicate self shadow) + cosill_type(i) = cos(slp_type_c(i))+tan(zen_rad)*sin(slp_type_c(i))*cos(asp_type_c(i)*PI/180) + IF (cosill_type(i)>1) cosill_type(i) = 1 + IF (cosill_type(i)<0) cosill_type(i) = 0 + + ! downscaling beam radiation + a_p = area_type_c(i) + IF (a_p.gt.1.0) a_p = 1 + IF (a_p.lt.0) a_p = 0 + beam_swrad_type(i) = sf_c*cosill_type(i)*opt_factor*a_p*beam_swrad_g + ENDDO + beam_swrad_c = sum(beam_swrad_type) + + ! downscaling diffuse radiation + svf = svf_c + IF (svf>1) svf = 1 + IF (svf<0) svf = 0 + diff_swrad_c = svf*diff_swrad_g + + ! downscaling reflected radiation + balb = alb + DO i = 1, num_type + tcf_type(i) = (1+cos(slp_type_c(i)))/2-svf + IF (tcf_type(i)<0) tcf_type(i) = 0 + + IF (isnan(alb)) THEN + refl_swrad_type(i) = -1.0e36 + ELSE + IF ((balb<0).or.(balb>1)) balb = 0 + refl_swrad_type(i) = balb*tcf_type(i)*(beam_swrad_c*coszen+(1-svf)*diff_swrad_c) + ENDIF + ENDDO + refl_swrad_c = sum(refl_swrad_type, mask = refl_swrad_type /= -1.0e36) + forc_swrad_c = beam_swrad_c+diff_swrad_c+refl_swrad_c + + ! But ensure that we don't depart too far from the atmospheric forcing value: + ! negative values of swrad are certainly bad, but small positive values might + ! also be bad. We can especially run into trouble due to the normalization: a + ! small swrad value in one column can lead to a big normalization factor, + ! leading to huge swrad values in other columns. + + forc_swrad_c = min(forc_swrad_c, & + forc_swrad_g * (1._r8 + shortwave_downscaling_limit)) + forc_swrad_c = max(forc_swrad_c, & + forc_swrad_g * (1._r8 - shortwave_downscaling_limit)) + ! for normalize + IF (forc_swrad_c==0.) forc_swrad_c = 0.0001 - END SUBROUTINE downscale_longwave_1c + END SUBROUTINE downscale_shortwave END MODULE MOD_ForcingDownscaling diff --git a/main/MOD_OrbCosazi.F90 b/main/MOD_OrbCosazi.F90 new file mode 100644 index 00000000..755d5b1e --- /dev/null +++ b/main/MOD_OrbCosazi.F90 @@ -0,0 +1,69 @@ +#include + +MODULE MOD_OrbCosazi + +!----------------------------------------------------------------------- + USE MOD_Precision + IMPLICIT NONE + SAVE + +! PUBLIC MEMBER FUNCTIONS: + PUBLIC :: orb_cosazi +!----------------------------------------------------------------------- + +CONTAINS + +!----------------------------------------------------------------------- + + FUNCTION orb_cosazi(calday, dlon, dlat, coszen) + +!----------------------------------------------------------------------- + USE MOD_Precision + IMPLICIT NONE + + REAL(r8), intent(in) :: calday !Julian cal day (1.xx to 365.xx) + REAL(r8), intent(in) :: dlat !Centered latitude (radians) + REAL(r8), intent(in) :: dlon !Centered longitude (radians) + REAL(r8), intent(in) :: coszen !cosine of sun zenith angle + REAL(r8) :: orb_cosazi !cosine of sun azimuth angle + + ! --- Local variables --- + REAL(r8) declin !Solar declination (radians) + REAL(r8) eccf !Earth-sun distance factor (ie. (1/r)**2) + REAL(r8) lambm !Lambda m, mean long of perihelion (rad) + REAL(r8) lmm !Intermediate argument involving lambm + REAL(r8) lamb !Lambda, the earths long of perihelion + REAL(r8) invrho !Inverse normalized sun/earth distance + REAL(r8) sinl !Sine of lmm + REAL(r8) pi !3.14159265358979323846... + REAL(r8), parameter :: & + dayspy=365.0, &!days per year + ve=80.5, &!Calday of vernal equinox assumes Jan 1 = calday 1 + eccen=1.672393084E-2, &!Eccentricity + obliqr=0.409214646, &!Earths obliquity in radians + lambm0=-3.2625366E-2, &!Mean long of perihelion at the vernal equinox (radians) + mvelpp=4.92251015 !moving vernal equinox longitude of + !perihelion plus pi (radians) + !------------------------------------------------------------------------------- + + pi = 4.*atan(1.) + lambm = lambm0 + (calday - ve)*2.*pi/dayspy + lmm = lambm - mvelpp + + sinl = sin(lmm) + lamb = lambm + eccen*(2.*sinl + eccen*(1.25*sin(2.*lmm) & + + eccen*((13.0/12.0)*sin(3.*lmm) - 0.25*sinl))) + invrho = (1. + eccen*cos(lamb - mvelpp)) / (1. - eccen*eccen) + + declin = asin(sin(obliqr)*sin(lamb)) + eccf = invrho*invrho + + orb_cosazi = (-1*cos(declin)*cos(calday*2.0*pi+dlon)- & + coszen*cos(dlat))/(sin(dlat)*sqrt(1-coszen*coszen)) + + IF (orb_cosazi<-1) orb_cosazi = -1 + IF (orb_cosazi>1) orb_cosazi = 1 + + END FUNCTION orb_cosazi + +END MODULE MOD_OrbCosazi diff --git a/main/MOD_OrbCoszen.F90 b/main/MOD_OrbCoszen.F90 index 7dc1093e..4ead0cc0 100644 --- a/main/MOD_OrbCoszen.F90 +++ b/main/MOD_OrbCoszen.F90 @@ -72,6 +72,9 @@ FUNCTION orb_coszen(calday,dlon,dlat) orb_coszen = sin(dlat)*sin(declin) & - cos(dlat)*cos(declin)*cos(calday*2.0*pi+dlon) + IF (orb_coszen<0) orb_coszen = 0 + IF (orb_coszen>1) orb_coszen = 1 + END FUNCTION orb_coszen END MODULE MOD_OrbCoszen diff --git a/main/MOD_Vars_1DForcing.F90 b/main/MOD_Vars_1DForcing.F90 index 9613a039..277095e6 100644 --- a/main/MOD_Vars_1DForcing.F90 +++ b/main/MOD_Vars_1DForcing.F90 @@ -30,6 +30,7 @@ MODULE MOD_Vars_1DForcing real(r8), allocatable :: forc_solsd (:) ! atm vis diffuse solar rad onto srf [W/m2] real(r8), allocatable :: forc_solld (:) ! atm nir diffuse solar rad onto srf [W/m2] real(r8), allocatable :: forc_frl (:) ! atmospheric infrared (longwave) radiation [W/m2] + real(r8), allocatable :: forc_swrad (:) ! atmospheric shortwave radiation [W/m2] real(r8), allocatable :: forc_hgt_u (:) ! observational height of wind [m] real(r8), allocatable :: forc_hgt_t (:) ! observational height of temperature [m] real(r8), allocatable :: forc_hgt_q (:) ! observational height of humidity [m] @@ -82,6 +83,7 @@ SUBROUTINE allocate_1D_Forcing allocate (forc_solsd (numpatch) ) ! atm vis diffuse solar rad onto srf [W/m2] allocate (forc_solld (numpatch) ) ! atm nir diffuse solar rad onto srf [W/m2] allocate (forc_frl (numpatch) ) ! atmospheric infrared (longwave) radiation [W/m2] + allocate (forc_swrad (numpatch) ) ! atmospheric shortwave radiation [W/m2] allocate (forc_hgt_u (numpatch) ) ! observational height of wind [m] allocate (forc_hgt_t (numpatch) ) ! observational height of temperature [m] allocate (forc_hgt_q (numpatch) ) ! observational height of humidity [m] @@ -133,6 +135,7 @@ SUBROUTINE deallocate_1D_Forcing () deallocate ( forc_solsd ) ! atm vis diffuse solar rad onto srf [W/m2] deallocate ( forc_solld ) ! atm nir diffuse solar rad onto srf [W/m2] deallocate ( forc_frl ) ! atmospheric infrared (longwave) radiation [W/m2] + deallocate ( forc_swrad ) ! atmospheric shortwave radiation [W/m2] deallocate ( forc_hgt_u ) ! observational height of wind [m] deallocate ( forc_hgt_t ) ! observational height of temperature [m] deallocate ( forc_hgt_q ) ! observational height of humidity [m] diff --git a/main/MOD_Vars_Global.F90 b/main/MOD_Vars_Global.F90 index 24a528be..d3a7bd8a 100644 --- a/main/MOD_Vars_Global.F90 +++ b/main/MOD_Vars_Global.F90 @@ -54,6 +54,11 @@ MODULE MOD_Vars_Global integer, parameter :: nl_roof = 10 integer, parameter :: nl_wall = 10 integer, parameter :: nvegwcs = 4 ! number of vegetation water potential nodes + + ! used for downscaling + integer, parameter :: num_type = 4 + integer, parameter :: num_zenith = 51 + integer, parameter :: num_azimuth = 36 ! bgc variables integer, parameter :: ndecomp_pools = 7 diff --git a/main/MOD_Vars_TimeInvariants.F90 b/main/MOD_Vars_TimeInvariants.F90 index f6f907c7..80e443fb 100644 --- a/main/MOD_Vars_TimeInvariants.F90 +++ b/main/MOD_Vars_TimeInvariants.F90 @@ -255,6 +255,14 @@ MODULE MOD_Vars_TimeInvariants real(r8) :: tcrit !critical temp. to determine rain or snow real(r8) :: wetwatmax !maximum wetland water (mm) + ! Used for downscaling + real(r8), allocatable :: svf_patches (:) ! sky view factor + real(r8), allocatable :: cur_patches (:) ! curvature + real(r8), allocatable :: sf_lut_patches (:,:,:) ! look up table of shadow factor of a patch + real(r8), allocatable :: asp_type_patches (:,:) ! topographic aspect of each character of one patch + real(r8), allocatable :: slp_type_patches (:,:) ! topographic slope of each character of one patch + real(r8), allocatable :: area_type_patches (:,:) ! area percentage of each character of one patch + ! PUBLIC MEMBER FUNCTIONS: PUBLIC :: allocate_TimeInvariants PUBLIC :: deallocate_TimeInvariants @@ -342,6 +350,14 @@ SUBROUTINE allocate_TimeInvariants () allocate (ibedrock (numpatch)) allocate (topoelv (numpatch)) allocate (topostd (numpatch)) + + ! Used for downscaling + allocate (svf_patches (numpatch)) + allocate (asp_type_patches (num_type,numpatch)) + allocate (slp_type_patches (num_type,numpatch)) + allocate (area_type_patches (num_type,numpatch)) + allocate (sf_lut_patches (num_azimuth,num_zenith,numpatch)) + allocate (cur_patches (numpatch)) ENDIF #if (defined LULC_IGBP_PFT || defined LULC_IGBP_PC) @@ -384,7 +400,7 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart) character(len=*), intent(in) :: dir_restart ! Local variables - character(len=256) :: file_restart, cyear + character(len=256) :: file_restart, cyear, lndname write(cyear,'(i4.4)') lc_year file_restart = trim(dir_restart) // '/const/' // trim(casename) //'_restart_const' // '_lc' // trim(cyear) // '.nc' @@ -468,6 +484,15 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart) CALL ncio_read_bcast_serial (file_restart, 'tcrit ', tcrit ) ! critical temp. to determine rain or snow CALL ncio_read_bcast_serial (file_restart, 'wetwatmax', wetwatmax) ! maximum wetland water (mm) + IF (DEF_USE_Forcing_Downscaling) THEN + CALL ncio_read_vector (file_restart, 'slp_type_patches', num_type, landpatch, slp_type_patches) + CALL ncio_read_vector (file_restart, 'svf_patches', landpatch, svf_patches) + CALL ncio_read_vector (file_restart, 'asp_type_patches', num_type, landpatch, asp_type_patches) + CALL ncio_read_vector (file_restart, 'area_type_patches', num_type, landpatch, area_type_patches) + CALL ncio_read_vector (file_restart, 'sf_lut_patches', num_azimuth, num_zenith, landpatch, sf_lut_patches) + CALL ncio_read_vector (file_restart, 'cur_patches', landpatch, cur_patches) + ENDIF + #if (defined LULC_IGBP_PFT || defined LULC_IGBP_PC) file_restart = trim(dir_restart) // '/const/' // trim(casename) //'_restart_pft_const' // '_lc' // trim(cyear) // '.nc' CALL READ_PFTimeInvariants (file_restart) @@ -546,6 +571,9 @@ SUBROUTINE WRITE_TimeInvariants (lc_year, casename, dir_restart) CALL ncio_define_dimension_vector (file_restart, landpatch, 'soilsnow', nl_soil-maxsnl) CALL ncio_define_dimension_vector (file_restart, landpatch, 'soil', nl_soil) CALL ncio_define_dimension_vector (file_restart, landpatch, 'lake', nl_lake) + CALL ncio_define_dimension_vector (file_restart, landpatch, 'type', num_type) + CALL ncio_define_dimension_vector (file_restart, landpatch, 'azi', num_azimuth) + CALL ncio_define_dimension_vector (file_restart, landpatch, 'zen', num_zenith) CALL ncio_write_vector (file_restart, 'patchclass', 'patch', landpatch, patchclass) ! CALL ncio_write_vector (file_restart, 'patchtype' , 'patch', landpatch, patchtype ) ! @@ -610,6 +638,15 @@ SUBROUTINE WRITE_TimeInvariants (lc_year, casename, dir_restart) CALL ncio_write_vector (file_restart, 'topoelv', 'patch', landpatch, topoelv) CALL ncio_write_vector (file_restart, 'topostd', 'patch', landpatch, topostd) + + IF (DEF_USE_Forcing_Downscaling) THEN + CALL ncio_write_vector (file_restart, 'svf_patches', 'patch', landpatch, svf_patches) + CALL ncio_write_vector (file_restart, 'cur_patches', 'patch', landpatch, cur_patches) + CALL ncio_write_vector (file_restart, 'slp_type_patches', 'type', num_type, 'patch', landpatch, slp_type_patches) + CALL ncio_write_vector (file_restart, 'asp_type_patches', 'type', num_type, 'patch', landpatch, asp_type_patches) + CALL ncio_write_vector (file_restart, 'area_type_patches', 'type', num_type, 'patch', landpatch, area_type_patches) + CALL ncio_write_vector (file_restart, 'sf_lut_patches', 'azi', num_azimuth, 'zen', num_zenith, 'patch', landpatch, sf_lut_patches) + ENDIF #ifdef USEMPI CALL mpi_barrier (p_comm_glb, p_err) @@ -662,6 +699,7 @@ END SUBROUTINE WRITE_TimeInvariants SUBROUTINE deallocate_TimeInvariants () + USE MOD_Namelist, only: DEF_USE_Forcing_Downscaling USE MOD_SPMD_Task USE MOD_LandPatch, only: numpatch @@ -736,6 +774,15 @@ SUBROUTINE deallocate_TimeInvariants () deallocate (topoelv ) deallocate (topostd ) + IF (DEF_USE_Forcing_Downscaling) THEN + deallocate(slp_type_patches ) + deallocate(svf_patches ) + deallocate(asp_type_patches ) + deallocate(area_type_patches ) + deallocate(sf_lut_patches ) + deallocate(cur_patches ) + ENDIF + ENDIF ENDIF @@ -758,7 +805,7 @@ SUBROUTINE check_TimeInvariants () USE MOD_SPMD_Task USE MOD_RangeCheck - USE MOD_Namelist, only : DEF_USE_BEDROCK + USE MOD_Namelist, only : DEF_USE_BEDROCK, DEF_USE_Forcing_Downscaling IMPLICIT NONE @@ -817,6 +864,15 @@ SUBROUTINE check_TimeInvariants () CALL check_vector_data ('topostd [m] ', topostd ) ! CALL check_vector_data ('BVIC [-] ', BVIC ) ! + IF (DEF_USE_Forcing_Downscaling) THEN + CALL check_vector_data ('slp_type_patches [rad] ' , slp_type_patches) ! slope + CALL check_vector_data ('svf_patches [-] ' , svf_patches) ! sky view factor + CALL check_vector_data ('asp_type_patches [rad] ' , asp_type_patches) ! aspect + CALL check_vector_data ('area_type_patches [-] ' , area_type_patches) ! area percent + CALL check_vector_data ('cur_patches [-]' , cur_patches ) + CALL check_vector_data ('sf_lut_patches [-] ' , sf_lut_patches) ! shadow mask + ENDIF + #ifdef USEMPI CALL mpi_barrier (p_comm_glb, p_err) #endif diff --git a/mkinidata/MOD_Initialize.F90 b/mkinidata/MOD_Initialize.F90 index d957747e..f013e01b 100644 --- a/mkinidata/MOD_Initialize.F90 +++ b/mkinidata/MOD_Initialize.F90 @@ -103,7 +103,7 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & character(len=256) :: fsoildat character(len=256) :: fsnowdat character(len=256) :: fcndat - character(len=256) :: ftopo + character(len=256) :: ftopo, lndname type(grid_type) :: gsoil type(grid_type) :: gsnow type(grid_type) :: gcn @@ -383,7 +383,32 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & ftopo = trim(dir_landdata)//'/topography/'//trim(cyear)//'/topostd_patches.nc' CALL ncio_read_vector (ftopo, 'topostd_patches', landpatch, topostd) #endif - +! ...................................... +! 1.5 Initialize topography factor data +! ...................................... +#ifdef SinglePoint + slp_type_patches(:,1) = SITE_slp_type + asp_type_patches(:,1) = SITE_asp_type + area_type_patches(:,1) = SITE_area_type + svf_patches(:) = SITE_svf + cur_patches(:) = SITE_cur + sf_lut_patches(:,:,1) = SITE_sf_lut +#else + IF (DEF_USE_Forcing_Downscaling) THEN + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/slp_type_patches.nc' ! slope + CALL ncio_read_vector (lndname, 'slp_type_patches', num_type, landpatch, slp_type_patches) + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/svf_patches.nc' ! sky view factor + CALL ncio_read_vector (lndname, 'svf_patches', landpatch, svf_patches) + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/asp_type_patches.nc' ! aspect + CALL ncio_read_vector (lndname, 'asp_type_patches', num_type, landpatch, asp_type_patches) + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/area_type_patches.nc' ! area percent + CALL ncio_read_vector (lndname, 'area_type_patches', num_type, landpatch, area_type_patches) + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/sf_lut_patches.nc' ! shadow mask + CALL ncio_read_vector (lndname, 'sf_lut_patches', num_azimuth, num_zenith, landpatch, sf_lut_patches) + lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/cur_patches.nc' ! curvature + CALL ncio_read_vector (lndname, 'cur_patches', landpatch, cur_patches) + ENDIF +#endif ! ................................ ! 1.6 Initialize TUNABLE constants ! ................................ diff --git a/mksrfdata/Aggregation_TopographyFactors.F90 b/mksrfdata/Aggregation_TopographyFactors.F90 new file mode 100644 index 00000000..6d53ab77 --- /dev/null +++ b/mksrfdata/Aggregation_TopographyFactors.F90 @@ -0,0 +1,423 @@ +#include + +SUBROUTINE Aggregation_TopographyFactors ( & + grid_topo_factor , dir_rawdata, dir_model_landdata, lc_year) + ! ---------------------------------------------------------------------- + ! Global topography-based factors data + ! + ! Created by Sisi Chen, Lu Li, 06/2024 + ! ---------------------------------------------------------------------- + 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 +#ifdef SrfdataDiag + USE MOD_Mesh, only : numelm + USE MOD_LandElm + USE MOD_SrfdataDiag +#endif + + IMPLICIT NONE + + ! arguments: + ! --------------------------------------------------------------- + INTEGER, intent(in) :: lc_year + TYPE(grid_type), intent(in) :: grid_topo_factor ! Grid structure for high resolution topography factors + CHARACTER(len=*), intent(in) :: dir_rawdata ! Direct of Rawdata + CHARACTER(len=*), intent(in) :: dir_model_landdata + + ! local variables: + ! --------------------------------------------------------------- + CHARACTER(len=256) :: landdir, lndname, cyear + CHARACTER(len=3) :: sdir + + TYPE (block_data_real8_2d) :: slp_grid ! slope + TYPE (block_data_real8_2d) :: asp_grid ! aspect + TYPE (block_data_real8_2d) :: svf_grid ! sky view factor + TYPE (block_data_real8_2d) :: cur_grid ! curvature + TYPE (block_data_real8_3d) :: tea_f_grid ! sine of terrain elevation angle at front of grid + TYPE (block_data_real8_3d) :: tea_b_grid ! sine of terrain elevation angle at back of grid + + ! patch + REAL(r8), allocatable :: svf_patches (:) + REAL(r8), allocatable :: cur_patches (:) + REAL(r8), allocatable :: tea_f_azi_patches (:,:) ! shape as (azimuth, patches) + REAL(r8), allocatable :: tea_b_azi_patches (:,:) + REAL(r8), allocatable :: sf_lut_patches (:,:,:) ! shape as (azimuth, zenith, patches) + + ! four defined types at all patches + REAL(r8), allocatable :: asp_type_patches (:,:) ! shape as (type, patches) + REAL(r8), allocatable :: slp_type_patches (:,:) + REAL(r8), allocatable :: area_type_patches (:,:) + + ! pixelsets + REAL(r8), allocatable :: slp_one (:) + REAL(r8), allocatable :: asp_one (:) + REAL(r8), allocatable :: svf_one (:) + REAL(r8), allocatable :: cur_one (:) + REAL(r8), allocatable :: area_one (:) + REAL(r8), allocatable :: sf_one (:) + REAL(r8), allocatable :: tea_f_azi_one (:,:) + REAL(r8), allocatable :: tea_b_azi_one (:,:) + REAL(r8), allocatable :: tea_f_one (:) + REAL(r8), allocatable :: tea_b_one (:) + LOGICAL , allocatable :: sf_mask_one (:) + LOGICAL , allocatable :: asp_mask_one (:) + LOGICAL , allocatable :: area_mask_one (:) + LOGICAL , allocatable :: slp_mask_one (:) + + ! pixelsets of four defined types at each patch + REAL(r8), allocatable :: asp_type_one (:,:) + REAL(r8), allocatable :: slp_type_one (:,:) + REAL(r8), allocatable :: area_type_one (:,:) + + REAL(r8) :: sum_area_one ! sum of pixel area of a patch + REAL(r8) :: zenith_angle(num_zenith) ! sine of sun zenith angle (divided by num_zenith part) + + ! local variables + INTEGER :: ipatch, i, ps, pe, type, a, z, count_pixels, num_pixels + +#ifdef SrfdataDiag + INTEGER :: typpatch(N_land_classification+1), ityp ! number of land classification +#endif + write(cyear,'(i4.4)') lc_year + landdir = trim(dir_model_landdata) // '/topography/' // trim(cyear) + +#ifdef USEMPI + CALL mpi_barrier (p_comm_glb, p_err) +#endif + IF (p_is_master) THEN + write(*,'(/, A)') 'Aggregate topography factor ...' + CALL system('mkdir -p ' // trim(adjustl(landdir))) + ENDIF +#ifdef USEMPI + CALL mpi_barrier (p_comm_glb, p_err) +#endif + +#ifdef SinglePoint + IF (USE_SITE_topography) THEN + RETURN + ENDIF +#endif + + ! ------------------------------------------------------------------- + ! read topography-based factor data + ! ------------------------------------------------------------------- + IF (p_is_io) THEN + lndname = trim(dir_rawdata)//"slope.nc" + CALL allocate_block_data (grid_topo_factor, slp_grid) + CALL ncio_read_block (lndname, 'slope', grid_topo_factor, slp_grid) + + lndname = trim(dir_rawdata)//"aspect.nc" + CALL allocate_block_data (grid_topo_factor, asp_grid) + CALL ncio_read_block (lndname, 'aspect', grid_topo_factor, asp_grid) + + lndname = trim(dir_rawdata)//"terrain_elev_angle_front.nc" + CALL allocate_block_data (grid_topo_factor, tea_f_grid, num_azimuth) + CALL ncio_read_block (lndname, 'tea_front', grid_topo_factor, num_azimuth, tea_f_grid) + + lndname = trim(dir_rawdata)//"terrain_elev_angle_back.nc" + CALL allocate_block_data (grid_topo_factor, tea_b_grid, num_azimuth) + CALL ncio_read_block (lndname, 'tea_back', grid_topo_factor, num_azimuth, tea_b_grid) + + lndname = trim(dir_rawdata)//"sky_view_factor.nc" + CALL allocate_block_data (grid_topo_factor, svf_grid) + CALL ncio_read_block (lndname, 'svf', grid_topo_factor, svf_grid) + + lndname = trim(dir_rawdata)//"curvature.nc" + CALL allocate_block_data (grid_topo_factor, cur_grid) + CALL ncio_read_block (lndname, 'curvature', grid_topo_factor, cur_grid) + + ! -------------------------------------------------------------------------- + ! aggregate the terrain factor data from the resolution of raw data to patch + ! -------------------------------------------------------------------------- +#ifdef USEMPI + ! mpi send + CALL aggregation_data_daemon ( grid_topo_factor, & + data_r8_2d_in1 = slp_grid, data_r8_2d_in2 = asp_grid, & + data_r8_2d_in3 = svf_grid, data_r8_2d_in4 = cur_grid, & + data_r8_3d_in1 = tea_f_grid, n1_r8_3d_in1 = num_azimuth, & + data_r8_3d_in2 = tea_b_grid, n1_r8_3d_in2 = num_azimuth) +#endif + ENDIF + + + IF (p_is_worker) THEN + ! allocate for output variables at patches + allocate (svf_patches (numpatch)) + allocate (cur_patches (numpatch)) + allocate (asp_type_patches (num_type, numpatch)) + allocate (slp_type_patches (num_type, numpatch)) + allocate (area_type_patches(num_type, numpatch)) + allocate (sf_lut_patches (num_azimuth, num_zenith, numpatch)) + ! generate sine of sun zenith angles at equal intervals + DO i = 1, num_zenith + zenith_angle(i) = pi/(2*num_zenith)*(i-1) + ENDDO + + ! aggregate loop + DO ipatch = 1, numpatch + CALL aggregation_request_data (landpatch, ipatch, grid_topo_factor, & + zip = USE_zip_for_aggregation, area = area_one, & + data_r8_2d_in1 = slp_grid, data_r8_2d_out1 = slp_one, & + data_r8_2d_in2 = asp_grid, data_r8_2d_out2 = asp_one, & + data_r8_2d_in3 = svf_grid, data_r8_2d_out3 = svf_one, & + data_r8_2d_in4 = cur_grid, data_r8_2d_out4 = cur_one, & + data_r8_3d_in1 = tea_f_grid, data_r8_3d_out1 = tea_f_azi_one, n1_r8_3d_in1 = num_azimuth, & + data_r8_3d_in2 = tea_b_grid, data_r8_3d_out2 = tea_b_azi_one, n1_r8_3d_in2 = num_azimuth) + + ! ------------------------------------------------------------------ + ! aggregate sky view factor, curvature at patches + ! ------------------------------------------------------------------ + IF (any(svf_one /= -9999.0)) THEN + svf_patches (ipatch) = & + sum(svf_one * area_one, mask = svf_one /= -9999.0) & + / sum(area_one, mask = svf_one /= -9999.0) + ELSE + svf_patches (ipatch) = -1.0e36 + ENDIF + + IF (any(cur_one /= -9999.0)) THEN + cur_patches (ipatch) = & + sum(cur_one * area_one, mask = cur_one /= -9999.0) & + / sum(area_one, mask = cur_one /= -9999.0) + ELSE + cur_patches (ipatch) = -1.0e36 + ENDIF + + ! ------------------------------------------------------------------------------ + ! aggregate look up table of shadow factor at patches + ! ------------------------------------------------------------------------------ + ! number of pixels + num_pixels = size(area_one) + + ! allocate pixel variables + allocate(tea_f_one (num_pixels)) + allocate(tea_b_one (num_pixels)) + allocate(sf_one (num_pixels)) + allocate(sf_mask_one (num_pixels)) + + ! sum of areas of one patch + sum_area_one = sum(area_one, mask = area_one>0) + + DO a = 1, num_azimuth + ! terrain elevation angle at each azimuth + tea_f_one(:) = tea_f_azi_one(a,:) + tea_b_one(:) = tea_b_azi_one(a,:) + + DO z = 1, num_zenith + ! count the pixels which are not missing value + count_pixels = 0 + + DO i = 1, num_pixels + IF ((isnan(tea_f_one(i))).or.(isnan(tea_b_one(i)))) THEN + sf_one(i) = 1 ! Not consider the effect of casting shadows + ELSE + IF (tea_f_one(i)>1) tea_f_one(i) = 1 + IF (tea_f_one(i)<-1) tea_f_one(i) = -1 + IF (tea_b_one(i)>1) tea_b_one(i) = 1 + IF (tea_b_one(i)<-1) tea_b_one(i) = -1 + tea_f_one(i) = asin(tea_f_one(i)) + tea_b_one(i) = asin(tea_b_one(i)) + + ! Compare the sun's altitude angle to the terrain's altitude angle + ! to determine the value of the shadow factor. + ! ----------------------------------------------------------------- + ! Sisi Chen, Lu Li, Yongjiu Dai, et al. Exploring Topography Downscaling + ! Methods for Hyper-Resolution Land Surface Modeling. + ! DOI: 10.22541/au.171403656.68476353/v1 + ! ----------------------------------------------------------------- + IF ((tea_b_one(i) /= -9999).and.(tea_f_one(i) /= -9999)) THEN + count_pixels = count_pixels+1 + + IF (pi*0.5 - zenith_angle(z) < tea_b_one(i)) THEN + sf_one(i) = 0 + ELSE IF (pi*0.5 - zenith_angle(z) > tea_f_one(i)) THEN + sf_one(i) = 1 + ELSE + IF (tea_f_one(i).eq.tea_b_one(i)) tea_f_one(i) = tea_b_one(i)+0.001 + sf_one(i) = (0.5*pi - zenith_angle(z) - tea_b_one(i))/(tea_f_one(i) - tea_b_one(i)) + ENDIF + + ENDIF + ENDIF + ENDDO + sf_mask_one(:) = sf_one(:) /= -9999 + sf_lut_patches(a,z,ipatch) = sum(sf_one(:), mask = sf_mask_one)/count_pixels + ENDDO + ENDDO + + ! deallocate + deallocate(tea_f_one) + deallocate(tea_b_one) + deallocate(sf_one) + deallocate(sf_mask_one) + + ! ----------------------------------------------------------------------------------------------- + ! aggregate slope and aspect at four defined types at patches + ! ----------------------------------------------------------------------------------------------- + ! allocate pixelsets variables + allocate(asp_type_one(1:num_type,1:num_pixels)) + allocate(slp_type_one(1:num_type,1:num_pixels)) + allocate(area_type_one(1:num_type,1:num_pixels)) + allocate(slp_mask_one(1:num_pixels)) + allocate(asp_mask_one(1:num_pixels)) + allocate(area_mask_one(1:num_pixels)) + + DO i = 1, num_type + asp_type_one(i,:) = -9999 + slp_type_one(i,:) = -9999 + area_type_one(i,:) = -9999 + ENDDO + + DO i = 1, num_pixels + ! Define the south slope, north slope, abrupt slope and gentle lope of target pixel + IF ((asp_one(i).ge.0 .and. asp_one(i).le.90) .or. (asp_one(i).ge.270 .and. asp_one(i).le.360).and.(slp_one(i).ge.15*pi/180)) THEN ! north abrupt slope + type = 1 + ELSE IF ((asp_one(i).ge.0 .and. asp_one(i).le.90) .or. (asp_one(i).ge.270 .and. asp_one(i).le.360).and.(slp_one(i)<15*pi/180)) THEN ! north gentle slope + type = 2 + ELSE IF ((asp_one(i).gt.90) .and. (asp_one(i).lt.270) .and. (slp_one(i).ge.15*pi/180)) THEN ! south abrupt slope + type = 3 + ELSE IF ((asp_one(i).gt.90) .and. (asp_one(i).lt.270) .and. (slp_one(i).lt.15*pi/180)) THEN ! south gentle slope + type = 4 + ELSE ! missing value=-9999 + cycle + END IF + + IF ((area_one(i)>0).and.(area_one(i)