Merge pull request #233 from derpycode/DEV_GMD23

Dev gmd23

genie-ecogem/data/input/3Diat4ZP_PiEu.eco

@@ -0,0 +1,41 @@
+
+ 01 02 03
+ \/ \/ \/
+
+-START-OF-DATA-
+ Diatom 2.00 1
+ Diatom 20.00 1
+ Diatom 200.00 1
+ Picoplankton 0.6 1
+ Picoplankton 2.00 1
+ Eukaryote 20.00 1
+ Eukaryote 200.00 1
+ Zooplankton 6.00 1
+ Zooplankton 20.00 1 
+ Zooplankton 200.00 1
+ Zooplankton 2000.00 1
+-END-OF-DATA-
+
+ /\ /\ /\
+ 01 02 03
+
+DATA FORMAT AND ORDER
+---------------------
+
+COLUMN #01: plankton functional type name
+COLUMN #02: plankton diameter (micrometers)
+COLUMN #03: number of randomised replicates
+
+INFO: TRACER ASSIGNMENT RULES
+-----------------------------
+Plankton functional type one of: Prochlorococcus
+ Synechococcus
+ Picoeukaryote
+ Picoplankton
+ Diatom
+ Coccolithophore
+ Diazotroph
+ Eukaryote
+ Phytoplankton
+ Zooplankton
+ Mixotroph

genie-ecogem/data/input/NoDiat4ZP_PiEu.eco

@@ -0,0 +1,38 @@
+
+ 01 02 03
+ \/ \/ \/
+
+-START-OF-DATA-
+ Picoplankton 0.6 1
+ Picoplankton 2.00 1
+ Eukaryote 20.00 1
+ Eukaryote 200.00 1
+ Zooplankton 6.00 1
+ Zooplankton 20.00 1 
+ Zooplankton 200.00 1
+ Zooplankton 2000.00 1
+-END-OF-DATA-
+
+ /\ /\ /\
+ 01 02 03
+
+DATA FORMAT AND ORDER
+---------------------
+
+COLUMN #01: plankton functional type name
+COLUMN #02: plankton diameter (micrometers)
+COLUMN #03: number of randomised replicates
+
+INFO: TRACER ASSIGNMENT RULES
+-----------------------------
+Plankton functional type one of: Prochlorococcus
+ Synechococcus
+ Picoeukaryote
+ Picoplankton
+ Diatom
+ Coccolithophore
+ Diazotroph
+ Eukaryote
+ Phytoplankton
+ Zooplankton
+ Mixotroph

genie-ecogem/src/fortran/ecogem.f90

