merged yunkebranch, resolved issue with c4 data (reversed lat), added…

… plotting example in Rmd
stineb · Feb 2, 2021 · fa5c806 · fa5c806
2 parents eda6b7d + 9715c2a
commit fa5c806
Show file tree

Hide file tree

Showing 15 changed files with 684 additions and 163 deletions.
diff --git a/Makefile b/Makefile
@@ -8,7 +8,7 @@
 # Set PROFILE to one of benilaptop, cx1, euler, pgi, or beniimac
 ##################################
 # change PROFILE to euler
-PROFILE=benilaptop
+PROFILE=beniimac
 
 ##################
 # pgf profile

diff --git a/linkdirs_sofun.py b/linkdirs_sofun.py
@@ -144,7 +144,7 @@
 
 	## c4 percentage
 	##--------------------------------------
-	src = dataroot + 'c4_still/final/c4_percentage.nc'
+	src = dataroot + 'c4_still/final/c4_percentage_revlat.nc'
 	dst = 'input/global/landcover'
 	if not os.path.isdir( dst ):
 		os.system( 'mkdir -p ' + dst )

diff --git a/params_std/params_gpp_pmodel.dat b/params_std/params_gpp_pmodel.dat
@@ -13,7 +13,7 @@ tau_acclim                              10.0000000
 
 kphio_Gr3                               0.09423773
 kphio_GN3                               0.09423773
-kphio_Gr4                               0.09423773
+kphio_Gr4                               1.00000000
 kphio_TND                               0.09423773
 kphio_TNE                               0.09423773
 kphio_TrD                               0.09423773

diff --git a/params_std/params_nimpl.dat b/params_std/params_nimpl.dat
@@ -1,5 +1,5 @@
 ###############################################################
-# PMODEL PARAMETERS
+# NIMPL PARAMETERS
 #--------------------------------------------------------------
 #                                      |
 # write value behind this position:    |->
@@ -25,14 +25,31 @@ vpd_alnpp                               -1.0192
 intersect_alnpp                         -9.2375
 
 #leaf C/N model
-lma_leafcn                              0.4470
-vcmax25_leafcn                          -0.1068
-intersect_leafcn                        1.5878
+lma_leafcn                              0.0161
+vcmax25_leafcn                          0.0041
+cmass_leafcn                            0.5000
+
+#nre model
+tg_nre                                  -0.0679
+vpd_nre                                 0.4217
+intersect_nre                           1.4541
 
 #constant ratio - derived from TRY database (N=2365 for root; N=13 for wood)
 root_cn                                 42
 wood_cn                                 97
 
+#Grassland models
+#BP/GPP model
+cue_grass                               0.4046
+
+#ANPP/NPP model
+tg_anpp_grass                           0.0385
+alpha_anpp_grass                        1.8720
+intersect_anpp_grass                    -2.5642
+
+#Constant of root C/N
+root_cn_grass                           46
+
 #                                      |
 # write value *behind* this position:  |->
 ## EOF 
