MODULE field_mod
  USE genmod
  IMPLICIT NONE
  
  INTEGER,PARAMETER :: field_T=1
  INTEGER,PARAMETER :: field_U=2
  INTEGER,PARAMETER :: field_Z=3

  INTEGER,PARAMETER :: type_integer=1
  INTEGER,PARAMETER :: type_real=2
  INTEGER,PARAMETER :: type_logical=3
    
  TYPE t_field
    CHARACTER(30)      :: name
    REAL(rstd),POINTER :: rval2d(:)
    REAL(rstd),POINTER :: rval3d(:,:)
    REAL(rstd),POINTER :: rval4d(:,:,:)

    INTEGER,POINTER :: ival2d(:)
    INTEGER,POINTER :: ival3d(:,:)
    INTEGER,POINTER :: ival4d(:,:,:)
    
    LOGICAL,POINTER :: lval2d(:)
    LOGICAL,POINTER :: lval3d(:,:)
    LOGICAL,POINTER :: lval4d(:,:,:)

    INTEGER :: ndim
    INTEGER :: field_type
    INTEGER :: data_type 
    INTEGER :: dim3
    INTEGER :: dim4
  END TYPE t_field   

  INTERFACE get_val
    MODULE PROCEDURE getval_r2d,getval_r3d,getval_r4d, &
                     getval_i2d,getval_i3d,getval_i4d, &
                     getval_l2d,getval_l3d,getval_l4d
  END INTERFACE
                   
  INTERFACE ASSIGNMENT(=)
    MODULE PROCEDURE getval_r2d,getval_r3d,getval_r4d, &
                     getval_i2d,getval_i3d,getval_i4d, &
                     getval_l2d,getval_l3d,getval_l4d 
  END INTERFACE

  PRIVATE :: allocate_field_, deallocate_field_

CONTAINS

  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name)
  USE domain_mod
  USE omp_para
    TYPE(t_field),POINTER :: field(:)
    INTEGER,INTENT(IN) :: field_type
    INTEGER,INTENT(IN) :: data_type
    INTEGER,OPTIONAL :: dim1,dim2
    CHARACTER(*), OPTIONAL :: name
!$OMP BARRIER
!$OMP MASTER
    ALLOCATE(field(ndomain))
!$OMP END MASTER
!$OMP BARRIER
    CALL allocate_field_(domain,assigned_domain, field,field_type,data_type,dim1,dim2,name)
!$OMP BARRIER
  END SUBROUTINE allocate_field

  SUBROUTINE allocate_fields(nfield,field,field_type,data_type,dim1,dim2,name)
  USE domain_mod
  USE omp_para
    INTEGER,INTENT(IN) :: nfield
    TYPE(t_field),POINTER :: field(:,:)
    INTEGER,INTENT(IN) :: field_type
    INTEGER,INTENT(IN) :: data_type
    INTEGER,OPTIONAL :: dim1,dim2
    CHARACTER(*), OPTIONAL :: name
    INTEGER :: i
!$OMP BARRIER
!$OMP MASTER
    ALLOCATE(field(ndomain,nfield))
!$OMP END MASTER
!$OMP BARRIER
    DO i=1,nfield
       CALL allocate_field_(domain,assigned_domain, field(:,i),field_type,data_type,dim1,dim2,name)
    END DO