@@ -398,6 +398,11 @@ subroutine ecogem( &
  elseif (io.eq.iChlo) then
  assimilated(io,:) = 0.0 ! don't assimilate chlorophyll
  endif
+ do jpred=1,npmax ! JDW - check do loop is needed
+ if (io.eq.iSili .and. pft(jpred).eq.'zooplankton') then ! .and. silicify(:).eq.0.0) then
+ assimilated(io,jpred) = 0.0 ! don't assimilate Si if not selected - Scott April 2019
+ endif
+ enddo ! Aaron Diatom 23
  unassimilated(:,:) = 1.0 - assimilated(:,:)
  enddo
  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -447,7 +452,8 @@ subroutine ecogem( &
  io=nut2quota(ii)
  dbiomassdt(io,:) = dbiomassdt(io,:) + up_inorg(ii ,:) * BioC(:)
  enddo
-
+ ! Take into account nitrogen fixation here scaled with dbiomass in phosphorus - Fanny Jun20
+ dbiomassdt(iNitr,:) = merge(dbiomassdt(iPhos,:)*40.0,dbiomassdt(iNitr,:),pft.eq.'diazotroph') ! Aaron Diatom 23
  if (eco_uptake_fluxes) then
  do io=1,iomax
  AP_uptake(io,:) = dbiomassdt(io,:) 
@@ -536,8 +542,8 @@ subroutine ecogem( &
 
  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  orgbGraz(:,:) = 0.0
-
- do io=1,iomax
+ do io=1,iomax !! Be aware this calculates also for Si "OM" - Fanny/Jamie - Jun19
+
  ! mortality
  dorgmatdt(io,1) = dorgmatdt(io,1) + sum(loc_biomass(io,:) * mortality(:) * beta_mort(:) ) ! fraction to DOM
  dorgmatdt(io,2) = dorgmatdt(io,2) + sum(loc_biomass(io,:) * mortality(:) * beta_mort_1(:)) ! fraction to POM
@@ -761,6 +767,7 @@ subroutine ecogem( &
  if (nquota) dum_egbg_sfcpart(is_PON ,:,:,:) = orgmat_flux(iNitr,2,:,:,:) / 1.0e3 / conv_m3_kg ! convert back to mol kg^{-1} s^{-1}
  if (pquota) dum_egbg_sfcpart(is_POP ,:,:,:) = orgmat_flux(iPhos,2,:,:,:) / 1.0e3 / conv_m3_kg ! convert back to mol kg^{-1} s^{-1}
  if (fquota) dum_egbg_sfcpart(is_POFe ,:,:,:) = orgmat_flux(iIron,2,:,:,:) / 1.0e3 / conv_m3_kg ! convert back to mol kg^{-1} s^{-1}
+ if (squota) dum_egbg_sfcpart(is_opal ,:,:,:) = (orgmat_flux(iSili,1,:,:,:) + orgmat_flux(iSili,2,:,:,:)) / 1.0e3 / conv_m3_kg ! convert back to mol kg^{-1} s^{-1} ! FMM/JDW June 2019, "DOSi + POSi" Aaron Diatom 23
  !ckc Isotopes particulate
  dum_egbg_sfcpart(is_POC_13C,:,:,:) = orgmatiso_flux(iCarb13C,2,:,:,:) / 1.0e3 / conv_m3_kg ! convert back to mol kg^{-1} s^{-1}
 
@@ -796,10 +803,30 @@ subroutine ecogem( &
  dum_egbg_sfcdiss(io_DIC_13C ,:,:,:) = dum_egbg_sfcdiss(io_DIC_13C,:,:,:) - 1.0 * dum_egbg_sfcpart(is_CaCO3_13C,:,:,:)
  endif
 
- ! set initial values for protected fraction of POM and CaCO3
- dum_egbg_sfcpart(is_POC_frac2 ,:,:,n_k) = par_bio_remin_POC_frac2
- dum_egbg_sfcpart(is_CaCO3_frac2,:,:,n_k) = par_bio_remin_CaCO3_frac2
+ ! set initial values for protected fraction of POM and CaCO3 Aaron Diatom 23
+ ! Added ballast parameterisation - Fanny, Aug20
+ ! NEED TO MAKE CARRYING COEF EXPLICIT (HARD CODING HERE BECAUSE WANT TO FIND A WAY TO DEFINE THEM ONLY ONE) - Fanny Aug20
+ ! Changed frac2 modification to be for imld:n_k rather than n_k
+ if (ctrl_bio_remin_POC_ballast_eco) then 
+ !dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k) = MERGE( &
+ ! & ( par_bio_remin_kc(:,:)*dum_egbg_sfcpart(is_CaCO3,:,:,imld:n_k) + &
+ ! & par_bio_remin_ko(:,:)*dum_egbg_sfcpart(is_opal,:,:,imld:n_k) + &
+ ! & par_bio_remin_kl(:,:)*dum_egbg_sfcpart(is_det,:,:,imld:n_k) ) &
+ ! & /dum_egbg_sfcpart(is_POC,:,:,imld:n_k) , &
+ ! & 0.0, dum_egbg_sfcpart(is_POC,:,:,imld:n_k) > const_real_nullsmall)
+ dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k) = MERGE( &
+ & ( 0.085*dum_egbg_sfcpart(is_CaCO3,:,:,imld:n_k) + &
+ & 0.025*dum_egbg_sfcpart(is_opal,:,:,imld:n_k) + &
+ & 0.0*dum_egbg_sfcpart(is_det,:,:,imld:n_k) ) &
+ & /dum_egbg_sfcpart(is_POC,:,:,imld:n_k) , &
+ & 0.0, dum_egbg_sfcpart(is_POC,:,:,imld:n_k) > const_real_nullsmall)
+ dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k) = MERGE( &
+ & 1.0,dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k), dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k) > 1.0)
+ else
+ dum_egbg_sfcpart(is_POC_frac2,:,:,imld:n_k) = par_bio_remin_POC_frac2
+ endif
 
+ dum_egbg_sfcpart(is_CaCO3_frac2,:,:,imld:n_k) = par_bio_remin_CaCO3_frac2
  ! ---------------------------------------------------------- !
  ! END
  ! ---------------------------------------------------------- !

genie-ecogem/src/fortran/ecogem_box.f90

