Feature/namelist control

Makefile

@@ -26,7 +26,7 @@ $(info DEBUG setting: '$(DEBUG)')
 FLFLAGS :=
 
 # Source objects
-OBJS := canopy_const_mod.o canopy_parm_mod.o canopy_utils_mod.o canopy_waf_mod.o canopy_wind_mod.o canopy_driver.o
+OBJS := canopy_files_mod.o canopy_const_mod.o canopy_parm_mod.o canopy_utils_mod.o canopy_waf_mod.o canopy_wind_mod.o canopy_readnml.o canopy_driver.o
 
 # Program name
 PROGRAM := canopy

README.md

@@ -5,8 +5,9 @@ Author(s):
 
 Patrick Campbell, Zachary Moon, and Wei-Ting Hung
 
-Compile and run canopy model:
+Compile, edit namelist, and run canopy model:
 - `make`
+- `namelist.canopy`
 - `./canopy`
 
 Canopy is parameterized by foliage distribution shape functions and parameters for different vegetation types (Katul et al., 2004; Massman et al., 2017)
@@ -22,15 +23,13 @@ Current Canopy-App components:
 
 2.  In-Canopy radiative transfer (i.e., photolysis rate adjusment) for air quality applications (in progress...)  (Makar et al., 2017)
 
-- `canopy_rad_mod.F90`
+- `canopy_prad_mod.F90`
 
 3.  In-Canopy vertical diffusion (i.e., eddy diffusivity adjustment) for air quality applications (in progress...) (Makar et al., 2017)
 
 - `canopy_vdiff_mod.F90`
 
 
-
-
 Citations:
 
 Katul, G.G., Mahrt, L., Poggi, D., and Sanz, C. (2004). One- and two-equation models for canopy turbulence. Boundary-Layer Meteorol. 113: 81–109. doi:10.1023/B:BOUN.0000037333.48760.e5.

canopy_driver.F90

@@ -16,18 +16,26 @@ program canopy_driver
 !-------------------------------------------------------------
        use canopy_utils_mod !utilities for canopy models
        use canopy_parm_mod  !main canopy parameters
+       use canopy_files_mod !main canopy input files
        use canopy_wind_mod  !main canopy wind model
        use canopy_waf_mod   !main Wind Adjustment Factor (WAF) model
 
       implicit none
         INTEGER, PARAMETER :: rk = SELECTED_REAL_KIND(15, 307)
-! !....this block defines geographic domain of inputs (SE CONUS, 0.1 degree resolution)
-        integer, parameter        ::    nlat=101        !length of x coordinate
-        integer, parameter        ::    nlon=201        !length of y coordinate
-!        integer, parameter        ::    nlat=1           !length of x coordinate --set for 1D check
-!        integer, parameter        ::    nlon=1           !length of y coordinate --set for 1D check
-        integer, parameter        ::    canlays=100     !Number of total above and below canopy layers
-        real(rk),    parameter    ::    href=10.0       !Reference Height above canopy @ 10 m  (m)
+
+! !....this block defines geographic domain of inputs (read from user namelist)
+        integer        ::    nlat        !length of x coordinate
+        integer        ::    nlon        !length of y coordinate
+        integer        ::    canlays     !Number of total above and below canopy layers
+        real           ::    canres      !Real value of canopy vertical resolution (m)
+        real           ::    href        !Reference Height above canopy @ 10 m  (m)
+        logical        ::    ifcanwind   !logical canopy wind/WAF option (default = .FALSE.)
+
+! !....this block gives assumed constant parameters for in-canopy conditions (read from user namelist)
+        real           ::    z0ghcm   ! ratio of ground roughness length to canopy top height
+                                                        !  (default case: Approx. currently does not account for understory
+                                                        !   variability)
+        real           ::   lamdars     ! Influence function associated with roughness sublayer (nondimensional)
 
 ! !....this block gives input canopy height and above reference conditions that should be passed (assume winds at href)
         real(rk)                  ::    lat             !latitude  (degrees)
@@ -43,13 +51,6 @@ program canopy_driver
         real(rk)                  ::    z0ref           !Input total/surface roughness length
         real(rk)                  ::    molref          !Input Monin-Obukhov Length
 
