BIO_EULER.f90 8.23 KB

!============================================================================
!
!                    EULERIAN CALCULATION OF PHYTOPLANKTON GROWTH
!
!       Description: In this routine, the light intensity at each depth is 
!                    calculated using the surface irradiance calculated in
!                    the subroutine DAYLIGHT.f90 and decreases exponentially
!                    with the depth. The phytoplankton growth is then computed
!                    in an Eulerian framework using Monod's equations if bio_model=1
!                    or Droop equations if bio_model=2.
!
!============================================================================

        subroutine bio_euler

        USE global_parameter                                            !shared parameters

!____________________________________________________________________________
!                       LOCAL PARAMETERS
!____________________________________________________________________________

        IMPLICIT NONE

        integer                         :: i,j,jj
        integer                         :: a
        double precision, allocatable   :: light_level(:)
        double precision, allocatable   :: light_bin(:)
        double precision, allocatable   :: dummy1(:),q_dummy(:),mu_dummy(:),nt(:),dt1(:)
        double precision                :: euler_av
        double precision                :: Pdmax=1
        double precision                :: Idark=50
        double precision, allocatable   :: Ndummy(:),Pl(:),Pd(:),CLQUOT(:)
        

        allocate(light_level(Nlevel))
        allocate(dummy1(Nlevel))
        allocate(light_bin(Nbin))
        allocate(q_dummy(Ncat))
        allocate(mu_dummy(Ncat))
        allocate(nt(Ncat))
        allocate(dt1(Ncat))
        allocate(Ndummy(Nbin))
        allocate(Pl(Nbin))
        allocate(Pd(Nbin))
       
        allocate(CLQUOT(Nbin))



!____________________________________________________________________________
!                                LIGHT
!____________________________________________________________________________


        light_level=rad*0.45*4.16*dexp(-k_bg*levels)                    !compute light intensity at each grid interfaces
                                                                        !
        a=0                                                             !
        do i=1,Nlevel                                                   !
        dummy1(Nlevel-a)=light_level(i)                                 !
        a=a+1                                                           !
        end do                                                          !
                                                                        ! 
        light_level=dummy1                                              !

        do i=1,Nbin                                                     !compute light intensity in each grid bin
        light_bin(i)=(light_level(i)+light_level(i+1))/2                !
        end do                                                          !
        


