MODULE create_limit_unstruct_mod
    PRIVATE
    INTEGER, PARAMETER                             :: lmdep=12

    PUBLIC create_limit_unstruct

CONTAINS


  SUBROUTINE create_limit_unstruct
   USE dimphy
  USE lmdz_xios
  USE ioipsl,             ONLY : ioget_year_len
  USE time_phylmdz_mod, ONLY : annee_ref
  USE indice_sol_mod
  USE phys_state_var_mod
  USE mod_phys_lmdz_para
  USE iniprint_mod_h
  IMPLICIT NONE
    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sic
    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sst
    REAL,    DIMENSION(klon,lmdep)                 :: rugos
    REAL,    DIMENSION(klon,lmdep)                 :: albedo
    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sic_mpi
    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sst_mpi
    REAL,    DIMENSION(klon_mpi,lmdep)             :: rugos_mpi
    REAL,    DIMENSION(klon_mpi,lmdep)             :: albedo_mpi
    INTEGER                                        :: ndays
    REAL                                           :: fi_ice(klon)
    REAL, ALLOCATABLE                              :: sic_year(:,:)
    REAL, ALLOCATABLE                              :: sst_year(:,:)
    REAL, ALLOCATABLE                              :: rugos_year(:,:)
    REAL, ALLOCATABLE                              :: albedo_year(:,:)
    REAL, ALLOCATABLE                              :: pctsrf_t(:,:,:)
    REAL, ALLOCATABLE                              :: phy_bil(:,:)
    REAL, ALLOCATABLE                              :: sst_year_mpi(:,:)
    REAL, ALLOCATABLE                              :: rugos_year_mpi(:,:)
    REAL, ALLOCATABLE                              :: albedo_year_mpi(:,:)
    REAL, ALLOCATABLE                              :: pctsrf_t_mpi(:,:,:)
    REAL, ALLOCATABLE                              :: phy_bil_mpi(:,:)
    INTEGER :: l,k
    INTEGER :: nbad
    INTEGER :: sic_time_axis_size
    INTEGER :: sst_time_axis_size
    CHARACTER(LEN=99)                  :: mess            ! error message


    ndays=ioget_year_len(annee_ref)

    IF (is_omp_master) CALL xios_get_axis_attr("time_sic",n_glo=sic_time_axis_size)
    CALL bcast_omp(sic_time_axis_size)
    ALLOCATE(sic_mpi(klon_mpi,sic_time_axis_size))
    ALLOCATE(sic(klon,sic_time_axis_size))


    IF (is_omp_master) CALL xios_get_axis_attr("time_sst",n_glo=sst_time_axis_size)
    CALL bcast_omp(sst_time_axis_size)
    ALLOCATE(sst_mpi(klon_mpi,sst_time_axis_size))
    ALLOCATE(sst(klon,sst_time_axis_size))

    IF (is_omp_master) THEN
      CALL xios_recv_field("sic_limit",sic_mpi)
      CALL xios_recv_field("sst_limit",sst_mpi)
      CALL xios_recv_field("rugos_limit",rugos_mpi)
      CALL xios_recv_field("albedo_limit",albedo_mpi)
    ENDIF
    CALL scatter_omp(sic_mpi,sic)
    CALL scatter_omp(sst_mpi,sst)
    CALL scatter_omp(rugos_mpi,rugos)
    CALL scatter_omp(albedo_mpi,albedo)

    ALLOCATE(sic_year(klon,ndays))
    ALLOCATE(sst_year(klon,ndays))
    ALLOCATE(rugos_year(klon,ndays))
    ALLOCATE(albedo_year(klon,ndays))
    ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
    ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0


! sic
    IF (sic_time_axis_size==lmdep) THEN
      CALL time_interpolation(ndays,sic,'gregorian',sic_year)
    ELSE IF (sic_time_axis_size==ndays) THEN
      sic_year=sic
    ELSE
      WRITE(mess,*) 'sic time axis is nor montly, nor daily. sic time interpolation ',&
                    'is requiered but is not currently managed'
      CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
    ENDIF

    sic_year(:,:)=sic_year(:,:)/100.  ! convert percent to fraction
    WHERE(sic_year(:,:)>1.0) sic_year(:,:)=1.0    ! Some fractions have some time large negative values
    WHERE(sic_year(:,:)<0.0) sic_year(:,:)=0.0    ! probably better to apply alse this filter before horizontal interpolation