!$OMP BARRIER
  END SUBROUTINE allocate_fields

  SUBROUTINE allocate_field_glo(field,field_type,data_type,dim1,dim2,name)
    USE domain_mod
    IMPLICIT NONE
    TYPE(t_field),POINTER :: field(:)
    INTEGER,INTENT(IN) :: field_type
    INTEGER,INTENT(IN) :: data_type
    INTEGER,OPTIONAL :: dim1,dim2
    CHARACTER(*), OPTIONAL :: name
    INTEGER :: ind
    INTEGER :: ii_size,jj_size
    LOGICAL :: assigned_dom(ndomain_glo)
    ! ONLY the master thread is allowed to call this routine

    ALLOCATE(field(ndomain_glo))    
    assigned_dom(:)=.TRUE.
    CALL allocate_field_(domain_glo,assigned_dom, field,field_type,data_type,dim1,dim2,name)  

  END SUBROUTINE allocate_field_glo

  SUBROUTINE allocate_field_(dom,assigned_dom, field,field_type,data_type,dim1,dim2,name) 
    USE domain_mod
    USE omp_para
    TYPE(t_domain), INTENT(IN) :: dom(:)
    LOGICAL, INTENT(IN) :: assigned_dom(:)
    TYPE(t_field) :: field(:)
    INTEGER,INTENT(IN) :: field_type
    INTEGER,INTENT(IN) :: data_type
    INTEGER,OPTIONAL :: dim1,dim2
    CHARACTER(*), OPTIONAL :: name
    INTEGER :: ind

    DO ind=1,SIZE(field)
      IF (.NOT. assigned_dom(ind)  .OR. .NOT. is_omp_level_master) CYCLE

      IF(PRESENT(name)) THEN
         field(ind)%name = name
      ELSE
         field(ind)%name = '(undefined)'
      END IF

      IF (PRESENT(dim2)) THEN
        field(ind)%ndim=4 
        field(ind)%dim4=dim2 
        field(ind)%dim3=dim1
      ELSE IF (PRESENT(dim1)) THEN
        field(ind)%ndim=3
        field(ind)%dim3=dim1
      ELSE
        field(ind)%ndim=2
      END IF
    
      field(ind)%data_type=data_type
      field(ind)%field_type=field_type

      CALL allocate_field_XvalY(field(ind), dom(ind))
     END DO
   
 END SUBROUTINE allocate_field_

  SUBROUTINE allocate_field_XvalY(field, dom)
    USE domain_mod, ONLY : t_domain
    USE grid_param
    TYPE(t_field) :: field
    TYPE(t_domain), INTENT(IN) :: dom
    INTEGER :: data_type, dim1, dim2, ij_size
    data_type = field%data_type
    dim1      = field%dim3
    dim2      = field%dim4
    SELECT CASE(grid_type)
    CASE(grid_ico)
       ij_size=dom%iim
       SELECT CASE(field%field_type)
       CASE(field_T)
          ij_size=ij_size*dom%jjm
       CASE(field_U)
          ij_size=3*ij_size*dom%jjm
       CASE(field_Z)
          ij_size=2*ij_size*dom%jjm
       END SELECT
       
       IF (field%ndim==4) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival4d(ij_size,dim1,dim2))
          IF (data_type==type_real)    ALLOCATE(field%rval4d(ij_size,dim1,dim2))
          IF (data_type==type_logical) ALLOCATE(field%lval4d(ij_size,dim1,dim2))
       ELSE IF (field%ndim==3) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival3d(ij_size,dim1))
          IF (data_type==type_real)    ALLOCATE(field%rval3d(ij_size,dim1))
          IF (data_type==type_logical) ALLOCATE(field%lval3d(ij_size,dim1))
       ELSE IF (field%ndim==2) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival2d(ij_size))
          IF (data_type==type_real)    ALLOCATE(field%rval2d(ij_size))
          IF (data_type==type_logical) ALLOCATE(field%lval2d(ij_size))
       ENDIF
       
    CASE(grid_unst)
       PRINT *, 'Allocating field ', field%name
       SELECT CASE(field%field_type)
       CASE(field_T)
          ij_size=primal_num
       CASE(field_U)
          ij_size=edge_num
       CASE(field_Z)
          ij_size=dual_num
       END SELECT

       PRINT *, 'Allocating field ', field%name
       PRINT *, '          with ij_size = ', ij_size
       
       IF (field%ndim==4) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival4d(dim1,ij_size,dim2))
          IF (data_type==type_real)    ALLOCATE(field%rval4d(dim1,ij_size,dim2))
          IF (data_type==type_logical) ALLOCATE(field%lval4d(dim1,ij_size,dim2))
       ELSE IF (field%ndim==3) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival3d(dim1,ij_size))
          IF (data_type==type_real)    ALLOCATE(field%rval3d(dim1,ij_size))
          IF (data_type==type_logical) ALLOCATE(field%lval3d(dim1,ij_size))
       ELSE IF (field%ndim==2) THEN
          IF (data_type==type_integer) ALLOCATE(field%ival2d(ij_size))
          IF (data_type==type_real)    ALLOCATE(field%rval2d(ij_size))
          IF (data_type==type_logical) ALLOCATE(field%lval2d(ij_size))
       ENDIF
       
    END SELECT
  END SUBROUTINE allocate_field_XvalY

  SUBROUTINE deallocate_field(field)
    USE domain_mod
    USE omp_para
    IMPLICIT NONE
    TYPE(t_field),POINTER :: field(:)
    !$OMP BARRIER
    CALL deallocate_field_(assigned_domain, field)
    !$OMP BARRIER
    !$OMP MASTER
    DEALLOCATE(field)
    !$OMP END MASTER
    !$OMP BARRIER
  END SUBROUTINE deallocate_field
  
  SUBROUTINE deallocate_fields(field)
    USE domain_mod
    USE omp_para
    IMPLICIT NONE
    TYPE(t_field),POINTER :: field(:,:)
    INTEGER :: i
    !$OMP BARRIER
    DO i=1,SIZE(field,2)
       CALL deallocate_field_(assigned_domain, field(:,i))
    END DO
    !$OMP BARRIER
    !$OMP MASTER
    DEALLOCATE(field)
    !$OMP END MASTER
    !$OMP BARRIER
  END SUBROUTINE deallocate_fields

  SUBROUTINE deallocate_field_(assigned_dom, field)
  USE omp_para
    IMPLICIT NONE
    LOGICAL, INTENT(IN) :: assigned_dom(:)
    TYPE(t_field) :: field(:)
    INTEGER :: data_type
    INTEGER :: ind
    DO ind=1,SIZE(field)
      IF (.NOT. assigned_dom(ind)  .OR. .NOT. is_omp_level_master) CYCLE

      data_type=field(ind)%data_type
        
      IF (field(ind)%ndim==4) THEN
        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d)
        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d)
        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d)
      ELSE IF (field(ind)%ndim==3) THEN
        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d)
        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d)
        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d)
      ELSE IF (field(ind)%ndim==2) THEN
        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d)
        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d)
        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d)
      ENDIF
      
   ENDDO
  END SUBROUTINE deallocate_field_

  SUBROUTINE deallocate_field_glo(field)
  USE domain_mod
  IMPLICIT NONE
    TYPE(t_field),POINTER :: field(:)
    INTEGER :: data_type
    INTEGER :: ind
    LOGICAL :: assigned_dom(ndomain_glo)
    ! ONLY the master thread is allowed to call this routine
    CALL deallocate_field_(assigned_dom, field)
    DEALLOCATE(field)       
  END SUBROUTINE deallocate_field_glo
    
  SUBROUTINE extract_slice(field_in, field_out, l)  
  USE domain_mod
  USE omp_para
  IMPLICIT NONE  
    TYPE(t_field) :: field_in(:)
    TYPE(t_field) :: field_out(:)
    INTEGER,INTENT(IN) :: l
    
    INTEGER :: ind
    INTEGER :: data_type

