Added LAD netcdf output and fixed up some things.

noaa-oar-arl · Feb 14, 2024 · 7122b6f · 7122b6f
1 parent ead1686
commit 7122b6f
Show file tree

Hide file tree

Showing 6 changed files with 53 additions and 5 deletions.
diff --git a/src/canopy_alloc.F90 b/src/canopy_alloc.F90
@@ -43,7 +43,8 @@ SUBROUTINE canopy_alloc
     if(.not.allocated(ppfd_sun))    allocate(ppfd_sun(modlays))
     if(.not.allocated(ppfd_shade))  allocate(ppfd_shade(modlays))
     if(.not.allocated(ppfd_ave))    allocate(ppfd_ave(modlays))
-    if(.not.allocated(lad))           allocate(lad(nlat*nlon,modlays))
+    if(.not.allocated(lad))         allocate(lad(nlat*nlon,modlays))
+    if(.not.allocated(lad_3d))      allocate(lad_3d(nlon,nlat,modlays))
 
 !-------------------------------------------------------------------------------
 ! Allocate arrays for Canopy Wind Outputs

diff --git a/src/canopy_calcs.F90 b/src/canopy_calcs.F90
@@ -174,6 +174,19 @@ SUBROUTINE canopy_calcs(nn)
                                 end if
                             end if
 
