source: branches/GRISLIv3/SOURCES/imposed_regions_rsl.f90 @ 483

!> \file imposed_regions_rsl.f90
!! Module permet d'imposer un scenario des variations locales, par regions - pas par pixel
!<

!> \namespace imposed_regions_rsl
!! Module qui permet d'imposer un scenario des variations locales, par regions - pas par pixel
!! \author A. Quiquet
!! \date January 2018
!! @note  used modules
!! @note  - use module3D_phy
!<

module imposed_regions_rsl

! 2 fichiers : 
!  i-  un fichier texte avec nb_reg colonnes et un certain nombre de lignes, comme fileforc dans forclim
!  ii- un fichier netcdf qui contient les regions, numerotees en partant de 0. La region 0 n'est pas touchee, niveau eustatique

  use geography, only:nx,ny
  use io_netcdf_grisli, only:read_ncdf_var

  implicit none


  integer :: nft   !<   number of lines in the rsl history file
  real,dimension(:),allocatable :: tdate         !< time for rsl forcing
  real,dimension(:,:),allocatable :: reg_rsltime !< rsl forcing for a specific region (time,region)

  integer nb_reg                                 !< number of regions
  integer,dimension(nx,ny) :: zregions              !< mask of the different regions

contains

!> SUBROUTINE: input_rsl
!! Routine permettant d'introduire des variations locales du niveau marin
!! en anomalie par rapport au niveau eustatique
!>
  subroutine input_rsl

    use module3D_phy, only:num_param,num_forc,num_rep_42
    use geography, only: dirforcage,dirnameinp
    implicit none

    character(len=8)   :: control   ! label to check clim. forc. file (filin) is usable
    character(len=150) :: filin     ! the file to be read
    character(len=150) :: filinreg  ! the file to be read
    character(len=100) :: filforc   ! the history of rsl, a row = a region
    character(len=100) :: filreg    ! the file of the different regions


    integer :: err                  ! recuperation erreur

    real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointe

    integer i

    namelist/imposed_regions_rsl/nb_reg,filreg,filforc


    rewind(num_param)
    read(num_param,imposed_regions_rsl)

    write(num_rep_42,'(A)')'!___________________________________________________________'
    write(num_rep_42,'(A)') '&imposed_regions_rsl                         ! nom du bloc '
    write(num_rep_42,*)
    write(num_rep_42,*) 'nb_reg = ', nb_reg
    write(num_rep_42,'(A,A,A)') 'filreg  = "',trim(filreg),'"'
    write(num_rep_42,'(A,A,A)') 'filforc = "',trim(filforc),'"'
    write(num_rep_42,*)'/'
    write(num_rep_42,*)

    filin=trim(dirforcage)//trim(filforc)

    open(num_forc,file=filin,status='old')

    read(num_forc,*) control,nft

    ! Determination of file size (line nb), allocation of perturbation array

    if (control.ne.'nb_lines') then
       write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
       write(6,*) 'le nb de lignes et le label de control nb_lines'
       stop
    endif

    if (.not.allocated(tdate)) then
       allocate(tdate(nft),stat=err)
       if (err/=0) then
          print *,"erreur a l'allocation du tableau Tdate dans imposed_regions_rsl",err
          stop 4
       end if
    end if

    if (.not.allocated(reg_rsltime)) then
       allocate(reg_rsltime(nb_reg,nft),stat=err)
       if (err/=0) then
          print *,"erreur a l'allocation du tableau reg_rsltime dans imposed_regions_rsl",err
          stop 4
       end if
    end if

    do i=1,nft
       read(num_forc,*) tdate(i),reg_rsltime(1:nb_reg,i)
    end do
    close(num_forc)

    filinreg=trim(dirnameinp)//trim(filreg)
    call read_ncdf_var('z',filinreg,tab)
    zregions(:,:) = tab(:,:)
    if (minval(zregions(:,:)).ne.0) then
       write(*,*) "Error in imposed_regions_rsl: should be numbered from 0 to nb_reg (min)"
    endif
    if (maxval(zregions(:,:)).ne.nb_reg) then
       write(*,*) "Error in imposed_regions_rsl: should be numbered from 0 to nb_reg (max)"
    endif

  end subroutine input_rsl


!> SUBROUTINE: rsl
!! Routine  
!>
  subroutine rsl
    
    use module3D_phy, only:sealevel,sealevel_2d,time

    implicit none
    
    integer i,j

    real,dimension(0:nb_reg) :: reg_rsl

    reg_rsl(0) = sealevel

    if(time.lt.tdate(1)) then
       reg_rsl(1:nb_reg)=reg_rsltime(:,1)

    else if (time.ge.tdate(nft)) then
       reg_rsl(1:nb_reg)=reg_rsltime(:,nft)

    else
       do i=1,nft-1
          if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general
             reg_rsl(1:nb_reg)=reg_rsltime(:,i)+(reg_rsltime(:,i+1)-reg_rsltime(:,i))*       &
                  (time-tdate(i))/(tdate(i+1)-tdate(i))
             goto 100
          endif
       end do
    endif
100     continue


    reg_rsl(0) = sealevel

    do j=1,ny
       do i=1,nx

          sealevel_2d(i,j)=reg_rsl(zregions(i,j))

       enddo
    enddo


  end subroutine rsl

end module imposed_regions_rsl