!$OMP BARRIER
    DO ind=1,ndomain
      data_type=field_in(ind)%data_type
      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
      
      IF (field_in(ind)%ndim==3 .AND. field_out(ind)%ndim==2) THEN  
        IF (data_type==type_integer)  field_out(ind)%ival2d=field_in(ind)%ival3d(:,l)
        IF (data_type==type_real)     field_out(ind)%rval2d=field_in(ind)%rval3d(:,l)
        IF (data_type==type_logical)  field_out(ind)%lval2d=field_in(ind)%lval3d(:,l)
      ELSE IF  (field_in(ind)%ndim==4 .AND. field_out(ind)%ndim==3) THEN
        IF (data_type==type_integer)  field_out(ind)%ival3d=field_in(ind)%ival4d(:,:,l)
        IF (data_type==type_real)     field_out(ind)%rval3d=field_in(ind)%rval4d(:,:,l)
        IF (data_type==type_logical)  field_out(ind)%lval3d=field_in(ind)%lval4d(:,:,l)
      ELSE
        PRINT *, 'extract_slice : cannot extract slice, dimension incompatible'
        STOP       
      ENDIF
   ENDDO 
!$OMP BARRIER    
  END  SUBROUTINE extract_slice  
  
  
  SUBROUTINE insert_slice(field_in, field_out, l)  
  USE domain_mod
  USE omp_para
  IMPLICIT NONE  
    TYPE(t_field) :: field_in(:)
    TYPE(t_field) :: field_out(:)
    INTEGER,INTENT(IN) :: l
    
    INTEGER :: ind
    INTEGER :: data_type