@@ -35,7 +35,7 @@ SUBROUTINE quota_status( &
  ! DEFINE LOCAL VARIABLES
  ! ---------------------------------------------------------- !
  integer :: io
- !
+  !
  ! *****************************************************************
  ! ******************** Evaluate Quota Status **********************
  ! *****************************************************************
@@ -80,6 +80,7 @@ SUBROUTINE quota_limitation( &
  ! DEFINE LOCAL VARIABLES
  ! ---------------------------------------------------------- !
  integer :: io
+ real,dimension(iomax,npmax) :: denom ! JDW Aaron Diatom 23 
  !
  ! *****************************************************************
  ! ******************** Evaluate Limitation ***********************
@@ -100,21 +101,30 @@ SUBROUTINE quota_limitation( &
  VLlimit(:) = 0.0 ! (npmax)
  qreg(:,:) = 0.0 ! (iomax,npmax)
  qreg_h(:,:) = 0.0 ! (iomax,npmax)
-
+ ! JDW: precalculate qmin-qmax to avoid divide by zero errorss
+ denom = merge(1.0/(qmax - qmin),0.0,qmax.gt.0.0) ! JDW: calculate 1.0/denominator and take care of instances of 1.0/0.0 Aaron Diatom 23
  ! Calculate quota limitation terms
  ! N and Si take linear form
- if (nquota) limit(iNitr,:) = (quota(iNitr,:) - qmin(iNitr,:)) / ( qmax(iNitr,:) - qmin(iNitr,:))
- if (squota) limit(iSili,:) = (quota(iSili,:) - qmin(iSili,:)) / ( qmax(iSili,:) - qmin(iSili,:))
+ ! Modified to take into account no N limitation by diazotrophs - Fanny Jun20
+ if (nquota) then
+ limit(iNitr,:) = (quota(iNitr,:) - qmin(iNitr,:)) * denom(iNitr,:)
+ limit(iNitr,:) = merge(1.0,limit(iNitr,:),pft.eq.'diazotroph')
+ endif
+ if (squota) limit(iSili,:) = (quota(iSili,:) - qmin(iSili,:)) * denom(iSili,:) ! JDW
+ !if (squota) limit(iSili,:) = (quota(iSili,:) - qmin(iSili,:)) / ( qmax(iSili,:) - qmin(iSili,:)) ! original Aaron Diatom 23
  ! P and Fe take normalised Droop form
  if (pquota) limit(iPhos,:) = (1.0 - qmin(iPhos,:)/quota(iPhos,:)) / (1.0 - qmin(iPhos,:)/qmax(iPhos,:) )
  if (fquota) limit(iIron,:) = (1.0 - qmin(iIron,:)/quota(iIron,:)) / (1.0 - qmin(iIron,:)/qmax(iIron,:) )
 
  ! Set Von Leibig limitation according to most limiting nutrient (excluding iCarb=1)
- VLlimit(:) = minval(limit(2:iomax,:),1)
+ ! VLlimit(:) = minval(limit(2:iomax,:),1) ! original
+ VLlimit(:) = minval(limit(2:max(iNitr,iPhos,iIron),:),1) ! JDW: calculate limitation for N,P,Fe only
+ VLlimit = merge(minval(limit(2:iomax,:),1),minval(limit(2:max(iNitr,iPhos,iIron),:),1),silicify.eq.1.0) ! JDW: in case of diatom reset taking into account SiO2 Aaron Diatom 23
 
  do io = 2,iomax ! skip carbon index; quota = X:C biomass ratio
  ! Calculate linear regulation term
- qreg(io,:) = (qmax(io,:) - quota(io,:)) / (qmax(io,:) - qmin(io,:) )
+ qreg(io,:) = (qmax(io,:) - quota(io,:)) * denom(io,:) ! JDW
+ !qreg(io,:) = (qmax(io,:) - quota(io,:)) / (qmax(io,:) - qmin(io,:) ) ! original Aaron Diatom 23
  ! Transform regulation term using hill number
  qreg_h(io,:) = qreg(io,:) ** hill
  enddo
@@ -295,7 +305,9 @@ SUBROUTINE photosynthesis( &
  if (useNO3) VCN(:) = VCN(:) + up_inorg(iNO3,:)
  if (useNO2) VCN(:) = VCN(:) + up_inorg(iNO2,:)
  if (useNH4) VCN(:) = VCN(:) + up_inorg(iNH4,:)
- elseif (pquota) then
+ ! To account for nitrogen fixation ! Fanny - June 2020 - Still need to check Moore et al (2002) + parameterise N:P_diazo (=40.0) - probably should scale with nitrogen fixation rate!!
+ VCN(:) = merge(up_inorg(iPO4,:) * 40.0,VCN(:),pft.eq.'diazotroph') ! Aaron Diatom 23
+ elseif (pquota) then
  VCN(:) = up_inorg(iPO4,:) * 16.0
  else
  print*,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"