-! !....this block gives assumed constant parameters for in-canopy conditions (Unavoidable Uncertainties!!!)
-!      Need to move to module canopy_parm_mod for all canopy/vegetation types
-        real(rk),    parameter    ::    z0ghcm=0.0025   ! ratio of ground roughness length to canopy top height
-                                                        !  (default case: Approx. currently does not account for understory
-                                                        !   variability)
-        real(rk),    parameter    ::   lamdars=1.25     ! Influence function associated with roughness sublayer (nondimensional)
-
 ! !....this block gives vegetion-type canopy dependent  parameters based on Katul et al. (2004)
         integer     ::    firetype         !1 = Above Canopy Fire; 0 = Below Canopy Fire
         real(rk)    ::    flameh         !Flame Height (m) 
@@ -61,31 +62,29 @@ program canopy_driver
 
 !Local variables        
         integer i,i0,loc
-        real(rk) :: zkcm       ( canlays )  ! in-canopy heights (m)
-        real(rk) :: resz       ( canlays )  ! canopy height resolution (m)
-        real(rk) :: ztothc     ( canlays )  ! z/h
-        real(rk) :: fainc      ( canlays )  ! incremental foliage shape function
-        real(rk) :: fafracz    ( canlays )  ! incremental fractional foliage shape function
-        real(rk) :: fafraczInt ( canlays )  ! integral of incremental fractional foliage shape function
-        real(rk) :: canBOT     ( canlays )  ! Canopy bottom wind reduction factors (nondimensional)
-        real(rk) :: canTOP     ( canlays )  ! Canopy top wind reduction factors (nondimensional)
+        real(rk), allocatable :: zkcm       ( : )  ! in-canopy heights (m)
+        real(rk), allocatable :: ztothc     ( : )  ! z/h
+        real(rk), allocatable :: fainc      ( : )  ! incremental foliage shape function
+        real(rk), allocatable :: fafracz    ( : )  ! incremental fractional foliage shape function
+        real(rk), allocatable :: fafraczInt ( : )  ! integral of incremental fractional foliage shape function
+        real(rk), allocatable :: canBOT     ( : )  ! Canopy bottom wind reduction factors (nondimensional)
+        real(rk), allocatable :: canTOP     ( : )  ! Canopy top wind reduction factors (nondimensional)
         real(rk) :: fatot                   ! integral of total fractional foliage shape function
-        real(rk) :: canWIND    ( canlays, nlat*nlon )  ! final mean canopy wind speeds (m/s)        
+        real(rk), allocatable :: canWIND    ( :, : )  ! final mean canopy wind speeds (m/s)        
 
         integer  ::    cansublays           ! number of sub-canopy layers
         integer  ::    canmidpoint          ! indice of the sub-canopy midpoint
         integer  ::    flamelays            ! number of flame layers
         integer  ::    midflamepoint        ! indice of the mid-flame point
-        real(rk) ::    waf     (nlat*nlon)  ! Calculated Wind Adjustment Factor
+        real(rk), allocatable ::    waf     (:)  ! Calculated Wind Adjustment Factor
 
 !     Test Generic 1D canopy data that should represent virtualized canopy
       TYPE :: profile_type
            integer     :: canlay       !profile layer for model
            real(rk)    :: zk           !profile heights for model (m)
-           real(rk)    :: dzk          !profile height increments (m)
       end TYPE profile_type
 
-      type(profile_type) :: profile( canlays )
+      type(profile_type), allocatable :: profile( : )
 
 !     Test Generic 2D met/sfc input variables that should be passed to canopy calculations
       TYPE :: variable_type
@@ -103,28 +102,55 @@ program canopy_driver
            real(rk)    :: mol          !Monin-Obukhov length
       end TYPE variable_type
 
