ADV_DIFF.f90 6.4 KB

!===============================================================================!
!                                                                               !
!			     EULERIAN ADVECTION-DIFFUSION                       !
!                                                                               !
!       Description: In this routine, advection and diffusion are computed      !
!                    using upwind sheme for following variables:                !
!                       - cell quota (for Eulerian model)                       !
!                       - phytoplankton concentration (for Eulerian model)      !
!                       - external nutrients (for both Lagrangian & Eulerian)   !
!                                                                               !
!                                                                               !
!===============================================================================!



        subroutine ADV_DIFF

        USE global_parameter                                    !shared parameters


!_______________________________________________________________________________
!                               LOCAL PARAMETERS
!_______________________________________________________________________________

        IMPLICIT NONE
        double precision, allocatable   :: DC(:)
        double precision, allocatable   :: DN_eul(:)
        double precision, allocatable   :: DN_lag(:)
        double precision, allocatable   :: Dq(:),DQQ(:,:)
        double precision, allocatable   :: q2(:),DNcell(:)
        integer                         :: i,j
        double precision                :: alpha
        double precision                :: beta
        double precision                :: teta
        double precision, allocatable   :: DCL(:)


        Allocate(DC(Nbin))
        Allocate(DN_eul(Nbin))
        Allocate(DN_lag(Nbin))
        Allocate(Dq(Nbin))
        Allocate(q2(Nbin))
         Allocate(DQQ(Nbin,Ncat))
        allocate(DNcell(Nbin))
        Allocate(DCL(Nbin))

        
        if(photo_model.eq.2) then

            do i=2,Nbin-1
        DCL(i)=w_s*(CL(i+1)-CL(i))/dz
        end do

        DCL(1)=w_s*CL(2)/dz
        DCL(Nbin)=-w_s*CL(Nbin)/dz

        CL=CL+DCL*time_step


         do i=2,Nbin-1
        DCL(i)=(K(i)*(CL(i-1)-CL(i))+K(i+1)*(CL(i+1)-CL(i)))/(dz**2) 
        end do

        DCL(1)=(K(2)*(CL(2)-CL(1)))/(dz**2)
        DCL(Nbin)=(K(Nlevel-1)*(CL(Nbin-1)-CL(Nbin)))/(dz**2)

        CL=CL+DCL*time_step

        end if 

        if(bio_model.le.2) then
!________________________________________________________________________________
!                ADVECTION OF CELL QUOTA FOR EULERIAN CALCULATIONS
!________________________________________________________________________________



        do i=2,Nbin-1
        DNcell(i)=w_s*(Ncell(i+1)-Ncell(i))/dz
        end do

        DNcell(1)=w_s*Ncell(2)/dz
        DNcell(Nbin)=-w_s*Ncell(Nbin)/dz

        Ncell=Ncell+DNcell*time_step


!________________________________________________________________________________
!                DIFFUSION OF CELL QUOTA FOR EULERIAN CALCULATIONS
!________________________________________________________________________________



        do i=2,Nbin-1
        DNcell(i)=(K(i)*(Ncell(i-1)-Ncell(i))+K(i+1)*(Ncell(i+1)-Ncell(i)))/(dz**2) 
        end do

        DNcell(1)=(K(2)*(Ncell(2)-Ncell(1)))/(dz**2) 
        DNcell(Nbin)=(K(Nlevel-1)*(Ncell(Nbin-1)-Ncell(Nbin)))/(dz**2) 

        Ncell=Ncell+DNcell*time_step




!________________________________________________________________________________
!                PHYTOPLANKTON ADVECTION  FOR EULERIAN CALCULATIONS
!________________________________________________________________________________


        do i=2,Nbin-1
        DC(i)=w_s*(euler(i+1)-euler(i))/dz
        end do

        DC(1)=w_s*euler(2)/dz
        DC(Nbin)=-w_s*euler(Nbin)/dz

        euler=euler+DC*time_step



!________________________________________________________________________________
!                PHYTOPLANKTON DIFFUSION  FOR EULERIAN CALCULATIONS
!________________________________________________________________________________


        do i=2,Nbin-1
        DC(i)=(K(i)*(euler(i-1)-euler(i))+K(i+1)*(euler(i+1)-euler(i)))/(dz**2)
        end do

        DC(1)=(K(2)*(euler(2)-euler(1)))/(dz**2)
        DC(Nbin)=(K(Nlevel-1)*(euler(Nbin-1)-euler(Nbin)))/(dz**2)

        euler=euler+DC*time_step

       
        else


         do j=1,Ncat
        do i=2,Nbin-1
        DQQ(i,j)=(K(i)*(q_dist(i-1,j)-q_dist(i,j))+K(i+1)*(q_dist(i+1,j)-q_dist(i,j)))/(dz**2)
        end do

        DQQ(1,j)=(K(2)*(q_dist(2,j)-q_dist(1,j)))/(dz**2)
        DQQ(Nbin,j)=(K(Nlevel-1)*(q_dist(Nbin-1,j)-q_dist(Nbin,j)))/(dz**2)

        q_dist(:,j)=q_dist(:,j)+DQQ(:,j)*time_step
        end do

        end if

!________________________________________________________________________________
!                NUTRIENT CONCENTRATION  FOR EULERIAN CALCULATIONS
!________________________________________________________________________________

        do i=2,Nbin-1
        DN_eul(i)=(K(i)*(N_euler(i-1)-N_euler(i))+K(i+1)*(N_euler(i+1)-N_euler(i)))/(dz**2)-Nut(i) 
        end do

        DN_eul(1)=(K(2)*(N_euler(2)-N_euler(1)))/(dz**2) -Nut(1)
        DN_eul(Nbin)=(K(Nlevel-1)*(N_euler(Nbin-1)-N_euler(Nbin)))/(dz**2)-Nut(Nbin)

        N_euler=N_euler+DN_eul*time_step
        N_euler(1:INT(2/dz))=N_euler(1:INT(2/dz))+bottom_flux*time_step

!________________________________________________________________________________
!                NUTRIENT CONCENTRATION   FOR LAGRANGIAN CALCULATIONS
!________________________________________________________________________________


        do i=2,Nbin-1
        DN_lag(i)=(K(i)*(N_lag(i-1)-N_lag(i))+K(i+1)*(N_lag(i+1)-N_lag(i)))/(dz**2) -Nut_taken(i)
        end do

        DN_lag(1)=(K(2)*(N_lag(2)-N_lag(1)))/(dz**2) -Nut_taken(1)
        DN_lag(Nbin)=(K(Nlevel-1)*(N_lag(Nbin-1)-N_lag(Nbin)))/(dz**2) -Nut_taken(Nbin)

        N_lag=N_lag+DN_lag*time_step
        N_lag(1:INT(2/dz))=N_lag(1:INT(2/dz))+bottom_flux*time_step
       
!___________________________________________________________________________________
!                       RELAXATION
!___________________________________________________________________________________

        if (nut_relax.eq.1) then

        N_lag=N_init
        N_euler=N_init

        end if
       
!==================================================================================
        end subroutine ADV_DIFF