!$OMP BARRIER
    DO ind=1,ndomain
      data_type=field_in(ind)%data_type
      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
      
      IF (field_in(ind)%ndim==2 .AND. field_out(ind)%ndim==3) THEN  
        IF (data_type==type_integer)  field_out(ind)%ival3d(:,l)=field_in(ind)%ival2d(:)
        IF (data_type==type_real)     field_out(ind)%rval3d(:,l)=field_in(ind)%rval2d(:)
        IF (data_type==type_logical)  field_out(ind)%lval3d(:,l)=field_in(ind)%lval2d(:)
      ELSE IF  (field_in(ind)%ndim==4 .AND. field_out(ind)%ndim==3) THEN
        IF (data_type==type_integer)  field_out(ind)%ival4d(:,:,l)=field_out(ind)%ival3d(:,:)
        IF (data_type==type_real)     field_out(ind)%rval4d(:,:,l)=field_out(ind)%rval3d(:,:)
        IF (data_type==type_logical)  field_out(ind)%lval4d(:,:,l)=field_out(ind)%lval3d(:,:)
      ELSE
        PRINT *, 'extract_slice : cannot insert slice, dimension incompatible'
        STOP       
      ENDIF
   ENDDO 
!$OMP BARRIER    
  
  END SUBROUTINE insert_slice
    
  SUBROUTINE getval_r2d(field_pt,field)
  IMPLICIT NONE  
    REAL(rstd),POINTER,INTENT(INOUT) :: field_pt(:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN
       PRINT *, 'get_val_r2d : bad pointer assignment with ' // TRIM(field%name) 
       STOP
    END IF
    field_pt=>field%rval2d
  END SUBROUTINE  getval_r2d

  SUBROUTINE getval_r3d(field_pt,field)
  IMPLICIT NONE  
    REAL(rstd),POINTER,INTENT(INOUT) :: field_pt(:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN
       PRINT *, 'get_val_r3d : bad pointer assignment with ' // TRIM(field%name) 
       STOP
!       CALL ABORT
    END IF
    field_pt=>field%rval3d
  END SUBROUTINE  getval_r3d

  SUBROUTINE getval_r4d(field_pt,field)
  IMPLICIT NONE  
    REAL(rstd),POINTER,INTENT(INOUT) :: field_pt(:,:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN
       PRINT *, 'get_val_r4d : bad pointer assignment with ' // TRIM(field%name)
       STOP
    END IF
    field_pt=>field%rval4d
  END SUBROUTINE  getval_r4d  

  
  SUBROUTINE getval_i2d(field_pt,field)
  IMPLICIT NONE  
    INTEGER,POINTER,INTENT(INOUT) :: field_pt(:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignment'        
    field_pt=>field%ival2d
  END SUBROUTINE  getval_i2d

  SUBROUTINE getval_i3d(field_pt,field)
  IMPLICIT NONE  
    INTEGER,POINTER,INTENT(INOUT) :: field_pt(:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignment'        
    field_pt=>field%ival3d
  END SUBROUTINE  getval_i3d

  SUBROUTINE getval_i4d(field_pt,field)
  IMPLICIT NONE  
    INTEGER,POINTER,INTENT(INOUT) :: field_pt(:,:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignment'        
    field_pt=>field%ival4d
  END SUBROUTINE  getval_i4d

  SUBROUTINE getval_l2d(field_pt,field)
  IMPLICIT NONE  
    LOGICAL,POINTER,INTENT(INOUT) :: field_pt(:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignment'        
    field_pt=>field%lval2d
  END SUBROUTINE  getval_l2d

  SUBROUTINE getval_l3d(field_pt,field)
  IMPLICIT NONE  
    LOGICAL,POINTER,INTENT(INOUT) :: field_pt(:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignment'        
    field_pt=>field%lval3d
  END SUBROUTINE  getval_l3d

  SUBROUTINE getval_l4d(field_pt,field)
  IMPLICIT NONE  
    LOGICAL,POINTER,INTENT(INOUT) :: field_pt(:,:,:)
    TYPE(t_field),INTENT(IN) :: field
    
    IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignment'        
    field_pt=>field%lval4d
  END SUBROUTINE  getval_l4d    

END MODULE field_mod