diff --git a/plots_sofun.Rmd b/plots_sofun.Rmd
@@ -0,0 +1,49 @@
+---
+title: "Visualise SOFUN"
+author: "Beni Stocker"
+date: "2/2/2021"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+library(ggplot2)
+library(dplyr)
+library(rbeni)
+```
+
+## Read output
+
+Annual
+
+```{r}
+nc_gpp <- read_nc_onefile("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.a.gpp.nc", varnam = "gpp")
+nc_pet <- read_nc_onefile("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.a.pet.nc", varnam = "pet")
+nc_alpha <- read_nc_onefile("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.a.alpha.nc", varnam = "alpha")
+```
+
+Daily: Take mean over daily values and read file.
+
+```{r}
+nc_wcont <- nc_timmean("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.d.wcont.nc")
+nc_vpd <- nc_timmean("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.d.vpd.nc")
+```
+
+## Global maps
+
+Annual
+
+```{r}
+plot_map3(nc_gpp, varnam = "gpp")
+plot_map3(nc_pet, varnam = "pet")
+plot_map3(nc_alpha, varnam = "alpha")
+plot_map3(nc_wcont, varnam = "wcont")
+plot_map3(nc_vpd, varnam = "vpd")
+```
+
+## Integrate globally
+
+```{r}
+df <- integrate_lonlat("~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.2000.a.gpp.nc", 
+                       "~/sofun/output_nc/global_FULL_MODIS-C006_MOD15A2_v3.3.fland.nc",
+                       varnam = "gpp")
+```
diff --git a/src/Makefile b/src/Makefile
@@ -33,7 +33,7 @@ DEMO_PMODEL_MODS=params_core.mod.f90 sofunutils.mod.f90 io_netcdf.mod.f90 interf
 PMODEL_SWBM_MODS=params_core.mod.f90 classdefs.mod.f90 sofunutils.mod.f90 io_netcdf.mod.f90 params_siml.mod.f90 params_domain.mod.f90 grid_siterun.mod.f90 params_soil.mod.f90 forcing_siterun_pmodel.mod.f90 interface_biosphere.mod.f90 tile_pmodel.mod.f90 waterbal_swbm.mod.f90 soiltemp_sitch.mod.f90 plant_pmodel.mod.f90 gpp_pmodel.mod.f90 vegdynamics_pmodel.mod.f90 biosphere_pmodel.mod.f90
 
 # List of source code files for reduced setup, executing SPLASH and P-MODEL for global simulations:
-GPMODEL_MODS=params_core.mod.f90 classdefs.mod.f90 sofunutils.mod.f90 io_netcdf.mod.f90 params_siml.mod.f90 params_domain_global.mod.f90 grid_global.mod.f90 params_soil_global.mod.f90 forcing_global_pmodel.mod.f90 interface_biosphere.mod.f90 plant_pmodel.mod.f90 tile_pmodel.mod.f90 waterbal_splash.mod.f90 soiltemp_sitch.mod.f90 gpp_pmodel.mod.f90 vegdynamics_pmodel.mod.f90 biosphere_pmodel.mod.f90
+GPMODEL_MODS=params_core.mod.f90 classdefs.mod.f90 sofunutils.mod.f90 io_netcdf.mod.f90 params_siml.mod.f90 params_domain_global.mod.f90 grid_global.mod.f90 params_soil_global.mod.f90 forcing_global_pmodel.mod.f90 interface_biosphere.mod.f90 plant_pmodel.mod.f90 tile_pmodel.mod.f90 waterbal_splash.mod.f90 soiltemp_sitch.mod.f90 gpp_pmodel.mod.f90 vegdynamics_pmodel.mod.f90 nuptake_impl.mod.f90 biosphere_pmodel.mod.f90
 
 # List of source code files for reduced setup, executing SPLASH for global simulations:
 GSPLASH_MODS=params_core.mod.f90 classdefs.mod.f90 sofunutils.mod.f90 io_netcdf.mod.f90 params_siml.mod.f90 params_domain_global.mod.f90 grid_global.mod.f90 params_soil_global.mod.f90 forcing_global_wmodel.mod.f90 interface_biosphere.mod.f90 tile_wmodel.mod.f90 waterbal_splash.mod.f90 soiltemp_sitch.mod.f90 biosphere_wmodel.mod.f90

diff --git a/src/biosphere_pmodel.mod.f90 b/src/biosphere_pmodel.mod.f90
@@ -10,6 +10,8 @@ module md_biosphere
   use md_tile, only: tile_type, tile_fluxes_type, initglobal_tile, initdaily_tile_fluxes, getpar_modl_tile, diag_daily, diag_annual, init_annual
   use md_interface, only: getout_daily_forcing, initoutput_forcing, initio_nc_forcing, writeout_nc_forcing
   use md_soiltemp, only: getout_daily_soiltemp, soiltemp, initoutput_soiltemp
+  use md_nuptake_impl, only: nuptake_impl, get_preds_nimpl, getpar_nimpl, initio_nc_nimpl, initoutput_nimpl, getout_annual_nimpl, writeout_nc_nimpl
+  ! use md_sofunutils, only: calc_patm
 
   implicit none
 
@@ -60,6 +62,7 @@ function biosphere_annual() result( out_biosphere )
       call getpar_modl_plant()
       call getpar_modl_waterbal()
       call getpar_modl_gpp()
+      call getpar_nimpl()
       if (verbose) print*, '... done'
 
       !----------------------------------------------------------------
@@ -71,6 +74,11 @@ function biosphere_annual() result( out_biosphere )
       call initglobal_tile(  tile(:,:),  size(interface%grid) )
       if (verbose) print*, '... done'
 
+      !----------------------------------------------------------------
+      ! Get nimpl predictor fields (only available inside nimpl module)
+      !----------------------------------------------------------------
+      call get_preds_nimpl( interface%domaininfo, interface%grid )
+
     endif 
 
     !----------------------------------------------------------------
@@ -81,6 +89,7 @@ function biosphere_annual() result( out_biosphere )
       call initio_nc_forcing()
       call initio_nc_gpp()
       call initio_nc_waterbal()
+      call initio_nc_nimpl()
       if (verbose) print*, '... done'
     end if
 
@@ -94,6 +103,7 @@ function biosphere_annual() result( out_biosphere )
       call initoutput_plant(    size(interface%grid) )
       call initoutput_forcing(  size(interface%grid) )
       call initoutput_soiltemp( size(interface%grid) )
+      call initoutput_nimpl(    size(interface%grid) )
       if (verbose) print*, '... done'
     end if
 
@@ -125,7 +135,7 @@ function biosphere_annual() result( out_biosphere )
         !----------------------------------------------------------------
         ! LOOP THROUGH MONTHS
         !----------------------------------------------------------------
-        doy=0   ! day of the year
+        doy=0
         monthloop: do moy=1,nmonth
 
           !----------------------------------------------------------------
@@ -254,7 +264,16 @@ function biosphere_annual() result( out_biosphere )
         !----------------------------------------------------------------
         ! annual diagnostics
         !----------------------------------------------------------------
+        if (verbose) print*,'calling diag_annual() ... '
         call diag_annual( tile(:,jpngr), tile_fluxes(:) )
+        if (verbose) print*,'... done'
+
+        !----------------------------------------------------------------
+        ! Statistical relationships with GPP to get N uptake
+        !----------------------------------------------------------------
+        if (verbose) print*,'calling nuptake_impl() ... '
+        call nuptake_impl( jpngr, interface%grid(jpngr)%dogridcell, tile(:,jpngr), tile_fluxes(:), interface%steering%init )
+        if (verbose) print*,'... done'
 
         !----------------------------------------------------------------
         ! collect annual output
@@ -263,6 +282,7 @@ function biosphere_annual() result( out_biosphere )
           if (verbose) print*,'calling getout_annual_() ... '
           ! call getout_annual_plant( tile(:)%plant(:), jpngr )
           call getout_annual_gpp( jpngr, tile_fluxes(:) )
+          call getout_annual_nimpl( jpngr, tile(:,jpngr), tile_fluxes(:) )
           if (verbose) print*,'... done'
         end if
 
@@ -283,6 +303,7 @@ function biosphere_annual() result( out_biosphere )
       call writeout_nc_forcing()
       call writeout_nc_gpp()
       call writeout_nc_waterbal()
+      call writeout_nc_nimpl()
       if (verbose) print*,'... done'
     end if
 

diff --git a/src/forcing_global_pmodel.mod.f90 b/src/forcing_global_pmodel.mod.f90
@@ -523,15 +523,16 @@ end function getfapar
 
   function get_fpc_grid( domaininfo, grid, params_siml ) result( fpc_grid_field )
     !////////////////////////////////////////////////////////////////
-    ! xxx add explanation here
+    ! Reads vegetation cover file to define fractional PFT coverage
+    ! for each PFT (currently distinguishing between C3 and C4)
     !----------------------------------------------------------------  
     use md_params_siml, only: paramstype_siml
     use md_params_core, only: npft, eps
 
     ! arguments
     type( domaininfo_type ), intent(in) :: domaininfo
     type( gridtype ), dimension(domaininfo%maxgrid), intent(in) :: grid
-    type( paramstype_siml ), intent(in) :: params_siml
+    type( paramstype_siml ), intent(inout) :: params_siml
 
     ! function return variable
     real, dimension(npft,domaininfo%maxgrid) :: fpc_grid_field
@@ -554,7 +555,7 @@ function get_fpc_grid( domaininfo, grid, params_siml ) result( fpc_grid_field )
     character(len=3), parameter :: latname = "lat"
     ! character(len=100), parameter :: dimname_pft = "z"
     character(len=100), parameter :: varname = "c4"
-    character(len=100), parameter :: filnam = "./input/global/landcover/c4_percentage.nc"
+    character(len=100), parameter :: filnam = "./input/global/landcover/c4_percentage_revlat.nc"
 
     !----------------------------------------------------------------  
     ! Get vegetation cover information from file
@@ -638,6 +639,12 @@ function get_fpc_grid( domaininfo, grid, params_siml ) result( fpc_grid_field )
     ! close NetCDF files
     call check( nf90_close( ncid ) )
 
+    ! override: activate both C3 and C4
+    params_siml%lGr3 = .true.
+    params_siml%lGr4 = .true.
+
+    if (npft /= 2) stop 'Aborting. Set npft = 2.'
+
     ! read from array to define land cover 'fpc_grid_field'
     n_noinfo = 0
     gridcellloop: do jpngr=1,domaininfo%maxgrid
@@ -646,58 +653,25 @@ function get_fpc_grid( domaininfo, grid, params_siml ) result( fpc_grid_field )
 
       if ( tmp==ncfillvalue .and. grid(jpngr)%dogridcell ) then
 
-        n_noinfo = n_noinfo + 1
-        ! fpc_grid_field(:,jpngr) = 1.0 / real( npft )
-
-      else        
-
-        fpc_grid_field(:,jpngr) = 0.0
-
-        if (npft==2) then
-          if (jpngr==1) print*,'GET_FPC_GRID: npft=2 ==> assuming distinction between C3 and C4 grasslands'
-
-          pft = 0
-
-          if ( params_siml%lGr3 ) then
-            ! xxx dirty: call all grass vegetation types 'Gr3'
-            pft = pft + 1
-
-            ! C3-PFT (can be grass or not grass)
-            fpc_grid_field(pft,jpngr) = 1.0 - tmp
-
-          end if
-
-          if ( params_siml%lGr4 ) then
-            ! xxx dirty: call all grass vegetation types 'Gr3'
-            pft = pft + 1
-
-            ! C4-PFT
-            fpc_grid_field(pft,jpngr) = tmp
-
-          end if
-
-        else if (npft==1) then
-
-          if (jpngr==1) print*,'GET_FPC_GRID: npft=1 ==> assuming all is C3'
+        ! For gridcells with missing information, assume 100% C3
+        ! C3-PFT (can be grass or not grass)
+        fpc_grid_field(1,jpngr) = 1.0
 
-          pft = 1
-          fpc_grid_field(pft,jpngr) = 1.0 !- tmp
+        ! C4-PFT
+        fpc_grid_field(2,jpngr) = 0.0
 
-        else
-
-          stop 'GET_FPC_GRID: only implemented for npft = 1 or 2.'
+      else        
 
-        end if
+        ! C3-PFT (can be grass or not grass)
+        fpc_grid_field(1,jpngr) = 1.0 - tmp
 
+        ! C4-PFT
+        fpc_grid_field(2,jpngr) = tmp
 
       end if
 
-      ! if ( abs(sum(fpc_grid_field(:,jpngr)) - 1.0)>eps .and. sum(fpc_grid_field(:,jpngr)) > 0.0 ) &
-      !   fpc_grid_field(:,jpngr) = fpc_grid_field(:,jpngr) / sum( fpc_grid_field(:,jpngr) )
-
     end do gridcellloop
 
-
     ! deallocate memory again (the problem is that climate input files are of unequal length in the record dimension)
     deallocate( vegtype_arr )
     deallocate( lon_arr )

diff --git a/src/gpp_pmodel.mod.f90 b/src/gpp_pmodel.mod.f90
@@ -322,9 +322,8 @@ subroutine gpp( tile, tile_fluxes, co2, climate, vegcover, grid, do_soilmstress,
       !----------------------------------------------------------------
       ! Vcmax25
       !----------------------------------------------------------------
-      tile(lu)%plant(pft)%vcmax25 = out_pmodel%vcmax25
+      ! tile(lu)%plant(pft)%vcmax25 = out_pmodel%vcmax25
       tile_fluxes(lu)%plant(pft)%vcmax25 = out_pmodel%vcmax25
-      !print*,'tile_fluxes(lu)%plant(pft)%vcmax25      ', tile_fluxes(lu)%plant(pft)%vcmax25
 
 
     end do pftloop
@@ -784,7 +783,11 @@ function pmodel( kphio, ppfd, co2, tc, vpd, patm, c4, method_optci, method_jmaxl
       jmax = 0.0
       jmax25 = 0.0
     else
-      fact_jmaxlim = vcmax * (ci + 2.0 * gammastar) / (kphio * ppfd * (ci + kmm))
+      if (kphio==0.0) then
+        fact_jmaxlim = 0.0
+      else
+        fact_jmaxlim = vcmax * (ci + 2.0 * gammastar) / (kphio * ppfd * (ci + kmm))
+      end if
       !print*,'fact_jmaxlim       ', fact_jmaxlim
       if (fact_jmaxlim >= 1 .or. fact_jmaxlim <= 0) then
         jmax = dummy
@@ -1753,6 +1756,11 @@ function calc_ftemp_kphio( dtemp, c4 ) result( ftemp )
 
     if (c4) then
       ftemp = (-0.008 + 0.00375 * dtemp - 0.58e-4 * dtemp**2) * 8.0 ! Based on calibrated values by Shirley
+      if (ftemp < 0.0) then
+        ftemp = 0.0
+      else
+        ftemp = ftemp
+      end if
     else
       ftemp = 0.352 + 0.022 * dtemp - 3.4e-4 * dtemp**2  ! Based on Bernacchi et al., 2003
     end if
@@ -2069,7 +2077,7 @@ subroutine getout_annual_gpp( jpngr, tile_fluxes )
 
       ! canopy level
       outagpp(jpngr)     = tile_fluxes(1)%canopy%agpp          ! take only LU = 1
-      outavcmax25(jpngr) = tile_fluxes(1)%canopy%avcmax25_mean ! take GPP-weighted mean, only LU = 1
+      outavcmax25(jpngr) = tile_fluxes(1)%canopy%avcmax25_max  ! take GPP-weighted mean, only LU = 1 !avcmax25_weightedmean
 
       ! pft level
       !outagpp(jpngr)     = tile_fluxes(1)%plant(pft)%agpp     ! canopy and pft level were the same!