+! ... calculate leaf area density profile from foliage shape function for output (m2/m3)
+                            do k=1, modlays
+                                if (zk(k) .gt. 0.0 .and. zk(k) .le. hcmref) then ! above ground level and at/below canopy top
+                                    if (k .lt. modlays)  then
+                                        lad_3d(i,j,k) = ((fafraczInt(k+1) - fafraczInt(k))*lairef)/modres
+                                    else
+                                        lad_3d(i,j,k) = lad_3d(i,j,modlays-1)
+                                    end if
+                                else
+                                    lad_3d(i,j,k) = 0.0_rk
+                                end if
+                            end do
+
 ! ... calculate zero-plane displacement height/hc and surface (soil+veg) roughness lengths/hc
                             call canopy_zpd(zhc(1:cansublays), fafraczInt(1:cansublays), &
                                 ubzref, z0ghc, lambdars, cdrag, pai, hcmref, hgtref, &
@@ -587,10 +600,14 @@ SUBROUTINE canopy_calcs(nn)
 
 ! ... calculate leaf area density profile from foliage shape function for output (m2/m3)
                         do k=1, modlays
-                            if (k .lt. modlays)  then
-                                lad(loc,k) = ((fafraczInt(k+1) - fafraczInt(k))*lairef)/modres
+                            if (zk(k) .gt. 0.0 .and. zk(k) .le. hcmref) then ! above ground level and at/below canopy top
+                                if (k .lt. modlays)  then
+                                    lad(loc,k) = ((fafraczInt(k+1) - fafraczInt(k))*lairef)/modres
+                                else
+                                    lad(loc,k) = lad(loc,modlays-1)
+                                end if
                             else
-                                lad(loc,k) = lad(loc,modlays-1)
+                                lad(loc,k) = 0.0_rk
                             end if
                         end do
 

diff --git a/src/canopy_canvars_mod.F90 b/src/canopy_canvars_mod.F90
@@ -43,6 +43,7 @@ MODULE canopy_canvars_mod
     real(rk), allocatable :: canBOT        ( : )          ! Canopy bottom wind reduction factors
     real(rk), allocatable :: canTOP        ( : )          ! Canopy top wind reduction factors
     real(rk), allocatable :: lad           ( : , : )      ! Leaf Area Density calculated from foliage shape function (m2/m3)
+    real(rk), allocatable :: lad_3d        ( : , : , : )  ! Leaf Area Density calculated from foliage shape function (m2/m3)
     real(rk), allocatable :: canWIND       ( : , : )      ! canopy wind speeds (m/s)
     real(rk), allocatable :: canWIND_3d    ( : , : , : )  ! canopy wind speeds -- 3D (m/s)
     real(rk), allocatable :: dx            ( : )          ! Model grid cell distance/resolution (m)
@@ -200,6 +201,7 @@ MODULE canopy_canvars_mod
 !-------------------------------------------------------------------------------
 
     TYPE(fld3ddata), ALLOCATABLE, TARGET :: fld3dxyzt ( : )
+    TYPE(fld3ddata), POINTER     :: c_lad
     TYPE(fld3ddata), POINTER     :: c_canwind
     TYPE(fld3ddata), POINTER     :: c_Kz
     TYPE(fld3ddata), POINTER     :: c_rjcf

diff --git a/src/canopy_dealloc.F90 b/src/canopy_dealloc.F90
@@ -45,6 +45,7 @@ SUBROUTINE canopy_dealloc
     if(allocated(ppfd_shade))  deallocate(ppfd_shade)
     if(allocated(ppfd_ave))    deallocate(ppfd_ave)
     if(allocated(lad))         deallocate(lad)
+    if(allocated(lad_3d))      deallocate(lad_3d)
 
 !-------------------------------------------------------------------------------
 ! Deallocate arrays for Canopy Wind

diff --git a/src/canopy_init.F90 b/src/canopy_init.F90
@@ -29,6 +29,7 @@ SUBROUTINE canopy_init
     if(allocated(ppfd_shade))   ppfd_shade(:)  = fillreal
     if(allocated(ppfd_ave))     ppfd_ave(:)    = fillreal
     if(allocated(lad))          lad(:,:)       = fillreal
+    if(allocated(lad_3d))       lad_3d(:,:,:)  = fillreal
 
 !-------------------------------------------------------------------------------
 ! Initialize arrays for Canopy Wind
@@ -41,6 +42,8 @@ SUBROUTINE canopy_init
         if(allocated(canWIND_3d))    canWIND_3d(:,:,:)    = fillreal
         if(allocated(dx))            dx(:)                = fillreal
         if(allocated(dx_2d))         dx_2d(:,:)           = fillreal
+        if(allocated(flameh))        flameh(:)            = fillreal
+        if(allocated(flameh_2d))     flameh_2d(:,:)       = fillreal
         if(allocated(waf))           waf(:)               = fillreal
         if(allocated(waf_2d))        waf_2d(:,:)          = fillreal
     end if

diff --git a/src/canopy_ncf_io_mod.F90 b/src/canopy_ncf_io_mod.F90
@@ -397,6 +397,21 @@ SUBROUTINE canopy_outncf_init
         !-------------------------------------------------------------------------------
         ! Time-varying 3d fields at cell centers.
         !-------------------------------------------------------------------------------
+        c_lad%fld = fillreal
+        c_lad%fldname = 'lad'
+        c_lad%long_name = 'leaf area density'
+        c_lad%units = 'm2 m-3'
+        c_lad%fillvalue = fillreal
+        c_lad%dimnames(1) = 'nlon'
+        c_lad%dimnames(2) = 'nlat'
+        c_lad%dimnames(3) = 'modlays'
+        c_lad%istart(1) = 1
+        c_lad%istart(2) = 1
+        c_lad%istart(3) = 1
+        c_lad%iend(1) = nlon
+        c_lad%iend(2) = nlat
+        c_lad%iend(3) = modlays
+
         if (ifcanwind .or. ifcanwaf) then
             c_canwind%fld = fillreal
             c_canwind%fldname = 'ws'
@@ -804,7 +819,9 @@ SUBROUTINE canopy_outncf_alloc
         ! Time-varying 2d fields at cell centers.
         !-------------------------------------------------------------------------------
 
-        nfld2dxyt = 1  ! canopy height
+        nfld2dxyt = 0
+
+        nfld2dxyt = nfld2dxyt +1  ! canopy height
 
         if (ifcanwind .or. ifcanwaf) then
             nfld2dxyt = nfld2dxyt + 1  !WAF
@@ -832,6 +849,8 @@ SUBROUTINE canopy_outncf_alloc
 
         nfld3dxyzt = 0
 
+        nfld3dxyzt = nfld3dxyzt + 1 !LAD
+
         if (ifcanwind .or. ifcanwaf) then
             nfld3dxyzt = nfld3dxyzt + 1 !CANWIND
         end if
@@ -873,6 +892,10 @@ SUBROUTINE canopy_outncf_alloc
         ENDDO
 
         set_index = 0
+
+        set_index = set_index + 1
+        c_lad    => fld3dxyzt( set_index )
+
         if (ifcanwind .or. ifcanwaf) then
             set_index = set_index + 1
             c_canwind    => fld3dxyzt( set_index )
@@ -1906,6 +1929,7 @@ SUBROUTINE canopy_write_ncf (OUTPREFX)
             !-------------------------------------------------------------------------------
             ! Time-varying 3d fields at cell centers.
             !-------------------------------------------------------------------------------
+            c_lad%fld  = lad_3d
             if (ifcanwind .or. ifcanwaf) then
                 c_canwind%fld  = canWIND_3d
             end if