! sst
    IF (sst_time_axis_size==lmdep) THEN
      CALL time_interpolation(ndays,sst,'gregorian',sst_year)
    ELSE IF (sst_time_axis_size==ndays) THEN
      sst_year=sst
    ELSE
      WRITE(mess,*)'sic time axis is nor montly, nor daily. sic time interpolation ',&
                   'is requiered but is not currently managed'
      CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
    ENDIF
    WHERE(sst_year(:,:)<271.38) sst_year(:,:)=271.38


! rugos
    DO l=1, lmdep
      WHERE(NINT(zmasq(:))/=1) rugos(:,l)=0.001
    ENDDO
    CALL time_interpolation(ndays,rugos,'360_day',rugos_year)

! albedo
    CALL time_interpolation(ndays,albedo,'360_day',albedo_year)


    DO k=1,ndays
      fi_ice=sic_year(:,k)
      WHERE(fi_ice>=1.0  ) fi_ice=1.0
      WHERE(fi_ice<EPSFRA) fi_ice=0.0
      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice

!!     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
!!        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
!!     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
!!        pctsrf_t(:,is_sic,k)=fi_ice(:)
!!     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
!     END IF
      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
      WHERE(1.0-zmasq<EPSFRA)
        pctsrf_t(:,is_sic,k)=0.0
        pctsrf_t(:,is_oce,k)=0.0
      ELSEWHERE
        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
          pctsrf_t(:,is_sic,k)=1.0-zmasq
          pctsrf_t(:,is_oce,k)=0.0
        ELSEWHERE
          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
             pctsrf_t(:,is_oce,k)=0.0
             pctsrf_t(:,is_sic,k)=1.0-zmasq
          END WHERE
        END WHERE
      END WHERE
      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
    END DO

    ALLOCATE(sst_year_mpi(klon_mpi,ndays))
    ALLOCATE(rugos_year_mpi(klon_mpi,ndays))
    ALLOCATE(albedo_year_mpi(klon_mpi,ndays))
    ALLOCATE(pctsrf_t_mpi(klon_mpi,nbsrf,ndays))
    ALLOCATE(phy_bil_mpi(klon_mpi,ndays))

    CALL gather_omp(pctsrf_t   , pctsrf_t_mpi)
    CALL gather_omp(sst_year   , sst_year_mpi)
    CALL gather_omp(phy_bil    , phy_bil_mpi)
    CALL gather_omp(albedo_year, albedo_year_mpi)
    CALL gather_omp(rugos_year , rugos_year_mpi)

    IF (is_omp_master) THEN
      CALL xios_send_field("foce_limout",pctsrf_t_mpi(:,is_oce,:))
      CALL xios_send_field("fsic_limout",pctsrf_t_mpi(:,is_sic,:))
      CALL xios_send_field("fter_limout",pctsrf_t_mpi(:,is_ter,:))
      CALL xios_send_field("flic_limout",pctsrf_t_mpi(:,is_lic,:))
      CALL xios_send_field("sst_limout", sst_year_mpi)
      CALL xios_send_field("bils_limout",phy_bil_mpi)
      CALL xios_send_field("alb_limout", albedo_year_mpi)
      CALL xios_send_field("rug_limout", rugos_year_mpi)
    ENDIF
  END SUBROUTINE create_limit_unstruct


  SUBROUTINE time_interpolation(ndays,field_in,calendar,field_out)
  USE pchsp_95_m, only: pchsp_95
  USE pchfe_95_m, only: pchfe_95
  USE arth_m, only: arth
  USE dimphy, ONLY : klon
  USE ioipsl,             ONLY : ioget_year_len
  USE time_phylmdz_mod, ONLY : annee_ref
  USE mod_phys_lmdz_para
  USE iniprint_mod_h
  IMPLICIT NONE

   INTEGER,         INTENT(IN)  :: ndays
   REAL,            INTENT(IN)  :: field_in(klon,lmdep)
   CHARACTER(LEN=*),INTENT(IN)  :: calendar
   REAL,            INTENT(OUT) :: field_out(klon,ndays)
 
   INTEGER :: ndays_in
   REAL    :: timeyear(lmdep)   
   REAL    :: yder(lmdep)   
   INTEGER :: ij,ierr, n_extrap
   LOGICAL :: skip

   CHARACTER (len = 50)         :: modname = 'create_limit_unstruct.time_interpolation'
   CHARACTER (len = 80)         :: abort_message

  
   IF (is_omp_master) ndays_in=year_len(annee_ref, calendar)
   CALL bcast_omp(ndays_in)
   IF (is_omp_master) timeyear=mid_months(annee_ref, calendar, lmdep)
   CALL bcast_omp(timeyear)
    
   n_extrap = 0
   skip=.FALSE.
   DO ij=1,klon
     yder = pchsp_95(timeyear, field_in(ij, :), ibeg=2, iend=2, vc_beg=0., vc_end=0.)
     CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr)
     if (ierr < 0) then
        abort_message='error in pchfe_95'
        CALL abort_physic(modname,abort_message,1)
     endif
     n_extrap = n_extrap + ierr
   END DO
   
   IF (n_extrap /= 0) then
     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
   ENDIF 
  
  
  END SUBROUTINE time_interpolation
  !-------------------------------------------------------------------------------
  !
  FUNCTION year_len(y,cal_in)
  !
  !-------------------------------------------------------------------------------
    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
    IMPLICIT NONE
  !-------------------------------------------------------------------------------
  ! Arguments:
    INTEGER                       :: year_len
    INTEGER,           INTENT(IN) :: y
    CHARACTER(LEN=*),  INTENT(IN) :: cal_in
  !-------------------------------------------------------------------------------
  ! Local variables:
    CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
  !-------------------------------------------------------------------------------
  !--- Getting the input calendar to reset at the end of the function
    CALL ioget_calendar(cal_out)
  
  !--- Unlocking calendar and setting it to wanted one
    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
  
  !--- Getting the number of days in this year
    year_len=ioget_year_len(y)
  
  !--- Back to original calendar
    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
  
  END FUNCTION year_len
  !
  !-------------------------------------------------------------------------------
  
  
  !-------------------------------------------------------------------------------
  !
  FUNCTION mid_months(y,cal_in,nm)
  !
  !-------------------------------------------------------------------------------
    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
    IMPLICIT NONE
  !-------------------------------------------------------------------------------
  ! Arguments:
    INTEGER,                INTENT(IN) :: y               ! year
    CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
    INTEGER,                INTENT(IN) :: nm              ! months/year number
    REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
  !-------------------------------------------------------------------------------
  ! Local variables:
    CHARACTER(LEN=99)                  :: mess            ! error message
    CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
    INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
    INTEGER                            :: m               ! months counter
    INTEGER                            :: nd              ! number of days
    INTEGER                            :: k
  !-------------------------------------------------------------------------------
    nd=year_len(y,cal_in)
  
    IF(nm==12) THEN
  
    !--- Getting the input calendar to reset at the end of the function
      CALL ioget_calendar(cal_out)
  
    !--- Unlocking calendar and setting it to wanted one
      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
  
    !--- Getting the length of each month
      DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
  
    !--- Back to original calendar
      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
  
    ELSE IF(MODULO(nd,nm)/=0) THEN
      WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
        nm,' months/year. Months number should divide days number.'
      CALL abort_physic('mid_months',TRIM(mess),1)
  
    ELSE
      mnth=(/(m,m=1,nm,nd/nm)/)
    END IF
  
  !--- Mid-months times
    mid_months(1)=0.5*REAL(mnth(1))
    DO k=2,nm
      mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
    END DO
  
  END FUNCTION mid_months
  

END MODULE create_limit_unstruct_mod