-      type(variable_type) :: variables(nlat*nlon)      
+      type(variable_type), allocatable :: variables( : )      
+
+!-------------------------------------------------------------------------------
+! Read user options from namelist.
+!-------------------------------------------------------------------------------
+
+      call  canopy_readnml(nlat,nlon,canlays,canres,href,z0ghcm,lamdars, &
+                           ifcanwind)
+      if (ifcanwind) then
+         write(*,*)  'Canopy wind/WAF option selected'
+      else
+         write(*,*)  'No option(s) selected'
+      end if
+
+!-------------------------------------------------------------------------------
+! Allocate necessary variables.
+!-------------------------------------------------------------------------------
+
+      if(.not.allocated(zkcm)) allocate(zkcm(canlays))
+      if(.not.allocated(ztothc)) allocate(ztothc(canlays))
+      if(.not.allocated(fainc)) allocate(fainc(canlays))
+      if(.not.allocated(fafracz)) allocate(fafracz(canlays))
+      if(.not.allocated(fafraczInt)) allocate(fafraczInt(canlays))
+      if(.not.allocated(canBOT)) allocate(canBOT(canlays))
+      if(.not.allocated(canTOP)) allocate(canTOP(canlays))
+      if(.not.allocated(canWIND)) allocate(canWIND(canlays,nlat*nlon))
+      if(.not.allocated(waf)) allocate(waf(nlat*nlon))
+      if(.not.allocated(profile)) allocate(profile(canlays))
+      if(.not.allocated(variables)) allocate(variables(nlat*nlon))
 
 
 ! ... read canopy profile data that should be passed
-      open(9,  file='input_profile.txt',  status='old')
+      open(9,  file=file_prof(1),  status='old')
       i0 = 0
       read(9,*,iostat=i0)  ! skip headline
       do i=1, canlays
         read(9, *) profile(i)
       end do
       close(9)
+
 ! ... read met/sfc input variables
-      open(8,  file='input_variables.txt',  status='old')
-!      open(8,  file='input_variables_1D.txt',  status='old') ! --set for 1D check
+      open(8,  file=file_vars(1),  status='old')
       i0 = 0
       read(8,*,iostat=i0)  ! skip headline
       do loc=1, nlat*nlon
         read(8, *) variables(loc)
       end do
       close(8)
 
-
       do loc=1, nlat*nlon
         lat      = variables(loc)%lat
         lon      = variables(loc)%lon
@@ -139,7 +165,7 @@ program canopy_driver
         z0ref    = variables(loc)%z0
         molref   = variables(loc)%mol
 
