BIO_LAGRANGE.f90 8.52 KB

!=============================================================================!
!                                                                             !
!                    LAGRANGIAN 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 Lagrangian framework using Monod's equations if    !
!                    bio_model=1 or Droop equations if bio_model=2.           !
!                                                                             !
!=============================================================================!


        subroutine light_growth

        use global_parameter                                                    !shared parameters


!____________________________________________________________________________
!                               LOCAL PARAMETERS
!____________________________________________________________________________

        IMPLICIT NONE

        integer                                 :: i
        integer                                 :: j
        integer                                 :: ii
        integer                                 :: a
        integer                                 :: div
        double precision, allocatable           :: dummy1(:)
        double precision, allocatable           :: light_level(:)
        double precision, allocatable           :: light_bin(:)
        double precision, allocatable           :: nb(:)
        double precision                        :: Pdmax=1
        double precision                        :: Idark=50      
        double precision, dimension(100000)     :: Fi
        double precision,allocatable            :: nb1(:)
        double precision,allocatable            :: Pl(:),Pd(:)
        double precision, allocatable           :: CLQUOT_lag(:)

        allocate(light_level(Nlevel))
        allocate(dummy1(Nlevel))
        allocate(light_bin(Nbin))
        allocate(nb(Nbin))
        allocate(nb1(Nbin))
       
        allocate(Pl(Nbpart))
        allocate(Pd(Nbpart))
        allocate(CLQUOT_lag(Nbpart))

        Fi(1:100000)=0
        div=0
        nb(1:Nbin)=0
        nb1=0d0


!___________________________________________________________________________
!                               LIGHT
!___________________________________________________________________________


        light_level=rad*0.45*4.16*dexp(-k_bg*levels)                            !compute light intensity at each grid interface
        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
!___________________________________________________________________________


        do i=1,nbpart

                qlag(i)=Ncell_lag(i)/weight(i)
                if (qlag(i).gt.q_max)then
                !write(*,*)qlag(i)
                Ncell_lag(i)=q_max*weight(i)
                qlag(i)=Ncell_lag(i)/weight(i)
                end if


                bin=1+floor(z_cell(i)*z_res)   
                                         !find the bin where the particle is.

                if (photo_model.eq.1) then
                uptlight(i)=Pdmax*(1d0-dexp(-light_bin(bin)/Idark))                     !compute P(I)

                else

                CLQUOT_lag(i)=CL_lag(i)/weight(i)
               ! uptlight(i)=Pdmax*(1d0-dexp(-light_bin(bin)/(Idark+1000*CLQUOT_lag(i))))
                 uptlight(i)=light_bin(bin)/(light_bin(bin)+(500*CLQUOT_lag(i))**2)


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

               

                end if
 

                if (bio_model.eq.1) then                                                ! Use Monod's equations

                Dphyto_lag(i)=(uptlight(i)*mu_max*(N_lag(bin)/(Kn+N_lag(bin))) -mortality)-mu_r                 !compute growth rate                 
                Nut_taken(bin)=Nut_taken(bin)+(uptlight(i)*mu_max*(N_lag(bin)/(Kn+N_lag(bin))))*q0*weight(i)    !nutrient taken
                weight(i)=weight(i)+Dphyto_lag(i)*time_step*weight(i)                            !compute growth for each particle
                
                 if (photo_model.eq.2) then
                 CL_lag(i)=CL_lag(i)+(1/tau*(Clquot_ACC(bin)-CLQUOT_lag(i))*weight(i)+Dphyto_lag(i)*CL_lag(i))*time_step
                end if

                else                                                                    ! Use Droop's equations
             

               ! do ii=1,Ncat
               ! if ((qlag(i).ge.q_edge(ii)).and.(qlag(i).le.q_edge(ii+1))) then
               !  Q_distrib(bin,ii)=Q_distrib(bin,ii)+1
               ! end if
               ! end do

 
                Uptake_lag(i)=Vmax*(N_lag(bin)/(Kn+N_lag(bin)))*(1-qlag(i)/q_max)       ! nutrient uptake
                Ncell_lag(i)=Ncell_lag(i)+(Uptake_lag(i)-kw*qlag(i))*time_step*weight(i)-mu_r*Ncell_lag(i)*time_step 
                Nut_taken(bin)=Nut_taken(bin)+(Uptake_lag(i)-kw*qlag(i))*weight(i)-mu_r*Ncell_lag(i)      ! nutrient taken


                if (qlag(i).ge.q0) then                                              
                Dphyto_lag(i)=uptlight(i)*mu_max*(1-q0/qlag(i)) -mu_r
               !Dphyto_lag(i)=uptlight(i)*mu_max*(qlag(i)/q_max) -mu_r
                else
                Dphyto_lag(i)=-mortality -mu_r
                end if
                ! do ii=1,Ncat
                !if ((Dphyto_lag(i).ge.mu_edge(ii)).and.(Dphyto_lag(i).le.mu_edge(ii+1))) then
               !  growth_distrib(bin,ii)=growth_distrib(bin,ii)+1
               ! end if
               ! end do

                weight(i)=weight(i)+Dphyto_lag(i)*time_step*weight(i)            !compute growth for each particle
                
        

                end if


               


                if (part_management.eq.1) then                                                  !particle management
                if(weight(i).gt.16*m0) then                                                         !if 4 divisions
                div=div+1                                                                       ! Number of particle to be divided
                Fi(div)=i                                                                       ! save which particle to be divided
                end if                                                                          !
                end if                                                                          !

                nb1(bin)=nb1(bin)+1
                nb(bin)=nb(bin)+weight(i)
                binned_weight(bin)=binned_weight(bin)+weight(i)
                binned_Qlag(bin)=binned_Qlag(bin) + qlag(i)*weight(i)
                binned_growth(bin)=binned_growth(bin) + Dphyto_lag(i)*weight(i)


        end do


        do i=1,Nbin
        Q_distrib(i,:)=Q_distrib(i,:)/nb1(i)
        growth_distrib(i,:)=growth_distrib(i,:)/nb1(i)
        binned_Qlag(i)=binned_Qlag(i)/nb(i)
        binned_growth(i)=binned_growth(i)/nb(i)
        end do

        Nut_taken=Nut_taken/dz

        

        if (part_management.eq.1) then                                                          !particle management
        CALL PART_MGMT(div,Fi)
        end if

!=================================================================================
        end subroutine light_growth