diff --git a/src/canopy_alloc.F90 b/src/canopy_alloc.F90 index 7fb91260..7e231ffd 100644 --- a/src/canopy_alloc.F90 +++ b/src/canopy_alloc.F90 @@ -24,6 +24,8 @@ SUBROUTINE canopy_alloc if(.not.allocated(variables_3d)) allocate(variables_3d(nlon,nlat,var3d_set)) if(.not.allocated(variables_1d)) allocate(variables_1d(var3d_set)) if(.not.allocated(variables_can)) allocate(variables_can(nlat*nlon)) + if(.not.allocated(pavdref)) allocate(pavdref(var3d_set)) + if(.not.allocated(levref)) allocate(levref(var3d_set)) if(.not.allocated(pavd_arr)) allocate(pavd_arr(var3d_set)) if(.not.allocated(lev_arr)) allocate(lev_arr(var3d_set)) end if diff --git a/src/canopy_calcs.F90 b/src/canopy_calcs.F90 index 50a0b8c5..4a75c9f7 100644 --- a/src/canopy_calcs.F90 +++ b/src/canopy_calcs.F90 @@ -35,18 +35,23 @@ SUBROUTINE canopy_calcs(nn) REAL(rk), save :: currentlai ! Current LAI [cm2/cm2] (saved from one timestep to the next) REAL(rk) :: tsteplai !Number of days between the past and current LAI !REAL(rk) :: tabovecanopy ! Above Canopy Temp (assigned = tmp2mref ), done in canopy_bioemi_mod.F90 - + REAL(rk) :: lat2d(nlon,nlat), lon2d(nlon,nlat), lat1d(nlon*nlat), lon1d(nlon*nlat) write(*,*) 'Calculating Canopy Parameters' write(*,*) '-------------------------------' + lat2d = variables_2d%lat + lon2d = variables_2d%lon + lat1d = variables%lat + lon1d = variables%lon + if (infmt_opt .eq. 0) then !Main input format is 2D NetCDF and output will be 2D NetCDf if (ifcanwind .or. ifcanwaf) then !only calculate if canopy wind or WAF option - call canopy_calcdx_2d(dx_opt, dx_set, nlat, nlon, variables_2d%lat, & - variables_2d%lon, dx_2d) + call canopy_calcdx_2d(dx_opt, dx_set, nlat, nlon, lat2d, & + lon2d, dx_2d) end if if (href_opt .eq. 0 ) then !setting entire array = href_set value from user NL @@ -151,11 +156,13 @@ SUBROUTINE canopy_calcs(nn) call canopy_foliage(modlays, zhc, zcanmax, sigmau, sigma1, & fafraczInt) else + pavdref = variables_3d(i,j,:)%pavd + levref = variables_1d%lev ! ... derive canopy/foliage distribution shape profile from interpolated GEDI PAVD profile - bottom up total in-canopy and fraction at z if (variables_2d(i,j)%lat .gt. (-1.0_rk*pavd_set) .and. & variables_2d(i,j)%lat .lt. pavd_set) then !use GEDI PAVD call canopy_pavd2fafrac(zcanmax, sigmau, sigma1, hcmref, zhc, & - variables_3d(i,j,:)%pavd, variables_1d%lev, fafraczInt) + pavdref, levref, fafraczInt) !check if there is observed canopy height but no PAVD profile if (hcmref .gt. 0.0 .and. maxval(fafraczInt) .le. 0.0) then !revert to prescribed shape profile call canopy_foliage(modlays, zhc, zcanmax, sigmau, sigma1, & @@ -423,8 +430,8 @@ SUBROUTINE canopy_calcs(nn) else if (infmt_opt .eq. 1) then !Main input format is 1D/2D text and output will be 1D/2D text if (ifcanwind .or. ifcanwaf) then !only calculate if canopy wind or WAF option - call canopy_calcdx(dx_opt, dx_set, nlat, nlon, variables%lat, & - variables%lon, dx) + call canopy_calcdx(dx_opt, dx_set, nlat, nlon, lat1d, & + lon1d, dx) end if if (href_opt .eq. 0 ) then !setting entire array = href_set value from user NL diff --git a/src/canopy_canmet_mod.F90 b/src/canopy_canmet_mod.F90 index 72977914..975049dd 100644 --- a/src/canopy_canmet_mod.F90 +++ b/src/canopy_canmet_mod.F90 @@ -123,7 +123,7 @@ MODULE canopy_canmet_mod ! real(rk) :: pavd01ref, pavd02ref, pavd03ref, pavd04ref, pavd05ref, & !Input canopy PAVD profile ! pavd06ref, pavd07ref, pavd08ref, pavd09ref, pavd10ref, & ! pavd11ref, pavd12ref, pavd13ref, pavd14ref - real(rk), allocatable :: pavd_arr ( : ) - real(rk), allocatable :: lev_arr ( : ) + real(rk), allocatable :: pavdref ( : ), pavd_arr ( : ) !plant area volume density (m2/m3) + real(rk), allocatable :: levref ( : ), lev_arr ( : ) !reference vertical levels with 3d input data END MODULE canopy_canmet_mod diff --git a/src/canopy_dealloc.F90 b/src/canopy_dealloc.F90 index 6d49afec..e28d9d2f 100644 --- a/src/canopy_dealloc.F90 +++ b/src/canopy_dealloc.F90 @@ -22,7 +22,11 @@ SUBROUTINE canopy_dealloc if(allocated(variables_1d)) deallocate(variables_1d) if(allocated(variables_2d)) deallocate(variables_2d) if(allocated(variables_3d)) deallocate(variables_3d) - if(allocated(variables_can)) deallocate(variables_can) + if(allocated(variables_can)) deallocate(variables_can) + if(allocated(pavdref)) deallocate(pavdref) + if(allocated(levref)) deallocate(levref) + if(allocated(pavd_arr)) deallocate(pavd_arr) + if(allocated(lev_arr)) deallocate(lev_arr) !------------------------------------------------------------------------------- ! Dellocate arrays for Canopy Distribution diff --git a/src/canopy_ncf_io_mod.F90 b/src/canopy_ncf_io_mod.F90 index 86c70b30..c11f5649 100644 --- a/src/canopy_ncf_io_mod.F90 +++ b/src/canopy_ncf_io_mod.F90 @@ -2066,6 +2066,7 @@ SUBROUTINE canopy_read_ncf(infile) ! Revised: 19 Dec 2022 Original version. (P.C. Campbell) !------------------------------------------------------------------------------- USE canopy_canopts_mod !main canopy option descriptions + use canopy_const_mod, ONLY: rk USE canopy_coord_mod USE canopy_canmet_mod USE netcdf @@ -2078,6 +2079,12 @@ SUBROUTINE canopy_read_ncf(infile) CHARACTER(LEN=32) :: date_init INTEGER :: cdfid,rcode,varid,it double precision :: rdtime + REAL(rk) :: variables_3d_real(nlon,nlat,var3d_set) + REAL(rk) :: variables_2d_real(nlon,nlat) + INTEGER :: variables_2d_int(nlon,nlat) + REAL(rk) :: variables_1d_lev_real(var3d_set) +! REAL(rk) :: variables_real(nlon*nlat) +! INTEGER :: variables_int(nlon*nlat) !------------------------------------------------------------------------------- ! Error, warning, and informational messages. @@ -2169,229 +2176,254 @@ SUBROUTINE canopy_read_ncf(infile) if (infmt_opt .eq. 0) then !Input format is 2D - CALL get_var_2d_real_cdf (cdfid, 'lat', variables_2d%lat, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'lat', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'lon', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%lat = variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%lat=reshape(variables_2d%lat,[size(variables_2d%lat)]) - CALL get_var_2d_real_cdf (cdfid, 'lon', variables_2d%lon, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'lon', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'lon', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%lon=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%lon=reshape(variables_2d%lon,[size(variables_2d%lon)]) !Canopy input met/sfc variables !Clumping index - CALL get_var_2d_real_cdf (cdfid, 'clu', variables_2d%clu, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'clu', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'clu', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%clu=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%clu=reshape(variables_2d%clu,[size(variables_2d%clu)]) !Cosine of solar zenith angle - CALL get_var_2d_real_cdf (cdfid, 'csz', variables_2d%csz, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'csz', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'csz', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%csz=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%csz=reshape(variables_2d%csz,[size(variables_2d%csz)]) !Forest Fraction - CALL get_var_2d_real_cdf (cdfid, 'ffrac', variables_2d%ffrac, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'ffrac', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'ffrac', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%ffrac=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%ffrac=reshape(variables_2d%ffrac,[size(variables_2d%ffrac)]) !Forest canopy height - CALL get_var_2d_real_cdf (cdfid, 'fh', variables_2d%fh, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'fh', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'fh', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%fh=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%fh=reshape(variables_2d%fh,[size(variables_2d%fh)]) !Fire Radiative Power - CALL get_var_2d_real_cdf (cdfid, 'frp', variables_2d%frp, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'frp', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'frp', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%frp=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%frp=reshape(variables_2d%frp,[size(variables_2d%frp)]) !Reference height above canopy - CALL get_var_2d_real_cdf (cdfid, 'href', variables_2d%href, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'href', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'href', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%href=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%href=reshape(variables_2d%href,[size(variables_2d%href)]) !Leaf Area Index - CALL get_var_2d_real_cdf (cdfid, 'lai', variables_2d%lai, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'lai', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'lai', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%lai=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%lai=reshape(variables_2d%lai,[size(variables_2d%lai)]) !Monin-Obukhov Length - CALL get_var_2d_real_cdf (cdfid, 'mol', variables_2d%mol, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'mol', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'mol', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%mol=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%mol=reshape(variables_2d%mol,[size(variables_2d%mol)]) !Friction velocity - CALL get_var_2d_real_cdf (cdfid, 'fricv', variables_2d%fricv, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'fricv', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'fricv', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%fricv=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%fricv=reshape(variables_2d%fricv,[size(variables_2d%fricv)]) !Reference U Wind Speed (at HREF) - CALL get_var_2d_real_cdf (cdfid, 'ugrd10m', variables_2d%ugrd10m, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'ugrd10m', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'ugrd10m', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%ugrd10m=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%ugrd10m=reshape(variables_2d%ugrd10m,[size(variables_2d%ugrd10m)]) !Reference V Wind Speed (at HREF) - CALL get_var_2d_real_cdf (cdfid, 'vgrd10m', variables_2d%vgrd10m, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'vgrd10m', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'vgrd10m', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%vgrd10m=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%vgrd10m=reshape(variables_2d%vgrd10m,[size(variables_2d%vgrd10m)]) !Surface (veg+soil) Roughness Length - CALL get_var_2d_real_cdf (cdfid, 'sfcr', variables_2d%sfcr, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'sfcr', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'sfcr', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%sfcr=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%sfcr=reshape(variables_2d%sfcr,[size(variables_2d%sfcr)]) !Vegetation Type - CALL get_var_2d_int_cdf (cdfid, 'vtype', variables_2d%vtype, it, rcode) + CALL get_var_2d_int_cdf (cdfid, 'vtype', variables_2d_int, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'vtype', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%vtype=variables_2d_int !Also reshape to 1D array for 1D calculation and output ! variables%vtype=reshape(variables_2d%vtype,[size(variables_2d%vtype)]) !Soil Type - CALL get_var_2d_int_cdf (cdfid, 'sotyp', variables_2d%sotyp, it, rcode) + CALL get_var_2d_int_cdf (cdfid, 'sotyp', variables_2d_int, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'sotyp', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%sotyp=variables_2d_int !Also reshape to 1D array for 1D calculation and output ! variables%sotyp=reshape(variables_2d%sotyp,[size(variables_2d%sotyp)]) !Surface pressure - CALL get_var_2d_real_cdf (cdfid, 'pressfc', variables_2d%pressfc, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'pressfc', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'pressfc', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%pressfc=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%pressfc=reshape(variables_2d%pressfc,[size(variables_2d%pressfc)]) !instantaneous surface downward shortwave flux - CALL get_var_2d_real_cdf (cdfid, 'dswrf', variables_2d%dswrf, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'dswrf', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'dswrf', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%dswrf=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%dswrf=reshape(variables_2d%dswrf,[size(variables_2d%dswrf)]) !instantaneous surface sensible heat net flux - CALL get_var_2d_real_cdf (cdfid, 'shtfl', variables_2d%shtfl, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'shtfl', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'shtfl', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%shtfl=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%shtfl=reshape(variables_2d%shtfl,[size(variables_2d%shtfl)]) !Surface temperature - CALL get_var_2d_real_cdf (cdfid, 'tmpsfc', variables_2d%tmpsfc, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'tmpsfc', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'tmpsfc', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%tmpsfc=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%tmpsfc=reshape(variables_2d%tmpsfc,[size(variables_2d%tmpsfc)]) !2-meter temperature - CALL get_var_2d_real_cdf (cdfid, 'tmp2m', variables_2d%tmp2m, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'tmp2m', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'tmp2m', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%tmp2m=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%tmp2m=reshape(variables_2d%tmp2m,[size(variables_2d%tmp2m)]) !2-meter specific humidity - CALL get_var_2d_real_cdf (cdfid, 'spfh2m', variables_2d%spfh2m, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'spfh2m', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'spfh2m', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%spfh2m=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%spfh2m=reshape(variables_2d%spfh2m,[size(variables_2d%spfh2m)]) !Height of planetary boundary layer - CALL get_var_2d_real_cdf (cdfid, 'hpbl', variables_2d%hpbl, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'hpbl', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'hpbl', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%hpbl=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%hpbl=reshape(variables_2d%hpbl,[size(variables_2d%hpbl)]) !Mass precipitation rate - CALL get_var_2d_real_cdf (cdfid, 'prate_ave', variables_2d%prate_ave, it, rcode) + CALL get_var_2d_real_cdf (cdfid, 'prate_ave', variables_2d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'prate_ave', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_2d%prate_ave=variables_2d_real !Also reshape to 1D array for 1D calculation and output ! variables%prate_ave=reshape(variables_2d%prate_ave,[size(variables_2d%prate_ave)]) !3D Input Level Profile if (var3d_opt .eq. 1) then - CALL get_var_1d_real_cdf (cdfid, 'lev', variables_1d%lev, it, rcode) + CALL get_var_1d_real_cdf (cdfid, 'lev', variables_1d_lev_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'lev', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_1d%lev=variables_1d_lev_real end if !GEDI 3D PAVD Profile if (pavd_opt .eq. 1 ) then @@ -2400,12 +2432,13 @@ SUBROUTINE canopy_read_ncf(infile) call exit(2) else !read PAVD !Forest Plant Area Volume Density (PAVD) Profile - CALL get_var_3d_real_cdf (cdfid, 'pavd', variables_3d%pavd, it, rcode) + CALL get_var_3d_real_cdf (cdfid, 'pavd', variables_3d_real, it, rcode) IF ( rcode /= nf90_noerr ) THEN WRITE (*,f9410) TRIM(pname), 'pavd', & TRIM(nf90_strerror(rcode)) CALL exit(2) ENDIF + variables_3d%pavd=variables_3d_real end if end if