-! ... call canopy parameters to get canopy and shape distribution parameters         
+! ... call canopy parameters to get canopy, fire info, and shape distribution parameters         
 
         call canopy_parm(lat, lon, vtyperef, hcm, ffracref, lairef, &
                          firetype, flameh, cdrag, &
@@ -148,13 +174,11 @@ program canopy_driver
 ! ... initialize canopy model and integrate to get fractional plant area distribution functions
         zkcm   = profile%zk
         ztothc = zkcm/hcm
-        resz   = profile%dzk
-        cansublays  = floor(hcm/resz(canlays))
+        cansublays  = floor(hcm/canres)
         canmidpoint = cansublays/2
-        flamelays    = floor(flameh/resz(canlays))
+        flamelays    = floor(flameh/canres)
         midflamepoint   = flamelays/2
 
-
 ! ... calculate canopy/foliage distribution shape profile - bottom up total in-canopy and fraction at z
         fainc = 0.0_rk  ! initialize
         fatot = 0.0_rk  ! initialize
@@ -166,42 +190,50 @@ program canopy_driver
           end if
         end do
         fatot = IntegrateTrapezoid(ztothc,fainc)
-
-! ... calculate plant distribution function and in-canopy wind speeds at z
+        
+! ... calculate plant distribution function
         fafracz    = 0.0_rk  ! initialize
         fafraczInt = 0.0_rk  ! initialize
         do i=1, canlays
           fafracz(i) = fainc(i)/fatot
           fafraczInt(i) = IntegrateTrapezoid(ztothc(1:i),fafracz(1:i))
-          call canopy_wind(hcm, zkcm(i), fafraczInt(i), ubzref, &
-                    z0ghcm, cdrag, pai, canBOT(i), canTOP(i), canWIND(i, loc))
         end do
 
+        if (ifcanwind) then
+! ... calculate in-canopy wind speeds at z
+           do i=1, canlays
+             call canopy_wind(hcm, zkcm(i), fafraczInt(i), ubzref, &
+                    z0ghcm, cdrag, pai, canBOT(i), canTOP(i), canWIND(i, loc))
+           end do
+
 ! ... calculate wind adjustment factor dependent on canopy winds and fire type
 
-        call canopy_waf(hcm, ztothc(1:cansublays), fafraczInt(1:cansublays), fafraczInt(1), & 
+           call canopy_waf(hcm, ztothc(1:cansublays), fafraczInt(1:cansublays), fafraczInt(1), & 
                         ubzref, z0ghcm, lamdars, cdrag, pai, href, flameh, firetype, & 
                         canBOT(midflamepoint), canTOP(midflamepoint), waf(loc))
-
-!            write(*,*)  'Wind Adjustment Factor:', waf(loc)
-      end do
-
-! ... save as text files for testing
-      open(10, file='output_canopy_wind.txt')
-      write(10, '(a30, f6.1, a2)') 'Reference height, h: ', href, 'm'
-      write(10, '(a30, i6)') 'Number of in-canopy layers: ', canlays
-      write(10, '(a8, a9, a12, a15)') 'Lat', 'Lon', 'Height (m)', 'WS (m/s)'
-      do loc=1, nlat*nlon
-        do i=1, canlays
-          write(10, '(f8.2, f9.2, f12.2, es15.7)')  variables(loc)%lat, variables(loc)%lon, profile(i)%zk, canWIND(i, loc)
-        end do
-      end do
+        end if
 
-      open(11, file='output_waf.txt')
-      write(11, '(a30, f6.1)') 'Reference height, h: ', href, 'm'
-      write(11, '(a8, a9, a19, a11)') 'Lat', 'Lon', 'Canopy height (m)', 'WAF'
-      do loc=1, nlat*nlon
-        write(11, '(f8.2, f9.2, f19.2, f11.7)')  variables(loc)%lat, variables(loc)%lon, variables(loc)%fh, waf(loc)
       end do
+
+      if (ifcanwind) then
+       write(*,*)  'Writing canopy wind/WAF output'
+! ... save as text files for testing
+       open(10, file='output_canopy_wind.txt')
+       write(10, '(a30, f6.1, a2)') 'Reference height, h: ', href, 'm'
+       write(10, '(a30, i6)') 'Number of in-canopy layers: ', canlays
+       write(10, '(a8, a9, a12, a15)') 'Lat', 'Lon', 'Height (m)', 'WS (m/s)'
+       do loc=1, nlat*nlon
+         do i=1, canlays
+           write(10, '(f8.2, f9.2, f12.2, es15.7)')  variables(loc)%lat, variables(loc)%lon, profile(i)%zk, canWIND(i, loc)
+         end do
+       end do
+
+       open(11, file='output_waf.txt')
+       write(11, '(a30, f6.1)') 'Reference height, h: ', href, 'm'
+       write(11, '(a8, a9, a19, a11)') 'Lat', 'Lon', 'Canopy height (m)', 'WAF'
+       do loc=1, nlat*nlon
+         write(11, '(f8.2, f9.2, f19.2, f11.7)')  variables(loc)%lat, variables(loc)%lon, variables(loc)%fh, waf(loc)
+       end do
+      end if
 
     end program canopy_driver

canopy_files_mod.F90

@@ -0,0 +1,17 @@
+
+MODULE canopy_files_mod
+
+!-------------------------------------------------------------------------------
+! Name:     Files
+! Purpose:  Contains FORTRAN units and file names.
+! Revised:  15 Jul 2022  Original version.  (P. C. Campbell)
+!-------------------------------------------------------------------------------
+
+  IMPLICIT NONE
+
+  INTEGER,            PARAMETER     :: max_mm     = 1
+  INTEGER,            PARAMETER     :: iutnml     =  8
+  CHARACTER(LEN=256)                :: file_prof ( max_mm ), file_vars ( max_mm )
+  CHARACTER(LEN=256), PARAMETER     :: file_nml   = 'namelist.canopy'
+
+END MODULE canopy_files_mod