!_____________________________________________________________________________
!                               GROWTH
!_____________________________________________________________________________
        Nut=0d0
        mu_dist=0d0
        do i=1,Nbin
        
        if (photo_model.eq.1) then
        growlight(i)=Pdmax*(1d0-dexp(-light_bin(i)/Idark))
        
        else
        
        CLQUOT=CL/euler
        
        !growlight(i)=Pdmax*(1d0-dexp(-light_bin(i)/(Idark+1000*CLQUOT(i))))
        growlight(i)=light_bin(i)/(light_bin(i)+(500*CLQUOT(i))**2)

        if (light_bin(i)/(rad_max*0.45*4.16).lt.0.01) then
        Clquot_ACC(i)=0.04
        else
        
        Clquot_ACC(i)=0.01+0.03*log(light_bin(i)/(rad_max*0.45*4.16))/log(0.01)
        end if

        end if
                             
        q(i)=Ncell(i)/euler(i)

        if(q(i).gt.q_max) then
        !write(*,*)q(i)
        Ncell(i)=q_max*euler(i)
        q(i)=Ncell(i)/euler(i)
        end if
        
        if (bio_model.eq.1) then                                                        ! Monod equations

        grow(i)=(growlight(i)*(mu_max*(N_euler(i))/(Kn+N_euler(i))) -mortality) -mu_r   ! growth rate
                   
        
        Nut(i)=(growlight(i)*mu_max*(N_euler(i)/(Kn+N_euler(i))))*q0*euler(i)
        euler(i)=euler(i)+grow(i)*time_step*euler(i)                                    !compute the growth with euler solver

        CL(i)=CL(i)+(1/tau*(Clquot_ACC(i)-CLQUOT(i))*euler(i)+grow(i)*CL(i))*time_step

        elseif (bio_model.eq.4) then
        
        
        CLQUOT_dumm=0
        do j=1,nCLcat
        if ((Clquot_ACC(i).gt.CLQUOTDIST(i,j).and.(j.ne.nCLcat).and.(CLQUOTDIST(i,j).gt.0d0)))then
        CLQUOT_dumm(j+1)=CLQUOT_dumm(j+1)+(1/tau*(Clquot_ACC(i)-CLQUOTDIST(i,j))/(CLquat(i,j+1)-CLquat(i,j)))*time_step
        elseif ((Clquot_ACC(i).lt.CLQUOTDIST(i,j).and.(j.ne.1).and.(CLQUOTDIST(i,j).gt.0d0)) then
        CLQUOT_dumm(j-1)=CLQUOT_dumm(j-1)+(1/tau*(Clquot_ACC(i)-CLQUOTDIST(i,j))/(CLquat(i,j)-CLquat(i,j-1)))*time_step
        end if
        end do        
        
        CL_to_mu=(light_bin(i)/(light_bin(i)+(500*CLcat)**2))*(mu_max*(N_euler(i))/(Kn+N_euler(i)))-mortality -mu_r 
        Nut(i)=sum(((light_bin(i)/(light_bin(i)+(500*CLcat)**2))*mu_max*(N_euler(i)/(Kn+N_euler(i))))*q0*CLQUOTDIST(i,:))

        CLQUOTDIST(i,:)=CLQUOTDIST(i,:)+CLQUOT_dumm
        CLQUOTDIST(i,:)=CLQUOTDIST(i,:)+CL_to_mu*CLQUOTDIST(i,:)*time_step

        euler(i)=sum(CLQUOTDIST(i,:))
        
       

        elseif (bio_model.eq.2) then                                                                            !Droop equations

        Uptake_eul(i)=Vmax*N_euler(i)/(Kn+N_euler(i))*(1-q(i)/q_max)

        Ncell(i)=Ncell(i)+(Uptake_eul(i)-kw*q(i))*time_step*euler(i) -mu_r*Ncell(i)*time_step                          
        Nut(i)=(Uptake_eul(i)-kw*q(i))*euler(i)-mu_r*Ncell(i)                                         !nutrient taken
         
        if (q(i).gt.q0) then                                            !compute the growth rate
        grow(i)=growlight(i)*mu_max*(1-q0/q(i))-mu_r    
        !grow(i)=growlight(i)*mu_max*(q(i)/q_max)-mu_r  
        else                                                                            !
        grow(i)=-mortality -mu_r                                                        !
        end if                                                                         ! 
        euler(i)=euler(i)+grow(i)*time_step*euler(i)                                    !compute the growth with euler solver

         

        elseif (bio_model.eq.3) then

          q_dummy=0
       nt=0
        dt1=0
        do j=1,Ncat
        uptake_dist(j)=Vmax*N_euler(i)/(Kn+N_euler(i))*(1-q_cat(j)/q_max)
        Nut(i)=Nut(i)+uptake_dist(j)*q_dist(i,j)-kw*q_cat(j)*q_dist(i,j)

        if (q_cat(j).gt.q0) then
        q_to_mu(j)=growlight(i)*mu_max*(1-q0/q_cat(j))-mu_r
        else
         q_to_mu(j)=-mortality -mu_r
        end if
       q_to_mu(1)=-mortality -mu_r

        if((uptake_dist(j).lt.q_to_mu(j)*q_cat(j)+kw*q_cat(j)).and.(j.ne.1).and.(q_dist(i,j).gt.0d0)) then
        dt1(j)=(q_cat(j-1)-q_cat(j))/(uptake_dist(j)-q_to_mu(j)*q_cat(j)-kw*q_cat(j))
        nt(j)=dt1(j)/time_step;
        q_dummy(j-1)=q_dummy(j-1)+q_dist(i,j)/nt(j)
        q_dummy(j)=q_dummy(j)-q_dist(i,j)/nt(j)
        elseif((uptake_dist(j).gt.q_to_mu(j)*q_cat(j)+kw*q_cat(j)).and.(j.ne.Ncat).and.(q_dist(i,j).gt.0d0)) then
         dt1(j)=(q_cat(j+1)-q_cat(j))/(uptake_dist(j)-q_to_mu(j)*q_cat(j)-kw*q_cat(j))
        nt(j)=dt1(j)/time_step;
        q_dummy(j+1)=q_dummy(j+1)+q_dist(i,j)/nt(j)
         q_dummy(j)=q_dummy(j)-q_dist(i,j)/nt(j)
        end if


        end do

       ! do j=1,Ncat
       ! do jj=1,Ncat
       ! if((q_to_mu(j).ge.mu_edge(jj)).and.(q_to_mu(j).le.mu_edge(jj+1))) then
       ! mu_dist(i,jj)=mu_dist(i,jj)+q_dist(i,j)
       ! end if
       ! end do
      !  end do
       ! write(*,*)euler(i),sum(q_dist(i,:))+ sum(q_dist(i,:)*q_to_mu*time_step)

        q_dist(i,:)=q_dist(i,:)+ q_dist(i,:)*q_to_mu*time_step
        q_dist(i,:)=q_dist(i,:)+q_dummy
      euler(i)=sum(q_dist(i,:))
        



        end if

        end do
            
!===============================================================================
end subroutine bio_euler